diff --git a/drivers/hrldas/PedoTransferSR2006Mod.F90 b/drivers/hrldas/PedoTransferSR2006Mod.F90
index 02090e82..261864cf 100644
--- a/drivers/hrldas/PedoTransferSR2006Mod.F90
+++ b/drivers/hrldas/PedoTransferSR2006Mod.F90
@@ -162,6 +162,9 @@ subroutine PedoTransferSR2006(NoahmpIO, noahmp, Sand, Clay, Orgm)
                   + sr2006_psi_e_b*psi_et         &
                   + sr2006_psi_e_c
     
+    theta_33    = max(10.0**-3.0,theta_33)  ! For numerical stability
+    theta_1500  = max(10.0**-5.0,theta_1500)  ! For numerical stability
+
     ! assign property values
     smcwlt = theta_1500
     smcref = theta_33
diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90
index 96a7105b..94451a74 100644
--- a/src/AtmosForcingMod.F90
+++ b/src/AtmosForcingMod.F90
@@ -151,21 +151,26 @@ subroutine ProcessAtmosForcing(noahmp)
        endif
     endif
 
-    ! wet-bulb scheme (Wang et al., 2019 GRL), C.He, 12/18/2020
+    ! wet-bulb scheme (Wang et al., 2019 GRL), C.He, 12/18/2020, R. Abolafia-Rosnezweig, 02/01/2024
     if ( OptRainSnowPartition == 5 ) then
-       TemperatureDegC = min( 50.0, max(-50.0,(TemperatureAirRefHeight-ConstFreezePoint)) )    ! Kelvin to degree Celsius with limit -50 to +50
-       if ( TemperatureAirRefHeight > ConstFreezePoint ) then
-          LatHeatVap = ConstLatHeatEvap
-       else
-          LatHeatVap = ConstLatHeatSublim
-       endif
-       PsychConst            = ConstHeatCapacAir * PressureAirRefHeight / (0.622 * LatHeatVap)
-       TemperatureWetBulb    = TemperatureDegC - 5.0    ! first guess wetbulb temperature
-       do LoopInd = 1, LoopNum
-          VapPresSat         = 610.8 * exp( (17.27*TemperatureWetBulb) / (237.3+TemperatureWetBulb) )
-          TemperatureWetBulb = TemperatureWetBulb - (VapPresSat - PressureVaporRefHeight) / PsychConst   ! Wang et al., 2019 GRL Eq.2
-       enddo
-       FrozenPrecipFrac      = 1.0 / (1.0 + 6.99e-5 * exp(2.0*(TemperatureWetBulb+3.97)))                ! Wang et al., 2019 GRL Eq. 1
+
+        if ( TemperatureAirRefHeight >= (ConstFreezePoint+10) ) then !avoid numerical errors when temperature is high
+            FrozenPrecipFrac = 0.0
+        else
+            TemperatureDegC = min( 50.0, max(-50.0,(TemperatureAirRefHeight-ConstFreezePoint)) )    ! Kelvin to degree Celsius with limit -50 to +50
+            if ( TemperatureAirRefHeight > ConstFreezePoint ) then
+                LatHeatVap = ConstLatHeatEvap
+            else
+                LatHeatVap = ConstLatHeatSublim
+            endif
+            PsychConst            = ConstHeatCapacAir * PressureAirRefHeight / (0.622 * LatHeatVap)
+            TemperatureWetBulb    = TemperatureDegC - 5.0    ! first guess wetbulb temperature
+            do LoopInd = 1, LoopNum
+                VapPresSat         = 610.8 * exp( (17.27*TemperatureWetBulb) / (237.3+TemperatureWetBulb) )
+                TemperatureWetBulb = TemperatureWetBulb - (VapPresSat - PressureVaporRefHeight) / PsychConst   ! Wang et al., 2019 GRL Eq.2
+            enddo
+            FrozenPrecipFrac      = 1.0 / (1.0 + 6.99e-5 * exp(2.0*(TemperatureWetBulb+3.97)))                ! Wang et al., 2019 GRL Eq. 1
+        endif
     endif
 
     ! rain-snow partitioning