Skip to content

Commit

Permalink
Merge pull request #2 from rgknox/regeneration
Browse files Browse the repository at this point in the history 
memory and restart updates to TRS
  • Loading branch information
@adamhb
adamhb authored Jul 19, 2023
2 parents da80dd9 + fa4d124 commit 37c44c8
Show file tree
Hide file tree
Showing 8 changed files with 302 additions and 228 deletions.
97 changes: 52 additions & 45 deletions biogeochem/EDPatchDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ module EDPatchDynamicsMod
use FatesConstantsMod , only : r8 => fates_r8
use FatesConstantsMod , only : itrue, ifalse
use FatesConstantsMod , only : t_water_freeze_k_1atm
use FatesConstantsMod , only : TRS_regeneration
use FatesPlantHydraulicsMod, only : InitHydrCohort
use FatesPlantHydraulicsMod, only : AccumulateMortalityWaterStorage
use FatesPlantHydraulicsMod, only : DeallocateHydrCohort
Expand All @@ -64,6 +65,7 @@ module EDPatchDynamicsMod
use EDLoggingMortalityMod, only : get_harvestable_carbon
use EDLoggingMortalityMod, only : get_harvest_debt
use EDParamsMod , only : fates_mortality_disturbance_fraction
use EDParamsMod , only : regeneration_model
use FatesAllometryMod , only : carea_allom
use FatesAllometryMod , only : set_root_fraction
use FatesConstantsMod , only : g_per_kg
Expand Down Expand Up @@ -661,17 +663,19 @@ subroutine spawn_patches( currentSite, bc_in)
! --------------------------------------------------------------------------
call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24)
call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa)
call new_patch%seedling_layer_par24%CopyFromDonor(currentPatch%seedling_layer_par24)
call new_patch%sdlng_mort_par%CopyFromDonor(currentPatch%sdlng_mort_par)
call new_patch%sdlng2sap_par%CopyFromDonor(currentPatch%sdlng2sap_par)

do pft = 1,maxpft
call new_patch%sdlng_emerg_smp(pft)%p%CopyFromDonor(currentPatch%sdlng_emerg_smp(pft)%p)
call new_patch%sdlng_mdd(pft)%p%CopyFromDonor(currentPatch%sdlng_mdd(pft)%p)
enddo
if ( regeneration_model == TRS_regeneration ) then
call new_patch%seedling_layer_par24%CopyFromDonor(currentPatch%seedling_layer_par24)
call new_patch%sdlng_mort_par%CopyFromDonor(currentPatch%sdlng_mort_par)
call new_patch%sdlng2sap_par%CopyFromDonor(currentPatch%sdlng2sap_par)
do pft = 1,numpft
call new_patch%sdlng_emerg_smp(pft)%p%CopyFromDonor(currentPatch%sdlng_emerg_smp(pft)%p)
call new_patch%sdlng_mdd(pft)%p%CopyFromDonor(currentPatch%sdlng_mdd(pft)%p)
enddo
end if

call new_patch%tveg_longterm%CopyFromDonor(currentPatch%tveg_longterm)

! --------------------------------------------------------------------------
! The newly formed patch from disturbance (new_patch), has now been given
! some litter from dead plants and pre-existing litter from the donor patches.
Expand Down Expand Up @@ -2139,22 +2143,25 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft)
call new_patch%tveg24%InitRMean(fixed_24hr,init_value=temp_init_veg,init_offset=real(hlm_current_tod,r8) )
allocate(new_patch%tveg_lpa)
call new_patch%tveg_lpa%InitRmean(ema_lpa,init_value=temp_init_veg)
allocate(new_patch%seedling_layer_par24)
call new_patch%seedling_layer_par24%InitRMean(fixed_24hr,init_value=init_seedling_par, init_offset=real(hlm_current_tod,r8))
allocate(new_patch%sdlng_mort_par)
call new_patch%sdlng_mort_par%InitRMean(ema_sdlng_mort_par,init_value=temp_init_veg)
allocate(new_patch%sdlng2sap_par)
call new_patch%sdlng2sap_par%InitRMean(ema_sdlng2sap_par,init_value=init_seedling_par)


