pkwater_equiv = pk_ice + freeh2o inconsistencies for values lt double…

… precision epsilon. solve by enforcing in PRMS (and in PWS)

autotest/test_prms_canopy.py

@@ -11,7 +11,7 @@
 # compare in memory (faster) or full output files? or both!
 do_compare_output_files = False
 do_compare_in_memory = True
-rtol = atol = 1.0e-5
+rtol = atol = 1e-13
 
 calc_methods = ("numpy", "numba", "fortran")
 params = ("params_sep", "params_one")

prms_src/prms5.2.1/prms/snowcomp.f90

@@ -1488,6 +1488,14 @@ INTEGER FUNCTION snorun()
         Basin_snowdepth = Basin_snowdepth + Pk_depth(i)*DBLE( Hru_area(i) )
         Basin_tcal = Basin_tcal + DBLE( Tcal(i)*Hru_area(i) )
 
+        ! Enforce consistency between prognostic and diagnostic vars when for these small values
+        ! causes issues with prms_canopy
+        if (pkwater_equiv(i) < DNEARZERO) then
+            pkwater_equiv(i) = 0.0D0
+            pk_ice(i) = 0.0D0
+            freeh2o(i) = 0.0D0
+        end if
+
       ENDDO
 
       ! Area normalize basin totals

pywatershed/hydrology/prms_canopy.f90

