diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index d69ffb254a..6c465b99b7 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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, & diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index a7a41c20ff..9149340074 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -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 @@ -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(:) @@ -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 @@ -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 @@ -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 + 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.' diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index d4c63b7ae7..faf91b1da6 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -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 @@ -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' @@ -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 diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index d750c55020..dcaf9cd972 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -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 diff --git a/parameter_files/fates_params_default.json b/parameter_files/fates_params_default.json index 1e799f8f1f..34bc923af4 100644 --- a/parameter_files/fates_params_default.json +++ b/parameter_files/fates_params_default.json @@ -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] },