Skip to content

site-level mass flux accounting part 2 #1473

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 3 commits into
base: main
Choose a base branch
from
Open

site-level mass flux accounting part 2 #1473

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0ba914c
Updated logic on restarts to assume newly calculated stocks should ma…
rgknox Sep 23, 2025
a7c739a
More consistency in site level mass flux accounting and restarting
rgknox Sep 23, 2025
2c9a15e
removed unecessary usage of fluxdiags%npp
rgknox Sep 23, 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
2 changes: 1 addition & 1 deletion main/ChecksBalancesMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,7 +302,7 @@ subroutine CheckIntegratedMassPools(site)

select case(element_list(el))
case(carbon12_element)
net_uptake = diag%npp + site_mass%net_root_uptake*area_inv
net_uptake = (site_mass%gpp_acc + site_mass%aresp_acc + site_mass%net_root_uptake)*area_inv
case(nitrogen_element)
net_uptake = site_mass%net_root_uptake*area_inv
case(phosphorus_element)
Expand Down
94 changes: 48 additions & 46 deletions main/EDMainMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out)
call ZeroBCOutCarbonFluxes(bc_out)

! Zero mass balance
call TotalBalanceCheck(currentSite, 0)
call TotalBalanceCheck(currentSite, 0, is_restarting=.false.)

! We do not allow phenology while in ST3 mode either, it is hypothetically
! possible to allow this, but we have not plugged in the litter fluxes
Expand Down Expand Up @@ -263,7 +263,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out)
currentPatch => currentPatch%younger
enddo

call TotalBalanceCheck(currentSite,1)
call TotalBalanceCheck(currentSite,1,is_restarting=.false.)

currentPatch => currentSite%oldest_patch
do while (associated(currentPatch))
Expand All @@ -286,7 +286,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out)

end if

call TotalBalanceCheck(currentSite,2)
call TotalBalanceCheck(currentSite,2,is_restarting=.false.)

!*********************************************************************************
! Patch dynamics sub-routines: fusion, new patch creation (spwaning), termination.
Expand All @@ -304,7 +304,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out)

call spawn_patches(currentSite, bc_in)

call TotalBalanceCheck(currentSite,3)
call TotalBalanceCheck(currentSite,3,is_restarting=.false.)

! fuse on the spawned patches.
call fuse_patches(currentSite, bc_in )
Expand All @@ -319,14 +319,14 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out)
end if

! SP has changes in leaf carbon but we don't expect them to be in balance.
call TotalBalanceCheck(currentSite,4)
call TotalBalanceCheck(currentSite,4,is_restarting=.false.)

! kill patches that are too small
call terminate_patches(currentSite, bc_in)
end if

! Final instantaneous mass balance check
call TotalBalanceCheck(currentSite,5)
call TotalBalanceCheck(currentSite,5,is_restarting=.false.)


end subroutine ed_ecosystem_dynamics
Expand Down Expand Up @@ -668,11 +668,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
currentSite%mass_balance(element_pos(carbon12_element))%net_root_uptake - &
currentCohort%daily_c_efflux*currentCohort%n

! Save NPP diagnostic for flux accounting [kg/m2/day]

currentSite%flux_diags%npp = currentSite%flux_diags%npp + &
currentCohort%npp_acc_hold/real( hlm_days_per_year,r8) * currentCohort%n * area_inv

! And simultaneously add the input fluxes to mass balance accounting
site_cmass%gpp_acc = site_cmass%gpp_acc + &
currentCohort%gpp_acc * currentCohort%n
Expand All @@ -682,8 +677,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
currentCohort%resp_excess_hold*currentCohort%n + &
currentCohort%resp_g_acc_hold*currentCohort%n/real( hlm_days_per_year,r8)



call currentCohort%prt%CheckMassConservation(ft,5)

! Update the leaf biophysical rates based on proportion of leaf
Expand Down Expand Up @@ -848,10 +841,7 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting )
call set_patchno(currentSite,.true.,1)
end if

! Pass site-level mass fluxes to output boundary conditions
! [kg/site/day] * [site/m2 day/sec] = [kgC/m2/s]
bc_out%gpp_site = site_cmass%gpp_acc * area_inv / sec_per_day
bc_out%ar_site = site_cmass%aresp_acc * area_inv / sec_per_day