allocate(new_patch%sdlng_mdd(maxpft))
allocate(new_patch%sdlng_emerg_smp(maxpft))

do pft = 1,maxpft
allocate(new_patch%sdlng_mdd(pft)%p)
call new_patch%sdlng_mdd(pft)%p%InitRMean(ema_sdlng_mdd, init_value=0.0_r8)
allocate(new_patch%sdlng_emerg_smp(pft)%p)
call new_patch%sdlng_emerg_smp(pft)%p%InitRMean(ema_sdlng_emerg_h2o,init_value=init_seedling_smp)
enddo
if ( regeneration_model == TRS_regeneration ) then
allocate(new_patch%seedling_layer_par24)
call new_patch%seedling_layer_par24%InitRMean(fixed_24hr,init_value=init_seedling_par, init_offset=real(hlm_current_tod,r8))
allocate(new_patch%sdlng_mort_par)
call new_patch%sdlng_mort_par%InitRMean(ema_sdlng_mort_par,init_value=temp_init_veg)
allocate(new_patch%sdlng2sap_par)
call new_patch%sdlng2sap_par%InitRMean(ema_sdlng2sap_par,init_value=init_seedling_par)
allocate(new_patch%sdlng_mdd(numpft))
allocate(new_patch%sdlng_emerg_smp(numpft))
do pft = 1,numpft
allocate(new_patch%sdlng_mdd(pft)%p)
call new_patch%sdlng_mdd(pft)%p%InitRMean(ema_sdlng_mdd, init_value=0.0_r8)
allocate(new_patch%sdlng_emerg_smp(pft)%p)
call new_patch%sdlng_emerg_smp(pft)%p%InitRMean(ema_sdlng_emerg_h2o,init_value=init_seedling_smp)
enddo
end if


allocate(new_patch%tveg_longterm)
call new_patch%tveg_longterm%InitRmean(ema_longterm,init_value=temp_init_veg)
Expand Down Expand Up @@ -2730,16 +2737,18 @@ subroutine fuse_2_patches(csite, dp, rp)
! Weighted mean of the running means
call rp%tveg24%FuseRMean(dp%tveg24,rp%area*inv_sum_area)
call rp%tveg_lpa%FuseRMean(dp%tveg_lpa,rp%area*inv_sum_area)
call rp%seedling_layer_par24%FuseRMean(dp%seedling_layer_par24,rp%area*inv_sum_area) !ahb
call rp%tveg_longterm%FuseRMean(dp%tveg_longterm,rp%area*inv_sum_area)

do pft = 1,maxpft
call rp%sdlng_emerg_smp(pft)%p%FuseRMean(dp%sdlng_emerg_smp(pft)%p,rp%area*inv_sum_area) !ahb
call rp%sdlng_mdd(pft)%p%FuseRMean(dp%sdlng_mdd(pft)%p,rp%area*inv_sum_area) !ahb
enddo

call rp%sdlng_mort_par%FuseRMean(dp%sdlng_mort_par,rp%area*inv_sum_area) !ahb
call rp%sdlng2sap_par%FuseRMean(dp%sdlng2sap_par,rp%area*inv_sum_area) !ahb
if ( regeneration_model == TRS_regeneration ) then
call rp%seedling_layer_par24%FuseRMean(dp%seedling_layer_par24,rp%area*inv_sum_area)
call rp%sdlng_mort_par%FuseRMean(dp%sdlng_mort_par,rp%area*inv_sum_area)
call rp%sdlng2sap_par%FuseRMean(dp%sdlng2sap_par,rp%area*inv_sum_area)
do pft = 1,numpft
call rp%sdlng_emerg_smp(pft)%p%FuseRMean(dp%sdlng_emerg_smp(pft)%p,rp%area*inv_sum_area)
call rp%sdlng_mdd(pft)%p%FuseRMean(dp%sdlng_mdd(pft)%p,rp%area*inv_sum_area)
enddo
end if

