From 63d270584263515bb76ec0c7b0c2afb3561cd77c Mon Sep 17 00:00:00 2001
From: KyleKlenk <kyle.c.klenk@gmail.com>
Date: Fri, 4 Aug 2023 13:59:33 -0600
Subject: [PATCH] Fix issue where GRUs were failing:

scalarCanopyLiq was not initailized properly. It now has the check in read_icond.f90
---
 .../actors/hru_actor/cpp_code/hru_actor.cpp   | 12 +++----
 .../hru_actor/fortran_code/hru_modelRun.f90   | 33 +++++--------------
 build/source/netcdf/read_icond.f90            |  9 +++--
 3 files changed, 22 insertions(+), 32 deletions(-)

diff --git a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
index 306d13a..9a66514 100644
--- a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
+++ b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
@@ -117,9 +117,7 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                 err = Run_HRU(self); // Simulate a Timestep
 
                 if (err != 0) {
-                    aout(self) << "Error: HRU_Actor - Run_HRU - HRU = " << self->state.indxHRU << 
-                        " - indxGRU = " << self->state.indxGRU << " - refGRU = "<< self->state.refGRU << std::endl;
-                    aout(self) << "Error = " << err << "\n";
+                    // We should have already printed the error to the screen if we get here
 
                     self->send(self->state.parent, hru_error::run_physics_unhandleable, self);
                     self->quit();
@@ -326,10 +324,12 @@ int Run_HRU(stateful_actor<hru_state>* self) {
         &self->state.dt_init, 
         &self->state.dt_init_factor,
         &self->state.err);
+
     if (self->state.err != 0) {
-        aout(self) << "Error: RunPhysics - HRU = " << self->state.indxHRU << 
-            " - indxGRU = " << self->state.indxGRU << " - refGRU = " << self->state.refGRU <<
-            " - Timestep = " << self->state.timestep << std::endl;
+        aout(self) << "\033[1;31mError: RunPhysics - HRU = " << self->state.indxHRU 
+                   << " - indxGRU = " << self->state.indxGRU 
+                   << " - refGRU = " << self->state.refGRU 
+                   << " - Timestep = " << self->state.timestep << "\033[0m" << std::endl;
         self->quit();
         return 20;
     }
diff --git a/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90 b/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
index 5268b47..63bcf61 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
@@ -223,11 +223,8 @@ subroutine runPhysics(&
                       notUsed_canopyDepth,            & ! intent(out): NOT USED: canopy depth (m)
                       notUsed_exposedVAI,             & ! intent(out): NOT USED: exposed vegetation area index (m2 m-2)
                       err,cmessage)                     ! intent(out): error control
-      if(err/=0)then
-        message=trim(message)//trim(cmessage)
-        print*, message
-        return
-      endif
+      if(err/=0)then;message=trim(message)//trim(cmessage); print*, char(27),'[33m',message,char(27),'[0m'; return; endif
+
     
       ! save the flag for computing the vegetation fluxes
       if(computeVegFluxFlag)      computeVegFlux = yes
@@ -278,11 +275,8 @@ subroutine runPhysics(&
 
   ! get height at bottom of each soil layer, negative downwards (used in Noah MP)
   allocate(zSoilReverseSign(nSoil),stat=err)
-  if(err/=0)then
-    message=trim(message)//'problem allocating space for zSoilReverseSign'
-    print*, message
-    err=20; return
-  endif
+  if(err/=0)then; message=trim(message)//'problem allocating space for zSoilReverseSign'; print*, char(27),'[33m',message,char(27),'[0m'; return; endif
+
   zSoilReverseSign(:) = -progStruct%var(iLookPROG%iLayerHeight)%dat(nSnow+1:nLayers)
  
   ! populate parameters in Noah-MP modules
@@ -297,11 +291,7 @@ subroutine runPhysics(&
 
   ! deallocate height at bottom of each soil layer(used in Noah MP)
   deallocate(zSoilReverseSign,stat=err)
-  if(err/=0)then
-    message=trim(message)//'problem deallocating space for zSoilReverseSign'
-    print*, message
-    err=20; return
-  endif
+  if(err/=0)then;message=trim(message)//'problem deallocating space for zSoilReverseSign'; print*, char(27),'[33m',message,char(27),'[0m'; return; endif
  
 
   ! overwrite the minimum resistance
@@ -328,12 +318,7 @@ subroutine runPhysics(&
         fluxStruct,         & ! data structure of model fluxes
         tmZoneOffsetFracDay,& ! time zone offset in fractional days
         err,cmessage)       ! error control
-  if(err/=0)then
-    err=20
-    message=trim(message)//trim(cmessage)
-    print*, message
-    return 
-  endif
+  if(err/=0)then;err=20; message=trim(message)//cmessage; print*, char(27),'[33m',message,char(27),'[0m'; return; endif
  
   ! initialize the number of flux calls
   diagStruct%var(iLookDIAG%numFluxCalls)%dat(1) = 0._dp
@@ -361,7 +346,7 @@ subroutine runPhysics(&
                   fluxStruct,         & ! intent(inout): model fluxes for a local HRU
                   ! error control
                   err,cmessage)       ! intent(out): error control
-  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif 
+  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); print*, char(27),'[33m',message,char(27),'[0m'; return; endif;
 
 
   !************************************* End of run_oneHRU *****************************************
@@ -395,7 +380,7 @@ subroutine runPhysics(&
   ! compute water balance for the basin aquifer
   if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then
     message=trim(message)//'multi_driver/bigBucket groundwater code not transferred from old code base yet'
-    err=20; return
+    err=20; print*, char(27),'[33m',message,char(27),'[0m'; return
   end if
 
   ! calculate total runoff depending on whether aquifer is connected
@@ -417,7 +402,7 @@ subroutine runPhysics(&
                   bvarStruct%var(iLookBVAR%averageInstantRunoff)%dat(1),           &  ! intent(out): instantaneous runoff (m s-1)
                   bvarStruct%var(iLookBVAR%averageRoutedRunoff)%dat(1),            &  ! intent(out): routed runoff (m s-1)
                   err,message)                                                                  ! intent(out): error control
-  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); print*, char(27),'[33m',message,char(27),'[0m'; return; endif;
   end associate
  
   !************************************* End of run_oneGRU *****************************************
diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90
index 07223ea..e3f3e01 100644
--- a/build/source/netcdf/read_icond.f90
+++ b/build/source/netcdf/read_icond.f90
@@ -236,7 +236,7 @@ subroutine read_icond(&
   do iVar = 1,size(prog_meta)
 
     ! skip variables that are computed later
-    if(prog_meta(iVar)%varName=='scalarCanopyWat'           .or. &
+    if(prog_meta(iVar)%varName=='scalarCanopyWat'          .or. &
       prog_meta(iVar)%varName=='spectralSnowAlbedoDiffuse' .or. &
       prog_meta(iVar)%varName=='scalarSurfaceTemp'         .or. &
       prog_meta(iVar)%varName=='mLayerVolFracWat'          .or. &
@@ -276,7 +276,12 @@ subroutine read_icond(&
 
     ! fix the snow albedo
     if(progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) < 0._dp)then
-    progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%var(iLookPARAM%albedoMax)%dat(1)
+      progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%var(iLookPARAM%albedoMax)%dat(1)
+    endif
+
+    ! make sure canopy water is positive
+    if(progData%var(iLookPROG%scalarCanopyliq)%dat(1) < 0.0001_rkind)then
+      progData%var(iLookPROG%scalarCanopyliq)%dat(1) = 0.0001_rkind
     endif
 
     ! initialize the spectral albedo
-- 
GitLab