Skip to content

Commit a00482e

Browse files
rgknoxrgknox
authored
Merge pull request #292 from ckoven/fates_ppm_cbalancefixes
bug fixes on prescribed physiology carbon balance checks
2 parents 94ef851 + 3c24978 commit a00482e

File tree

8 files changed

+95
-67
lines changed

8 files changed

+95
-67
lines changed

‎biogeochem/EDCohortDynamicsMod.F90‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -863,7 +863,6 @@ subroutine fuse_cohorts(patchptr, bc_in)
863863
currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn
864864
currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn
865865
currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn
866-
currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn
867866
currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn
868867

869868
! logging mortality, Yi Xu
@@ -1254,7 +1253,6 @@ subroutine copy_cohort( currentCohort,copyc )
12541253
! Mortality diagnostics
12551254
n%cmort = o%cmort
12561255
n%bmort = o%bmort
1257-
n%imort = o%imort
12581256
n%fmort = o%fmort
12591257
n%hmort = o%hmort
12601258

‎biogeochem/EDPatchDynamicsMod.F90‎

Lines changed: 29 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,10 @@ module EDPatchDynamicsMod
3131
use EDLoggingMortalityMod, only : logging_litter_fluxes
3232
use EDLoggingMortalityMod, only : logging_time
3333
use EDParamsMod , only : fates_mortality_disturbance_fraction
34+
use FatesConstantsMod , only : g_per_kg
35+
use FatesConstantsMod , only : ha_per_m2
36+
use FatesConstantsMod , only : days_per_sec
37+
use FatesConstantsMod , only : years_per_day
3438

3539

3640
! CIME globals
@@ -117,7 +121,6 @@ subroutine disturbance_rates( site_in)
117121
currentCohort%cmort = cmort
118122
currentCohort%bmort = bmort
119123
currentCohort%hmort = hmort
120-
currentCohort%imort = 0.0_r8 ! Impact mortality is always zero except in new patches
121124
currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed
122125

123126
call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, &
@@ -199,7 +202,6 @@ subroutine disturbance_rates( site_in)
199202
currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction)
200203
currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction)
201204
currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction)
202-
! currentCohort%imort will likely exist with logging
203205
end if
204206
currentCohort => currentCohort%taller
205207
enddo !currentCohort
@@ -399,7 +401,6 @@ subroutine spawn_patches( currentSite, bc_in)
399401
nc%hmort = nan
400402
nc%bmort = nan
401403
nc%fmort = nan
402-
nc%imort = nan
403404
nc%lmort_logging = nan
404405
nc%lmort_collateral = nan
405406
nc%lmort_infra = nan
@@ -414,7 +415,18 @@ subroutine spawn_patches( currentSite, bc_in)
414415
! The number density per square are doesn't change, but since the patch is smaller
415416
! and cohort counts are absolute, reduce this number.
416417
nc%n = currentCohort%n * patch_site_areadis/currentPatch%area
417-
418+
419+
! because the mortality rate due to impact for the cohorts which had been in the understory and are now in the newly-
420+
! disturbed patch is very high, passing the imort directly to history results in large numerical errors, on account
421+
! of the sharply reduced number densities. so instead pass this info via a site-level diagnostic variable before reducing
422+
! the number density.
423+
currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = &
424+
currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + &
425+
nc%n * ED_val_understorey_death / hlm_freq_day
426+
currentSite%imort_carbonflux = currentSite%imort_carbonflux + &
427+
(nc%n * ED_val_understorey_death / hlm_freq_day ) * &
428+
currentCohort%b * g_per_kg * days_per_sec * years_per_day * ha_per_m2
429+
418430
! Step 2: Apply survivor ship function based on the understory death fraction
419431
! remaining of understory plants of those that are knocked over by the overstorey trees dying...
420432
nc%n = nc%n * (1.0_r8 - ED_val_understorey_death)
@@ -428,7 +440,6 @@ subroutine spawn_patches( currentSite, bc_in)
428440
! so with the number density must come the effective mortality rates.
429441

430442
nc%fmort = 0.0_r8 ! Should had also been zero in the donor
431-
nc%imort = ED_val_understorey_death/hlm_freq_day ! This was zero in the donor
432443
nc%cmort = currentCohort%cmort
433444
nc%hmort = currentCohort%hmort
434445
nc%bmort = currentCohort%bmort
@@ -443,6 +454,7 @@ subroutine spawn_patches( currentSite, bc_in)
443454
! Besides, the current and newly created patch sum to unity
444455

445456
currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)
457+
446458
else
447459
! grass is not killed by mortality disturbance events. Just move it into the new patch area.
448460
! Just split the grass into the existing and new patch structures
@@ -452,7 +464,6 @@ subroutine spawn_patches( currentSite, bc_in)
452464
currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)
453465

454466
nc%fmort = 0.0_r8
455-
nc%imort = 0.0_r8
456467
nc%cmort = currentCohort%cmort
457468
nc%hmort = currentCohort%hmort
458469
nc%bmort = currentCohort%bmort
@@ -478,7 +489,6 @@ subroutine spawn_patches( currentSite, bc_in)
478489
nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort)
479490

