Merge pull request #4 from rgknox/fire_mortality

fire mortality -- restart diagnostics

Author: jkshuman

biogeochem/EDCohortDynamicsMod.F90

@@ -21,6 +21,7 @@ module EDCohortDynamicsMod
   use EDTypesMod            , only : min_npm2, min_nppatch
   use EDTypesMod            , only : min_n_safemath
   use EDTypesMod            , only : nlevleaf
+  use EDTypesMod            , only : ican_upper
   use FatesInterfaceMod      , only : hlm_use_planthydro
   use FatesInterfaceMod      , only : hlm_parteh_mode
   use FatesPlantHydraulicsMod, only : FuseCohortHydraulics
@@ -606,12 +607,19 @@ subroutine terminate_cohorts( currentSite, currentPatch, level )
           ! preserve a record of the to-be-terminated cohort for mortality accounting
           levcan = currentCohort%canopy_layer
 
-          currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) = &
-               currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) + currentCohort%n
-          !
-          currentSite%termination_carbonflux(levcan) = currentSite%termination_carbonflux(levcan) + &
-                currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
-          
+          if(levcan==ican_upper) then
+             currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) = &
+                   currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) + currentCohort%n
+ 
+             currentSite%term_carbonflux_canopy = currentSite%term_carbonflux_canopy + &
+                   currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
+          else
+             currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) = &
+                   currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) + currentCohort%n
+ 
+             currentSite%term_carbonflux_ustory = currentSite%term_carbonflux_ustory + &
+                   currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
+          end if
 
           !put the litter from the terminated cohorts straight into the fragmenting pools
           if (currentCohort%n.gt.0.0_r8) then

biogeochem/EDPatchDynamicsMod.F90

@@ -17,6 +17,7 @@ module EDPatchDynamicsMod
   use EDTypesMod           , only : dtype_ifall
   use EDTypesMod           , only : dtype_ilog
   use EDTypesMod           , only : dtype_ifire
+  use EDTypesMod           , only : ican_upper
   use FatesInterfaceMod    , only : hlm_use_planthydro
   use FatesInterfaceMod    , only : hlm_numSWb
   use FatesInterfaceMod    , only : bc_in_type
@@ -538,20 +539,35 @@ subroutine spawn_patches( currentSite, bc_in)
 
                 levcan = currentCohort%canopy_layer 
 
-                ! before changing number densities, track total rate of trees that died due to fire, as well as from each fire mortality term
-                currentSite%fmort_rate(currentCohort%size_class, currentCohort%pft, levcan) = &
-                     currentSite%fmort_rate(currentCohort%size_class, currentCohort%pft, levcan) + &
-                     nc%n * currentCohort%fire_mort / hlm_freq_day
+                if(levcan==ican_upper) then
+                   
+                   ! before changing number densities, track total rate of trees that died 
+                   ! due to fire, as well as from each fire mortality term
+                   currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = &
+                         currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + &
+                         nc%n * currentCohort%fire_mort / hlm_freq_day
+
+                   currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + &
+                         (nc%n * currentCohort%fire_mort) * &
+                         total_c * g_per_kg * days_per_sec * ha_per_m2
+                   
+                else
+                   currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = &
+                         currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + &
+                         nc%n * currentCohort%fire_mort / hlm_freq_day
+
+                   currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + &
+                         (nc%n * currentCohort%fire_mort) * &
+                         total_c * g_per_kg * days_per_sec * ha_per_m2
+                end if
+
                 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
                 currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = &
                      currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + &
                      nc%n * currentCohort%crownfire_mort / hlm_freq_day
-                currentSite%fmort_carbonflux(levcan) = currentSite%fmort_carbonflux(levcan) + &
-                     (nc%n * currentCohort%fire_mort) * &
-                     total_c * g_per_kg * days_per_sec * ha_per_m2
-
+               
                 ! loss of individual from fire in new patch.
                 nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) 
 

main/EDInitMod.F90

@@ -70,11 +70,13 @@ subroutine init_site_vars( site_in )
     ! !LOCAL VARIABLES:
     !----------------------------------------------------------------------
     !
-    allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,1:nclmax))
+    allocate(site_in%term_nindivs_canopy(1:nlevsclass,1:numpft))
+    allocate(site_in%term_nindivs_ustory(1:nlevsclass,1:numpft))
     allocate(site_in%demotion_rate(1:nlevsclass))
     allocate(site_in%promotion_rate(1:nlevsclass))
     allocate(site_in%imort_rate(1:nlevsclass,1:numpft))
