Merge pull request #124 from NCAR/develop

Sync with Develop branch for bug fix

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
 

src/AtmosForcingMod.F90

@@ -154,7 +154,7 @@ subroutine ProcessAtmosForcing(noahmp)
     ! wet-bulb scheme (Wang et al., 2019 GRL), C.He, 12/18/2020, R. Abolafia-Rosnezweig, 02/01/2024
     if ( OptRainSnowPartition == 5 ) then
 
-        if ( TemperatureAirRefHeight >= (ConstFreezePoint+10) ) then !avoid numerical errors when temperature is high
+        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

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