Sync latest develop branch bug fix to LIS noahmp code

README.md

@@ -28,6 +28,8 @@ Technical documentation freely available at http://dx.doi.org/10.5065/ew8g-yr95
 
 **Noah-MP version 5.0 model description paper**:  He, C., Valayamkunnath, P., Barlage, M., Chen, F., Gochis, D., Cabell, R., Schneider, T., Rasmussen, R., Niu, G.-Y., Yang, Z.-L., Niyogi, D., and Ek, M.: Modernizing the open-source community Noah with multi-parameterization options (Noah-MP) land surface model (version 5.0) with enhanced modularity, interoperability, and applicability, Geosci. Model Dev., 16, 5131–5151, https://doi.org/10.5194/gmd-16-5131-2023, 2023.
 
+**Noah-MP development future priority paper**: He, C., Chen, F., Barlage, M., Yang, Z.-L., Wegiel, J. W., Niu, G.-Y., Gochis, D., Mocko, D. M., Abolafia-Rosenzweig, R., Zhang, Z., Lin, T.-S., Valayamkunnath, P., Ek, M., and Niyogi, D. (2023): Enhancing the community Noah-MP land model capabilities for Earth sciences and applications, Bull. Amer. Meteor. Soc., E2023–E2029, https://doi.org/10.1175/BAMS-D-23-0249.1
+
 
 ## Noah-MP GitHub structure
 

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

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+8) ) 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

src/PhenologyMainMod.F90

@@ -111,7 +111,7 @@ subroutine PhenologyMain (noahmp)
     ThicknessCanBury      = min(max(SnowDepth-HeightCanopyBot,0.0), (HeightCanopyTop-HeightCanopyBot))
     CanopyFracSnowBury    = ThicknessCanBury / max(1.0e-06, (HeightCanopyTop-HeightCanopyBot))           ! snow buried fraction
     if ( (HeightCanopyTop > 0.0) .and. (HeightCanopyTop <= 1.0) ) then                                   ! MB: change to 1.0 & 0.2 to reflect changes to HeightCanopyTop in MPTABLE
-       SnowDepthVegBury   = HeightCanopyTop * exp(-SnowDepth / 0.2)
+       SnowDepthVegBury   = HeightCanopyTop * exp(-min(SnowDepth,10.0) / 0.2)     ! CH: add limit to avoid numerical floating issue
        CanopyFracSnowBury = min(SnowDepth, SnowDepthVegBury) / SnowDepthVegBury
     endif
 

src/ResistanceAboveCanopyChen97Mod.F90

@@ -62,9 +62,9 @@ subroutine ResistanceAboveCanopyChen97(noahmp, IterationInd)
 ! local statement functions
     ! LECH'S surface functions
     PSLMU(ZZ) = -0.96 * log(1.0 - 4.5 * ZZ)
-    PSLMS(ZZ) = ZZ * RRIC - 2.076 * (1.0 - 1.0/(ZZ + 1.0))
+    PSLMS(ZZ) = ZZ / RFC - 2.076 * (1.0 - 1.0/(ZZ + 1.0))
     PSLHU(ZZ) = -0.96 * log(1.0 - 4.5 * ZZ)
-    PSLHS(ZZ) = ZZ * RFAC - 2.076 * (1.0 - 1.0/(ZZ + 1.0))
+    PSLHS(ZZ) = ZZ * RFAC - 2.076 * (1.0 - exp(-1.2 * ZZ))
     ! PAULSON'S surface functions
     PSPMU(XX) = -2.0*log( (XX+1.0)*0.5 ) - log( (XX*XX+1.0)*0.5 ) + 2.0*atan(XX) - PIHF
     PSPMS(YY) = 5.0 * YY

src/ResistanceBareGroundChen97Mod.F90

@@ -63,9 +63,10 @@ subroutine ResistanceBareGroundChen97(noahmp, IndIter)
 ! local statement functions
     ! LECH'S surface functions
     PSLMU(ZZ) = -0.96 * log(1.0 - 4.5 * ZZ)
-    PSLMS(ZZ) = ZZ * RRIC - 2.076 * (1.0 - 1.0/(ZZ + 1.0))
+    PSLMS(ZZ) = ZZ / RFC - 2.076 * (1.0 - 1.0/(ZZ + 1.0))
     PSLHU(ZZ) = -0.96 * log(1.0 - 4.5 * ZZ)
-    PSLHS(ZZ) = ZZ * RFAC - 2.076 * (1.0 - 1.0/(ZZ + 1.0))
+    PSLHS(ZZ) = ZZ * RFAC - 2.076 * (1.0 - exp(-1.2 * ZZ))
+
     ! PAULSON'S surface functions
     PSPMU(XX) = -2.0*log( (XX+1.0)*0.5 ) - log( (XX*XX+1.0)*0.5 ) + 2.0*atan(XX) - PIHF
     PSPMS(YY) = 5.0 * YY