Skip to content

Split initial density into to parameters with a switch on HLM side #1550

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
mvdebolskiy wants to merge 8 commits into
base: main
Choose a base branch
from
+76 −31
Open

Split initial density into to parameters with a switch on HLM side #1550

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
7164f22
split init_density parameter and make it deoend on the switch
mvdebolskiy Apr 3, 2026
a78990c
add protection against bad param values
mvdebolskiy Apr 3, 2026
cc15828
address review
mvdebolskiy Apr 27, 2026
44ffff1
Update main/FatesInterfaceMod.F90
mvdebolskiy Apr 30, 2026
06a887f
revert init_seed
mvdebolskiy May 4, 2026
92f414d
minor format change to json
glemieux May 6, 2026
8937fe8
fix param order
mvdebolskiy May 6, 2026
83cd29c
Merge tag 'sci.1.92.4_api.45.0.0' into add-init-dbh-nocomp
glemieux May 8, 2026
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
18 changes: 6 additions & 12 deletions main/EDInitMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@ module EDInitMod
use FatesInterfaceTypesMod , only : nlevcoage
use FatesInterfaceTypesMod , only : nlevdamage
use FatesInterfaceTypesMod , only : hlm_use_nocomp
use FatesInterfaceTypesMod , only : hlm_use_dbh_init
use FatesInterfaceTypesMod , only : nlevage
use FatesAllometryMod , only : h2d_allom
use FatesAllometryMod , only : h_allom
Expand Down Expand Up @@ -1217,7 +1218,7 @@ subroutine init_cohorts(site_in, patch_in, bc_in)
! if any pfts are starting with a non-recruitment size then the whole site
! needs the inventory type of spread
do pft = 1, numpft
if (EDPftvarcon_inst%initd(pft) < 0.0_r8) then
if (hlm_use_dbh_init .eq. itrue) then
site_in%spread = init_spread_inventory
end if
end do
Expand Down Expand Up @@ -1308,13 +1309,6 @@ subroutine init_cohorts(site_in, patch_in, bc_in)
end select phen_select
end if if_spmode

! If EDPftvarcon_inst%initd is positive, then it is interpreted as
! initial recruit density (stems/m2)
! If EDPftvarcon_inss%initd is negative, then it is interpreted as
! the initial DBH of the plant, and that the canopy is closed
! at one layer. However, this is only possible in nocomp since
! each PFT has its own patch area, and we don't have to assume
! the differences in coverage.

if_fullfates: if (hlm_use_nocomp .eq. ifalse) then

Expand All @@ -1328,8 +1322,8 @@ subroutine init_cohorts(site_in, patch_in, bc_in)

else ! We are in a nocomp simulation

! interpret as initial density and calculate diameter
if_init_dens: if (EDPftvarcon_inst%initd(pft) > nearzero) then
! Use initial density and calculate diameter
if_init_dens: if (hlm_use_dbh_init .eq. ifalse) then

cohort_n = EDPftvarcon_inst%initd(pft)*patch_in%area
! in nocomp mode we only have one PFT per patch
Expand All @@ -1346,9 +1340,9 @@ subroutine init_cohorts(site_in, patch_in, bc_in)
! calculate the plant diameter from height
call h2d_allom(height, pft, dbh)

else ! interpret as initial diameter and calculate density
else ! Use initial diameter and calculate density

dbh = abs(EDPftvarcon_inst%initd(pft))
dbh = EDPftvarcon_inst%initdbh(pft)