@@ -23,7 +23,8 @@ pure subroutine calc_canopy( &
         net_ppt, &  ! in
         net_rain, &  ! in
         net_snow, &  ! in
-        pkwater_ante, &  ! in
+        pk_ice_prev, &  ! in
+        freeh2o_prev, &  ! in
         potet, &  ! in
         potet_sublim, &  ! in
         snow_intcp, &  ! in
@@ -71,7 +72,7 @@ pure subroutine calc_canopy( &
         real(kind=8), intent(in), dimension(nhru) :: &
             covden_sum, covden_win, hru_intcpstor, hru_intcpevap, hru_ppt, &
             hru_rain, hru_snow, intcp_changeover, intcp_evap, intcp_stor,&
-            net_ppt, net_rain, net_snow, pkwater_ante, potet, &
+            net_ppt, net_rain, net_snow, pk_ice_prev, freeh2o_prev, potet, &
             potet_sublim, snow_intcp, srain_intcp, transp_on, wrain_intcp
 
         ! Output vectors
@@ -178,7 +179,7 @@ pure subroutine calc_canopy( &
                         else if (cov_type(ii) == GRASSES) then
                             ! if there is no snowpack and no snowfall, then apparently,
                             ! grasses can intercept rain.
-                            if ((pkwater_ante(ii) < DNEARZERO) &
+                            if ((pk_ice_prev(ii) + freeh2o_prev(ii) < DNEARZERO) &
                                 .and. (netsnow < NEARZERO)) then
                                 intcpstor_in = intcpstor
                                 netrain_in = netrain

pywatershed/hydrology/prms_canopy.py

@@ -20,7 +20,7 @@
 
 # set constants (may need .value for enum to be used in > comparisons)
 NEARZERO = 1.0e-6
-DNEARZERO = np.finfo(float).eps  # EPSILON(0.0D0)
+DNEARZERO = np.finfo(np.float64).eps  # EPSILON(0.0D0)
 BARESOIL = CovType.BARESOIL.value
 GRASSES = CovType.GRASSES.value
 LAND = HruType.LAND.value
@@ -38,7 +38,8 @@ class PRMSCanopy(ConservativeProcess):
         control: a Control object
         discretization: a discretization of class Parameters
         parameters: a parameter object of class Parameters
-        pkwater_ante: Previous snowpack water equivalent on each HRU
+        pk_ice_prev: Previous snowpack ice on each HRU
+        freeh2o_prev: Previous snowpack free water on each HRU
         transp_on: Flag indicating whether transpiration is occurring
             (0=no;1=yes)
         hru_ppt: Precipitation on each HRU
@@ -56,7 +57,8 @@ def __init__(
         control: Control,
         discretization: Parameters,
         parameters: Parameters,
-        pkwater_ante: adaptable,
+        pk_ice_prev: adaptable,
+        freeh2o_prev: adaptable,
         transp_on: adaptable,
         hru_ppt: adaptable,
         hru_rain: adaptable,
@@ -103,7 +105,8 @@ def get_parameters() -> tuple:
     @staticmethod
     def get_inputs() -> tuple:
         return (
-            "pkwater_ante",
+            "pk_ice_prev",
+            "freeh2o_prev",
             "transp_on",
             "hru_ppt",
             "hru_rain",
@@ -229,7 +232,8 @@ def _init_calc_method(self):
             #         nb.float64[:],  # net_ppt
             #         nb.float64[:],  # net_rain
             #         nb.float64[:],  # net_snow
-            #         nb.float64[:],  # self.pkwater_ante
+            #         nb.float64[:],  # self.pk_ice_prev
+            #         nb.float64[:],  # self.freeh2o_prev
             #         nb.float64[:],  # potet
             #         nb.float64[:],  # potet_sublim
             #         nb.float64[:],  # snow_intcp
@@ -316,7 +320,8 @@ def _calculate(self, time_length):
                 net_ppt=self.net_ppt,
                 net_rain=self.net_rain,
                 net_snow=self.net_snow,
-                pkwater_ante=self.pkwater_ante,
+                pk_ice_prev=self.pk_ice_prev,
+                freeh2o_prev=self.freeh2o_prev,
                 potet=self.potet,
                 potet_sublim=self.potet_sublim,
                 snow_intcp=self.snow_intcp,
@@ -366,7 +371,8 @@ def _calculate(self, time_length):
                 self.net_ppt,
                 self.net_rain,
                 self.net_snow,
-                self.pkwater_ante,
+                self.pk_ice_prev,
+                self.freeh2o_prev,
                 self.potet,
                 self.potet_sublim,
                 self.snow_intcp,
@@ -412,7 +418,8 @@ def _calculate_numpy(
         net_ppt,
         net_rain,
         net_snow,
-        pkwater_ante,
+        pk_ice_prev,
+        freeh2o_prev,
         potet,
         potet_sublim,
         snow_intcp,
@@ -513,9 +520,8 @@ def _calculate_numpy(
                             # can intercept rain.
                             # IF ( pkwater_ante(i)<DNEARZERO .AND. netsnow<NEARZERO ) THEN
                             if (
-                                pkwater_ante[i] < DNEARZERO
-                                and netsnow < NEARZERO
-                            ):
+                                pk_ice_prev[i] + freeh2o_prev[i]
+                            ) < DNEARZERO and netsnow < NEARZERO:
                                 intcpstor, netrain = intercept(
                                     hru_rain[i],
                                     stor_max_rain,

pywatershed/hydrology/prms_canopy.pyf

@@ -6,7 +6,7 @@ python module prms_canopy_f ! in
         module canopy ! in :prms_canopy_f:prms_canopy.f90
             real(kind=8), parameter,optional :: zero=0.0
             real(kind=8), parameter,optional :: one=1.0
-            subroutine calc_canopy(nhru,cov_type,covden_sum,covden_win,hru_intcpstor,hru_intcpevap,hru_ppt,hru_rain,hru_snow,intcp_changeover,intcp_evap,intcp_stor,intcp_transp_on,net_ppt,net_rain,net_snow,pkwater_ante,potet,potet_sublim,snow_intcp,srain_intcp,transp_on,wrain_intcp,time_length,hru_type,nearzero,dnearzero,baresoil,grasses,land,lake,rain,snow,off,active,intcp_evap_out,intcp_form_out,intcp_stor_out,net_rain_out,net_snow_out,net_ppt_out,hru_intcpstor_out,hru_intcpevap_out,intcp_changeover_out,intcp_transp_on_out) ! in :prms_canopy_f:prms_canopy.f90:canopy
+            subroutine calc_canopy(nhru,cov_type,covden_sum,covden_win,hru_intcpstor,hru_intcpevap,hru_ppt,hru_rain,hru_snow,intcp_changeover,intcp_evap,intcp_stor,intcp_transp_on,net_ppt,net_rain,net_snow,pk_ice_prev, freeh2o_prev,potet,potet_sublim,snow_intcp,srain_intcp,transp_on,wrain_intcp,time_length,hru_type,nearzero,dnearzero,baresoil,grasses,land,lake,rain,snow,off,active,intcp_evap_out,intcp_form_out,intcp_stor_out,net_rain_out,net_snow_out,net_ppt_out,hru_intcpstor_out,hru_intcpevap_out,intcp_changeover_out,intcp_transp_on_out) ! in :prms_canopy_f:prms_canopy.f90:canopy
                 !integer(kind=4) intent(in) :: nhru
                 integer, optional,intent(in),check(len(cov_type)>=nhru),depend(cov_type) :: nhru=len(cov_type)                
                 integer(kind=8) dimension(nhru),intent(in) :: cov_type
@@ -24,7 +24,8 @@ python module prms_canopy_f ! in
                 real(kind=8) dimension(nhru),intent(in) :: net_ppt
                 real(kind=8) dimension(nhru),intent(in) :: net_rain
                 real(kind=8) dimension(nhru),intent(in) :: net_snow
-                real(kind=8) dimension(nhru),intent(in) :: pkwater_ante
+                real(kind=8) dimension(nhru),intent(in) :: pk_ice_prev
+                real(kind=8) dimension(nhru),intent(in) :: freeh2o_prev
                 real(kind=8) dimension(nhru),intent(in) :: potet
                 real(kind=8) dimension(nhru),intent(in) :: potet_sublim
                 real(kind=8) dimension(nhru),intent(in) :: snow_intcp

pywatershed/hydrology/prms_snow.py

@@ -1058,6 +1058,12 @@ def _calculate_numpy(
 
         # << end of space loop and previous if
 
+        mask_pkwater_lt_eps = pkwater_equiv < epsilon64
+        if mask_pkwater_lt_eps.any():
+            pkwater_equiv = np.where(mask_pkwater_lt_eps, zero, pkwater_equiv)
+            pk_ice = np.where(mask_pkwater_lt_eps, zero, pk_ice)
+            freeh2o = np.where(mask_pkwater_lt_eps, zero, freeh2o)
+
         freeh2o_change[:] = freeh2o - freeh2o_prev
         pk_ice_change[:] = pk_ice - pk_ice_prev
 
@@ -3085,8 +3091,8 @@ def _calc_snowevap(
             snow_evap = pkwater_equiv  # [inches]
             pkwater_equiv = zero  # [inches]
             pk_ice = zero  # [inches]
-            pk_def = zero  # [cal/cm^2]
             freeh2o = zero  # [inches]
+            pk_def = zero  # [cal/cm^2]
             pk_temp = zero  # [degrees C]
 
         else:
@@ -3175,6 +3181,8 @@ def _calc_snowevap(
     @staticmethod
     def set_snow_zero():
         pkwater_equiv = zero
+        pk_ice = zero
+        freeh2o = zero
         pk_depth = zero
         pss = zero
         snsv = zero
@@ -3186,8 +3194,6 @@ def set_snow_zero():
         snowcov_area = zero
         pk_def = zero
         pk_temp = zero
-        pk_ice = zero
-        freeh2o = zero
         snowcov_areasv = zero
         ai = zero
         frac_swe = zero