480491
nc%fmort = currentCohort%fire_mort/hlm_freq_day
481-
nc%imort = 0.0_r8
482492

483493
nc%cmort = currentCohort%cmort
484494
nc%hmort = currentCohort%hmort
@@ -508,7 +518,6 @@ subroutine spawn_patches( currentSite, bc_in)
508518
nc%hmort = nan
509519
nc%bmort = nan
510520
nc%fmort = nan
511-
nc%imort = nan
512521
nc%lmort_logging = nan
513522
nc%lmort_collateral = nan
514523
nc%lmort_infra = nan
@@ -526,6 +535,17 @@ subroutine spawn_patches( currentSite, bc_in)
526535
! The number density per square are doesn't change, but since the patch is smaller
527536
! and cohort counts are absolute, reduce this number.
528537
nc%n = currentCohort%n * patch_site_areadis/currentPatch%area
538+
539+
! because the mortality rate due to impact for the cohorts which had been in the understory and are now in the newly-
540+
! disturbed patch is very high, passing the imort directly to history results in large numerical errors, on account
541+
! of the sharply reduced number densities. so instead pass this info via a site-level diagnostic variable before reducing
542+
! the number density.
543+
currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = &
544+
currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + &
545+
nc%n * ED_val_understorey_death / hlm_freq_day
546+
currentSite%imort_carbonflux = currentSite%imort_carbonflux + &
547+
(nc%n * ED_val_understorey_death / hlm_freq_day ) * &
548+
currentCohort%b * g_per_kg * days_per_sec * years_per_day * ha_per_m2
529549

530550
! Step 2: Apply survivor ship function based on the understory death fraction
531551

@@ -540,7 +560,6 @@ subroutine spawn_patches( currentSite, bc_in)
540560

541561

542562
nc%fmort = 0.0_r8
543-
nc%imort = ED_val_understorey_death/hlm_freq_day
544563
nc%cmort = currentCohort%cmort
545564
nc%hmort = currentCohort%hmort
546565
nc%bmort = currentCohort%bmort
@@ -560,7 +579,6 @@ subroutine spawn_patches( currentSite, bc_in)
560579

561580
! No grass impact mortality imposed on the newly created patch
562581
nc%fmort = 0.0_r8
563-
nc%imort = 0.0_r8
564582
nc%cmort = currentCohort%cmort
565583
nc%hmort = currentCohort%hmort
566584
nc%bmort = currentCohort%bmort
@@ -992,7 +1010,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat
9921010
!currentCohort%dmort = mortality_rates(currentCohort)
9931011
!the disturbance calculations are done with the previous n, c_area and d_mort. So it's probably &
9941012
!not right to recalcualte dmort here.
995-
canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day)
1013+
canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day * fates_mortality_disturbance_fraction)
9961014

9971015
canopy_mortality_woody_litter = canopy_mortality_woody_litter + &
9981016
canopy_dead*(currentCohort%bdead+currentCohort%bsw)

‎biogeochem/EDPhysiologyMod.F90‎

Lines changed: 21 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -843,8 +843,10 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
843843
if (hlm_use_ed_prescribed_phys .eq. itrue) then
844844
if (currentCohort%canopy_layer .eq. 1) then
845845
currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(currentCohort%pft) * currentCohort%c_area / currentCohort%n
846+
currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year ! add these for balance checking purposes
846847
else
847848
currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(currentCohort%pft) * currentCohort%c_area / currentCohort%n
849+
currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year ! add these for balance checking purposes
848850
endif
849851
endif
850852

@@ -1077,22 +1079,27 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
10771079
temp_cohort%bstore = EDPftvarcon_inst%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) &
10781080
+ EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))
10791081

1080-
if (hlm_use_ed_prescribed_phys .eq. ifalse) then
1082+
if (hlm_use_ed_prescribed_phys .eq. ifalse .or. EDPftvarcon_inst%prescribed_recruitment(ft) .lt. 0. ) then
10811083
temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day &
10821084
/ (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore)
10831085
else
10841086
! prescribed recruitment rates. number per sq. meter per year
10851087
temp_cohort%n = currentPatch%area * EDPftvarcon_inst%prescribed_recruitment(ft) * hlm_freq_day
1088+
! modify the carbon balance accumulators to take into account the different way of defining recruitment
1089+
! add prescribed rates as an input C flux, and the recruitment that would have otherwise occured as an output flux
1090+
! (since the carbon associated with them effectively vanishes)
1091+
currentSite%flux_in = currentSite%flux_in + temp_cohort%n * (temp_cohort%bstore + temp_cohort%balive + temp_cohort%bdead)
1092+
currentSite%flux_out = currentSite%flux_out + currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day
10861093
endif
10871094

10881095
temp_cohort%laimemory = 0.0_r8
10891096
if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == 1)then
1090-
temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
1091-
EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
1097+
temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
1098+
EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
10921099
endif
10931100
if (EDPftvarcon_inst%stress_decid(temp_cohort%pft) == 1.and.currentSite%dstatus == 1)then
1094-
temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
1095-
EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
1101+
temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
1102+
EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
10961103
endif
10971104

