Skip to content

Add some biophysical variables by land-use type #1407

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Add some biophysical variables by land-use type #1407

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
baadc32
added biophys vars by land use: LH, SH, LW, SW, TV, TSA
ckoven May 9, 2025
093b3b2
merged main and resolved conflicts
ckoven May 12, 2025
125c01b
moved FATES_GPP_LU hist var to the hifrq landuse call
ckoven May 9, 2025
deed888
updates per code review comments
ckoven Oct 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
224 changes: 212 additions & 12 deletions main/FatesHistoryInterfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ module FatesHistoryInterfaceMod
use FatesConstantsMod , only : mg_per_kg
use FatesConstantsMod , only : pi_const
use FatesConstantsMod , only : nearzero
use FatesConstantsMod , only : rsnbl_math_prec
use FatesConstantsMod , only : t_water_freeze_k_1atm
use FatesConstantsMod , only : n_term_mort_types
use FatesConstantsMod , only : i_term_mort_type_cstarv
Expand Down Expand Up @@ -138,6 +139,7 @@ module FatesHistoryInterfaceMod
use FatesSizeAgeTypeIndicesMod, only : get_cdamagesize_class_index
use FatesSizeAgeTypeIndicesMod, only : get_cdamagesizepft_class_index
use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index
use FatesInterfaceTypesMod , only : hlm_use_luh

implicit none
private ! By default everything is private
Expand Down Expand Up @@ -362,6 +364,12 @@ module FatesHistoryInterfaceMod
integer :: ih_burnedarea_si_landuse
integer :: ih_gpp_si_landuse
integer :: ih_npp_si_landuse
integer :: ih_tveg_si_landuse
integer :: ih_tsa_si_landuse
integer :: ih_sw_abs_si_landuse
integer :: ih_lw_net_si_landuse
integer :: ih_shflux_si_landuse
integer :: ih_lhflux_si_landuse

! land use by land use variables
integer :: ih_disturbance_rate_si_lulu
Expand Down Expand Up @@ -838,6 +846,7 @@ module FatesHistoryInterfaceMod
procedure :: update_history_hifrq_sitelevel
procedure :: update_history_hifrq_subsite
procedure :: update_history_hifrq_subsite_ageclass
procedure :: update_history_hifrq_landuse
procedure :: update_history_hydraulics
procedure :: update_history_nutrflux

Expand Down Expand Up @@ -5007,6 +5016,9 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,bc_out,dt_tstep)

if(hlm_hist_level_hifrq>0) then
call update_history_hifrq_sitelevel(this,nc,nsites,sites,bc_in,bc_out,dt_tstep)
if (hlm_use_luh .eq. itrue) then
call update_history_hifrq_landuse(this,nc,nsites,sites,bc_in,dt_tstep)
end if
if(hlm_hist_level_hifrq>1) then
call update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tstep)
call update_history_hifrq_subsite_ageclass(this,nsites,sites,dt_tstep)
Expand Down Expand Up @@ -5240,6 +5252,161 @@ end subroutine update_history_hifrq_sitelevel

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

subroutine update_history_hifrq_landuse(this,nc,nsites,sites,bc_in,dt_tstep)

!
! Arguments
class(fates_history_interface_type) :: this
integer , intent(in) :: nc ! clump index
integer , intent(in) :: nsites
type(ed_site_type) , intent(inout), target :: sites(nsites)
type(bc_in_type) , intent(in) :: bc_in(nsites)
real(r8) , intent(in) :: dt_tstep

! Locals
integer :: s ! The local site index
integer :: io_si ! The site index of the IO array

real(r8) :: landuse_statevector(n_landuse_cats)
real(r8) :: canopy_area_bylanduse(n_landuse_cats)
integer :: i_lu
logical :: foundbaregroundpatch
Copy link
Contributor

Choose a reason for hiding this comment

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

foundbaregroundpatch appears to be set but unused.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It is used within the if clause on line https://github.com/NGEET/fates/pull/1407/files#diff-03feca0f1df8d3791c0d8f7b6055a0368d5241716d19d6a97c3d227fd6c5ad23R5337


type(fates_patch_type),pointer :: cpatch
type(fates_cohort_type),pointer :: ccohort
real(r8) :: dt_tstep_inv ! Time step in frequency units (/s)
real(r8) :: vegarea_per_patcharea ! temporary weighting variable (unitless)