if(hlm_use_sp.eq.ifalse .and. (.not.is_restarting))then
call canopy_spread(currentSite)
Expand All @@ -860,13 +850,13 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting )
site_cmass%aresp_acc = 0._r8
end if

call TotalBalanceCheck(currentSite,6)
call TotalBalanceCheck(currentSite,6,is_restarting=is_restarting)

if(hlm_use_sp.eq.ifalse .and. (.not.is_restarting) )then
call canopy_structure(currentSite, bc_in)
endif

call TotalBalanceCheck(currentSite,final_check_id)
call TotalBalanceCheck(currentSite,final_check_id,is_restarting=is_restarting)

! Update recruit L2FRs based on new canopy position
call SetRecruitL2FR(currentSite)
Expand Down Expand Up @@ -927,26 +917,34 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting )
! Set boundary condition to HLM for carbon loss to atm from fires and grazing
! [kgC/ha/day]*[ha/m2]*[day/s] = [kg/m2/s]

bc_out%fire_closs_to_atm_si = site_cmass%burn_flux_to_atm * ha_per_m2 * days_per_sec
bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * ha_per_m2 * days_per_sec


bc_out%fire_closs_to_atm_si = site_cmass%burn_flux_to_atm * area_inv * days_per_sec
bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * area_inv * days_per_sec
bc_out%gpp_site = site_cmass%gpp_acc * area_inv * days_per_sec
bc_out%ar_site = site_cmass%aresp_acc * area_inv * days_per_sec

end subroutine ed_update_site

!-------------------------------------------------------------------------------!

subroutine TotalBalanceCheck (currentSite, call_index )
subroutine TotalBalanceCheck (currentSite, call_index, is_restarting )

!
! !DESCRIPTION:
! This routine looks at the mass flux in and out of the FATES and compares it to
! the change in total stocks (states).
! Fluxes in are NPP. Fluxes out are decay of CWD and litter into SOM pools.
! Note: If the model is restarting, it is assumed that the mass stocks
! that were saved in the restart file are the "old" stocks, and they
! should equal the stocks that are currently in the sites. However,
! the fluxes on a restart are saved so that they can inform the HLM
! on the next day, so we have to modify our mass balance check to ignore
! fluxes on restarts.
!
! !ARGUMENTS:
type(ed_site_type) , intent(inout) :: currentSite
integer , intent(in) :: call_index
logical , intent(in) :: is_restarting

!
! !LOCAL VARIABLES:
type(site_massbal_type),pointer :: site_mass
Expand Down Expand Up @@ -987,32 +985,35 @@ subroutine TotalBalanceCheck (currentSite, call_index )

change_in_stock = 0.0_r8


! Loop through the number of elements in the system

do el = 1, num_elements
do_elem_loop: do el = 1, num_elements

site_mass => currentSite%mass_balance(el)

call SiteMassStock(currentSite,el,total_stock,biomass_stock,litter_stock,seed_stock)

change_in_stock = total_stock - site_mass%old_stock

flux_in = site_mass%seed_in + &
site_mass%net_root_uptake + &
site_mass%gpp_acc + &
site_mass%flux_generic_in + &
site_mass%patch_resize_err

flux_out = sum(site_mass%wood_product_harvest(:)) + &
sum(site_mass%wood_product_landusechange(:)) + &
site_mass%burn_flux_to_atm + &
site_mass%seed_out + &
site_mass%flux_generic_out + &
site_mass%frag_out + &
site_mass%aresp_acc + &
site_mass%herbivory_flux_out

if(is_restarting) then
flux_in = 0._r8
flux_out = 0._r8
else
flux_in = site_mass%seed_in + &
site_mass%net_root_uptake + &
site_mass%gpp_acc + &
site_mass%flux_generic_in + &
site_mass%patch_resize_err

flux_out = sum(site_mass%wood_product_harvest(:)) + &
sum(site_mass%wood_product_landusechange(:)) + &
site_mass%burn_flux_to_atm + &
site_mass%seed_out + &
site_mass%flux_generic_out + &
site_mass%frag_out + &
site_mass%aresp_acc + &
site_mass%herbivory_flux_out
end if

