Merge pull request #248 from jmccreight/bug_snow_mass_err

Bug snow mass err

autotest/test_prms_canopy.py

@@ -9,9 +9,9 @@
 from utils_compare import compare_in_memory, compare_netcdfs
 
 # 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
+do_compare_output_files = True
+do_compare_in_memory = False
+rtol = atol = 1e-12
 
 calc_methods = ("numpy", "numba", "fortran")
 params = ("params_sep", "params_one")

prms_src/prms5.2.1/prms/prms_constants.f90

@@ -7,9 +7,11 @@ MODULE PRMS_CONSTANTS
     !! Define real precision and range
     integer, parameter :: dp = REAL64
     !! Define double precision and range
-    real(REAL32), parameter :: CLOSEZERO = EPSILON(0.0)
+    ! https://en.wikipedia.org/wiki/Machine_epsilon
+    ! use values slightly larger than the informal definition
+    real(REAL32), parameter :: CLOSEZERO = 1.20e-7  ! EPSILON(0.0)
     real(REAL32), parameter :: NEARZERO = 1.0E-6
-    real(REAL64), parameter :: DNEARZERO = EPSILON(0.0D0)
+    real(REAL64), parameter :: DNEARZERO = 2.23e-16  ! EPSILON(0.0D0)
 
     integer, parameter :: MAXFILE_LENGTH = 256
     integer, parameter :: MAXLINE_LENGTH = 256

prms_src/prms5.2.1/prms/snowcomp.f90

@@ -1439,7 +1439,7 @@ INTEGER FUNCTION snorun()
 
 
 ! LAST check to clear out all arrays if packwater is gone
-        IF ( Pkwater_equiv(i)<=0.0D0 ) THEN
+        IF ( Pkwater_equiv(i) <= DNEARZERO ) THEN
           IF ( Print_debug>DEBUG_less ) THEN
             IF ( Pkwater_equiv(i)<-DNEARZERO ) &
      &           PRINT *, 'Snowpack problem, pkwater_equiv negative, HRU:', i, ' value:', Pkwater_equiv(i)
@@ -1458,8 +1458,8 @@ INTEGER FUNCTION snorun()
           Snowcov_area(i) = 0.0
           Pk_def(i) = 0.0
           Pk_temp(i) = 0.0
-          Pk_ice(i) = 0.0
-          Freeh2o(i) = 0.0
+          Pk_ice(i) = 0.0D0
+          Freeh2o(i) = 0.0D0
           Snowcov_areasv(i) = 0.0 ! rsr, not in original code
           Ai(i) = 0.0D0
           Frac_swe(i) = 0.0