associate( hio_tveg_si_landuse => this%hvars(ih_tveg_si_landuse)%r82d,&
hio_gpp_si_landuse => this%hvars(ih_gpp_si_landuse)%r82d, &
hio_tsa_si_landuse => this%hvars(ih_tsa_si_landuse)%r82d,&
hio_sw_abs_si_landuse => this%hvars(ih_sw_abs_si_landuse)%r82d,&
hio_lw_net_si_landuse => this%hvars(ih_lw_net_si_landuse)%r82d,&
hio_shflux_si_landuse => this%hvars(ih_shflux_si_landuse)%r82d,&
hio_lhflux_si_landuse => this%hvars(ih_lhflux_si_landuse)%r82d)

do_sites: do s = 1,nsites

io_si = sites(s)%h_gid

dt_tstep_inv = 1.0_r8/dt_tstep

! biophysical properties that are indexed by land use
landuse_statevector(:) = sites(s)%get_current_landuse_statevector() * AREA

! get the total canopy area for each land use type
canopy_area_bylanduse(:) = 0._r8
cpatch => sites(s)%oldest_patch
do while(associated(cpatch))
if ( cpatch%land_use_label .ne. nocomp_bareground_land ) then
canopy_area_bylanduse(cpatch%land_use_label) = canopy_area_bylanduse(cpatch%land_use_label) + &
cpatch%total_canopy_area
endif
cpatch => cpatch%younger
end do

cpatch => sites(s)%oldest_patch
do while(associated(cpatch))
if (cpatch%total_canopy_area .gt. rsnbl_math_prec) then
! for TVEG, since it is only defined on vegetated area of vegetated patches, normalize by the total vegetated area
hio_tveg_si_landuse(io_si,cpatch%land_use_label) = hio_tveg_si_landuse(io_si,cpatch%land_use_label) + &
bc_in(s)%t_veg_pa(cpatch%patchno) * cpatch%total_canopy_area/canopy_area_bylanduse(cpatch%land_use_label)

! for the rest of these, first weight by the vegetated area of each patch over the total patch area for each land use type
vegarea_per_patcharea = cpatch%total_canopy_area/landuse_statevector(cpatch%land_use_label)

hio_tsa_si_landuse(io_si,cpatch%land_use_label) = hio_tsa_si_landuse(io_si,cpatch%land_use_label) + &
bc_in(s)%t2m_pa(cpatch%patchno) * vegarea_per_patcharea

hio_sw_abs_si_landuse(io_si,cpatch%land_use_label) = hio_sw_abs_si_landuse(io_si,cpatch%land_use_label) + &
bc_in(s)%swabs_pa(cpatch%patchno) * vegarea_per_patcharea

hio_lw_net_si_landuse(io_si,cpatch%land_use_label) = hio_lw_net_si_landuse(io_si,cpatch%land_use_label) + &
bc_in(s)%netlw_pa(cpatch%patchno) * vegarea_per_patcharea

hio_shflux_si_landuse(io_si,cpatch%land_use_label) = hio_shflux_si_landuse(io_si,cpatch%land_use_label) + &
bc_in(s)%shflux_pa(cpatch%patchno) * vegarea_per_patcharea

hio_lhflux_si_landuse(io_si,cpatch%land_use_label) = hio_lhflux_si_landuse(io_si,cpatch%land_use_label) + &
bc_in(s)%lhflux_pa(cpatch%patchno) * vegarea_per_patcharea
endif
cpatch => cpatch%younger
end do

! for all the land-use indexed variables, except for TVEG, also add in the component for the unvegetated area of each land use
foundbaregroundpatch = .false.
cpatch => sites(s)%oldest_patch
do while(associated(cpatch))
if (cpatch%land_use_label .eq. nocomp_bareground_land .and. .not. foundbaregroundpatch) then
foundbaregroundpatch = .true.
do i_lu = 1, n_landuse_cats
if ( landuse_statevector(i_lu) .gt. rsnbl_math_prec ) then
Copy link
Contributor

Choose a reason for hiding this comment

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

To reduce indentation and thus readability, this if ( landuse_statevector(i_lu) .gt. rsnbl_math_prec ) then ... could be replaced by cycling on its inverse—if ( landuse_statevector(i_lu) .le. rsnbl_math_prec ) then cycle.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That's interesting -- I guess for me I prefer the extra indentation as it is easier to see the code block for which the if statement applies. for some reason I find cycle statements harder to parse visually.

hio_tsa_si_landuse(io_si,i_lu) = hio_tsa_si_landuse(io_si,i_lu) + &
bc_in(s)%t2m_pa(cpatch%patchno) * &
(landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu)

hio_sw_abs_si_landuse(io_si,i_lu) = hio_sw_abs_si_landuse(io_si,i_lu) + &
bc_in(s)%swabs_pa(cpatch%patchno) * &
(landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu)