call rp%tveg_longterm%FuseRMean(dp%tveg_longterm,rp%area*inv_sum_area)

rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area
rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area
Expand Down Expand Up @@ -3122,19 +3131,17 @@ subroutine dealloc_patch(cpatch)
endif

! Deallocate any running means
deallocate(cpatch%seedling_layer_par24)
deallocate(cpatch%sdlng_mort_par)
deallocate(cpatch%sdlng2sap_par)

do pft = 1, maxpft
deallocate(cpatch%sdlng_mdd(pft)%p)
deallocate(cpatch%sdlng_emerg_smp(pft)%p)
enddo


deallocate(cpatch%sdlng_mdd)
deallocate(cpatch%sdlng_emerg_smp)

if ( regeneration_model == TRS_regeneration ) then
deallocate(cpatch%seedling_layer_par24)
deallocate(cpatch%sdlng_mort_par)
deallocate(cpatch%sdlng2sap_par)
do pft = 1, numpft
deallocate(cpatch%sdlng_mdd(pft)%p)
deallocate(cpatch%sdlng_emerg_smp(pft)%p)
enddo
deallocate(cpatch%sdlng_mdd)
deallocate(cpatch%sdlng_emerg_smp)
end if

deallocate(cpatch%tveg24, stat=istat, errmsg=smsg)
if (istat/=0) then
Expand Down
86 changes: 43 additions & 43 deletions biogeochem/EDPhysiologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2312,61 +2312,61 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, bc_in, currentPatch )
! that times the ratio of (hypothetical) seed mass to recruit biomass

!==============================================================================================
do pft = 1,numpft
do pft = 1,numpft

! If the TRS's seedling dynamics is switched off, then we use FATES's default approach
! to germination
if ( regeneration_model == default_regeneration .or. &
if_tfs_or_def: if ( regeneration_model == default_regeneration .or. &
regeneration_model == TRS_no_seedling_dyn .or. &
prt_params%allom_dbh_maxheight(pft) < min_max_dbh_for_trees ) then

litt%seed_germ_in(pft) = min(litt%seed(pft) * EDPftvarcon_inst%germination_rate(pft), &
max_germination)*years_per_day
! If TRS seedling dynamics is switched on we calculate seedling emergence (i.e. germination)
! as a pft-specific function of understory light and soil moisture.
litt%seed_germ_in(pft) = min(litt%seed(pft) * EDPftvarcon_inst%germination_rate(pft), &
max_germination)*years_per_day

! If TRS seedling dynamics is switched on we calculate seedling emergence (i.e. germination)
! as a pft-specific function of understory light and soil moisture.
else if ( regeneration_model == TRS_regeneration .and. &
prt_params%allom_dbh_maxheight(pft) > min_max_dbh_for_trees ) then

! Step 1. Calculate how germination rate is modified by understory light
! This applies to photoblastic germinators (e.g. many tropical pioneers)
prt_params%allom_dbh_maxheight(pft) > min_max_dbh_for_trees ) then

! Calculate mean PAR at the seedling layer (MJ m-2 day-1) over the prior 24 hours
seedling_layer_par = currentPatch%seedling_layer_par24%GetMean() * sec_per_day * megajoules_per_joule
! Step 1. Calculate how germination rate is modified by understory light
! This applies to photoblastic germinators (e.g. many tropical pioneers)

! Calculate the photoblastic germination rate modifier (Eq. 3 Hanbury-Brown et al., 2022)
photoblastic_germ_modifier = seedling_layer_par / &
(seedling_layer_par + EDPftvarcon_inst%par_crit_germ(pft))
! Calculate mean PAR at the seedling layer (MJ m-2 day-1) over the prior 24 hours
seedling_layer_par = currentPatch%seedling_layer_par24%GetMean() * sec_per_day * megajoules_per_joule

