From e9b85a03dfcea81f32a9f5ce1c6ceac3e2f63456 Mon Sep 17 00:00:00 2001 From: tslin2 Date: Wed, 17 Jan 2024 17:13:47 -0700 Subject: [PATCH 01/13] merge --- src/CanopyHydrologyMod.F90 | 5 +++++ src/ResistanceCanopyStomataBallBerryMod.F90 | 5 +++++ src/ResistanceCanopyStomataJarvisMod.F90 | 5 +++++ 3 files changed, 15 insertions(+) diff --git a/src/CanopyHydrologyMod.F90 b/src/CanopyHydrologyMod.F90 index 24fab3b4..f41c2919 100644 --- a/src/CanopyHydrologyMod.F90 +++ b/src/CanopyHydrologyMod.F90 @@ -34,6 +34,7 @@ subroutine CanopyHydrology(noahmp) 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 CanopyLiqWater => noahmp%water%state%CanopyLiqWater ,& ! inout, intercepted canopy liquid water [mm] CanopyIce => noahmp%water%state%CanopyIce ,& ! inout, intercepted canopy ice [mm] TemperatureCanopy => noahmp%energy%state%TemperatureCanopy ,& ! inout, vegetation temperature [K] @@ -68,7 +69,11 @@ subroutine CanopyHydrology(noahmp) ! canopy liquid water ! maximum canopy intercepted water +<<<<<<< Updated upstream CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) +======= + CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) +>>>>>>> Stashed changes ! canopy evaporation, transpiration, and dew if ( FlagFrozenCanopy .eqv. .false. ) then ! Barlage: change to FlagFrozenCanopy diff --git a/src/ResistanceCanopyStomataBallBerryMod.F90 b/src/ResistanceCanopyStomataBallBerryMod.F90 index d479bec0..c453b676 100644 --- a/src/ResistanceCanopyStomataBallBerryMod.F90 +++ b/src/ResistanceCanopyStomataBallBerryMod.F90 @@ -79,6 +79,7 @@ subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade) TemperatureCanopy => noahmp%energy%state%TemperatureCanopy ,& ! in, vegetation temperature [K] VapPresSatCanopy => noahmp%energy%state%VapPresSatCanopy ,& ! in, canopy saturation vapor pressure at TV [Pa] PressureVaporCanAir => noahmp%energy%state%PressureVaporCanAir ,& ! in, canopy air vapor pressure [Pa] + VegFrac => noahmp%energy%state%VegFrac ,& ! in, greeness vegetation fraction 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] @@ -101,7 +102,11 @@ subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade) ResistanceStomataTmp = 1.0 / ConductanceLeafMin * CF PhotosynLeafTmp = 0.0 if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit / max(VegFrac,1.0e-6) ! Sunlit case +<<<<<<< Updated upstream if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade / max(VegFrac,1.0e-6) ! Shaded case +======= + if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade / max(VegFrac,1.0e-6) ! Shaded case +>>>>>>> Stashed changes ! only compute when there is radiation absorption if ( RadPhotoActAbsTmp > 0.0 ) then diff --git a/src/ResistanceCanopyStomataJarvisMod.F90 b/src/ResistanceCanopyStomataJarvisMod.F90 index 39388bd1..0e60f8b1 100644 --- a/src/ResistanceCanopyStomataJarvisMod.F90 +++ b/src/ResistanceCanopyStomataJarvisMod.F90 @@ -68,8 +68,13 @@ subroutine ResistanceCanopyStomataJarvis(noahmp, IndexShade) ResistanceTemp = 0.0 ResistanceVapDef = 0.0 ResistanceStomataTmp = 0.0 +<<<<<<< Updated upstream 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 +>>>>>>> Stashed changes ! compute MixingRatioTmp and MixingRatioSat SpecHumidityTmp = 0.622 * PressureVaporCanAir / (PressureAirRefHeight - 0.378*PressureVaporCanAir) ! specific humidity From 074a66eadfc128d54e1fec71e82544a32b4b8fc9 Mon Sep 17 00:00:00 2001 From: tslin2 Date: Wed, 17 Jan 2024 17:15:37 -0700 Subject: [PATCH 02/13] negative theta_1500 in PedoTransferSR2006Mod.F90 --- drivers/hrldas/PedoTransferSR2006Mod.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/drivers/hrldas/PedoTransferSR2006Mod.F90 b/drivers/hrldas/PedoTransferSR2006Mod.F90 index 02090e82..848c9074 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**-5.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 From 53767f48f322594b32e82e2b5ff173d5e41d21dd Mon Sep 17 00:00:00 2001 From: tslin2 Date: Wed, 17 Jan 2024 17:19:43 -0700 Subject: [PATCH 03/13] update --- src/CanopyHydrologyMod.F90 | 1 - src/ResistanceCanopyStomataBallBerryMod.F90 | 5 ----- src/ResistanceCanopyStomataJarvisMod.F90 | 6 +----- 3 files changed, 1 insertion(+), 11 deletions(-) diff --git a/src/CanopyHydrologyMod.F90 b/src/CanopyHydrologyMod.F90 index f41c2919..31b50385 100644 --- a/src/CanopyHydrologyMod.F90 +++ b/src/CanopyHydrologyMod.F90 @@ -34,7 +34,6 @@ subroutine CanopyHydrology(noahmp) 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 CanopyLiqWater => noahmp%water%state%CanopyLiqWater ,& ! inout, intercepted canopy liquid water [mm] CanopyIce => noahmp%water%state%CanopyIce ,& ! inout, intercepted canopy ice [mm] TemperatureCanopy => noahmp%energy%state%TemperatureCanopy ,& ! inout, vegetation temperature [K] diff --git a/src/ResistanceCanopyStomataBallBerryMod.F90 b/src/ResistanceCanopyStomataBallBerryMod.F90 index c453b676..b455adcd 100644 --- a/src/ResistanceCanopyStomataBallBerryMod.F90 +++ b/src/ResistanceCanopyStomataBallBerryMod.F90 @@ -83,7 +83,6 @@ 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] @@ -102,11 +101,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 -<<<<<<< Updated upstream if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade / max(VegFrac,1.0e-6) ! Shaded case -======= - if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade / max(VegFrac,1.0e-6) ! Shaded case ->>>>>>> Stashed changes ! only compute when there is radiation absorption if ( RadPhotoActAbsTmp > 0.0 ) then diff --git a/src/ResistanceCanopyStomataJarvisMod.F90 b/src/ResistanceCanopyStomataJarvisMod.F90 index 0e60f8b1..845a9b1d 100644 --- a/src/ResistanceCanopyStomataJarvisMod.F90 +++ b/src/ResistanceCanopyStomataJarvisMod.F90 @@ -68,13 +68,9 @@ subroutine ResistanceCanopyStomataJarvis(noahmp, IndexShade) ResistanceTemp = 0.0 ResistanceVapDef = 0.0 ResistanceStomataTmp = 0.0 -<<<<<<< Updated upstream + 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 ->>>>>>> Stashed changes ! compute MixingRatioTmp and MixingRatioSat SpecHumidityTmp = 0.622 * PressureVaporCanAir / (PressureAirRefHeight - 0.378*PressureVaporCanAir) ! specific humidity From 75cd9a497ad2064f33fba687b0f7a32acb440b66 Mon Sep 17 00:00:00 2001 From: tslin2 Date: Wed, 17 Jan 2024 17:21:30 -0700 Subject: [PATCH 04/13] update theta --- src/CanopyHydrologyMod.F90 | 4 ---- src/ResistanceCanopyStomataBallBerryMod.F90 | 2 +- src/ResistanceCanopyStomataJarvisMod.F90 | 1 - 3 files changed, 1 insertion(+), 6 deletions(-) diff --git a/src/CanopyHydrologyMod.F90 b/src/CanopyHydrologyMod.F90 index 31b50385..24fab3b4 100644 --- a/src/CanopyHydrologyMod.F90 +++ b/src/CanopyHydrologyMod.F90 @@ -68,11 +68,7 @@ subroutine CanopyHydrology(noahmp) ! canopy liquid water ! maximum canopy intercepted water -<<<<<<< Updated upstream CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) -======= - CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) ->>>>>>> Stashed changes ! canopy evaporation, transpiration, and dew if ( FlagFrozenCanopy .eqv. .false. ) then ! Barlage: change to FlagFrozenCanopy diff --git a/src/ResistanceCanopyStomataBallBerryMod.F90 b/src/ResistanceCanopyStomataBallBerryMod.F90 index b455adcd..d479bec0 100644 --- a/src/ResistanceCanopyStomataBallBerryMod.F90 +++ b/src/ResistanceCanopyStomataBallBerryMod.F90 @@ -79,10 +79,10 @@ subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade) TemperatureCanopy => noahmp%energy%state%TemperatureCanopy ,& ! in, vegetation temperature [K] VapPresSatCanopy => noahmp%energy%state%VapPresSatCanopy ,& ! in, canopy saturation vapor pressure at TV [Pa] PressureVaporCanAir => noahmp%energy%state%PressureVaporCanAir ,& ! in, canopy air vapor pressure [Pa] - VegFrac => noahmp%energy%state%VegFrac ,& ! in, greeness vegetation fraction 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] diff --git a/src/ResistanceCanopyStomataJarvisMod.F90 b/src/ResistanceCanopyStomataJarvisMod.F90 index 845a9b1d..39388bd1 100644 --- a/src/ResistanceCanopyStomataJarvisMod.F90 +++ b/src/ResistanceCanopyStomataJarvisMod.F90 @@ -68,7 +68,6 @@ 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 From 2ba002a647a6468e1152f88ccaf4d5f5f814b6b2 Mon Sep 17 00:00:00 2001 From: tslin2 Date: Wed, 17 Jan 2024 17:33:46 -0700 Subject: [PATCH 05/13] update for value --- drivers/hrldas/PedoTransferSR2006Mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/drivers/hrldas/PedoTransferSR2006Mod.F90 b/drivers/hrldas/PedoTransferSR2006Mod.F90 index 848c9074..261864cf 100644 --- a/drivers/hrldas/PedoTransferSR2006Mod.F90 +++ b/drivers/hrldas/PedoTransferSR2006Mod.F90 @@ -162,7 +162,7 @@ subroutine PedoTransferSR2006(NoahmpIO, noahmp, Sand, Clay, Orgm) + sr2006_psi_e_b*psi_et & + sr2006_psi_e_c - theta_33 = max(10.0**-5.0,theta_33) ! For numerical stability + 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 From c27b8689b8714ae4235a0dca5b8ec927b011ddca Mon Sep 17 00:00:00 2001 From: Ronnie Abolafia-Rosenzweig <36647844+RAbolafiaRosenzweig@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:44:38 -0700 Subject: [PATCH 06/13] Add files via upload Updated to avoid numerical errors associated with low FrozenPrecipFrac for OptRainSnowPartition = 5 --- src/AtmosForcingMod.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index 96a7105b..112c6d81 100644 --- a/src/AtmosForcingMod.F90 +++ b/src/AtmosForcingMod.F90 @@ -165,7 +165,13 @@ subroutine ProcessAtmosForcing(noahmp) 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 structure created by Ronnie Abolafia-Rosenzweig (Feb 1, 2024) to impose FrozenPrecipFrac=0 on high-temp instances to avoid numerical issues + if ( TemperatureWetBulb >= (ConstFreezePoint+5)) then + FrozenPrecipFrac = 0 + else + 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 From 8569cc6340a9f223a6a4e374feea8be65cb26b48 Mon Sep 17 00:00:00 2001 From: Ronnie Abolafia-Rosenzweig <36647844+RAbolafiaRosenzweig@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:54:43 -0700 Subject: [PATCH 07/13] Add files via upload AtmosForcingMod.F90 updated to avoid numerical errors associated with warm wet bulb temperature for OptRainSnowPartition = 5 --- src/AtmosForcingMod.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index 112c6d81..aa3f13c5 100644 --- a/src/AtmosForcingMod.F90 +++ b/src/AtmosForcingMod.F90 @@ -166,8 +166,7 @@ subroutine ProcessAtmosForcing(noahmp) TemperatureWetBulb = TemperatureWetBulb - (VapPresSat - PressureVaporRefHeight) / PsychConst ! Wang et al., 2019 GRL Eq.2 enddo - ! If structure created by Ronnie Abolafia-Rosenzweig (Feb 1, 2024) to impose FrozenPrecipFrac=0 on high-temp instances to avoid numerical issues - if ( TemperatureWetBulb >= (ConstFreezePoint+5)) then + if ( TemperatureWetBulb >= (ConstFreezePoint+5)) then !added by RAR 02/2024 FrozenPrecipFrac = 0 else FrozenPrecipFrac = 1.0 / (1.0 + 6.99e-5 * exp(2.0*(TemperatureWetBulb+3.97))) ! Wang et al., 2019 GRL Eq. 1 From 7e302f636307c88c9b7722bf836b4e24c011b08e Mon Sep 17 00:00:00 2001 From: Ronnie Abolafia-Rosenzweig <36647844+RAbolafiaRosenzweig@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:56:57 -0700 Subject: [PATCH 08/13] Add files via upload Update to OptRainSnowPartition == 5 to avoid numerical errors in instances of warm temperatures by forcing FrozenPrecipFrac = 0 --- src/AtmosForcingMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index aa3f13c5..a55d434c 100644 --- a/src/AtmosForcingMod.F90 +++ b/src/AtmosForcingMod.F90 @@ -151,7 +151,7 @@ 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 @@ -166,7 +166,7 @@ subroutine ProcessAtmosForcing(noahmp) TemperatureWetBulb = TemperatureWetBulb - (VapPresSat - PressureVaporRefHeight) / PsychConst ! Wang et al., 2019 GRL Eq.2 enddo - if ( TemperatureWetBulb >= (ConstFreezePoint+5)) then !added by RAR 02/2024 + if ( TemperatureWetBulb >= (ConstFreezePoint+5)) then !avoid numerical errors when temperature is high FrozenPrecipFrac = 0 else FrozenPrecipFrac = 1.0 / (1.0 + 6.99e-5 * exp(2.0*(TemperatureWetBulb+3.97))) ! Wang et al., 2019 GRL Eq. 1 From 50c3657877455606be662aa9bba07d7261a294f8 Mon Sep 17 00:00:00 2001 From: Ronnie Abolafia-Rosenzweig <36647844+RAbolafiaRosenzweig@users.noreply.github.com> Date: Thu, 1 Feb 2024 14:18:53 -0700 Subject: [PATCH 09/13] Add files via upload if ( TemperatureWetBulb >= 5) then FrozenPrecipFrac = 0 to avoid crashing in instances of high temperatures for OptRainSnowPartition == 5 --- src/AtmosForcingMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index a55d434c..b87c77f3 100644 --- a/src/AtmosForcingMod.F90 +++ b/src/AtmosForcingMod.F90 @@ -166,7 +166,7 @@ subroutine ProcessAtmosForcing(noahmp) TemperatureWetBulb = TemperatureWetBulb - (VapPresSat - PressureVaporRefHeight) / PsychConst ! Wang et al., 2019 GRL Eq.2 enddo - if ( TemperatureWetBulb >= (ConstFreezePoint+5)) then !avoid numerical errors when temperature is high + if ( TemperatureWetBulb >= 5) then !avoid numerical errors when temperature is high FrozenPrecipFrac = 0 else FrozenPrecipFrac = 1.0 / (1.0 + 6.99e-5 * exp(2.0*(TemperatureWetBulb+3.97))) ! Wang et al., 2019 GRL Eq. 1 From 8138618b8b49836fe212f4ccc514ddc48e493e0f Mon Sep 17 00:00:00 2001 From: Ronnie Abolafia-Rosenzweig <36647844+RAbolafiaRosenzweig@users.noreply.github.com> Date: Fri, 2 Feb 2024 08:40:34 -0700 Subject: [PATCH 10/13] Add files via upload FrozenPrecipFrac = 0.0 constraint imposed for instances of TemperatureAirRefHeight > 10C to avoid numerical errors in warm instances for OptRainSnowPartition == 5. --- src/AtmosForcingMod.F90 | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index b87c77f3..94451a74 100644 --- a/src/AtmosForcingMod.F90 +++ b/src/AtmosForcingMod.F90 @@ -153,22 +153,22 @@ 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 - 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 - - if ( TemperatureWetBulb >= 5) then !avoid numerical errors when temperature is high - FrozenPrecipFrac = 0 + + 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 From 1ac7515db1e159a363cdd314c3898e0e56ca0987 Mon Sep 17 00:00:00 2001 From: cenlinhe Date: Wed, 17 Apr 2024 14:52:14 -0600 Subject: [PATCH 11/13] bug fix for snow to avoid numerical floating error --- src/AtmosForcingMod.F90 | 2 +- src/PhenologyMainMod.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AtmosForcingMod.F90 b/src/AtmosForcingMod.F90 index 94451a74..9209a2dd 100644 --- a/src/AtmosForcingMod.F90 +++ b/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 diff --git a/src/PhenologyMainMod.F90 b/src/PhenologyMainMod.F90 index e2319350..e455a032 100644 --- a/src/PhenologyMainMod.F90 +++ b/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 From ca2219246cebdcd9157f327c4527c9294ad6852c Mon Sep 17 00:00:00 2001 From: cenlinhe Date: Thu, 18 Apr 2024 10:47:45 -0600 Subject: [PATCH 12/13] bug fix for LECH SURFACE FUNCTIONS in Chen97 --- src/ResistanceAboveCanopyChen97Mod.F90 | 4 ++-- src/ResistanceBareGroundChen97Mod.F90 | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/ResistanceAboveCanopyChen97Mod.F90 b/src/ResistanceAboveCanopyChen97Mod.F90 index 1020acb2..a9e58e5f 100644 --- a/src/ResistanceAboveCanopyChen97Mod.F90 +++ b/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 diff --git a/src/ResistanceBareGroundChen97Mod.F90 b/src/ResistanceBareGroundChen97Mod.F90 index f3510ce0..db7c3174 100644 --- a/src/ResistanceBareGroundChen97Mod.F90 +++ b/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 From ce8636494ff9cd7030586dbf637c20a24743e807 Mon Sep 17 00:00:00 2001 From: Cenlin_He Date: Thu, 18 Apr 2024 14:54:11 -0600 Subject: [PATCH 13/13] Update README.md --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index de0a3fb9..3090b930 100644 --- a/README.md +++ b/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