-    allocate(site_in%fmort_rate(1:nlevsclass,1:numpft,1:nclmax))
+    allocate(site_in%fmort_rate_canopy(1:nlevsclass,1:numpft))
+    allocate(site_in%fmort_rate_ustory(1:nlevsclass,1:numpft))
     allocate(site_in%fmort_rate_cambial(1:nlevsclass,1:numpft))
     allocate(site_in%fmort_rate_crown(1:nlevsclass,1:numpft))
     allocate(site_in%growthflux_fusion(1:nlevsclass,1:numpft))
@@ -127,13 +129,17 @@ subroutine zero_site( site_in )
     site_in%fates_to_bgc_last_ts = 0.0_r8
 
     ! termination and recruitment info
-    site_in%terminated_nindivs(:,:,:) = 0._r8
-    site_in%termination_carbonflux(:) = 0._r8
+    site_in%term_nindivs_canopy(:,:) = 0._r8
+    site_in%term_nindivs_ustory(:,:) = 0._r8
+    site_in%term_carbonflux_canopy = 0._r8
+    site_in%term_carbonflux_ustory = 0._r8
     site_in%recruitment_rate(:) = 0._r8
     site_in%imort_rate(:,:) = 0._r8
     site_in%imort_carbonflux = 0._r8
-    site_in%fmort_rate(:,:,:) = 0._r8
-    site_in%fmort_carbonflux(:) = 0._r8
+    site_in%fmort_rate_canopy(:,:) = 0._r8
+    site_in%fmort_rate_ustory(:,:) = 0._r8
+    site_in%fmort_carbonflux_canopy = 0._r8
+    site_in%fmort_carbonflux_ustory = 0._r8
     site_in%fmort_rate_cambial(:,:) = 0._r8
     site_in%fmort_rate_crown(:,:) = 0._r8
 

main/EDTypesMod.F90

@@ -605,26 +605,40 @@ module EDTypesMod
 
      ! PLANT HYDRAULICS
      type(ed_site_hydr_type), pointer :: si_hydr
-        
+     
+     ! DIAGNOSTICS
      ! TERMINATION, RECRUITMENT, DEMOTION, and DISTURBANCE
-
-     real(r8), allocatable :: terminated_nindivs(:,:,:)          ! number of individuals that were in cohorts which were terminated this timestep, on size x pft x canopy array. 
-     real(r8) :: termination_carbonflux(nclmax)                  ! carbon flux from live to dead pools associated with termination mortality, per canopy level
-     real(r8) :: recruitment_rate(1:maxpft)                      ! number of individuals that were recruited into new cohorts
-     real(r8), allocatable :: demotion_rate(:)                   ! rate of individuals demoted from canopy to understory per FATES timestep
+     
+     real(r8), allocatable :: term_nindivs_canopy(:,:) ! number of canopy individuals that were in cohorts which 
+                                                       ! were terminated this timestep, on size x pft
+     real(r8), allocatable :: term_nindivs_ustory(:,:) ! number of understory individuals that were in cohorts which 
+                                                       ! were terminated this timestep, on size x pft
+     real(r8) :: term_carbonflux_canopy                ! carbon flux from live to dead pools associated 
+                                                       ! with termination mortality, per canopy level
+     real(r8) :: term_carbonflux_ustory                ! carbon flux from live to dead pools associated 
+                                                       ! with termination mortality, per canopy level
+     
+     real(r8) :: recruitment_rate(1:maxpft)            ! number of individuals that were recruited into new cohorts
+     real(r8), allocatable :: demotion_rate(:)         ! rate of individuals demoted from canopy to understory per FATES timestep
      real(r8) :: demotion_carbonflux                             ! biomass of demoted individuals from canopy to understory [kgC/ha/day]
      real(r8), allocatable :: promotion_rate(:)                  ! rate of individuals promoted from understory to canopy per FATES timestep
      real(r8) :: promotion_carbonflux                            ! biomass of promoted individuals from understory to canopy [kgC/ha/day]
      real(r8), allocatable :: imort_rate(:,:)                    ! rate of individuals killed due to impact mortality per year.  on size x pft array
      real(r8) :: imort_carbonflux                                ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day]
 