! Step 2. Calculate how germination rate is modified by soil moisture in the rooting zone of
! the seedlings. This is a pft-specific running mean based on pft-specific seedling rooting
! depth.
! Calculate the photoblastic germination rate modifier (Eq. 3 Hanbury-Brown et al., 2022)
photoblastic_germ_modifier = seedling_layer_par / &
(seedling_layer_par + EDPftvarcon_inst%par_crit_germ(pft))

! Get running mean of soil matric potential (mm of H2O suction) at the seedling rooting depth
! This running mean based on pft-specific seedling rooting depth.
seedling_layer_smp = currentPatch%sdlng_emerg_smp(pft)%p%GetMean()
! Step 2. Calculate how germination rate is modified by soil moisture in the rooting zone of
! the seedlings. This is a pft-specific running mean based on pft-specific seedling rooting
! depth.

! Calculate a soil wetness index (1 / -soil matric pontential (MPa) ) used by the TRS
! to calculate seedling mortality from moisture stress.
wetness_index = 1.0_r8 / (seedling_layer_smp * (-1.0_r8) * mpa_per_mm_suction)
! Get running mean of soil matric potential (mm of H2O suction) at the seedling rooting depth
! This running mean based on pft-specific seedling rooting depth.
seedling_layer_smp = currentPatch%sdlng_emerg_smp(pft)%p%GetMean()

! Step 3. Calculate the seedling emergence rate based on soil moisture and germination
! rate modifier (Step 1). See Eq. 4 of Hanbury-Brown et al., 2022

! If SMP is below a pft-specific value, then no germination occurs
if ( seedling_layer_smp .GE. EDPftvarcon_inst%seedling_psi_emerg(pft) ) then
seedling_emerg_rate = photoblastic_germ_modifier * EDPftvarcon_inst%a_emerg(pft) * &
wetness_index**EDPftvarcon_inst%b_emerg(pft)
else

seedling_emerg_rate = 0.0_r8

end if ! End soil-moisture based seedling emergence rate

! Step 4. Calculate the amount of carbon germinating out of the seed bank
litt%seed_germ_in(pft) = litt%seed(pft) * seedling_emerg_rate

end if !End use TRS with seedling dynamics
! Calculate a soil wetness index (1 / -soil matric pontential (MPa) ) used by the TRS
! to calculate seedling mortality from moisture stress.
wetness_index = 1.0_r8 / (seedling_layer_smp * (-1.0_r8) * mpa_per_mm_suction)

! Step 3. Calculate the seedling emergence rate based on soil moisture and germination
! rate modifier (Step 1). See Eq. 4 of Hanbury-Brown et al., 2022

! If SMP is below a pft-specific value, then no germination occurs
if ( seedling_layer_smp .GE. EDPftvarcon_inst%seedling_psi_emerg(pft) ) then
seedling_emerg_rate = photoblastic_germ_modifier * EDPftvarcon_inst%a_emerg(pft) * &
wetness_index**EDPftvarcon_inst%b_emerg(pft)
else

seedling_emerg_rate = 0.0_r8

end if ! End soil-moisture based seedling emergence rate

! Step 4. Calculate the amount of carbon germinating out of the seed bank
litt%seed_germ_in(pft) = litt%seed(pft) * seedling_emerg_rate

end if if_tfs_or_def

!set the germination only under the growing season...c.xu

Expand Down
21 changes: 15 additions & 6 deletions biogeophys/FatesPlantRespPhotosynthMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2101,7 +2101,8 @@ subroutine LeafLayerMaintenanceRespiration_Ryan_1991(lnc_top, &
! Locals
real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s)
real(r8) :: lmr25top ! canopy top leaf maint resp rate at 25C for this pft (umol CO2/m**2/s)

integer :: c3c4_path_index ! Index for which photosynthetic pathway

