diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index f509aec000..111ac0d5f3 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -541,6 +541,7 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch) integer :: nlevsoil ! number of soil layers integer :: ilyr ! soil layer loop counter integer :: dcmpy ! decomposability index + real(r8) :: net_seed_available ! Available seed after inputs and decay do el = 1, num_elements @@ -550,15 +551,20 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch) ! ----------------------------------------------------------------------------------- do pft = 1,numpft - litt%seed(pft) = litt%seed(pft) + & - litt%seed_in_local(pft) + & - litt%seed_in_extern(pft) - & - litt%seed_decay(pft) - & - litt%seed_germ_in(pft) - - ! Note that the recruitment scheme will use seed_germ - ! for its construction costs. - litt%seed_germ(pft) = litt%seed_germ(pft) + & + + ! Calculate net available seed after inputs and decay + net_seed_available = litt%seed(pft) + litt%seed_in_local(pft) + & + litt%seed_in_extern(pft) - litt%seed_decay(pft) + + ! Cap germination at available seed to prevent negative seed pool + litt%seed_germ_in(pft) = min(litt%seed_germ_in(pft), net_seed_available) + + ! Update pools + litt%seed(pft) = net_seed_available - litt%seed_germ_in(pft) + + ! Note that the recruitment scheme will use seed_germ + ! for its construction costs. + litt%seed_germ(pft) = litt%seed_germ(pft) + & litt%seed_germ_in(pft) - & litt%seed_germ_decay(pft) @@ -2253,6 +2259,7 @@ subroutine SeedDecay( litt , currentPatch ) real(r8) :: seedling_light_mort_rate ! daily seedling mortality rate from light stress real(r8) :: seedling_h2o_mort_rate ! daily seedling mortality rate from moisture stress real(r8) :: seedling_mdds ! moisture deficit days accumulated in the seedling layer + real(r8) :: total_seedling_mort_rate ! summed seedling mortality rates !---------------------------------------------------------------------- @@ -2317,25 +2324,42 @@ subroutine SeedDecay( litt , currentPatch ) ! Get the current seedling moisture deficit days (tracked as a pft-specific exponential ! average) - seedling_mdds = currentPatch%sdlng_mdd(pft)%p%GetMean() + seedling_mdds = currentPatch%sdlng_mdd(pft)%p%GetMean() ! Calculate seedling mortality as a function of moisture deficit days (mdd) ! If the seedling mmd value is below a critical threshold then moisture-based mortality is zero + if (seedling_mdds < EDPftvarcon_inst%seedling_mdd_crit(pft)) then seedling_h2o_mort_rate = 0.0_r8 else seedling_h2o_mort_rate = EDPftvarcon_inst%seedling_h2o_mort_a(pft) * seedling_mdds**2 + & EDPftvarcon_inst%seedling_h2o_mort_b(pft) * seedling_mdds + & EDPftvarcon_inst%seedling_h2o_mort_c(pft) + ! Cap h2o mortality rate at 1, quadratic can produce values >1 for large moisture + ! deficit days, if rate exceeds 1 write to fates log + + if (seedling_h2o_mort_rate > 1.0_r8) then + write(fates_log(), *) 'TRS seedling_h2o_mort_rate > 1', & + 'pft=', pft, 'uncapped=', seedling_h2o_mort_rate, 'mdds=', seedling_mdds + end if + seedling_h2o_mort_rate = min(1.0_r8, seedling_h2o_mort_rate) end if ! mdd threshold check ! Step 3. Sum modes of mortality (including background mortality) and send dead seedlings ! to litter - litt%seed_germ_decay(pft) = (litt%seed_germ(pft) * seedling_light_mort_rate) + & - (litt%seed_germ(pft) * seedling_h2o_mort_rate) + & - (litt%seed_germ(pft) * EDPftvarcon_inst%background_seedling_mort(pft) & - * years_per_day) - + ! Cap total mortality rate at 1 to prevent negative seedling pool + total_seedling_mort_rate = seedling_light_mort_rate + & + seedling_h2o_mort_rate + & + (EDPftvarcon_inst%background_seedling_mort(pft) * years_per_day) + if (total_seedling_mort_rate > 1.0_r8) then + write(fates_log(), *) 'TRS total seedling mortality rate > 1', & + 'pft=', pft, 'uncapped=', total_seedling_mort_rate, & + 'seedling_light_mort_rate=', seedling_light_mort_rate, & + 'seedling_h2o_mort_rate=', seedling_h2o_mort_rate + end if + litt%seed_germ_decay(pft) = litt%seed_germ(pft) * & + min(1.0_r8, total_seedling_mort_rate) + else litt%seed_germ_decay(pft) = litt%seed_germ(pft) * & @@ -2440,6 +2464,15 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, currentPatch ) 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) + ! Cap emergence rate at 1, rate can exceed 1 for some parameter combinations, + ! leading to negative seed bank, write to log if this is the case + if (seedling_emerg_rate > 1.0_r8) then + write(fates_log(),*) 'TRS seedling_emerg_rate > 1', & + 'pft=', pft, 'uncapped=', seedling_emerg_rate, & + 'photoblastic_germ_modifier=', photoblastic_germ_modifier, & + 'wetness_index=', wetness_index + end if + seedling_emerg_rate = min(1.0_r8, seedling_emerg_rate ) else seedling_emerg_rate = 0.0_r8 @@ -2522,6 +2555,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in) real(r8) :: stem_drop_fraction ! real(r8) :: fnrt_drop_fraction ! real(r8) :: sdlng2sap_par ! running mean of PAR at the seedling layer [MJ/m2/day] + real(r8) :: sdlng2sap_rate ! seedling-to-sapling transition rate factor (unitless) real(r8) :: seedling_layer_smp ! soil matric potential at seedling rooting depth [mm H2O suction] integer, parameter :: recruitstatus = 1 ! whether the newly created cohorts are recruited or initialized integer :: ilayer_seedling_root ! the soil layer at seedling rooting depth @@ -2667,11 +2701,22 @@ subroutine recruitment(currentSite, currentPatch, bc_in) sdlng2sap_par = currentPatch%sdlng2sap_par%GetMean()* & sec_per_day*megajoules_per_joule + sdlng2sap_rate = EDPftvarcon_inst%seedling_light_rec_a(ft)* & + sdlng2sap_par**EDPftvarcon_inst%seedling_light_rec_b(ft) + + ! If the seedling to sapling transition rate exceeds 1, + ! cap at 1 (prevent seed_germ from going negative) and write to fates log + + if (sdlng2sap_rate > 1.0_r8) then + write(fates_log(), *) 'TRS seedling-to-sapling transition rate > 1', & + 'pft=', ft, 'uncapped=', sdlng2sap_rate, & + 'sdlng2sap_par=', sdlng2sap_par + end if + mass_avail = currentPatch%area* & currentPatch%litter(el)%seed_germ(ft)* & - EDPftvarcon_inst%seedling_light_rec_a(ft)* & - sdlng2sap_par**EDPftvarcon_inst%seedling_light_rec_b(ft) - + min(1.0_r8, sdlng2sap_rate) + ! If soil moisture is below pft-specific seedling moisture stress threshold the ! recruitment does not occur. ilayer_seedling_root = minloc(abs(bc_in%z_sisl(:) - & @@ -2780,6 +2825,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in) currentPatch%litter(el)%seed_germ(ft) = & currentPatch%litter(el)%seed_germ(ft) - cohort_n / currentPatch%area * & (m_struct + m_leaf + m_fnrt + m_sapw + m_store + m_repro) + end if end do diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 2f358cda63..557fcecaa3 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -666,7 +666,7 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft) call this%seedling_layer_par24%InitRMean(fixed_24hr, & init_value=init_seedling_par, init_offset=real(current_tod, r8)) call this%sdlng_mort_par%InitRMean(ema_sdlng_mort_par, & - init_value=temp_init_veg) + init_value=init_seedling_par) call this%sdlng2sap_par%InitRMean(ema_sdlng2sap_par, & init_value=init_seedling_par) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 6982dd700b..13422ac3e4 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -2335,8 +2335,12 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) endif ! Update the seedling layer smp and mdd running means - call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp) - call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd) + ! Only update after first model day to avoid adding + ! initialization spike in smp to running means + if (hlm_model_day > 2.0_r8) then + call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd) + call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp) + endif enddo !end pft loop @@ -2410,8 +2414,8 @@ subroutine SeedlingParPatch(cpatch, & ! Start with the assumption that there is a single canopy layer seedling_par_high = atm_par_dir+atm_par_dif - par_high_frac = 1._r8-cpatch%total_canopy_area - par_low_frac = cpatch%total_canopy_area + par_high_frac = 1._r8 - (cpatch%total_canopy_area / cpatch%area) + par_low_frac = cpatch%total_canopy_area / cpatch%area ! Work up through the canopy layers from the bottom layer do cl = cpatch%NCL_p,max(1,cpatch%NCL_p-1),-1