Skip to content
76 changes: 54 additions & 22 deletions biogeochem/EDCanopyStructureMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ module EDCanopyStructureMod
use FatesInterfaceTypesMod , only : hlm_use_planthydro
use FatesInterfaceTypesMod , only : hlm_use_cohort_age_tracking
use FatesInterfaceTypesMod , only : hlm_use_sp
use FatesInterfaceTypesMod , only : hlm_use_interstitial_bareground
use FatesInterfaceTypesMod , only : numpft
use FatesInterfaceTypesMod, only : bc_in_type
use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
Expand Down Expand Up @@ -1328,7 +1329,7 @@ end subroutine leaf_area_profile

! ======================================================================================

subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out,bc_in)

! ----------------------------------------------------------------------------------
! The purpose of this routine is to package output boundary conditions related
Expand All @@ -1345,22 +1346,25 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
type(ed_site_type), intent(inout), target :: sites(nsites)
integer, intent(in) :: fcolumn(nsites)
type(bc_out_type), intent(inout) :: bc_out(nsites)
type(bc_in_type), intent(in) :: bc_in(nsites)

! Locals
type (fates_cohort_type) , pointer :: currentCohort
integer :: s, ifp, c, p
type (fates_patch_type) , pointer :: currentPatch
real(r8) :: bare_frac_area
real(r8) :: local_patch_fraction
real(r8) :: total_patch_area
real(r8) :: total_canopy_area
real(r8) :: total_patch_leaf_stem_area
real(r8) :: weight ! Weighting for cohort variables in patch

real(r8) :: weight ! Weighting for cohort variables in patch
real(r8) :: weighting_area ! Area to normalize against depending on insterstitial bareground handling
real(r8) :: bareground_fraction ! fraction of the patch that is bareground

do s = 1,nsites

total_patch_area = 0._r8
total_canopy_area = 0._r8
bc_out(s)%canopy_fraction_pa(:) = 0._r8
local_patch_fraction = 0._r8
bc_out(s)%patch_fraction(:) = 0._r8
bc_out(s)%dleaf_pa(:) = 0._r8
bc_out(s)%z0m_pa(:) = 0._r8
bc_out(s)%displa_pa(:) = 0._r8
Expand Down Expand Up @@ -1389,24 +1393,46 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)

bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8))

! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic
! Roughness length and displacement height used by the HLM
! use total LAI + SAI to weight the leaft characteristic dimension
! Define the characteristic roughness length and displacement height used by the HLM
! Nominally, FATES uses canopy-only crown area weighting for all cohorts in the patch
! when the interstitial bareground is being collapsed into the HLM bareground patch.
! When we account for the bareground as part of the patch, the full area is accounted
! for in the calculation.
! If there is no canopy and we are considering the interstitial bareground we do
! something else.
! Use total LAI + SAI to weight the leaft characteristic dimension
! Avoid this if running in satellite phenology mode
! ----------------------------------------------------------------------------

! Set the are to be used for weighting against the cohort area
! If we are not using the interstitial bareground paradigm, simply
! use the patch area. Otherwise use the total canopy area.
weighting_area = currentPatch%total_canopy_area
if (hlm_use_interstitial_bareground .eq. ifalse) then
weighting_area = currentPatch%area
end if

if (currentPatch%total_canopy_area > nearzero) then
currentCohort => currentPatch%shortest
do while(associated(currentCohort))
if (currentCohort%canopy_layer .eq. 1) then
weight = min(1.0_r8,currentCohort%c_area/currentPatch%total_canopy_area)
weight = min(1.0_r8,currentCohort%c_area/weighting_area)
bc_out(s)%z0m_pa(ifp) = bc_out(s)%z0m_pa(ifp) + &
EDPftvarcon_inst%z0mr(currentCohort%pft) * currentCohort%height * weight
Comment thread
glemieux marked this conversation as resolved.
bc_out(s)%displa_pa(ifp) = bc_out(s)%displa_pa(ifp) + &
EDPftvarcon_inst%displar(currentCohort%pft) * currentCohort%height * weight
endif
currentCohort => currentCohort%taller
end do

! If we are not using the interstitial bareground, we need to account for the bareground area
! fraction in calculating the roughness length. Here we mimic ZengWang 2007. This should
! eventually be replaced with a more robust approach that more accurately accounts for the
! the bareground fraction of the patch.
if (hlm_use_interstitial_bareground .eq. ifalse) then
bareground_fraction = (currentPatch%area - currentPatch%total_canopy_area) / currentPatch%area
bc_out(s)%z0m_pa(ifp) = exp(log(bc_out(s)%z0m_pa(ifp)) + bareground_fraction * log(bc_in(s)%z0mg))
end if

! for lai, scale to total LAI + SAI in patch. first add up all the LAI and SAI in the patch
total_patch_leaf_stem_area = 0._r8
Expand All @@ -1433,9 +1459,9 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
endif
else
! if no canopy, then use dummy values (first PFT) of aerodynamic properties
bc_out(s)%z0m_pa(ifp) = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp)
bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp)
bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1)
bc_out(s)%z0m_pa(ifp) = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp)
bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp)
bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1)
Comment thread
glemieux marked this conversation as resolved.
endif
! -----------------------------------------------------------------------------

Expand All @@ -1445,19 +1471,25 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants
! currentPatch%area/AREA is the fraction of the soil covered by this patch.