! Parameter
real(r8), parameter :: lmrha = 46390._r8 ! activation energy for lmr (J/mol)
real(r8), parameter :: lmrhd = 150650._r8 ! deactivation energy for lmr (J/mol)
Expand All @@ -2117,7 +2118,10 @@ subroutine LeafLayerMaintenanceRespiration_Ryan_1991(lnc_top, &
! ----------------------------------------------------------------------------------
lmr25 = lmr25top * nscaler

if (nint(EDPftvarcon_inst%c3psn(ft)) /= 1) then
! photosynthetic pathway: 0. = c4, 1. = c3
c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft))

if (c3c4_path_index == c3_path_index) then
! temperature sensitivity of C3 plants
lmr = lmr25 * ft1_f(veg_tempk, lmrha) * &
fth_f(veg_tempk, lmrhd, lmrse, lmrc)
Expand Down Expand Up @@ -2249,7 +2253,7 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
! (umol electrons/m**2/s)
real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve
! (C4 plants) at 25C

integer :: c3c4_path_index ! Index for which photosynthetic pathway

! Parameters
! ---------------------------------------------------------------------------------
Expand Down Expand Up @@ -2312,14 +2316,19 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler

! Adjust for temperature
vcmax = vcmax25 * ft1_f(veg_tempk, vcmaxha) * fth_f(veg_tempk, vcmaxhd, vcmaxse, vcmaxc)
jmax = jmax25 * ft1_f(veg_tempk, jmaxha) * fth_f(veg_tempk, jmaxhd, jmaxse, jmaxc)
! photosynthetic pathway: 0. = c4, 1. = c3
c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft))

if (nint(EDPftvarcon_inst%c3psn(ft)) /= 1) then
if (c3c4_path_index == c3_path_index) then
vcmax = vcmax25 * ft1_f(veg_tempk, vcmaxha) * fth_f(veg_tempk, vcmaxhd, vcmaxse, vcmaxc)
else
vcmax = vcmax25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8)
vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-veg_tempk ) ))
vcmax = vcmax / (1._r8 + exp( 0.3_r8*(veg_tempk-(tfrz+40._r8)) ))
end if

jmax = jmax25 * ft1_f(veg_tempk, jmaxha) * fth_f(veg_tempk, jmaxhd, jmaxse, jmaxc)

!q10 response of product limited psn.
co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8)
end if
Expand Down
13 changes: 8 additions & 5 deletions main/EDInitMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,6 +274,10 @@ subroutine zero_site( site_in )
call site_in%flux_diags(el)%ZeroFluxDiags()
end do

! This will be initialized in FatesSoilBGCFluxMod:PrepCH4BCs()
! It checks to see if the value is below -9000. If it is,
! it will assume the first value of the smoother is set
site_in%ema_npp = -9999.9_r8

! termination and recruitment info
site_in%term_nindivs_canopy(:,:) = 0._r8
Expand Down Expand Up @@ -847,6 +851,10 @@ subroutine init_cohorts( site_in, patch_in, bc_in)
temp_cohort%canopy_trim = 1.0_r8
temp_cohort%crowndamage = 1 ! Assume no damage to begin with

! Retrieve drop fraction of non-leaf tissues for phenology initialisation
fnrt_drop_fraction = prt_params%phen_fnrt_drop_fraction(pft)
stem_drop_fraction = prt_params%phen_stem_drop_fraction(pft)


! Initialise phenology variables.
spmode_case: select case (hlm_use_sp)
Expand Down Expand Up @@ -913,11 +921,6 @@ subroutine init_cohorts( site_in, patch_in, bc_in)
temp_cohort%n = temp_cohort%n * sum(site_in%use_this_pft)
endif

! Retrieve drop fraction of non-leaf tissues for phenology initialisation
fnrt_drop_fraction = prt_params%phen_fnrt_drop_fraction(temp_cohort%pft)
stem_drop_fraction = prt_params%phen_stem_drop_fraction(temp_cohort%pft)



! h,dbh,leafc,n from SP values or from small initial size.
if(hlm_use_sp.eq.itrue)then
Expand Down
Loading

0 comments on commit 37c44c8

Please sign in to comment.