net_flux = flux_in - flux_out
error = abs(net_flux - change_in_stock)

Expand Down Expand Up @@ -1110,12 +1111,13 @@ subroutine TotalBalanceCheck (currentSite, call_index )

! This is the last check of the sequence, where we update our total
! error check and the final fates stock
if(call_index == final_check_id) then
site_mass%old_stock = total_stock
site_mass%err_fates = net_flux - change_in_stock
if(call_index == final_check_id .and. .not.is_restarting) then
site_mass%old_stock = total_stock
site_mass%err_fates = net_flux - change_in_stock
end if

end do
end do do_elem_loop

end if ! not SP mode
end subroutine TotalBalanceCheck

Expand Down
21 changes: 6 additions & 15 deletions main/EDTypesMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -217,18 +217,15 @@ module EDTypesMod
type, public :: site_fluxdiags_type


! This is for all diagnostics that are uniform over all elements (C,N,P)
! These are site level flux diagnostics that are not used
! in mass balance checks. We use these structures
! to inform the history output. These values are not
! zero'd when dynamics are completed. These values
! are zero'd on cold-starts, and on restarts prior to the read

! This is for all diagnostics that are uniform over all elements (C,N,P)
type(elem_diag_type), pointer :: elem(:)

! This variable is slated as to-do, but the fluxdiags type needs
! to be refactored first. Currently this type is allocated
! by chemical species (ie C, N or P). GPP is C, but not N or P (RGK 0524)
! Previous day GPP [kgC/m2/year], partitioned by size x pft
!real(r8),allocatable :: gpp_prev_scpf(:)

real(r8) :: npp ! kg m-2 day-1

! Nutrient Flux Diagnostics

real(r8) :: resp_excess ! plant carbon respired due to carbon overflow
Expand Down Expand Up @@ -678,7 +675,6 @@ subroutine ZeroFluxDiags(this)

end do

this%npp = 0._r8
this%resp_excess = 0._r8
this%nh4_uptake = 0._r8
this%no3_uptake = 0._r8
Expand All @@ -694,11 +690,6 @@ subroutine ZeroFluxDiags(this)
this%p_uptake_scpf(:) = 0._r8
this%p_efflux_scpf(:) = 0._r8

! We don't zero gpp_prev_scpf because this is not
! incremented like others, it is assigned at the end
! of the daily history write process


return
end subroutine ZeroFluxDiags

Expand Down
42 changes: 35 additions & 7 deletions main/FatesRestartInterfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,8 @@ module FatesRestartInterfaceMod
integer :: ir_landuse_config_si
integer :: ir_gpp_acc_si
integer :: ir_aresp_acc_si
integer :: ir_herbivory_flux_out_si
integer :: ir_burn_flux_to_atm_si

integer :: ir_ncohort_pa
integer :: ir_canopy_layer_co
Expand Down Expand Up @@ -1199,6 +1201,15 @@ subroutine define_restart_vars(this, initialize_variables)
units='kg/ha', veclength=num_elements, flushval = flushzero, &
hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_errfates_mbal)

call this%RegisterCohortVector(symbol_base='herbivory_flux_out', vtype=site_r8, &
long_name_base='Mass flux of herbivory losses at the site level', &
units='kg/ha/day', veclength=num_elements, flushval = flushzero, &
hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_herbivory_flux_out_si)

call this%RegisterCohortVector(symbol_base='burn_flux_to_atm', vtype=site_r8, &
long_name_base='Mass flux of burn loss to the atmosphere at site level', &
units='kg/ha/day', veclength=num_elements, flushval = flushzero, &
hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_burn_flux_to_atm_si)

