Repository: EC-USGS/pywatershed (language: python)

Issue:
'PRMSSoilzone' object has no attribute '_adjust_parameters'
I'm running my NHM CONUS input files through the latest pywatershed develop version.  I get the exception "AttributeError: 'PRMSSoilzone' object has no attribute '_adjust_parameters'" at pywatershed\hydrology\prms_soilzone.py:319 in _set_initial_conditions

I forgot to say that when I ran the much smaller Willamette model I didn't have this problem.



Discussion / hints:

Fix patch:
diff --git a/bin/prms_linux_gfort_dbl_prec b/bin/prms_linux_gfort_dbl_prec
index 00b2ab53..df8801e2 100755
Binary files a/bin/prms_linux_gfort_dbl_prec and b/bin/prms_linux_gfort_dbl_prec differ
diff --git a/bin/prms_mac_intel_gfort_dbl_prec b/bin/prms_mac_intel_gfort_dbl_prec
index 9abb32c4..23db984e 100755
Binary files a/bin/prms_mac_intel_gfort_dbl_prec and b/bin/prms_mac_intel_gfort_dbl_prec differ
diff --git a/bin/prms_mac_m1_ifort_dbl_prec b/bin/prms_mac_m1_ifort_dbl_prec
index a9e1ee79..3410b374 100755
Binary files a/bin/prms_mac_m1_ifort_dbl_prec and b/bin/prms_mac_m1_ifort_dbl_prec differ
diff --git a/bin/prms_win_gfort_dbl_prec.exe b/bin/prms_win_gfort_dbl_prec.exe
index 9ac58623..f6b139a3 100644
Binary files a/bin/prms_win_gfort_dbl_prec.exe and b/bin/prms_win_gfort_dbl_prec.exe differ
diff --git a/prms_src/prms5.2.1/prms/prms_constants.f90 b/prms_src/prms5.2.1/prms/prms_constants.f90
index 152f1dc5..e5f57421 100644
--- a/prms_src/prms5.2.1/prms/prms_constants.f90
+++ b/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
diff --git a/prms_src/prms5.2.1/prms/snowcomp.f90 b/prms_src/prms5.2.1/prms/snowcomp.f90
index ef2bd151..1a3b9eb7 100644
--- a/prms_src/prms5.2.1/prms/snowcomp.f90
+++ b/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
diff --git a/pywatershed/atmosphere/prms_atmosphere.py b/pywatershed/atmosphere/prms_atmosphere.py
index e6e7be6a..ea065964 100644
--- a/pywatershed/atmosphere/prms_atmosphere.py
+++ b/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]
diff --git a/pywatershed/atmosphere/prms_solar_geometry.py b/pywatershed/atmosphere/prms_solar_geometry.py
index acf3f05e..693c0cc0 100644
--- a/pywatershed/atmosphere/prms_solar_geometry.py
+++ b/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]):
diff --git a/pywatershed/constants.py b/pywatershed/constants.py
index de351298..3c574293 100644
--- a/pywatershed/constants.py
+++ b/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
 
diff --git a/pywatershed/hydrology/prms_canopy.f90 b/pywatershed/hydrology/prms_canopy.f90
index 82cdc7bf..7047363e 100644
--- a/pywatershed/hydrology/prms_canopy.f90
+++ b/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
diff --git a/pywatershed/hydrology/prms_canopy.py b/pywatershed/hydrology/prms_canopy.py
index f8cdbc81..91e2fe08 100644
--- a/pywatershed/hydrology/prms_canopy.py
+++ b/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]
diff --git a/pywatershed/hydrology/prms_canopy.pyf b/pywatershed/hydrology/prms_canopy.pyf
index c6d7f6aa..df9b66b8 100644
--- a/pywatershed/hydrology/prms_canopy.pyf
+++ b/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
diff --git a/pywatershed/hydrology/prms_channel.py b/pywatershed/hydrology/prms_channel.py
index 1ff0330c..1dc6eb9e 100644
--- a/pywatershed/hydrology/prms_channel.py
+++ b/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
diff --git a/pywatershed/hydrology/prms_runoff.py b/pywatershed/hydrology/prms_runoff.py
index ec2c1766..73f7b9d3 100644
--- a/pywatershed/hydrology/prms_runoff.py
+++ b/pywatershed/hydrology/prms_runoff.py
@@ -9,17 +9,14 @@
 from ..base.control import Control
 from ..constants import (
     HruType,
-    epsilon32,
-    epsilon64,
+    dnearzero,
     nan,
+    nearzero,
     numba_num_threads,
     zero,
 )
 from ..parameters import Parameters
 