hio_lw_net_si_landuse(io_si,i_lu) = hio_lw_net_si_landuse(io_si,i_lu) + &
bc_in(s)%netlw_pa(cpatch%patchno) * &
(landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu)

hio_shflux_si_landuse(io_si,i_lu) = hio_shflux_si_landuse(io_si,i_lu) + &
bc_in(s)%shflux_pa(cpatch%patchno) * &
(landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu)

hio_lhflux_si_landuse(io_si,i_lu) = hio_lhflux_si_landuse(io_si,i_lu) + &
bc_in(s)%lhflux_pa(cpatch%patchno) * &
(landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu)
end if
end do
end if
cpatch => cpatch%younger
end do

! instead of leaving the values for unoccupied areas as zero, set as missing values
do i_lu = 1, n_landuse_cats

! if a given land use type is not present, set the value as missing
if ( landuse_statevector(i_lu) .le. rsnbl_math_prec ) then
hio_tsa_si_landuse(io_si,i_lu) = hlm_hio_ignore_val
hio_sw_abs_si_landuse(io_si,i_lu) = hlm_hio_ignore_val
hio_lw_net_si_landuse(io_si,i_lu) = hlm_hio_ignore_val
hio_shflux_si_landuse(io_si,i_lu) = hlm_hio_ignore_val
hio_lhflux_si_landuse(io_si,i_lu) = hlm_hio_ignore_val
end if

! for tveg, ignore if there is no vegetation present on any patches of a given land use type
if ( canopy_area_bylanduse(i_lu) .le. rsnbl_math_prec ) then
hio_tveg_si_landuse(io_si,i_lu) = hlm_hio_ignore_val
end if

end do

! for GPP by land use, we need to loop over both patches and cohorts
cpatch => sites(s)%oldest_patch
do while(associated(cpatch))
Copy link
Contributor

Choose a reason for hiding this comment

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

Why is this do loop separate from the others? Is it because GPP is the only one with the ccohort%n * dt_tstep_inv normalization factor? A comment would be useful to explain this.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I made it separate because it is the only loop in this subroutine that loops over both patches and cohorts. adding comment to that effect.

ccohort => cpatch%shortest
do while(associated(ccohort))
if ( (.not. ccohort%isnew) .and. (cpatch%land_use_label .gt. nocomp_bareground_land) ) then
hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) &
+ ccohort%gpp_tstep * ccohort%n * dt_tstep_inv
end if
ccohort => ccohort%taller
end do
cpatch => cpatch%younger
end do

end do do_sites

end associate
return
end subroutine update_history_hifrq_landuse

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

subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tstep)

! ---------------------------------------------------------------------------------
Expand Down Expand Up @@ -5298,7 +5465,6 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tst
hio_froot_mr_understory_si_scls => this%hvars(ih_froot_mr_understory_si_scls)%r82d, &
hio_resp_g_understory_si_scls => this%hvars(ih_resp_g_understory_si_scls)%r82d, &
hio_resp_m_understory_si_scls => this%hvars(ih_resp_m_understory_si_scls)%r82d, &
hio_gpp_si_landuse => this%hvars(ih_gpp_si_landuse)%r82d, &
hio_parsun_z_si_cnlf => this%hvars(ih_parsun_z_si_cnlf)%r82d, &
hio_parsha_z_si_cnlf => this%hvars(ih_parsha_z_si_cnlf)%r82d, &
hio_ts_net_uptake_si_cnlf => this%hvars(ih_ts_net_uptake_si_cnlf)%r82d, &
Expand Down Expand Up @@ -5387,11 +5553,6 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tst
hio_ar_frootm_si_scpf(io_si,scpf) = hio_ar_frootm_si_scpf(io_si,scpf) + &
ccohort%froot_mr * n_perm2

if (cpatch%land_use_label .gt. nocomp_bareground_land) then
hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) &
+ ccohort%gpp_tstep * ccohort%n * dt_tstep_inv
end if

! accumulate fluxes on canopy- and understory- separated fluxes
if (ccohort%canopy_layer .eq. 1) then

Expand Down Expand Up @@ -8743,6 +8904,51 @@ subroutine define_history_vars(this, initialize_variables)
avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_tveg_si )

if (hlm_use_luh .eq. itrue) then
! biophysics variables that are indexed by land use type
call this%set_history_var(vname='FATES_TVEG_LU', units='degrees Kelvin', &
long='fates instantaneous mean vegetation temperature by land use type', &
use_default='active', &
avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_tveg_si_landuse )