! If we are accounting for the interstitial bareground within the HLM bareground patch
! we use the projected total canopy area fraction as the output boundary condition.
! Otherwise we simply use the fraction of the area that the patch contributes to the given site.
if (hlm_use_interstitial_bareground == itrue) then
local_patch_fraction = min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)
else
local_patch_fraction = 1.0_r8
end if

if(currentPatch%area.gt.0.0_r8)then
bc_out(s)%canopy_fraction_pa(ifp) = &
min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA)
bc_out(s)%patch_fraction(ifp) = local_patch_fraction*(currentPatch%area/AREA)
else
bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8
bc_out(s)%patch_fraction(ifp) = 0.0_r8
endif

bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * &
bare_frac_area = (1.0_r8 - local_patch_fraction) * &
(currentPatch%area/AREA)

total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area

@glemieux glemieux Jun 14, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that this is dead and duplicate code being removed. total_canopy_area wasn't being used downstream.


total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp)
total_patch_area = total_patch_area + bc_out(s)%patch_fraction(ifp) + bare_frac_area

bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label

Expand All @@ -1484,7 +1516,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8
end if

else ! nocomp or SP, and currentPatch%nocomp_pft_label .eq. 0
else ! nocomp or SP, (i.e.currentPatch%nocomp_pft_label .eq. 0

total_patch_area = total_patch_area + currentPatch%area/AREA

Expand All @@ -1511,7 +1543,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
do while(associated(currentPatch))
ifp = currentPatch%patchno
if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! for vegetated patches only
bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area
bc_out(s)%patch_fraction(ifp) = bc_out(s)%patch_fraction(ifp)/total_patch_area
endif ! veg patch
currentPatch => currentPatch%younger
end do
Expand Down
11 changes: 9 additions & 2 deletions main/FatesInterfaceMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -411,7 +411,7 @@ subroutine zero_bcs(fates,s)
fates%bc_out(s)%dleaf_pa(:) = 0.0_r8
fates%bc_out(s)%nocomp_pft_label_pa(:) = 0

fates%bc_out(s)%canopy_fraction_pa(:) = 0.0_r8
fates%bc_out(s)%patch_fraction(:) = 0.0_r8
fates%bc_out(s)%frac_veg_nosno_alb_pa(:) = 0.0_r8

if (hlm_use_planthydro.eq.itrue) then
Expand Down Expand Up @@ -744,7 +744,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
allocate(bc_out%displa_pa(maxpatch_total))
allocate(bc_out%z0m_pa(maxpatch_total))

allocate(bc_out%canopy_fraction_pa(maxpatch_total))
allocate(bc_out%patch_fraction(maxpatch_total))
allocate(bc_out%frac_veg_nosno_alb_pa(maxpatch_total))

allocate(bc_out%nocomp_pft_label_pa(maxpatch_total))
Expand Down Expand Up @@ -1574,6 +1574,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
hlm_use_fixed_biogeog = unset_int
hlm_use_nocomp = unset_int
hlm_use_sp = unset_int
hlm_use_interstitial_bareground = unset_int
hlm_use_inventory_init = unset_int
hlm_use_dbh_init = unset_int
hlm_inventory_ctrl_file = 'unset'
Expand Down Expand Up @@ -2050,6 +2051,12 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
write(fates_log(),*) 'Transfering hlm_use_sp= ',ival,' to FATES'
end if

case('use_fates_interstitial_bareground')
hlm_use_interstitial_bareground = ival
if (fates_global_verbose()) then
write(fates_log(),*) 'Transfering hlm_use_interstitial_bareground= ',ival,' to FATES'
end if

case('use_planthydro')
hlm_use_planthydro = ival
if (fates_global_verbose()) then
Expand Down
16 changes: 12 additions & 4 deletions main/FatesInterfaceTypesMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -225,6 +225,9 @@ module FatesInterfaceTypesMod
integer, public :: hlm_use_sp ! Flag to use FATES satellite phenology (LAI) mode
! 1 = TRUE, 0 = FALSE

integer, public :: hlm_use_interstitial_bareground ! Flag to use count the FATES interstitial bareground area
! with the HLM bareground patch
! 1 = TRUE, 0 = FALSE

! Flag specifying what types of history fields to allocate and prepare
! The "_dynam" refers to history fields that can be updated on the dynamics (daily) step
Expand Down Expand Up @@ -503,6 +506,9 @@ module FatesInterfaceTypesMod

! soil temperature (Kelvin)
real(r8), allocatable :: t_soisno_sl(:)

! surface roughness length, site bareground (m)
real(r8) :: z0mg

! Canopy Radiation Boundaries
! ---------------------------------------------------------------------------------
Expand Down Expand Up @@ -766,10 +772,12 @@ module FatesInterfaceTypesMod
real(r8), allocatable :: displa_pa(:) ! displacement height [m]
real(r8), allocatable :: dleaf_pa(:) ! leaf characteristic dimension/width/diameter [m]

real(r8), allocatable :: canopy_fraction_pa(:) ! Area fraction of each patch in the site
! Use most likely for weighting
! This is currently the projected canopy
! area of each patch [0-1]
real(r8), allocatable :: patch_fraction(:) ! Area fraction of each patch in the site
! Use most likely for weighting
! This is either the projected canopy
! area of each patch [0-1] or the actual
! are fraction that each patch contributes
! to the site

real(r8), allocatable :: frac_veg_nosno_alb_pa(:) ! This is not really a fraction
! this is actually binary based on if any
Expand Down