! Time integrated mass balance accounting [kg/m2]
call this%RegisterCohortVector(symbol_base='fates_liveveg_intflux', vtype=site_r8, &
Expand Down Expand Up @@ -2543,7 +2554,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
do i_lu_donor = 1, n_landuse_cats
do i_lu_receiver = 1, n_landuse_cats
do i_dist = 1, n_dist_types
rio_disturbance_rates_siluludi(io_idx_si_luludi) = sites(s)%disturbance_rates(i_dist,i_lu_donor, i_lu_receiver)
rio_disturbance_rates_siluludi(io_idx_si_luludi) = &
sites(s)%disturbance_rates(i_dist,i_lu_donor, i_lu_receiver)
io_idx_si_luludi = io_idx_si_luludi + 1
end do
end do
Expand All @@ -2557,19 +2569,30 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
io_idx_si_scpf = io_idx_co_1st

do i_cwd=1,ncwd
this%rvars(ir_cwdagin_flxdg+el-1)%r81d(io_idx_si_cwd) = sites(s)%flux_diags%elem(el)%cwd_ag_input(i_cwd)
this%rvars(ir_cwdbgin_flxdg+el-1)%r81d(io_idx_si_cwd) = sites(s)%flux_diags%elem(el)%cwd_bg_input(i_cwd)
this%rvars(ir_cwdagin_flxdg+el-1)%r81d(io_idx_si_cwd) = &
sites(s)%flux_diags%elem(el)%cwd_ag_input(i_cwd)
this%rvars(ir_cwdbgin_flxdg+el-1)%r81d(io_idx_si_cwd) = &
sites(s)%flux_diags%elem(el)%cwd_bg_input(i_cwd)
io_idx_si_cwd = io_idx_si_cwd + 1
end do

do i_pft=1,numpft
this%rvars(ir_leaflittin_flxdg+el-1)%r81d(io_idx_si_pft) = sites(s)%flux_diags%elem(el)%surf_fine_litter_input(i_pft)
this%rvars(ir_rootlittin_flxdg+el-1)%r81d(io_idx_si_pft) = sites(s)%flux_diags%elem(el)%root_litter_input(i_pft)
this%rvars(ir_woodprod_harvest_mbal+el-1)%r81d(io_idx_si_pft) = sites(s)%mass_balance(el)%wood_product_harvest(i_pft)
this%rvars(ir_woodprod_landusechange_mbal+el-1)%r81d(io_idx_si_pft) = sites(s)%mass_balance(el)%wood_product_landusechange(i_pft)
this%rvars(ir_leaflittin_flxdg+el-1)%r81d(io_idx_si_pft) = &
sites(s)%flux_diags%elem(el)%surf_fine_litter_input(i_pft)
this%rvars(ir_rootlittin_flxdg+el-1)%r81d(io_idx_si_pft) = &
sites(s)%flux_diags%elem(el)%root_litter_input(i_pft)
this%rvars(ir_woodprod_harvest_mbal+el-1)%r81d(io_idx_si_pft) = &
sites(s)%mass_balance(el)%wood_product_harvest(i_pft)
this%rvars(ir_woodprod_landusechange_mbal+el-1)%r81d(io_idx_si_pft) = &
sites(s)%mass_balance(el)%wood_product_landusechange(i_pft)
io_idx_si_pft = io_idx_si_pft + 1
end do

this%rvars(ir_herbivory_flux_out_si+el-1)%r81d(io_idx_si) = &
sites(s)%mass_balance(el)%herbivory_flux_out
this%rvars(ir_burn_flux_to_atm_si+el-1)%r81d(io_idx_si) = &
sites(s)%mass_balance(el)%burn_flux_to_atm

this%rvars(ir_oldstock_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%old_stock
this%rvars(ir_errfates_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%err_fates

Expand Down Expand Up @@ -3607,6 +3630,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
io_idx_si_pft = io_idx_si_pft + 1
end do

sites(s)%mass_balance(el)%herbivory_flux_out = &
this%rvars(ir_herbivory_flux_out_si+el-1)%r81d(io_idx_si)
sites(s)%mass_balance(el)%burn_flux_to_atm = &
this%rvars(ir_burn_flux_to_atm_si+el-1)%r81d(io_idx_si)

sites(s)%mass_balance(el)%old_stock = this%rvars(ir_oldstock_mbal+el-1)%r81d(io_idx_si)
sites(s)%mass_balance(el)%err_fates = this%rvars(ir_errfates_mbal+el-1)%r81d(io_idx_si)

Expand Down