Develop

README.md

@@ -1,13 +1,12 @@
-![image](https://user-images.githubusercontent.com/43385564/236596277-45d492e7-79fd-4006-83e2-3044106afd94.png)
+![noahmp_logo_update](https://github.com/NCAR/noahmp/assets/43385564/1fb47fc2-99bd-4360-9ed0-6d5656c29626)
 
 
 [![DOI](https://zenodo.org/badge/236657733.svg)](https://zenodo.org/badge/latestdoi/236657733)
 
 
-# Noah-MP Community Model Repository
+# Noah-MP<sup>®</sup> Community Model Repository
 
-
-Noah-MP is a widely-used state-of-the-art land surface model used in many research and operational weather/climate models (e.g., HRLDAS, WRF, MPAS, WRF-Hydro/NWM, NOAA/UFS, NASA/LIS, etc.).
+Noah-MP<sup>®</sup> is a widely-used state-of-the-art land surface model used in many research and operational weather/climate models (e.g., HRLDAS, WRF, MPAS, WRF-Hydro/NWM, NOAA/UFS, NASA/LIS, etc.).
 
 This is the official Noah-MP land surface model unified repository for code downloading and contribution. Noah-MP is a community open-source model developed with the contributions from the entire scientific community. For development, maintenance, and release of the community Noah-MP GitHub code, please contact: Cenlin He (cenlinhe@ucar.edu) and Fei Chen (feichen@ucar.edu).
 
@@ -27,7 +26,7 @@ Technical documentation freely available at http://dx.doi.org/10.5065/ew8g-yr95
 
 **Original Noah-MP model description paper**:   Niu, G. Y., Yang, Z. L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., ... & Xia, Y. (2011). The community Noah land surface model with multiparameterization options (Noah‐MP): 1. Model description and evaluation with local‐scale measurements. Journal of Geophysical Research: Atmospheres, 116(D12).
 
-**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-MP land surface model (version 5.0) with enhanced modularity, interoperability, and applicability, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2023-675, 2023.
+**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 GitHub structure

drivers/hrldas/ConfigVarInTransferMod.F90

@@ -146,11 +146,12 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
 
     ! treatment for urban point
     if ( (NoahmpIO%IVGTYP(I,J) == NoahmpIO%ISURBAN_TABLE) .or. (NoahmpIO%IVGTYP(I,J) > NoahmpIO%URBTYPE_beg) ) then
-       noahmp%config%domain%FlagUrban = .true. 
-       if(NoahmpIO%SF_URBAN_PHYSICS == 0 ) then
-           noahmp%config%domain%VegType = NoahmpIO%ISURBAN_TABLE
+       if ( NoahmpIO%SF_URBAN_PHYSICS == 0 ) then
+           noahmp%config%domain%VegType = NoahmpIO%ISURBAN_TABLE  ! treat as bulk urban point
+           noahmp%config%domain%FlagUrban = .true.
        else
-           noahmp%config%domain%VegType = NoahmpIO%NATURAL_TABLE  ! set urban vegetation type based on table natural
+           noahmp%config%domain%VegType = NoahmpIO%NATURAL_TABLE  ! set rural vegetation type based on table natural
+                                                                  ! urban is handled by explicit urban scheme outside Noah-MP
            NoahmpIO%GVFMAX(I,J)         = 0.96 * 100.0            ! unit: %
        endif         
     endif

drivers/hrldas/EnergyVarOutTransferMod.F90

@@ -137,13 +137,11 @@ subroutine EnergyVarOutTransfer(noahmp, NoahmpIO)
     NoahmpIO%CHB2XY  (I,J) = noahmp%energy%state%ExchCoeffSh2mBare
     NoahmpIO%Q2MVXY  (I,J) = noahmp%energy%state%SpecHumidity2mVeg /(1.0-noahmp%energy%state%SpecHumidity2mVeg)  ! spec humidity to mixing ratio
     NoahmpIO%Q2MBXY  (I,J) = noahmp%energy%state%SpecHumidity2mBare/(1.0-noahmp%energy%state%SpecHumidity2mBare)
+    NoahmpIO%ALBEDO  (I,J) = noahmp%energy%state%AlbedoSfc
     NoahmpIO%IRRSPLH (I,J) = NoahmpIO%IRRSPLH(I,J) + &
                              (noahmp%energy%flux%HeatLatentIrriEvap * noahmp%config%domain%MainTimeStep)
     NoahmpIO%TSLB    (I,1:NumSoilLayer,J)       = noahmp%energy%state%TemperatureSoilSnow(1:NumSoilLayer)
     NoahmpIO%TSNOXY  (I,-NumSnowLayerMax+1:0,J) = noahmp%energy%state%TemperatureSoilSnow(-NumSnowLayerMax+1:0)
-    if ( noahmp%energy%state%AlbedoSfc > -999 ) then
-       NoahmpIO%ALBEDO(I,J) = noahmp%energy%state%AlbedoSfc
-    endif
 
     !SNICAR
     if (noahmp%config%nmlist%OptSnowAlbedo == 3 )then

drivers/hrldas/NoahmpInitMainMod.F90

@@ -187,8 +187,13 @@ subroutine NoahmpInitMain(NoahmpIO)
                 NoahmpIO%LAI(I,J)      = 0.0
                 NoahmpIO%LAI(I,J)      = max(NoahmpIO%LAI(I,J), 0.05)                       ! at least start with 0.05 for arbitrary initialization (v3.7)
                 NoahmpIO%XSAIXY(I,J)   = max(0.1*NoahmpIO%LAI(I,J), 0.05)                   ! MB: arbitrarily initialize SAI using input LAI (v3.7)
-                NoahmpIO%LFMASSXY(I,J) = NoahmpIO%LAI(I,J) * 1000.0 / &
-                                         max(NoahmpIO%SLA_TABLE(NoahmpIO%IVGTYP(I,J)),1.0)  ! use LAI to initialize (v3.7)
+                if ( urbanpt_flag .eqv. .true. ) then
+                   NoahmpIO%LFMASSXY(I,J) = NoahmpIO%LAI(I,J) * 1000.0 / &
+                                            max(NoahmpIO%SLA_TABLE(NoahmpIO%NATURAL_TABLE),1.0)! use LAI to initialize (v3.7)
+                else
+                   NoahmpIO%LFMASSXY(I,J) = NoahmpIO%LAI(I,J) * 1000.0 / &
+                                            max(NoahmpIO%SLA_TABLE(NoahmpIO%IVGTYP(I,J)),1.0)  ! use LAI to initialize (v3.7)
+                endif
                 NoahmpIO%STMASSXY(I,J) = NoahmpIO%XSAIXY(I,J) * 1000.0 / 3.0                ! use SAI to initialize (v3.7)
                 NoahmpIO%RTMASSXY(I,J) = 500.0                                              ! these are all arbitrary and probably should be
                 NoahmpIO%WOODXY(I,J)   = 500.0                                              ! in the table or read from initialization

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

parameters/NoahmpTable.TBL

@@ -47,7 +47,7 @@
  ISCROP       =  2  ! crop land type in USGS
  EBLFOREST    = 13  ! evergreen broadleaf forest land type in USGS
  NATURAL      =  5  ! natural vegation type in urban pixel in USGS
- URBTYPE_beg  = 40  ! land type number above which are urban (e.g., LCZ)
+ URBTYPE_beg  = 50  ! land type number above which are urban (e.g., LCZ)
  LCZ_1        = 51  ! urban local climate zone (LCZ) type 1: compact highrise
  LCZ_2        = 52  ! urban local climate zone (LCZ) type 2: compact midrise
  LCZ_3        = 53  ! urban local climate zone (LCZ) type 3: compact lowrise
@@ -243,7 +243,7 @@
  ISCROP       = 12  ! crop land type in MODIS
  EBLFOREST    =  2  ! evergreen broadleaf forest land type in MODIS
  NATURAL      = 14  ! natural vegation type in urban pixel in MODIS
- URBTYPE_beg  = 40  ! land type number above which are urban (e.g., LCZ)
+ URBTYPE_beg  = 50  ! land type number above which are urban (e.g., LCZ)
  LCZ_1        = 51  ! urban local climate zone (LCZ) type 1: compact highrise
  LCZ_2        = 52  ! urban local climate zone (LCZ) type 2: compact midrise
  LCZ_3        = 53  ! urban local climate zone (LCZ) type 3: compact lowrise

src/AtmosForcingMod.F90

@@ -166,21 +166,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

src/CanopyHydrologyMod.F90

@@ -31,6 +31,7 @@ subroutine CanopyHydrology(noahmp)
               LeafAreaIndEff    => noahmp%energy%state%LeafAreaIndEff    ,& ! in,    leaf area index, after burying by snow
               StemAreaIndEff    => noahmp%energy%state%StemAreaIndEff    ,& ! in,    stem area index, after burying by snow
               FlagFrozenCanopy  => noahmp%energy%state%FlagFrozenCanopy  ,& ! in,    used to define latent heat pathway
+              VegFrac           => noahmp%energy%state%VegFrac           ,& ! in,    greeness vegetation fraction
               SnowfallDensity   => noahmp%water%state%SnowfallDensity    ,& ! in,    bulk density of snowfall [kg/m3]
               CanopyLiqHoldCap  => noahmp%water%param%CanopyLiqHoldCap   ,& ! in,    maximum intercepted liquid water per unit veg area index [mm]
               VegFrac           => noahmp%energy%state%VegFrac           ,& ! in,    greeness vegetation fraction
@@ -68,7 +69,7 @@ subroutine CanopyHydrology(noahmp)
 
     ! canopy liquid water
     ! maximum canopy intercepted water
-    CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff)
+    CanopyLiqWaterMax =  VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff)
 
     ! canopy evaporation, transpiration, and dew
     if ( FlagFrozenCanopy .eqv. .false. ) then    ! Barlage: change to FlagFrozenCanopy

src/EnergyMainGlacierMod.F90

@@ -163,7 +163,7 @@ subroutine EnergyMainGlacier(noahmp)
     if ( RadSwDownRefHeight > 0.0 ) then
        AlbedoSfc = RadSwReflSfc / RadSwDownRefHeight
     else
-       AlbedoSfc = -999.9
+       AlbedoSfc = undefined_real
     endif
 
     end associate

src/EnergyMainMod.F90

@@ -340,7 +340,7 @@ subroutine EnergyMain(noahmp)
     if ( RadSwDownRefHeight > 0.0 ) then
        AlbedoSfc = RadSwReflSfc / RadSwDownRefHeight
     else
-       AlbedoSfc = -999.9
+       AlbedoSfc = undefined_real
     endif
     end associate
 

src/ResistanceCanopyStomataBallBerryMod.F90

@@ -83,6 +83,7 @@ subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade)
               PressureAtmosO2         => noahmp%energy%state%PressureAtmosO2         ,& ! in,  atmospheric o2 pressure [Pa]
               PressureAtmosCO2        => noahmp%energy%state%PressureAtmosCO2        ,& ! in,  atmospheric co2 pressure [Pa]
               ResistanceLeafBoundary  => noahmp%energy%state%ResistanceLeafBoundary  ,& ! in,  leaf boundary layer resistance [s/m]
+              VegFrac                 => noahmp%energy%state%VegFrac                 ,& ! in,  greeness vegetation fraction
               RadPhotoActAbsSunlit    => noahmp%energy%flux%RadPhotoActAbsSunlit     ,& ! in,  average absorbed par for sunlit leaves [W/m2]
               RadPhotoActAbsShade     => noahmp%energy%flux%RadPhotoActAbsShade      ,& ! in,  average absorbed par for shaded leaves [W/m2]
               ResistanceStomataSunlit => noahmp%energy%state%ResistanceStomataSunlit ,& ! out, sunlit leaf stomatal resistance [s/m]
@@ -101,7 +102,7 @@ subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade)
     ResistanceStomataTmp = 1.0 / ConductanceLeafMin * CF
     PhotosynLeafTmp      = 0.0
     if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit / max(VegFrac,1.0e-6)  ! Sunlit case
