bug fix for VegFrac scaling, issue #91-92

NCAR · Sep 19, 2023 · a7dd399 · a7dd399
1 parent e0d2064
commit a7dd399
Show file tree

Hide file tree

Showing 4 changed files with 11 additions and 8 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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