10981105
cohortstatus = currentSite%status
@@ -1101,17 +1108,16 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
11011108
endif
11021109

11031110
if (temp_cohort%n > 0.0_r8 )then
1104-
if ( DEBUG ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort '
1105-
call create_cohort(currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
1106-
temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore, &
1107-
temp_cohort%laimemory, cohortstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, &
1108-
bc_in)
1111+
if ( DEBUG ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort '
1112+
call create_cohort(currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
1113+
temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore, &
1114+
temp_cohort%laimemory, cohortstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, &
1115+
bc_in)
11091116

1110-
! keep track of how many individuals were recruited for passing to history
1111-
currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
1117+
! keep track of how many individuals were recruited for passing to history
1118+
currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
11121119

11131120
endif
1114-
11151121
enddo !pft loop
11161122

11171123
deallocate(temp_cohort) ! delete temporary cohort
@@ -1172,7 +1178,7 @@ subroutine CWD_Input( currentSite, currentPatch)
11721178
SF_val_CWD_frac(c) * currentCohort%n/currentPatch%area *(1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft))
11731179
enddo
11741180

1175-
if (currentCohort%canopy_layer > 1)then
1181+
!if (currentCohort%canopy_layer > 1)then
11761182

11771183
! ================================================
11781184
! Litter fluxes for understorey mortality. KgC/m2/year
@@ -1271,7 +1277,7 @@ subroutine CWD_Input( currentSite, currentPatch)
12711277
currentSite%resources_management%delta_individual + &
12721278
(dead_n_dlogging+dead_n_ilogging) * hlm_freq_day * currentPatch%area
12731279

1274-
endif !canopy layer
1280+
!endif !canopy layer
12751281

12761282
currentCohort => currentCohort%taller
12771283
enddo ! end loop over cohorts

‎main/EDInitMod.F90‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@ subroutine init_site_vars( site_in )
6464
allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,2))
6565
allocate(site_in%demotion_rate(1:nlevsclass))
6666
allocate(site_in%promotion_rate(1:nlevsclass))
67+
allocate(site_in%imort_rate(1:nlevsclass,1:numpft))
6768
!
6869
end subroutine init_site_vars
6970

@@ -115,6 +116,8 @@ subroutine zero_site( site_in )
115116
site_in%terminated_nindivs(:,:,:) = 0._r8
116117
site_in%termination_carbonflux(:) = 0._r8
117118
site_in%recruitment_rate(:) = 0._r8
119+
site_in%imort_rate(:,:) = 0._r8
120+
site_in%imort_carbonflux = 0._r8
118121

119122
! demotion/promotion info
120123
site_in%demotion_rate(:) = 0._r8

‎main/EDMainMod.F90‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -609,7 +609,6 @@ subroutine bypass_dynamics(currentSite)
609609
currentCohort%bmort = 0.0_r8
610610
currentCohort%hmort = 0.0_r8
611611
currentCohort%cmort = 0.0_r8
612-
currentCohort%imort = 0.0_r8
613612
currentCohort%fmort = 0.0_r8
614613

615614
currentCohort => currentCohort%taller

‎main/EDTypesMod.F90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -226,7 +226,6 @@ module EDTypesMod
226226
real(r8) :: bmort ! background mortality rate n/year
227227
real(r8) :: cmort ! carbon starvation mortality rate n/year
228228
real(r8) :: hmort ! hydraulic failure mortality rate n/year
229-
real(r8) :: imort ! mortality from impacts by others n/year
230229
real(r8) :: fmort ! fire mortality n/year
231230

232231
! Logging Mortality Rate
@@ -540,6 +539,8 @@ module EDTypesMod
540539
real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day]
541540
real(r8), allocatable :: promotion_rate(:) ! rate of individuals promoted from understory to canopy per FATES timestep
542541
real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day]
542+
real(r8), allocatable :: imort_rate(:,:) ! rate of individuals killed due to impact mortality per year. on size x pft array
543+
real(r8) :: imort_carbonflux ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day]
543544

544545
! some diagnostic-only (i.e. not resolved by ODE solver) flux of carbon to CWD and litter pools from termination and canopy mortality
545546
real(r8) :: CWD_AG_diagnostic_input_carbonflux(1:ncwd) ! diagnostic flux to AG CWD [kg C / m2 / yr]
@@ -743,7 +744,6 @@ subroutine dump_cohort(ccohort)
743744
write(fates_log(),*) 'co%woody_turnover = ', ccohort%woody_turnover
744745
write(fates_log(),*) 'co%cmort = ', ccohort%cmort
745746
write(fates_log(),*) 'co%bmort = ', ccohort%bmort
746-
write(fates_log(),*) 'co%imort = ', ccohort%imort
747747
write(fates_log(),*) 'co%fmort = ', ccohort%fmort
748748
write(fates_log(),*) 'co%hmort = ', ccohort%hmort
749749
write(fates_log(),*) 'co%isnew = ', ccohort%isnew

0 commit comments

Comments
 (0)