call this%set_history_var(vname='FATES_TSA_LU', units='degrees Kelvin', &
long='fates instantaneous mean near-surface (2m) air temperature by land use type', &
use_default='active', &
avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_tsa_si_landuse )

call this%set_history_var(vname='FATES_SWABS_LU', units='W m-2', &
long='fates absorbed shortwave radiation by land use type', &
use_default='active', &
avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_sw_abs_si_landuse )

call this%set_history_var(vname='FATES_NETLW_LU', units='W m-2', &
long='fates net longwave flux by land use type', &
use_default='active', &
avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_lw_net_si_landuse )

call this%set_history_var(vname='FATES_SHFLUX_LU', units='W m-2', &
long='fates sensible heat flux by land use type', &
use_default='active', &
avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_shflux_si_landuse )

call this%set_history_var(vname='FATES_LHFLUX_LU', units='W m-2', &
long='fates latent heat flux by land use type', &
use_default='active', &
avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
ivar=ivar, initialize=initialize_variables, index = ih_lhflux_si_landuse )

call this%set_history_var(vname='FATES_GPP_LU', units='kg m-2 s-1', &
long='gross primary productivity by land use type in kg carbon per m2 per second', &
use_default='inactive', avgflag='A', vtype=site_landuse_r8, &
hlms='CLM:ALM', upfreq=group_hifr_simple, ivar=ivar, initialize=initialize_variables, &
index = ih_gpp_si_landuse)
endif

call this%set_history_var(vname='FATES_VIS_RAD_ERROR', units='-', &
long='mean two-stream solver error for VIS', use_default='active', &
avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, &
Expand Down Expand Up @@ -8902,12 +9108,6 @@ subroutine define_history_vars(this, initialize_variables)
hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, &
index = ih_gpp_si_age)

call this%set_history_var(vname='FATES_GPP_LU', units='kg m-2 s-1', &
long='gross primary productivity by land use type in kg carbon per m2 per second', &
use_default='inactive', avgflag='A', vtype=site_landuse_r8, &
hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, &
index = ih_gpp_si_landuse)

call this%set_history_var(vname='FATES_RDARK_USTORY_SZ', &
units = 'kg m-2 s-1', &
long='dark respiration for understory plants in kg carbon per m2 per second by size', &
Expand Down
12 changes: 12 additions & 0 deletions main/FatesInterfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -321,6 +321,12 @@ subroutine zero_bcs(fates,s)
fates%bc_in(s)%hksat_sisl(:) = 0.0_r8
end if

! Diagnostic quantities for outputting FATES patch-resolved
fates%bc_in(s)%lhflux_pa(:) = 0._r8
fates%bc_in(s)%shflux_pa(:) = 0._r8
fates%bc_in(s)%swabs_pa(:) = 0._r8
fates%bc_in(s)%netlw_pa(:) = 0._r8
fates%bc_in(s)%t2m_pa(:) = 0._r8

! Output boundaries
fates%bc_out(s)%active_suction_sl(:) = .false.
Expand Down Expand Up @@ -552,6 +558,12 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats,
end if

! Land use
! Diagnostic quantities for outputting FATES patch-resolved
allocate(bc_in%lhflux_pa(0:maxpatch_total))
allocate(bc_in%shflux_pa(0:maxpatch_total))
allocate(bc_in%swabs_pa(0:maxpatch_total))
allocate(bc_in%netlw_pa(0:maxpatch_total))
allocate(bc_in%t2m_pa(0:maxpatch_total))

! harvest flag denote data from hlm,
! while the logging flag signifies only that logging is occurring (which could just be FATES logging)
Expand Down
7 changes: 7 additions & 0 deletions main/FatesInterfaceTypesMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -593,6 +593,13 @@ module FatesInterfaceTypesMod
real(r8),allocatable :: h2o_liq_sisl(:) ! Liquid water mass in each layer (kg/m2)
real(r8) :: smpmin_si ! restriction for min of soil potential (mm)

! Diagnostic quantities for outputting FATES patch-resolved
real(r8),allocatable :: lhflux_pa(:) ! latent heat flux
real(r8),allocatable :: shflux_pa(:) ! sensible heat flux
real(r8),allocatable :: swabs_pa(:) ! shortwave absorbed radiation
real(r8),allocatable :: netlw_pa(:) ! longwave net radiation
real(r8),allocatable :: t2m_pa(:) ! 2m air temeprature

! Land use
! ---------------------------------------------------------------------------------
real(r8),allocatable :: hlm_harvest_rates(:) ! annual harvest rate per cat from hlm for a site
Expand Down