diff --git a/README.md b/README.md index 993659b3..de0a3fb9 100644 --- a/README.md +++ b/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® 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® 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 diff --git a/drivers/hrldas/ConfigVarInTransferMod.F90 b/drivers/hrldas/ConfigVarInTransferMod.F90 index d29396a9..0218b97c 100644 --- a/drivers/hrldas/ConfigVarInTransferMod.F90 +++ b/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 diff --git a/drivers/hrldas/EnergyVarOutTransferMod.F90 b/drivers/hrldas/EnergyVarOutTransferMod.F90 index 21c855e1..1bb72971 100644 --- a/drivers/hrldas/EnergyVarOutTransferMod.F90 +++ b/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 diff --git a/drivers/hrldas/NoahmpInitMainMod.F90 b/drivers/hrldas/NoahmpInitMainMod.F90 index f753391d..efdd9f3e 100644 --- a/drivers/hrldas/NoahmpInitMainMod.F90 +++ b/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 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/parameters/NoahmpTable.TBL b/parameters/NoahmpTable.TBL index a49f5c96..2ff46d21 100644 --- a/parameters/NoahmpTable.TBL +++ b/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 diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index f4ee49cb..488689d5 100644 --- a/src/AtmosForcingMod.F90 +++ b/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 diff --git a/src/CanopyHydrologyMod.F90 b/src/CanopyHydrologyMod.F90 index 6bf57453..caae6b86 100644 --- a/src/CanopyHydrologyMod.F90 +++ b/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 diff --git a/src/EnergyMainGlacierMod.F90 b/src/EnergyMainGlacierMod.F90 index 25958bcd..3fc0bf07 100644 --- a/src/EnergyMainGlacierMod.F90 +++ b/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 diff --git a/src/EnergyMainMod.F90 b/src/EnergyMainMod.F90 index 42c74f14..b040014d 100644 --- a/src/EnergyMainMod.F90 +++ b/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 diff --git a/src/ResistanceCanopyStomataBallBerryMod.F90 b/src/ResistanceCanopyStomataBallBerryMod.F90 index 81c7a80d..e094deaf 100644 --- a/src/ResistanceCanopyStomataBallBerryMod.F90 +++ b/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 diff --git a/src/ResistanceCanopyStomataJarvisMod.F90 b/src/ResistanceCanopyStomataJarvisMod.F90 index c4ba68d2..8daf2a18 100644 --- a/src/ResistanceCanopyStomataJarvisMod.F90 +++ b/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 diff --git a/src/SnowLayerCombineMod.F90 b/src/SnowLayerCombineMod.F90 index cf40b1f9..bd1950d2 100644 --- a/src/SnowLayerCombineMod.F90 +++ b/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) diff --git a/utility/Machine.F90 b/utility/Machine.F90 index 30dc26a8..f5d8f2b2 100644 --- a/utility/Machine.F90 +++ b/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