Merge pull request #926 from JessicaNeedham/JessicaNeedham-IO-bagw-vars

Add mortality and productivity aboveground carbon fluxes as history variables

biogeochem/EDCohortDynamicsMod.F90

@@ -899,6 +899,12 @@ subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)
             currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
    end if
 
+   currentSite%term_abg_flux(currentCohort%size_class, currentCohort%pft) = &
+        currentSite%term_abg_flux(currentCohort%size_class, currentCohort%pft) + &
+        currentCohort%n * ( (struct_c+sapw_c+store_c) * prt_params%allom_agb_frac(currentCohort%pft) + &
+        leaf_c )
+
+
    ! put the litter from the terminated cohorts
    ! straight into the fragmenting pools
 

biogeochem/EDPatchDynamicsMod.F90

@@ -759,6 +759,14 @@ subroutine spawn_patches( currentSite, bc_in)
                                        (nc%n * ED_val_understorey_death / hlm_freq_day ) * &
                                        total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2
 
+                                  currentSite%imort_abg_flux(currentCohort%size_class, currentCohort%pft) = &
+                                       currentSite%imort_abg_flux(currentCohort%size_class, currentCohort%pft) + &
+                                       (nc%n * ED_val_understorey_death / hlm_freq_day ) * &
+                                       ( (sapw_c + struct_c + store_c) * prt_params%allom_agb_frac(currentCohort%pft) + &
+                                       leaf_c ) * &
+                                       g_per_kg * days_per_sec * years_per_day * ha_per_m2
+
+
                                   ! Step 2:  Apply survivor ship function based on the understory death fraction
                                   ! remaining of understory plants of those that are knocked over
                                   ! by the overstorey trees dying...
@@ -849,6 +857,14 @@ subroutine spawn_patches( currentSite, bc_in)
                                     total_c * g_per_kg * days_per_sec * ha_per_m2
                             end if
 
+                            currentSite%fmort_abg_flux(currentCohort%size_class, currentCohort%pft) = &
+                                 currentSite%fmort_abg_flux(currentCohort%size_class, currentCohort%pft) + &
+                                 (nc%n * currentCohort%fire_mort) * &
+                                 ( (sapw_c + struct_c + store_c) * prt_params%allom_agb_frac(currentCohort%pft) + &
+                                 leaf_c ) * &
+                                 g_per_kg * days_per_sec * ha_per_m2
+
+
                             currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = &
                                  currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + &
                                  nc%n * currentCohort%cambial_mort / hlm_freq_day

main/EDInitMod.F90

@@ -162,6 +162,10 @@ subroutine init_site_vars( site_in, bc_in, bc_out )
     allocate(site_in%fmort_carbonflux_canopy(1:numpft))
     allocate(site_in%fmort_carbonflux_ustory(1:numpft))
 
+    allocate(site_in%term_abg_flux(1:nlevsclass,1:numpft))
+    allocate(site_in%imort_abg_flux(1:nlevsclass,1:numpft))
+    allocate(site_in%fmort_abg_flux(1:nlevsclass,1:numpft))
+
     site_in%nlevsoil   = bc_in%nlevsoil
     allocate(site_in%rootfrac_scr(site_in%nlevsoil))
     allocate(site_in%zi_soil(0:site_in%nlevsoil))
@@ -278,6 +282,10 @@ subroutine zero_site( site_in )
     site_in%fmort_carbonflux_ustory(:) = 0._r8
     site_in%fmort_rate_cambial(:,:) = 0._r8
     site_in%fmort_rate_crown(:,:) = 0._r8
+    site_in%term_abg_flux(:,:) = 0._r8
+    site_in%imort_abg_flux(:,:) = 0._r8
+    site_in%fmort_abg_flux(:,:) = 0._r8
+
 
     ! fusoin-induced growth flux of individuals
     site_in%growthflux_fusion(:,:) = 0._r8

main/EDTypesMod.F90

@@ -847,6 +847,11 @@ module EDTypesMod
      real(r8), allocatable :: fmort_carbonflux_canopy(:) ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec]
      real(r8), allocatable :: fmort_carbonflux_ustory(:) ! biomass of understory indivs killed due to fire per year [gC/m2/sec] 
 
+     real(r8), allocatable :: term_abg_flux(:,:)          ! aboveground biomass lost due to termination mortality x size x pft
+     real(r8), allocatable :: imort_abg_flux(:,:)         ! aboveground biomass lost due to impact mortality x size x pft
+     real(r8), allocatable :: fmort_abg_flux(:,:)         ! aboveground biomass lost due to fire mortality x size x pft
+
+
      real(r8) :: demotion_carbonflux                     ! biomass of demoted individuals from canopy to understory [kgC/ha/day]
      real(r8) :: promotion_carbonflux                    ! biomass of promoted individuals from understory to canopy [kgC/ha/day]
      real(r8) :: recruitment_rate(1:maxpft)              ! number of individuals that were recruited into new cohorts