@@ -2625,9 +2625,9 @@ SUBROUTINE snowevap(Potet_sublim, Potet, Snowcov_area, Snow_evap, &
           IF ( Print_debug>DEBUG_less ) THEN
             IF ( Pkwater_equiv<-DNEARZERO ) &
      &           PRINT *, 'snowpack issue, negative pkwater_equiv in snowevap', Pkwater_equiv
-            Pkwater_equiv = 0.0D0  ! JLM: this is INSIDE a debug statement? will change the answers
             !  is this in the originial source
           ENDIF
+          Pkwater_equiv = 0.0D0  ! JLM: this is INSIDE a debug statement? will change the answers
         ENDIF
 
         Snow_evap = 0.0

pywatershed/atmosphere/prms_atmosphere.py

@@ -8,7 +8,7 @@
 
 from ..base.adapter import adaptable
 from ..base.control import Control
-from ..constants import epsilon, inch2cm, nan, one, zero
+from ..constants import inch2cm, nan, nearzero, one, zero
 from ..parameters import Parameters
 from ..utils.time_utils import datetime_doy, datetime_month
 from .solar_constants import solf
@@ -379,7 +379,7 @@ def adjust_precip(self):
         # Calculate the mix everywhere, then set the precip/rain/snow amounts
         # from the conditions.
         tdiff = self.tmaxf.data - self.tminf.data
-        tdiff = np.where(tdiff < epsilon, 0.000100, tdiff)
+        tdiff = np.where(tdiff < nearzero, 1.0e-4, tdiff)
         self.prmx.data[:] = (
             (self.tmaxf.data - self.tmax_allsnow[month_ind]) / tdiff
         ) * self.adjmix_rain[month_ind]

pywatershed/atmosphere/prms_solar_geometry.py

@@ -8,13 +8,11 @@
 from pywatershed.utils.netcdf_utils import NetCdfWrite
 
 from ..base.control import Control
-from ..constants import epsilon32, nan, one, zero
+from ..constants import dnearzero, nan, one, zero
 from ..parameters import Parameters
 from ..utils.prms5util import load_soltab_debug
 from .solar_constants import ndoy, pi, pi_12, r1, solar_declination, two_pi
 
-epsilon = epsilon32
-
 doy = np.arange(ndoy) + 1
 
 # Dimensions
@@ -226,10 +224,7 @@ def compute_soltab(
 
         # d1 is the denominator of equation 12, Lee, 1963
         d1 = sl_cos * x0_cos - sl_sin * np.sin(x0) * aspects_cos
-        eps_d1 = epsilon
-        wh_d1_lt_eps = np.where(d1 < eps_d1)
-        if len(wh_d1_lt_eps[0]) > 0:
-            d1[wh_d1_lt_eps] = eps_d1
+        d1 = np.where(np.abs(d1) < dnearzero, dnearzero, d1)
 
         # x2 is the difference in longitude between the location of
         # the HRU and the equivalent horizontal surface expressed in angle hour
@@ -307,14 +302,15 @@ def compute_soltab(
             sunh[wh_t6_lt_t1] = (t3 - t2 + t1 - t6)[wh_t6_lt_t1] * pi_12
 
         # The first condition checked
-        wh_sl_zero = np.where(tile_space_to_time(np.abs(sl)) < epsilon)
-        if len(wh_sl_zero[0]):
-            solt[wh_sl_zero] = func3(np.zeros(nhru), x0, t1, t0)[wh_sl_zero]
-            sunh[wh_sl_zero] = (t1 - t0)[wh_sl_zero] * pi_12
-
-        wh_sunh_lt_zero = np.where(sunh < epsilon)
-        if len(wh_sunh_lt_zero[0]):
-            sunh[wh_sunh_lt_zero] = zero
+        mask_sl_lt_dnearzero = tile_space_to_time(np.abs(sl)) < dnearzero
+        solt = np.where(
+            mask_sl_lt_dnearzero, func3(np.zeros(nhru), x0, t1, t0), solt
+        )
+
+        sunh = np.where(mask_sl_lt_dnearzero, (t1 - t0) * pi_12, sunh)
+
+        mask_sunh_lt_dnearzero = sunh < dnearzero
+        sunh = np.where(mask_sunh_lt_dnearzero, zero, sunh)
 
         wh_solt_lt_zero = np.where(solt < zero)
         if len(wh_solt_lt_zero[0]):

pywatershed/constants.py

@@ -26,8 +26,15 @@
 nan = np.nan
 
 epsilon = np.finfo(zero).eps
-epsilon64 = epsilon
-epsilon32 = np.finfo(zero.astype("float32")).eps
+# https://en.wikipedia.org/wiki/Machine_epsilon
+# use values slightly larger than the informal definition
+epsilon64 = 2.23e-16  # epsilon
+epsilon32 = 1.20e-07  # np.finfo(zero.astype("float32")).eps
+
+# These are PRMS conventions, should not be used elsewhere
+nearzero = 1.0e-6
+dnearzero = epsilon64
+closezero = epsilon32
 
 fill_value_f4 = 9.96921e36
 

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

@@ -7,7 +7,15 @@
 from ..base.adapter import adaptable
 from ..base.conservative_process import ConservativeProcess
 from ..base.control import Control
-from ..constants import CovType, HruType, nan, numba_num_threads, zero
+from ..constants import (
+    CovType,
+    HruType,
+    dnearzero,
+    nan,
+    nearzero,
+    numba_num_threads,
+    zero,
+)
 from ..parameters import Parameters
 
 try:
@@ -19,8 +27,6 @@
     has_prmscanopy_f = False
 
 # set constants (may need .value for enum to be used in > comparisons)
-NEARZERO = 1.0e-6
-DNEARZERO = np.finfo(float).eps  # EPSILON(0.0D0)
 BARESOIL = CovType.BARESOIL.value
 GRASSES = CovType.GRASSES.value
 LAND = HruType.LAND.value
@@ -38,7 +44,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 +63,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 +111,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 +238,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
@@ -238,8 +248,8 @@ def _init_calc_method(self):
             #         nb.float64[:],  # wrain_intcp
             #         nb.float64,  # np.float64(time_length),
             #         nb.int64[:],  # self._hru_type
-            #         nb.float64,  # NEARZERO
-            #         nb.float64,  # DNEARZERO,
+            #         nb.float64,  # nearzero
+            #         nb.float64,  # dnearzero,
             #         nb.int32,  # BARESOIL,
             #         nb.int32,  # GRASSES,
             #         nb.int32,  # LAND,
@@ -316,7 +326,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,
@@ -325,8 +336,8 @@ def _calculate(self, time_length):
                 wrain_intcp=self.wrain_intcp,
                 time_length=time_length,
                 hru_type=self._hru_type,
-                NEARZERO=np.float64(NEARZERO),
-                DNEARZERO=np.float64(DNEARZERO),
+                nearzero=nearzero,
+                dnearzero=dnearzero,
                 BARESOIL=np.int32(BARESOIL),
                 GRASSES=np.int32(GRASSES),
                 LAND=np.int32(LAND),
@@ -366,7 +377,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,
@@ -375,8 +387,8 @@ def _calculate(self, time_length):
                 self.wrain_intcp,
                 time_length,
                 self._hru_type,
-                np.float64(NEARZERO),
-                np.float64(DNEARZERO),
+                nearzero,
+                dnearzero,
                 np.int32(BARESOIL),
                 np.int32(GRASSES),
                 np.int32(LAND),
@@ -412,7 +424,8 @@ def _calculate_numpy(
         net_ppt,
         net_rain,
         net_snow,
-        pkwater_ante,
+        pk_ice_prev,
+        freeh2o_prev,
         potet,
         potet_sublim,
         snow_intcp,
@@ -421,8 +434,8 @@ def _calculate_numpy(
         wrain_intcp,
         time_length,
         hru_type,
-        NEARZERO,
-        DNEARZERO,
+        nearzero,
+        dnearzero,
         BARESOIL,
         GRASSES,
         LAND,
@@ -511,11 +524,10 @@ def _calculate_numpy(
                         elif cov_type[i] == GRASSES:
                             # if there is no snowpack and no snowfall, then apparently, grasses
                             # can intercept rain.
-                            # IF ( pkwater_ante(i)<DNEARZERO .AND. netsnow<NEARZERO ) THEN
+                            # 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,
@@ -535,7 +547,7 @@ def _calculate_numpy(
                             intcpstor,
                             netsnow,
                         )
-                        if netsnow < NEARZERO:
+                        if netsnow < nearzero:
                             netrain = netrain + netsnow
                             netsnow = 0.0
                             # todo: deal with newsnow and pptmix?
@@ -548,7 +560,7 @@ def _calculate_numpy(
 
             # if precipitation assume no evaporation or sublimation
             if intcpstor > 0.0:
-                if hru_ppt[i] < NEARZERO:
+                if hru_ppt[i] < nearzero:
                     epan_coef = 1.0
                     evrn = potet[i] / epan_coef
                     evsn = potet[i] * potet_sublim[i]

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_channel.py

@@ -250,7 +250,7 @@ def _initialize_channel_data(self) -> None:
             # not in prms6
             self.seg_slope = np.where(mask_too_flat, 1.0e-4, self.seg_slope)
 
-        # initialize Kcoef to 24.0 for segments with zero velocities
+        # JDH: initialize Kcoef to 24.0 for segments with zero velocities
         # this is different from PRMS, which relied on divide by zero resulting
         # in a value of infinity that when evaluated relative to a maximum
         # desired Kcoef value of 24 would be reset to 24. This approach is