From a7dd399972fca4d477c43654d384287ad62114a2 Mon Sep 17 00:00:00 2001 From: cenlinhe Date: Tue, 19 Sep 2023 16:51:14 -0600 Subject: [PATCH] bug fix for VegFrac scaling, issue #91-92 --- src/CanopyHydrologyMod.F90 | 5 +++-- src/CanopyWaterInterceptMod.F90 | 4 ++-- src/ResistanceCanopyStomataBallBerryMod.F90 | 5 +++-- src/ResistanceCanopyStomataJarvisMod.F90 | 5 +++-- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/CanopyHydrologyMod.F90 b/src/CanopyHydrologyMod.F90 index 264da43..24fab3b 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] CanopyLiqWater => noahmp%water%state%CanopyLiqWater ,& ! inout, intercepted canopy liquid water [mm] @@ -67,7 +68,7 @@ subroutine CanopyHydrology(noahmp) ! canopy liquid water ! maximum canopy intercepted water - CanopyLiqWaterMax = CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) + CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) ! canopy evaporation, transpiration, and dew if ( FlagFrozenCanopy .eqv. .false. ) then ! Barlage: change to FlagFrozenCanopy @@ -92,7 +93,7 @@ subroutine CanopyHydrology(noahmp) ! canopy ice ! maximum canopy intercepted ice - CanopyIceMax = 6.6 * (0.27 + 46.0/SnowfallDensity) * (LeafAreaIndEff + StemAreaIndEff) + CanopyIceMax = VegFrac * 6.6 * (0.27 + 46.0/SnowfallDensity) * (LeafAreaIndEff + StemAreaIndEff) ! canopy sublimation and frost SublimCanopyIce = min( CanopyIce/MainTimeStep, SublimCanopyIce ) diff --git a/src/CanopyWaterInterceptMod.F90 b/src/CanopyWaterInterceptMod.F90 index c235b07..274d0c2 100644 --- a/src/CanopyWaterInterceptMod.F90 +++ b/src/CanopyWaterInterceptMod.F90 @@ -79,7 +79,7 @@ subroutine CanopyWaterIntercept(noahmp) ! ----------------------- canopy liquid water ------------------------------ ! maximum canopy water - CanopyLiqWaterMax = CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) + CanopyLiqWaterMax = VegFrac * CanopyLiqHoldCap * (LeafAreaIndEff + StemAreaIndEff) ! average rain interception and throughfall if ( (LeafAreaIndEff+StemAreaIndEff) > 0.0 ) then @@ -102,7 +102,7 @@ subroutine CanopyWaterIntercept(noahmp) ! ----------------------- canopy ice ------------------------------ ! maximum canopy ice - CanopyIceMax = 6.6 * (0.27 + 46.0/SnowfallDensity) * (LeafAreaIndEff + StemAreaIndEff) + CanopyIceMax = VegFrac * 6.6 * (0.27 + 46.0/SnowfallDensity) * (LeafAreaIndEff + StemAreaIndEff) ! average snow interception and throughfall if ( (LeafAreaIndEff+StemAreaIndEff) > 0.0 ) then diff --git a/src/ResistanceCanopyStomataBallBerryMod.F90 b/src/ResistanceCanopyStomataBallBerryMod.F90 index 72b84b5..d479bec 100644 --- a/src/ResistanceCanopyStomataBallBerryMod.F90 +++ b/src/ResistanceCanopyStomataBallBerryMod.F90 @@ -82,6 +82,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] @@ -99,8 +100,8 @@ subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade) CF = PressureAirRefHeight / (8.314 * TemperatureAirRefHeight) * 1.0e06 ! unit conversion factor ResistanceStomataTmp = 1.0 / ConductanceLeafMin * CF PhotosynLeafTmp = 0.0 - if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit ! Sunlit case - if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade ! 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 ! only compute when there is radiation absorption if ( RadPhotoActAbsTmp > 0.0 ) then diff --git a/src/ResistanceCanopyStomataJarvisMod.F90 b/src/ResistanceCanopyStomataJarvisMod.F90 index 86575de..39388bd 100644 --- a/src/ResistanceCanopyStomataJarvisMod.F90 +++ b/src/ResistanceCanopyStomataJarvisMod.F90 @@ -53,6 +53,7 @@ subroutine ResistanceCanopyStomataJarvis(noahmp, IndexShade) VaporPresDeficitFac => noahmp%energy%param%VaporPresDeficitFac ,& ! in, Parameter used in vapor pressure deficit function TemperatureCanopy => noahmp%energy%state%TemperatureCanopy ,& ! in, vegetation temperature [K] PressureVaporCanAir => noahmp%energy%state%PressureVaporCanAir ,& ! in, canopy air vapor pressure [Pa] + 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] @@ -67,8 +68,8 @@ subroutine ResistanceCanopyStomataJarvis(noahmp, IndexShade) ResistanceTemp = 0.0 ResistanceVapDef = 0.0 ResistanceStomataTmp = 0.0 - if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit ! Sunlit case - if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade ! 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