main/FatesHistoryInterfaceMod.F90

@@ -447,6 +447,9 @@ module FatesHistoryInterfaceMod
   integer :: ih_crownfiremort_si_scpf
   integer :: ih_cambialfiremort_si_scpf
 
+  integer :: ih_abg_mortality_cflux_si_scpf
+  integer :: ih_abg_productivity_cflux_si_scpf
+
   integer :: ih_m10_si_capf
   integer :: ih_nplant_si_capf
 
@@ -2339,6 +2342,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
                hio_crownfiremort_si_scpf     => this%hvars(ih_crownfiremort_si_scpf)%r82d, &
                hio_cambialfiremort_si_scpf   => this%hvars(ih_cambialfiremort_si_scpf)%r82d, &
 
+               hio_abg_mortality_cflux_si_scpf    => this%hvars(ih_abg_mortality_cflux_si_scpf)%r82d, &
+               hio_abg_productivity_cflux_si_scpf => this%hvars(ih_abg_productivity_cflux_si_scpf)%r82d, &
+
                hio_fire_c_to_atm_si  => this%hvars(ih_fire_c_to_atm_si)%r81d, &
                hio_burn_flux_elem    => this%hvars(ih_burn_flux_elem)%r82d, &
 
@@ -3083,6 +3089,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
                hio_npp_stor_si_scpf(io_si,scpf) = hio_npp_stor_si_scpf(io_si,scpf) + &
                   store_m_net_alloc*n_perm2 / days_per_year / sec_per_day
 
+
                ! Woody State Variables (basal area growth increment)
                if ( prt_params%woody(ft) == itrue) then
 
@@ -3133,6 +3140,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
                      ccohort%asmort*ccohort%n / m2_per_ha
                end if
 
+
+
+
                hio_m1_si_scls(io_si,scls) = hio_m1_si_scls(io_si,scls) + ccohort%bmort*ccohort%n / m2_per_ha
                hio_m2_si_scls(io_si,scls) = hio_m2_si_scls(io_si,scls) + ccohort%hmort*ccohort%n / m2_per_ha
                hio_m3_si_scls(io_si,scls) = hio_m3_si_scls(io_si,scls) + ccohort%cmort*ccohort%n / m2_per_ha
@@ -3227,12 +3237,30 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
                     (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_m * &
                     ccohort%n * ha_per_m2
 
+
                hio_hydraulicmortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_hydraulicmortality_carbonflux_si_pft(io_si,ccohort%pft) + &
                     ccohort%hmort * total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2
 
                hio_cstarvmortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_cstarvmortality_carbonflux_si_pft(io_si,ccohort%pft) + &
                     ccohort%cmort * total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2
 
+               ! Aboveground mortality
+               hio_abg_mortality_cflux_si_scpf(io_si,scpf) = hio_abg_mortality_cflux_si_scpf(io_si,scpf) + &
+                    (ccohort%bmort + ccohort%hmort + ccohort%cmort + &
+                    ccohort%frmort + ccohort%smort + ccohort%asmort) * &
+                    ( (sapw_m + struct_m + store_m ) * prt_params%allom_agb_frac(ccohort%pft) + &
+                    leaf_m ) * ccohort%n * days_per_sec * years_per_day * ha_per_m2 + &
+                    (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * &
+                    ( (sapw_m + struct_m + store_m ) * prt_params%allom_agb_frac(ccohort%pft) + &
+                    leaf_m ) * ccohort%n * ha_per_m2
+
+               ! Aboveground woody productivity
+               hio_abg_productivity_cflux_si_scpf(io_si,scpf) = hio_abg_productivity_cflux_si_scpf(io_si,scpf) + &
+                    ( (sapw_m_net_alloc + struct_m_net_alloc + store_m_net_alloc) * prt_params%allom_agb_frac(ccohort%pft) + &
+                    leaf_m_net_alloc ) * n_perm2 / &
+                    days_per_year / sec_per_day
+
+
                ! number density by size and biomass
                hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + &
                   total_m * ccohort%n * prt_params%allom_agb_frac(ccohort%pft) * AREA_INV
@@ -3654,6 +3682,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
          cpatch => cpatch%younger
       end do patchloop !patch loop
 
+
       ! divide basal-area-weighted height by basal area to get mean
       if ( sum(hio_ba_si_scpf(io_si,:)) .gt. nearzero ) then
          hio_ba_weighted_height_si(io_si) = hio_ba_weighted_height_si(io_si) / sum(hio_ba_si_scpf(io_si,:))
@@ -3792,11 +3821,23 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
          hio_mortality_carbonflux_si_pft(io_si,i_pft) = hio_mortality_carbonflux_si_pft(io_si,i_pft) + &
               (sites(s)%fmort_carbonflux_canopy(i_pft) + &
               sites(s)%fmort_carbonflux_ustory(i_pft) + &
-              sites(s)%imort_carbonflux(i_pft) ) / g_per_kg !cdk
-
+              sites(s)%imort_carbonflux(i_pft) ) / g_per_kg  ! cdk
+   
          hio_firemortality_carbonflux_si_pft(io_si,i_pft) = sites(s)%fmort_carbonflux_canopy(i_pft) / g_per_kg
       end do
 
+      ! add imort and fmort to aboveground woody mortality 
+      do i_pft = 1, numpft
+         do i_scls = 1,nlevsclass
+            i_scpf = (i_pft-1)*nlevsclass + i_scls
+            hio_abg_mortality_cflux_si_scpf(io_si,i_scpf) = hio_abg_mortality_cflux_si_scpf(io_si,i_scpf) + &
+                 (sites(s)%fmort_abg_flux(i_scls,i_pft) / g_per_kg ) + &
+                 (sites(s)%imort_abg_flux(i_scls,i_pft) / g_per_kg) + &
+                 (sites(s)%term_abg_flux(i_scls,i_pft)  * days_per_sec * ha_per_m2 ) ! jfn
+         end do
+      end do
+
+
       if(hlm_use_tree_damage .eq. itrue) then
 
          do i_pft = 1, numpft
@@ -3842,6 +3883,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in)
       sites(s)%fmort_rate_cambial(:,:) = 0._r8
       sites(s)%fmort_rate_crown(:,:) = 0._r8
       sites(s)%growthflux_fusion(:,:) = 0._r8