-     real(r8), allocatable :: fmort_rate(:,:,:)                  ! rate of individuals killed due to fire mortality per year.  on size x pft x can-layer array 
-                                                                 ! (1:nlevsclass,1:numpft,1:nclmax)
-     real(r8) :: fmort_carbonflux(nclmax)                        ! biomass of individuals killed due to fire mortality per year. [gC/m2/sec]
-     real(r8), allocatable :: fmort_rate_cambial(:,:)            ! rate of individuals killed due to fire mortality from cambial damage per year.  on size x pft array
-     real(r8), allocatable :: fmort_rate_crown(:,:)              ! rate of individuals killed due to fire mortality from crown damage per year.  on size x pft array
-
-     real(r8), allocatable :: growthflux_fusion(:,:)             ! rate of individuals moving into a given size class bin due to fusion in a given day. on size x pft array 
+     real(r8), allocatable :: fmort_rate_canopy(:,:)             ! rate of canopy individuals killed due to fire mortality per year.  
+                                                                 ! on size x pft array  (1:nlevsclass,1:numpft)
+     real(r8), allocatable :: fmort_rate_ustory(:,:)             ! rate of understory individuals killed due to fire mortality per year.  
+                                                                 ! on size x pft array  (1:nlevsclass,1:numpft)
+     real(r8) :: fmort_carbonflux_canopy                         ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec]
+     real(r8) :: fmort_carbonflux_ustory                         ! biomass of understory indivs killed due to fire per year [gC/m2/sec] 
+     real(r8), allocatable :: fmort_rate_cambial(:,:)            ! rate of individuals killed due to fire mortality 
+                                                                 ! from cambial damage per year.  on size x pft array
+     real(r8), allocatable :: fmort_rate_crown(:,:)              ! rate of individuals killed due to fire mortality 
+                                                                 ! from crown damage per year.  on size x pft array
+
+     real(r8), allocatable :: growthflux_fusion(:,:)             ! rate of individuals moving into a given size class bin 
+                                                                 ! due to fusion in a given day. on size x pft array 
 
 
      ! some diagnostic-only (i.e. not resolved by ODE solver) flux of carbon to CWD and litter pools from termination and canopy mortality

main/FatesHistoryInterfaceMod.F90

@@ -2123,24 +2123,27 @@ subroutine update_history_dyn(this,nc,nsites,sites)
                i_scpf = (i_pft-1)*nlevsclass + i_scls
                !
                ! termination mortality. sum of canopy and understory indices
-               hio_m6_si_scpf(io_si,i_scpf) = sum(sites(s)%terminated_nindivs(i_scls,i_pft,1:nclmax)) * days_per_year
+               hio_m6_si_scpf(io_si,i_scpf) = (sites(s)%term_nindivs_canopy(i_scls,i_pft) + &
+                                               sites(s)%term_nindivs_ustory(i_scls,i_pft)) * days_per_year
 
                hio_m6_si_scls(io_si,i_scls) = hio_m6_si_scls(io_si,i_scls) + &
-                     sum(sites(s)%terminated_nindivs(i_scls,i_pft,1:nclmax)) * days_per_year
+                     (sites(s)%term_nindivs_canopy(i_scls,i_pft) + &
+                      sites(s)%term_nindivs_ustory(i_scls,i_pft)) * days_per_year
+                     
 
                !
                ! add termination mortality to canopy and understory mortality
                hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + &
-                    sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper) * days_per_year
+                    sites(s)%term_nindivs_canopy(i_scls,i_pft) * days_per_year
 
                hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + &
-                    sum(sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper+1:nclmax)) * days_per_year
+                    sites(s)%term_nindivs_ustory(i_scls,i_pft) * days_per_year
 
                hio_mortality_canopy_si_scpf(io_si,i_scpf) = hio_mortality_canopy_si_scpf(io_si,i_scpf) + &
-                    sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper) * days_per_year
+                     sites(s)%term_nindivs_canopy(i_scls,i_pft) * days_per_year
 
                hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + &
-                    sum(sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper+1:nclmax)) * days_per_year
+                     sites(s)%term_nindivs_ustory(i_scls,i_pft) * days_per_year
 
                !
                ! imort on its own
@@ -2158,44 +2161,45 @@ subroutine update_history_dyn(this,nc,nsites,sites)
                iscag = i_scls ! since imort is by definition something that only happens in newly disturbed patches, treat as such
                hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + &
                     sites(s)%imort_rate(i_scls, i_pft)
-               !
 
                ! fire mortality from the site-level diagnostic rates
-               hio_m5_si_scpf(io_si,i_scpf) = sum(sites(s)%fmort_rate(i_scls, i_pft,1:nclmax))
-               hio_m5_si_scls(io_si,i_scls) = hio_m5_si_scls(io_si,i_scls) + sum(sites(s)%fmort_rate(i_scls, i_pft,1:nclmax))
+               hio_m5_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_canopy(i_scls, i_pft) + &
+                     sites(s)%fmort_rate_ustory(i_scls, i_pft)
+               hio_m5_si_scls(io_si,i_scls) = hio_m5_si_scls(io_si,i_scls) + &
+                     sites(s)%fmort_rate_canopy(i_scls, i_pft) +  sites(s)%fmort_rate_ustory(i_scls, i_pft)
                !
                hio_crownfiremort_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_crown(i_scls, i_pft)
                hio_cambialfiremort_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_cambial(i_scls, i_pft)
                !
                ! fire components of overall canopy and understory mortality
                hio_mortality_canopy_si_scpf(io_si,i_scpf) = hio_mortality_canopy_si_scpf(io_si,i_scpf) + &