-NEARZERO = 1.0e-6
-DNEARZERO = np.finfo(float).eps  # EPSILON(0.0D0)
-
 RAIN = 0
 SNOW = 1
 
@@ -329,7 +326,7 @@ def dprst_init(self):
                     open_vol_r = (
                         self.dprst_vol_open[i] / self.dprst_vol_open_max[i]
                     )
-                    if open_vol_r < NEARZERO:
+                    if open_vol_r < nearzero:
                         frac_op_ar = 0.0
                     elif open_vol_r > 1.0:
                         frac_op_ar = 1.0
@@ -348,7 +345,7 @@ def dprst_init(self):
                     clos_vol_r = (
                         self.dprst_vol_clos[i] / self.dprst_vol_clos_max[i]
                     )
-                    if clos_vol_r < NEARZERO:
+                    if clos_vol_r < nearzero:
                         frac_cl_ar = 0.0
                     elif clos_vol_r > 1.0:
                         frac_cl_ar = 1.0
@@ -890,10 +887,10 @@ def compute_infil(
         # procedure is used.  If there is no snowpack and no precip,
         # then check for melt from last of snowpack.  If rain/snow mix
         # with no antecedent snowpack, compute snowmelt portion of runoff.
-        cond3 = snowmelt < NEARZERO
-        cond4 = pkwater_equiv < epsilon64
-        cond1 = net_ppt > zero
-        cond6 = net_snow < NEARZERO
+        # cond3 = snowmelt < nearzero
+        cond4 = pkwater_equiv < dnearzero
+        # cond1 = net_ppt > zero
+        cond6 = net_snow < nearzero
 
         if snowmelt > 0.0:
             # if not cond3:  # this results in less precision accuracy
@@ -902,7 +899,7 @@ def compute_infil(
             avail_water = avail_water + snowmelt
             infil = infil + snowmelt
             if hru_flag == 1:
-                if (pkwater_equiv > 0.0) or (net_rain < NEARZERO):
+                if (pkwater_equiv > 0.0) or (net_rain < nearzero):
                     # Pervious area computations
                     infil, srp = check_capacity(
                         soil_moist_prev,
@@ -929,11 +926,11 @@ def compute_infil(
                         srp,
                     )
 
-        # elif pkwater_equiv < DNEARZERO:
+        # elif pkwater_equiv < Dnearzero:
         elif cond4:
             # If no snowmelt and no snowpack but there was net snow then
             # snowpack was small and was lost to sublimation.
-            # if net_snow < NEARZERO and net_rain > 0.0:
+            # if net_snow < nearzero and net_rain > 0.0:
             if cond6 and through_rain > 0.0:
                 ##  if net_snow < 1.0e-6 and net_rain > 0.0:
                 # cond3 & cond4 & cond6 & cond1
@@ -1038,15 +1035,15 @@ def dprst_comp(
         # of snowpack.
         # If rain/snow mix with no antecedent snowpack, compute snowmelt
         # portion of runoff.
-        if snowmelt:
+        if snowmelt > zero:
             inflow = inflow + snowmelt
 
         # !******There was no snowmelt but a snowpack may exist.  If there is
         # !******no snowpack then check for rain on a snowfree HRU.
-        elif pkwater_equiv < DNEARZERO:
+        elif pkwater_equiv < dnearzero:
             # !      If no snowmelt and no snowpack but there was net snow then
             # !      snowpack was small and was lost to sublimation.
-            if net_snow < NEARZERO and net_rain > 0.0:
+            if net_snow < nearzero and net_rain > 0.0:
                 inflow = inflow + net_rain
 
         dprst_add_water_use = OFF
@@ -1080,7 +1077,7 @@ def dprst_comp(
             srp = srp - dprst_srp / perv_frac
             if srp < zero:
                 srp = zero
-                # if srp < -NEARZERO: this should raise an error
+                # if srp < -nearzero: this should raise an error
 
         if sri > 0.0:
             tmp = sri * imperv_frac * sro_to_dprst_imperv * hru_area
@@ -1094,8 +1091,8 @@ def dprst_comp(
                 dprst_vol_clos = dprst_vol_clos + dprst_sri_clos
             sri = sri - dprst_sri / imperv_frac
             if sri < 0.0:
-                if sri < -NEARZERO:
-                    sri = 0.0
+                # if sri < -nearzero: warn??
+                sri = 0.0
 
         # <<<
         dprst_insroff_hru = dprst_srp + dprst_sri
@@ -1103,7 +1100,7 @@ def dprst_comp(
         dprst_area_open = 0.0
         if dprst_vol_open > 0.0:
             open_vol_r = dprst_vol_open / dprst_vol_open_max
-            if open_vol_r < NEARZERO:
+            if open_vol_r < nearzero:
                 frac_op_ar = 0.0
             elif open_vol_r > 1.0:
                 frac_op_ar = 1.0
@@ -1117,7 +1114,7 @@ def dprst_comp(
             dprst_area_clos = 0.0
             if dprst_vol_clos > 0.0:
                 clos_vol_r = dprst_vol_clos / dprst_vol_clos_max
-                if clos_vol_r < NEARZERO:
+                if clos_vol_r < nearzero:
                     frac_cl_ar = 0.0
                 elif clos_vol_r > 1.0:
                     frac_cl_ar = 1.0
@@ -1126,8 +1123,8 @@ def dprst_comp(
                 dprst_area_clos = dprst_area_clos_max * frac_cl_ar
                 if dprst_area_clos > dprst_area_clos_max:
                     dprst_area_clos = dprst_area_clos_max
-                if dprst_area_clos < NEARZERO:
-                    dprst_area_clos = 0.0
+                # if dprst_area_clos < nearzero:
+                #     dprst_area_clos = 0.0
 
         #
         # evaporate water from depressions based on snowcov_area
@@ -1186,7 +1183,7 @@ def dprst_comp(
                 dprst_vol_open = 0.0
 
         if dprst_area_clos_max > 0.0:
-            if dprst_area_clos > NEARZERO:
+            if dprst_area_clos > nearzero:
                 seep_clos = dprst_vol_clos * dprst_seep_rate_clos
                 dprst_vol_clos = dprst_vol_clos - seep_clos
                 if dprst_vol_clos < 0.0:
diff --git a/pywatershed/hydrology/prms_snow.py b/pywatershed/hydrology/prms_snow.py
index 3e797e45..5eaa3432 100644
--- a/pywatershed/hydrology/prms_snow.py
+++ b/pywatershed/hydrology/prms_snow.py
@@ -9,10 +9,11 @@
 from ..base.control import Control
 from ..constants import (
     HruType,
-    epsilon32,
-    epsilon64,
+    closezero,
+    dnearzero,
     inch2cm,
     nan,
+    nearzero,
     numba_num_threads,
     one,
     zero,
@@ -222,9 +223,7 @@ def get_init_values() -> dict:
             "pk_precip": zero,
             "pk_temp": zero,
             "pksv": zero,
-            "pkwater_ante": nan,
             "pkwater_equiv": nan,
-            "pkwater_equiv_change": nan,
             "pptmix_nopack": False,
             "pss": nan,
             "pst": nan,
@@ -249,15 +248,7 @@ def get_mass_budget_terms():
                 "snowmelt",
                 "through_rain",
             ],
-            "storage_changes": [
-                "freeh2o_change",
-                "pk_ice_change"
-                # eventaully pakwater_equiv_change should be removed
-                # entirely from the code and the metadata.
-                # But there are significant gaps between the above and
-                # pkwater_equiv
-                # "pkwater_equiv_change",
-            ],
+            "storage_changes": ["freeh2o_change", "pk_ice_change"],
         }
 
     @staticmethod
@@ -321,62 +312,78 @@ def _set_initial_conditions(self):
 
         if True:
             # For now there is no restart capability. we'll use the following
-            # line when there is
+            # line above when there is
             # if self.control.options["restart"] in [0, 2, 3]:
 
             # The super().__init__ already set_initial_conditions using its
             # set_initial_conditions
             # Below Im just following PRMS6, will reconcile later with the
             # super (may be redundant).
-            vars_init = [
-                "albedo",
-                "iasw",
-                "int_alb",
-                "iso",
-                "lso",
-                "lst",
-                "mso",
-                "pk_def",
-                "pk_temp",
-                "pksv",
-                "salb",
-                "scrv",
-                "slst",
-                "snowcov_areasv",
-                "snsv",
-            ]
-            for vv in vars_init:
-                self._initialize_var(vv)
-
-            pkweq_gt_zero = self.pkwater_equiv > zero
-            wh_pkweq_gt_zero = np.where(pkweq_gt_zero)
-            self.pk_depth[wh_pkweq_gt_zero] = (
-                self.pkwater_equiv[wh_pkweq_gt_zero]
-                * self.deninv[wh_pkweq_gt_zero]
+            # vars_init = [
+            #     "albedo",
+            #     "iasw",
+            #     "int_alb",
+            #     "iso",
+            #     "lso",
+            #     "lst",
+            #     "mso",
+            #     "pk_def",
+            #     "pk_temp",
+            #     "pksv",
+            #     "salb",
+            #     "scrv",
+            #     "slst",
+            #     "snowcov_areasv",
+            #     "snsv",
+            # ]
+            # for vv in vars_init:
+            #     self._initialize_var(vv)
+
+            mask_pkweq_gt_zero = self.pkwater_equiv > zero
+
+            self.pk_depth = np.where(
+                mask_pkweq_gt_zero,
+                self.pkwater_equiv * self.deninv,
+                self.pk_depth,
             )
-            self.pk_den[wh_pkweq_gt_zero] = (
-                self.pkwater_equiv[wh_pkweq_gt_zero]
-                / self.pk_depth[wh_pkweq_gt_zero]
+
+            self.pk_den = np.where(
+                mask_pkweq_gt_zero,
+                self.pkwater_equiv / self.pk_depth,
+                self.pk_den,
             )
-            self.pk_ice[wh_pkweq_gt_zero] = self.pkwater_equiv[
-                wh_pkweq_gt_zero
-            ]
-            self.freeh2o[wh_pkweq_gt_zero] = (
-                self.pk_ice[wh_pkweq_gt_zero]
-                * self.freeh2o_cap[wh_pkweq_gt_zero]
+
+            self.pk_ice = np.where(
+                mask_pkweq_gt_zero, self.pkwater_equiv, self.pk_ice
             )
-            self.ai[wh_pkweq_gt_zero] = self.pkwater_equiv[
-                wh_pkweq_gt_zero
-            ]  # inches
 
-            ai_gt_snarea_thresh = self.ai > self.snarea_thresh
-            wh_pkweq_gt_zero_and_ai_gt_snth = np.where(
-                pkweq_gt_zero & ai_gt_snarea_thresh
+            self.freeh2o = np.where(
+                mask_pkweq_gt_zero,
+                self.pk_ice * self.freeh2o_cap,
+                self.freeh2o,
+            )
+
+            self.ai = np.where(
+                mask_pkweq_gt_zero,
+                self.pkwater_equiv,
+                self.ai,
+            )
+
+            mask_ai_gt_snarea_thresh = self.ai > self.snarea_thresh
+            self.ai = np.where(
+                mask_pkweq_gt_zero & mask_ai_gt_snarea_thresh,
+                self.snarea_thresh,
+                self.ai,
+            )
+
+            mask_ai_gt_zero = self.ai > dnearzero
+            self.frac_swe = np.where(
+                mask_ai_gt_zero,
+                np.minimum(self.pkwater_equiv / self.ai, 1),
+                self.frac_swe,
             )
-            self.ai[wh_pkweq_gt_zero_and_ai_gt_snth] = self.snarea_thresh[
-                wh_pkweq_gt_zero_and_ai_gt_snth
-            ]
 
+            wh_pkweq_gt_zero = np.where(mask_pkweq_gt_zero)
             for ww in range(len(wh_pkweq_gt_zero[0])):
                 self.snowcov_area[wh_pkweq_gt_zero[ww]] = self.sca_deplcrv(
                     self.snarea_curve_2d[
@@ -447,7 +454,6 @@ def _init_calc_method(self):
         return
 
     def _advance_variables(self) -> None:
-        self.pkwater_ante[:] = self.pkwater_equiv
         self.freeh2o_prev[:] = self.freeh2o
         self.pk_ice_prev[:] = self.pk_ice
         return
@@ -475,7 +481,6 @@ def _calculate(self, simulation_time):
             self.pk_temp[:],
             self.pksv[:],
             self.pkwater_equiv[:],
-            self.pkwater_equiv_change[:],
             self.pptmix_nopack[:],
             self.pss[:],
             self.pst[:],
@@ -551,9 +556,7 @@ def _calculate(self, simulation_time):
             pk_precip=self.pk_precip,
             pk_temp=self.pk_temp,
             pksv=self.pksv,
-            pkwater_ante=self.pkwater_ante,
             pkwater_equiv=self.pkwater_equiv,
-            pkwater_equiv_change=self.pkwater_equiv_change,
             potet=self.potet,
             potet_sublim=self.potet_sublim,
             pptmix=self.pptmix,
@@ -652,9 +655,7 @@ def _calculate_numpy(
         pk_precip,
         pk_temp,
         pksv,
-        pkwater_ante,
         pkwater_equiv,
-        pkwater_equiv_change,
         potet,
         potet_sublim,
         pptmix,
@@ -764,14 +765,14 @@ def _calculate_numpy(
             #     print(f"pk_ice 0 : {pk_ice[dbgind]}")
             #     print(f"tcal 0 : {tcal[dbgind]}")
 
-            if (pkwater_equiv[jj] < epsilon64) and (newsnow[jj] == 0):
+            if (pkwater_equiv[jj] < dnearzero) and (newsnow[jj] == 0):
                 # Skip the HRU if there is no snowpack and no new snow
                 # Reset to be sure it is zero if snowpack melted on last
                 # timestep.
                 snowcov_area[jj] = zero
                 continue
 
-            if newsnow[jj] and (pkwater_equiv[jj] < epsilon64):
+            if newsnow[jj] and (pkwater_equiv[jj] < dnearzero):
                 snowcov_area[jj] = one
 
             # <<
@@ -1000,7 +1001,7 @@ def _calculate_numpy(
                 # <<
                 elif pkwater_equiv[jj] < zero:
                     # if verbose:
-                    #     if pkwater_equiv[jj] < (-1 * epsilon64):
+                    #     if pkwater_equiv[jj] < (-1 * dnearzero):
                     #         print(
                     #             f"snowpack issue 3, negative pkwater_equiv, "
                     #             f"HRU: {jj}, value: {pkwater_equiv[jj]}"
@@ -1043,9 +1044,9 @@ def _calculate_numpy(
 
             # <<<<  The if starting between steps 3 & 4 ends here
             # LAST check to clear out all arrays if packwater is gone
-            if pkwater_equiv[jj] <= zero:
+            if pkwater_equiv[jj] <= dnearzero:
                 # if verbose:
-                #     if pkwater_equiv[jj] < -epsilon64:
+                #     if pkwater_equiv[jj] < -dnearzero:
                 #         print(
                 #             "Snowpack problem, pkwater_equiv negative, "
                 #             f"HRU: {jj}, value: {pkwater_equiv[jj]}"
@@ -1078,18 +1079,18 @@ def _calculate_numpy(
 
         # << end of space loop and previous if
 
-        pkwater_equiv_change[:] = pkwater_equiv - pkwater_ante
         freeh2o_change[:] = freeh2o - freeh2o_prev
         pk_ice_change[:] = pk_ice - pk_ice_prev
 
-        nearzero = 1.0e-6
         cond1 = net_ppt > zero
         cond2 = pptmix_nopack != 0
         cond3 = snowmelt < nearzero
-        cond4 = pkwater_equiv < epsilon32
+        cond4 = pkwater_equiv < dnearzero
         cond5 = snow_evap < nearzero
         cond6 = net_snow < nearzero
+        cond7 = snow_evap > (-1 * (pk_ice_change + freeh2o_change))
         # reverse order from the if statements
+
         through_rain[:] = np.where(
             cond1 & cond3 & cond4 & cond6, net_rain, zero
         )
@@ -1098,6 +1099,16 @@ def _calculate_numpy(
         )
         through_rain[:] = np.where(cond1 & cond2, net_rain, through_rain)
 
+        # This condition does not exist in PRMS as far as I can tell
+        # but is necessary for mass balance
+        # This is when it rains on snow (no new snow) and then snow_evap
+        # consumes the pack during the timestep.
+        through_rain[:] = np.where(
+            cond1 & cond6 & cond7,
+            zero,
+            through_rain,
+        )
+
         return (
             ai,
             albedo,
@@ -1120,7 +1131,6 @@ def _calculate_numpy(
             pk_temp,
             pksv,
             pkwater_equiv,
-            pkwater_equiv_change,
             pptmix_nopack,
             pss,
             pst,
@@ -1230,7 +1240,7 @@ def _calc_ppt_to_pack(
             train = (tmaxc + tmax_allsnow_c_current) * 0.5  # [degrees C]
 
             # Temperatures will differ, depending on the presence of existing
-            # snowpack. should this be epsilon32?
+            # snowpack. should this be closezero?
             if pkwater_equiv > zero:
                 # If there is a snowpack, snow temperature is halfway between
                 # the minimum daily temperature and maximum temperature for
@@ -1248,7 +1258,7 @@ def _calc_ppt_to_pack(
             # If there is any rain, the rain temperature is the average
             # temperature.
             train = tavgc  # [degrees C]
-            if train < epsilon32:
+            if train < closezero:
                 # If average temperature is close to freezing, the rain
                 # temperature is halfway between the maximum daily temperature
                 # and maximum temperature
@@ -1312,7 +1322,7 @@ def _calc_ppt_to_pack(
                     # rain to bring the snowpack to isothermal at 0 degC or
                     # not 3 options below (if-then, elseif, else).
 
-                    if abs(net_rain - pndz) < epsilon32:
+                    if abs(net_rain - pndz) < closezero:
                         # (1.1.1) Exactly enough rain to bring pack to
                         # isothermal...
                         # Heat deficit and temperature of the snowpack go to 0.
@@ -1704,9 +1714,9 @@ def _calc_calin(
                     pss = pkwater_equiv  # [inches]
 
         # <<<
-        else:  # abs(dif) < epsilon64:
+        else:  # abs(dif) < dnearzero:
             # JLM: The test had been equality with zero, changed to
-            #      less than epsilon64.
+            #      less than dnearzero.
             # (2) Just enough heat to overcome heat deficit
             # Set temperature and heat deficit to zero. the pack is "ripe"
             pk_temp = zero  # [degrees C]
@@ -1744,7 +1754,7 @@ def _calc_caloss(
         # Loss of heat is handled differently if there is liquid water in the
         # snowpack or not.
         # 2 options below (if-then, else)
-        if freeh2o < epsilon32:
+        if freeh2o < closezero:
             # (1) No free water exists in pack
             # Heat deficit increases because snow is colder than pack
             # (minus a negative number = plus).
@@ -1768,8 +1778,8 @@ def _calc_caloss(
                 # (2) Only part of free water freezes
                 # The calories absorbed by the new snow freezes some of the
                 # free water (increase in ice, decrease in free water).
-                pk_ice = pk_ice + (-cal / 203.2)  # [inches]
-                freeh2o = freeh2o - (-cal / 203.2)  # [inches]
+                pk_ice = pk_ice - (cal / 203.2)  # [inches]
+                freeh2o = freeh2o + (cal / 203.2)  # [inches]
                 return (
                     freeh2o,
                     pk_def,
@@ -1868,7 +1878,7 @@ def _calc_snowcov(
         # <
         # Calculate the ratio of the current packwater equivalent to the
         # maximum packwater equivalent for the given snowpack.
-        if ai > epsilon64:
+        if ai > dnearzero:
             frac_swe = pkwater_equiv / ai  # [fraction]
             frac_swe = min(one, frac_swe)
         else:
@@ -2403,7 +2413,7 @@ def _calc_step_4(
         calc_snowbal,
         canopy_covden,
         albedo,
-        cecn_coef,  # [control.current_month
+        cecn_coef,  # control.current_month
         cov_type,
         deninv,
         den_max,
@@ -2935,7 +2945,7 @@ def _calc_snowbal(
                 # snowpack. Right now cal isn't used for anything outside this
                 # subroutine, but care should be taken if it is.
         # <<
-        elif qcond < epsilon32:
+        elif qcond < closezero:
             # (2)  There is no heat conduction, qcond = zero
             # If the pack temperature is isothermal at 0 degC, then apply any
             # incoming radiation, condensation (latent heat), and convection
@@ -3082,7 +3092,7 @@ def _calc_snowevap(
         # for evaporation, and if the potential evapotation is enough to
         # completely deplete the snow pack or not.
         # 3 options below (if-then, elseif, else)
-        if ez < epsilon32:
+        if ez < closezero:
             # (1) There is no potential for evaporation...
             snow_evap = 0.0  # [inches]
 
@@ -3095,8 +3105,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:
@@ -3135,7 +3145,7 @@ def _calc_snowevap(
 
             if pkwater_equiv < zero:
                 # if verbose:
-                #     if pkwater_equiv < -epsilon64:
+                #     if pkwater_equiv < -dnearzero:
                 #         print(
                 #             "snowpack issue, negative pkwater_equiv in "
                 #             f"snowevap: {pkwater_equiv}"
@@ -3143,7 +3153,7 @@ def _calc_snowevap(
                 # #  <<
                 # zero of pkwater_equiv is inside a debug statement that
                 # is dosent seem to be triggered in test runs
-                # pkwater_equiv = zero
+                pkwater_equiv = zero
                 pass
 
             # <
@@ -3160,7 +3170,7 @@ def _calc_snowevap(
 
                 if pkwater_equiv < zero:
                     # if verbose:
-                    #     if pkwater_equiv < -epsilon64:
+                    #     if pkwater_equiv < -dnearzero:
                     #         print(
                     #             "snowpack issue 2, negative pkwater_equiv in"
                     #             f"snowevap: {pkwater_equiv}"
@@ -3185,6 +3195,8 @@ def _calc_snowevap(
     @staticmethod
     def set_snow_zero():
         pkwater_equiv = zero
+        pk_ice = zero
+        freeh2o = zero
         pk_depth = zero
         pss = zero
         snsv = zero
@@ -3196,8 +3208,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
diff --git a/pywatershed/hydrology/prms_soilzone.py b/pywatershed/hydrology/prms_soilzone.py
index f0e4182a..4c98f1a8 100644
--- a/pywatershed/hydrology/prms_soilzone.py
+++ b/pywatershed/hydrology/prms_soilzone.py
@@ -11,8 +11,8 @@
     ETType,
     HruType,
     SoilType,
-    epsilon,
     nan,
+    nearzero,
     numba_num_threads,
     one,
     zero,
@@ -89,6 +89,9 @@ def __init__(
         self._set_inputs(locals())
         self._set_options(locals())
 
+        # This uses options
+        self._initialize_soilzone_data()
+
         self._set_budget()
         self._init_calc_method()
 
@@ -214,6 +217,10 @@ def get_mass_budget_terms():
         }
 
     def _set_initial_conditions(self):
+        # this is called in the super before options are set on self
+        pass
+
+    def _initialize_soilzone_data(self):
         # Derived parameters
         # JLM: is this awkward here?
         # JLM: it's definitely awkward to edit a parameter. maybe
@@ -918,7 +925,7 @@ def _calculate_numpy(
 
                 # PRMSIV Step 9
                 # Compute slow contribution to interflow, if any
-                if slow_stor[hh] > epsilon:
+                if slow_stor[hh] > zero:  # single precision?
                     (
                         slow_stor[hh],
                         slow_flow[hh],
@@ -935,7 +942,7 @@ def _calculate_numpy(
                 slow_stor[hh] = availh2o
 
             # <
-            if (slow_stor[hh] > epsilon) and (ssr2gw_rate[hh] > zero):
+            if (slow_stor[hh] > zero) and (ssr2gw_rate[hh] > zero):
                 (
                     ssr_to_gw[hh],
                     slow_stor[hh],
@@ -981,7 +988,7 @@ def _calculate_numpy(
                 dunnianflw_gvr = topfr  # ?? is this right
 
             # <
-            pervactet = zero
+            perv_actet[hh] = zero
 
             # Compute actual evapotranspiration
             if soil_moist[hh] > zero:
@@ -991,7 +998,7 @@ def _calculate_numpy(
                     avail_potet,
                     potet_rechr[hh],
                     potet_lower[hh],
-                    pervactet,
+                    perv_actet[hh],
                 ) = compute_szactet(
                     transp_on[hh],
                     cov_type[hh],
@@ -1007,9 +1014,8 @@ def _calculate_numpy(
                 )
 
             # <
-            hru_actet[hh] = hru_actet[hh] + pervactet * hru_frac_perv[hh]
+            hru_actet[hh] = hru_actet[hh] + perv_actet[hh] * hru_frac_perv[hh]
             avail_potet = potet[hh] - hru_actet[hh]
-            perv_actet[hh] = pervactet
             soil_lower[hh] = soil_moist[hh] - soil_rechr[hh]
 
             if hru_type[hh] == HruType.LAND.value:
@@ -1277,7 +1283,7 @@ def _compute_szactet(
         #   et_type=3    - transpiration plus evaporation
         #   et_type=1    - default... JLM: which is ??
 
-        if avail_potet < epsilon:
+        if avail_potet < nearzero:
             et_type = ETType.ET_DEFAULT  # 1
             avail_potet = zero
 
@@ -1391,5 +1397,5 @@ def _compute_szactet(
             avail_potet,
             potet_rechr,
             potet_lower,
-            et,
+            et,  # -> perv_actet
         )


Test patch:
diff --git a/autotest/test_prms_canopy.py b/autotest/test_prms_canopy.py
index 21fe3d1c..57b95bf6 100644
--- a/autotest/test_prms_canopy.py
+++ b/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")
diff --git a/test_data/generate/prms_diagnostic_variables.py b/test_data/generate/prms_diagnostic_variables.py
index 60ee3416..43da2926 100644
--- a/test_data/generate/prms_diagnostic_variables.py
+++ b/test_data/generate/prms_diagnostic_variables.py
@@ -1,7 +1,7 @@
 import pathlib as pl
 
 import pywatershed as pws
-from pywatershed.constants import epsilon32, nan, zero
+from pywatershed.constants import dnearzero, nearzero, nan, zero
 
 import numpy as np
 import xarray as xr
@@ -193,6 +193,8 @@ def diagnose_final_vars_to_nc(
             "pptmix_nopack",
             "snowmelt",
             "pkwater_equiv",
+            "pk_ice_change",
+            "freeh2o_change",
             "snow_evap",
             "net_snow",
             "net_rain",
@@ -203,14 +205,15 @@ def diagnose_final_vars_to_nc(
             data_file = data_dir / f"{vv}.nc"
             data[vv] = xr.open_dataarray(data_file)
 
-        nearzero = 1.0e-6
-
         cond1 = data["net_ppt"] > zero
         cond2 = data["pptmix_nopack"] != 0
         cond3 = data["snowmelt"] < nearzero
-        cond4 = data["pkwater_equiv"] < epsilon32
+        cond4 = data["pkwater_equiv"] < dnearzero
         cond5 = data["snow_evap"] < nearzero
         cond6 = data["net_snow"] < nearzero
+        cond7 = data["snow_evap"] > -1 * (
+            data["pk_ice_change"] + data["freeh2o_change"]
+        )
 
         through_rain = data["net_rain"] * zero
         # these are in reverse order
@@ -224,6 +227,16 @@ def diagnose_final_vars_to_nc(
             cond1 & cond2, data["net_rain"], through_rain
         )
 
+        # This condition does not exist in PRMS as far as I can tell
+        # but is necessary for mass balance
+        # This is when it rains on snow (no new snow) and then snow_evap
+        # consumes the pack during the timestep.
+        through_rain[:] = np.where(
+            cond1 & cond6 & cond7,
+            zero,
+            through_rain,
+        )
+
         through_rain.to_dataset(name="through_rain").to_netcdf(out_file)
         through_rain.close()