+      sites(s)%fmort_abg_flux(:,:) = 0._r8
+      sites(s)%imort_abg_flux(:,:) = 0._r8
+      sites(s)%term_abg_flux(:,:) = 0._r8
 
       sites(s)%imort_rate_damage(:,:,:) = 0.0_r8
       sites(s)%term_nindivs_canopy_damage(:,:,:) = 0.0_r8
@@ -6610,6 +6654,18 @@ subroutine define_history_vars(this, initialize_variables)
          upfreq=1, ivar=ivar, initialize=initialize_variables,                 &
          index=ih_cstarvmortality_carbonflux_si_pft)
 
+    call this%set_history_var(vname='FATES_ABOVEGROUND_MORT_SZPF', units='kg m-2 s-1',    &
+         long='Aboveground flux of carbon from AGB to necromass due to mortality', &
+         use_default='inactive', avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', &
+         upfreq=1, ivar=ivar, initialize=initialize_variables,                 &
+         index=ih_abg_mortality_cflux_si_scpf)
+
+    call this%set_history_var(vname='FATES_ABOVEGROUND_PROD_SZPF', units='kg m-2 s-1',    &
+         long='Aboveground carbon productivity', &
+         use_default='inactive', avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', &
+         upfreq=1, ivar=ivar, initialize=initialize_variables,                 &
+         index=ih_abg_productivity_cflux_si_scpf)
+
     call this%set_history_var(vname='MORTALITY_CROWNAREA_CANOPY',              &
           units = 'm2/ha/year',                                                &
           long='Crown area of canopy trees that died',                         &

main/FatesRestartInterfaceMod.F90

@@ -228,6 +228,9 @@ module FatesRestartInterfaceMod
   integer :: ir_imortcflux_sipft
   integer :: ir_fmortcflux_cano_sipft
   integer :: ir_fmortcflux_usto_sipft
+  integer :: ir_abg_term_flux_siscpf
+  integer :: ir_abg_imort_flux_siscpf
+  integer :: ir_abg_fmort_flux_siscpf
 
   integer :: ir_cwdagin_flxdg
   integer :: ir_cwdbgin_flxdg
@@ -1288,7 +1291,22 @@ subroutine define_restart_vars(this, initialize_variables)
    call this%set_restart_var(vname='fates_termcflux_ustory', vtype=cohort_r8, &
          long_name='fates diagnostic term carbon flux understory', &
          units='', flushval = flushzero, &
-         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index =   ir_termcflux_usto_sipft )
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_usto_sipft )
+
+   call this%set_restart_var(vname='fates_abg_term_flux', vtype=cohort_r8, &
+         long_name='fates aboveground biomass loss from termination mortality', &
+         units='', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_abg_term_flux_siscpf )
+
+   call this%set_restart_var(vname='fates_abg_imort_flux', vtype=cohort_r8, &
+         long_name='fates aboveground biomass loss from impact mortality', &
+         units='', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_abg_imort_flux_siscpf )
+
+   call this%set_restart_var(vname='fates_abg_fmort_flux', vtype=cohort_r8, &
+         long_name='fates aboveground biomass loss from fire mortality', &
+         units='', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_abg_fmort_flux_siscpf )
 
    call this%set_restart_var(vname='fates_democflux', vtype=site_r8, &
          long_name='fates diagnostic demotion carbon flux', &
@@ -1955,6 +1973,10 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
            rio_imortcflux_sipft        => this%rvars(ir_imortcflux_sipft)%r81d, &
            rio_fmortcflux_cano_sipft   => this%rvars(ir_fmortcflux_cano_sipft)%r81d, &
            rio_fmortcflux_usto_sipft   => this%rvars(ir_fmortcflux_usto_sipft)%r81d, &
+           rio_abg_imort_flux_siscpf => this%rvars(ir_abg_imort_flux_siscpf)%r81d, &
+           rio_abg_fmort_flux_siscpf => this%rvars(ir_abg_fmort_flux_siscpf)%r81d, &
+           rio_abg_term_flux_siscpf  => this%rvars(ir_abg_term_flux_siscpf)%r81d, &
+
            rio_imortrate_sicdpf        => this%rvars(ir_imortrate_sicdpf)%r81d, &
            rio_imortcflux_sicdsc       => this%rvars(ir_imortcflux_sicdsc)%r81d, &
            rio_termcflux_cano_sicdsc   => this%rvars(ir_termcflux_cano_sicdsc)%r81d, &
@@ -2335,6 +2357,10 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
                 rio_termnindiv_usto_siscpf(io_idx_si_scpf) = sites(s)%term_nindivs_ustory(i_scls,i_pft)
                 rio_growflx_fusion_siscpf(io_idx_si_scpf)  = sites(s)%growthflux_fusion(i_scls, i_pft)
 
+                rio_abg_term_flux_siscpf(io_idx_si_scpf) = sites(s)%term_abg_flux(i_scls, i_pft)
+                rio_abg_imort_flux_siscpf(io_idx_si_scpf) = sites(s)%imort_abg_flux(i_scls, i_pft)
+                rio_abg_fmort_flux_siscpf(io_idx_si_scpf) = sites(s)%fmort_abg_flux(i_scls, i_pft)
+
                 io_idx_si_scpf = io_idx_si_scpf + 1
              end do
 
@@ -2853,7 +2879,10 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
           rio_crownarea_usto_damage_si=> this%rvars(ir_crownarea_usto_si)%r81d, &
           rio_imortcflux_sipft        => this%rvars(ir_imortcflux_sipft)%r81d, &
           rio_fmortcflux_cano_sipft   => this%rvars(ir_fmortcflux_cano_sipft)%r81d, &
-          rio_fmortcflux_usto_sipft   => this%rvars(ir_fmortcflux_usto_sipft)%r81d)
+          rio_fmortcflux_usto_sipft   => this%rvars(ir_fmortcflux_usto_sipft)%r81d, &
+          rio_abg_term_flux_siscpf   => this%rvars(ir_abg_term_flux_siscpf)%r81d, &
+          rio_abg_imort_flux_siscpf  => this%rvars(ir_abg_imort_flux_siscpf)%r81d, &
+          rio_abg_fmort_flux_siscpf  => this%rvars(ir_abg_fmort_flux_siscpf)%r81d )
 
 
        totalcohorts = 0
@@ -3262,6 +3291,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
                 sites(s)%term_nindivs_canopy(i_scls,i_pft) = rio_termnindiv_cano_siscpf(io_idx_si_scpf)
                 sites(s)%term_nindivs_ustory(i_scls,i_pft) = rio_termnindiv_usto_siscpf(io_idx_si_scpf)
                 sites(s)%growthflux_fusion(i_scls, i_pft)  = rio_growflx_fusion_siscpf(io_idx_si_scpf)
+
+                sites(s)%term_abg_flux(i_scls,i_pft) = rio_abg_term_flux_siscpf(io_idx_si_scpf)
+                sites(s)%imort_abg_flux(i_scls,i_pft) = rio_abg_imort_flux_siscpf(io_idx_si_scpf)
+                sites(s)%fmort_abg_flux(i_scls,i_pft) = rio_abg_fmort_flux_siscpf(io_idx_si_scpf)
+
                 io_idx_si_scpf = io_idx_si_scpf + 1
              end do