! calculate crown area of a single plant
call carea_allom(dbh, 1.0_r8, site_in%spread, pft, crown_damage, &
Expand Down
60 changes: 42 additions & 18 deletions main/EDPftvarcon.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,10 @@ module EDPftvarcon
real(r8), allocatable :: displar(:) ! ratio of displacement height to canopy top height
real(r8), allocatable :: bark_scaler(:) ! scaler from dbh to bark thickness. For fire model.
real(r8), allocatable :: crown_kill(:) ! scaler on fire death. For fire model.
real(r8), allocatable :: initd(:) ! initial seedling density [/m2] (positive values)
! or -dbh [cm] (negative values)
real(r8), allocatable :: initd(:) ! initial seedling density [/m2]
real(r8), allocatable :: initdbh(:) ! initial seedling dbh [cm]
! alternative to initd for nocomp
! used only when use_dbh_init and use_nocomp are both true
real(r8), allocatable :: init_seed(:) ! Initial seed bank [kg/m2]
! For SP: this is unused
! For Nocomp: This only applies the seed from the
Expand Down Expand Up @@ -349,6 +351,10 @@ subroutine TransferParamsPFT(pstruct)
allocate(EDPftvarcon_inst%initd(numpft))
EDPftvarcon_inst%initd(:) = param_p%r_data_1d(:)

param_p => pstruct%GetParamFromName('fates_recruit_init_dbh')
allocate(EDPftvarcon_inst%initdbh(numpft))
EDPftvarcon_inst%initdbh(:) = param_p%r_data_1d(:)

param_p => pstruct%GetParamFromName('fates_recruit_init_seed')
allocate(EDPftvarcon_inst%init_seed(numpft))
EDPftvarcon_inst%init_seed(:) = param_p%r_data_1d(:)
Expand Down Expand Up @@ -851,6 +857,7 @@ subroutine FatesReportPFTParams(is_master)
write(fates_log(),fmt0) 'bark_scaler = ',EDPftvarcon_inst%bark_scaler
write(fates_log(),fmt0) 'crown_kill = ',EDPftvarcon_inst%crown_kill
write(fates_log(),fmt0) 'initd = ',EDPftvarcon_inst%initd
write(fates_log(),fmt0) 'initdbh = ',EDPftvarcon_inst%initdbh
write(fates_log(),fmt0) 'init_seed = ',EDPftvarcon_inst%init_seed
write(fates_log(),fmt0) 'seed_suppl = ',EDPftvarcon_inst%seed_suppl
write(fates_log(),fmt0) 'lf_flab = ',EDPftvarcon_inst%lf_flab
Expand Down Expand Up @@ -957,6 +964,7 @@ subroutine FatesCheckParams(is_master)
use FatesInterfaceTypesMod, only : hlm_use_fixed_biogeog,hlm_use_sp, hlm_name
use FatesInterfaceTypesMod, only : hlm_use_inventory_init
use FatesInterfaceTypesMod, only : hlm_use_nocomp
use FatesInterfaceTypesMod, only : hlm_use_dbh_init
use EDParamsMod , only : max_nocomp_pfts_by_landuse, maxpatches_by_landuse
use FatesConstantsMod , only : n_landuse_cats

Expand Down Expand Up @@ -1259,30 +1267,46 @@ subroutine FatesCheckParams(is_master)
! Check that in initial density is not equal to zero in a cold-start run
!-----------------------------------------------------------------------------------

if ( EDPftvarcon_inst%initd(ipft) < -nearzero ) then
write(fates_log(),*) ' Option to use negative values for initial density'
write(fates_log(),*) ' has been depricated. Set use_fates_dbh_init to true'
write(fates_log(),*) ' to initialize with dbh in nocomp mode. The corresponding'
write(fates_log(),*) ' is called fates_recruit_init_dbh'
write(fates_log(),*) ' Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if

if ( hlm_use_inventory_init == ifalse .and. &
abs( EDPftvarcon_inst%initd(ipft) ) < nearzero ) then
EDPftvarcon_inst%initd(ipft) < nearzero ) then

write(fates_log(),*) ' In a cold start run initial density cannot be zero.'
write(fates_log(),*) ' For a bare ground run set to initial recruit density.'
write(fates_log(),*) ' In a cold start run initial density cannot be zero'
write(fates_log(),*) ' without inventory initialization.'
write(fates_log(),*) ' If no-comp is on it is possible to initialize with larger '
write(fates_log(),*) ' plants by setting fates_recruit_init_density to a negative number'
write(fates_log(),*) ' which will be interpreted as (absolute) initial dbh. '
write(fates_log(),*) ' plants by setting use_fates_dbh_init to true'
write(fates_log(),*) ' so that initial dbh parameter is used instead of density. '
write(fates_log(),*) ' Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))

end if

if ( hlm_use_nocomp .eq. ifalse .and. EDPftvarcon_inst%initd(ipft) < -nearzero ) then
write(fates_log(),*) ' When not in a noncomp configuration, FATES does not'
write(fates_log(),*) ' know how to interpret a negative %initd (number density)'
write(fates_log(),*) ' on a cold-start. In nocomp with a negative, we assume the absolute'
write(fates_log(),*) ' value of the %inidt parameter is the initial plant size.'
write(fates_log(),*) ' And since the fractional area of each PFT is known from'
write(fates_log(),*) ' the surface file, we can derive a number density from this'
write(fates_log(),*) ' However, we do not have a hypothesis to do this in full FATES.'
write(fates_log(),*) ' Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if



if (hlm_use_dbh_init .eq. itrue) then
if (EDPftvarcon_inst%initdbh(ipft) < nearzero) then
write(fates_log(),*) ' You are running in nocomp mode using initial dbh instead of density.'
write(fates_log(),*) ' In a cold start run initial dbh cannot be less or equal zero.'
write(fates_log(),*) ' Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))
endif
if (EDPftvarcon_inst%init_seed(ipft) < nearzero) then
Copy link
Copy Markdown
Contributor

@JessicaNeedham JessicaNeedham Apr 24, 2026