-    if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade / max(VegFrac,1.0e-6)  ! Shaded case
+    if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade  / max(VegFrac,1.0e-6)  ! Shaded case
 
     ! only compute when there is radiation absorption
     if ( RadPhotoActAbsTmp > 0.0 ) then

src/ResistanceCanopyStomataJarvisMod.F90

@@ -68,8 +68,9 @@ subroutine ResistanceCanopyStomataJarvis(noahmp, IndexShade)
     ResistanceTemp       = 0.0
     ResistanceVapDef     = 0.0
     ResistanceStomataTmp = 0.0
-    if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit / max(VegFrac,1.0e-6)  ! Sunlit case
-    if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade / max(VegFrac,1.0e-6)  ! Shaded case
+    if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit / max(VegFrac,1.0e-6) ! Sunlit case
+    if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade  / max(VegFrac,1.0e-6) ! Shaded case
+
 
     ! compute MixingRatioTmp and MixingRatioSat
     SpecHumidityTmp = 0.622 * PressureVaporCanAir / (PressureAirRefHeight - 0.378*PressureVaporCanAir) ! specific humidity

src/SnowLayerCombineMod.F90

@@ -85,8 +85,8 @@ subroutine SnowLayerCombine(noahmp)
              endif
 
           else
+             if ( NumSnowLayerNeg < -1 ) then    ! MB/KM: change to NumSnowLayerNeg
 !             if ( NumSnowLayerOld < -1 ) then    ! MB/KM: change to NumSnowLayerNeg !samlin