-                    sites(s)%fmort_rate(i_scls, i_pft,ican_upper)
+                    sites(s)%fmort_rate_canopy(i_scls, i_pft)
                hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + &
-                    sites(s)%fmort_rate(i_scls, i_pft,ican_upper)
+                    sites(s)%fmort_rate_canopy(i_scls, i_pft)
 
                ! the fire mortality rates for each layer are total dead, since the usable
                ! output will then normalize by the counts, we are allowed to sum over layers
                hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + &
-                     sum(sites(s)%fmort_rate(i_scls, i_pft,ican_upper+1:nclmax))
+                     sites(s)%fmort_rate_ustory(i_scls, i_pft)
 
                hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + &
-                     sum(sites(s)%fmort_rate(i_scls, i_pft,ican_upper+1:nclmax))
+                     sites(s)%fmort_rate_ustory(i_scls, i_pft)
 
                !
                ! carbon flux associated with mortality of trees dying by fire
                hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + &
-                    sites(s)%fmort_carbonflux(ican_upper)
-
+                     sites(s)%fmort_carbonflux_canopy
+               
                hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + &
-                     sum(sites(s)%fmort_carbonflux(ican_upper+1:nclmax))
-
+                     sites(s)%fmort_carbonflux_ustory
+               
                !
                ! for scag variables, also treat as happening in the newly-disurbed patch
 
                hio_mortality_canopy_si_scag(io_si,iscag) = hio_mortality_canopy_si_scag(io_si,iscag) + &
-                    sites(s)%fmort_rate(i_scls, i_pft,ican_upper)
+                    sites(s)%fmort_rate_canopy(i_scls, i_pft)
                hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + &
-                    sum(sites(s)%fmort_rate(i_scls, i_pft,ican_upper+1:nclmax))
+                    sites(s)%fmort_rate_ustory(i_scls, i_pft)
 
                ! while in this loop, pass the fusion-induced growth rate flux to history
                hio_growthflux_fusion_si_scpf(io_si,i_scpf) = hio_growthflux_fusion_si_scpf(io_si,i_scpf) + &
@@ -2208,11 +2212,14 @@ subroutine update_history_dyn(this,nc,nsites,sites)
          hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + &
               sites(s)%imort_carbonflux
          !
-         sites(s)%terminated_nindivs(:,:,:) = 0._r8
+         sites(s)%term_nindivs_canopy(:,:) = 0._r8
+         sites(s)%term_nindivs_ustory(:,:) = 0._r8
          sites(s)%imort_carbonflux = 0._r8
          sites(s)%imort_rate(:,:) = 0._r8
-         sites(s)%fmort_rate(:,:,:) = 0._r8
-         sites(s)%fmort_carbonflux(:) = 0._r8
+         sites(s)%fmort_rate_canopy(:,:) = 0._r8
+         sites(s)%fmort_rate_ustory(:,:) = 0._r8
+         sites(s)%fmort_carbonflux_canopy = 0._r8
+         sites(s)%fmort_carbonflux_ustory = 0._r8
          sites(s)%fmort_rate_cambial(:,:) = 0._r8
          sites(s)%fmort_rate_crown(:,:) = 0._r8
          sites(s)%growthflux_fusion(:,:) = 0._r8
@@ -2254,13 +2261,14 @@ subroutine update_history_dyn(this,nc,nsites,sites)
          ! mortality-associated carbon fluxes
 
          hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + &
-              sites(s)%termination_carbonflux(ican_upper) * g_per_kg * days_per_sec * ha_per_m2
-
+               sites(s)%term_carbonflux_canopy * g_per_kg * days_per_sec * ha_per_m2
+         
          hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + &
-              sum(sites(s)%termination_carbonflux(ican_upper+1:nclmax)) * g_per_kg * days_per_sec * ha_per_m2
+               sites(s)%term_carbonflux_ustory * g_per_kg * days_per_sec * ha_per_m2
 
          ! and zero the site-level termination carbon flux variable
-         sites(s)%termination_carbonflux(:) = 0._r8
+         sites(s)%term_carbonflux_canopy = 0._r8
+         sites(s)%term_carbonflux_ustory = 0._r8
          !
          ! add the site-level disturbance-associated cwd and litter input fluxes to thir respective flux fields
          do i_cwd = 1, ncwd