Choose a reason for hiding this comment

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

Why do we need an extra check for init_seed being positive if we are in dbh initialisation mode? There's another check right below that should catch it.

Copy link
Copy Markdown
Contributor Author

@mvdebolskiy mvdebolskiy Apr 24, 2026

Choose a reason for hiding this comment

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

You can not initiate the model with 0 seed and dbh. It will crash within first few days because of the canopy. Basically you are starting with fat grown tries with no carbon?

Copy link
Copy Markdown
Contributor

@JessicaNeedham JessicaNeedham Apr 24, 2026

Choose a reason for hiding this comment

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

That makes sense. I guess I was wondering if this is specific to dbh initialisation? I'm assuming the same applies if you have 0 seeds and 0 number density? So is the check here not sufficient? https://github.com/mvdebolskiy/fates/blob/a78990cfcd80856c12ec9a25d91f9df8bb409490/main/EDPftvarcon.F90#L1302

Copy link
Copy Markdown
Contributor Author

@mvdebolskiy mvdebolskiy Apr 25, 2026

Choose a reason for hiding this comment

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

With 0 seeds non-zero density is fine.

write(fates_log(),*) ' You are running in nocomp mode using initial dbh instead of density.'
write(fates_log(),*) ' In a cold start run initial seed cannot equal zero '
write(fates_log(),*) ' Otherwize mass-balance will not be conserved.'
write(fates_log(),*) ' Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))
endif
endif

if ( EDPftvarcon_inst%init_seed(ipft) < 0._r8) then
write(fates_log(),*) ' Initial seed pool fates_init_seed can not be negative.'
Expand Down
15 changes: 15 additions & 0 deletions main/FatesInterfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1575,6 +1575,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
hlm_use_nocomp = unset_int
hlm_use_sp = unset_int
hlm_use_inventory_init = unset_int
hlm_use_dbh_init = unset_int
hlm_inventory_ctrl_file = 'unset'
hlm_hist_level_dynam = unset_int
hlm_hist_level_hifrq = unset_int
Expand Down Expand Up @@ -1683,6 +1684,14 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
call endrun(msg=errMsg(sourcefile, __LINE__))
end if

if ( .not.((hlm_use_dbh_init.eq.itrue).or.(hlm_use_dbh_init.eq.ifalse)) ) then
write(fates_log(), *) 'The Fates dbh init flag must be 0 or 1, exiting'
call endrun(msg=errMsg(sourcefile, __LINE__))
elseif ((hlm_use_dbh_init .eq. itrue) .and. (hlm_use_nocomp .eq. ifalse)) then
write(fates_log(), *) 'The Fates dbh init can only be ON in NOCOMP mode, exiting'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if

if(hlm_ivis .ne. ivis) then
write(fates_log(), *) 'FATES assumption about the index of visible shortwave'
write(fates_log(), *) 'radiation is different from the HLM, exiting'
Expand Down Expand Up @@ -2191,6 +2200,12 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
write(fates_log(),*) 'Transfering hlm_use_inventory_init= ',ival,' to FATES'
end if

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

case('hist_hifrq_dimlevel')
hlm_hist_level_hifrq = ival
if (fates_global_verbose()) then
Expand Down
5 changes: 5 additions & 0 deletions main/FatesInterfaceTypesMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,6 +211,11 @@ module FatesInterfaceTypesMod
! This need only be defined when
! hlm_use_inventory_init = 1

integer, public :: hlm_use_dbh_init ! Flag to use fates_recruit_init_dbh
! instead of fates_recruit_init_density
! only works in nocomp mode
! 1 = TRUE, 0 = FALSE

integer, public :: hlm_use_fixed_biogeog ! Flag to use FATES fixed biogeography mode
! 1 = TRUE, 0 = FALSE

Expand Down
9 changes: 8 additions & 1 deletion parameter_files/fates_params_default.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1269,10 +1269,17 @@
"units": "m",
"data": [1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 0.2, 0.2, 0.2, 0.8, 0.8, 0.11, 0.2, 0.2]
},
"fates_recruit_init_dbh": {
"dtype": "float",
"dims": ["fates_pft"],
"long_name": "initial seedling dbh for a cold-start near-bare-ground simulation. Can only to be used in nocomp mode. To use it set use_dbh_init to true",
"units": "cm",
"data": [ 1, 1, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.1, 0.1, 0.1]
},
"fates_recruit_init_density": {
"dtype": "float",
"dims": ["fates_pft"],
"long_name": "initial seedling density for a cold-start near-bare-ground simulation. If negative sets initial tree dbh - only to be used in nocomp mode",
"long_name": "initial seedling density for a cold-start near-bare-ground simulation.",
"units": "stems/m2",
"data": [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.16, 0.2, 0.2, 0.2, 0.2]
},
Expand Down