-             if ( NumSnowLayerNeg < -1 ) then    ! MB/KM: change to NumSnowLayerNeg !samlin
                 SnowLiqWater(J-1)           = SnowLiqWater(J-1) + SnowLiqWater(J)
                 SnowIce(J-1)                = SnowIce(J-1) + SnowIce(J)
                 ThicknessSnowSoilLayer(J-1) = ThicknessSnowSoilLayer(J-1) + ThicknessSnowSoilLayer(J)

utility/Machine.F90

@@ -18,9 +18,9 @@ module Machine
   integer, public, parameter :: kind_noahmp = 4 ! single precision
 #endif
 
-  integer,                public, parameter :: undefined_int  = huge(1)    ! undefined integer for variable initialization
-  real(kind=kind_noahmp), public, parameter :: undefined_real = huge(1.0)  ! undefined real for variable initializatin
-  integer,                public, parameter :: undefined_int_neg  = -999   ! undefined integer negative for variable initialization
-  real(kind=kind_noahmp), public, parameter :: undefined_real_neg = -999.0 ! undefined real negative for variable initializatin
+  integer,                public, parameter :: undefined_int  = -9999       ! undefined integer for variable initialization
+  real(kind=kind_noahmp), public, parameter :: undefined_real = -9999.0     ! undefined real for variable initializatin
+  integer,                public, parameter :: undefined_int_neg  = -9999   ! undefined integer negative for variable initialization
+  real(kind=kind_noahmp), public, parameter :: undefined_real_neg = -9999.0 ! undefined real negative for variable initializatin
 
 end module Machine