Merge pull request #431 from jkshuman/fire_mortality

Fire mortality diagnostics update

Author: rgknox

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
@@ -602,19 +603,23 @@ subroutine terminate_cohorts( currentSite, currentPatch, level )
       endif    !  if (.not.currentCohort%isnew .and. level == 2) then
 
       if (terminate == 1) then 
+         
           ! preserve a record of the to-be-terminated cohort for mortality accounting
-          if (currentCohort%canopy_layer .eq. 1) then
-             levcan = 1
-          else
-             levcan = 2
-          endif
+          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
@@ -927,7 +932,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                                    currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn
                                    currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn
                                    currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn
-                                   currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn
                                    currentCohort%frmort = (currentCohort%n*currentCohort%frmort + nextc%n*nextc%frmort)/newn
 
                                    ! logging mortality, Yi Xu
@@ -1303,7 +1307,6 @@ subroutine copy_cohort( currentCohort,copyc )
     ! Mortality diagnostics
     n%cmort = o%cmort
     n%bmort = o%bmort
-    n%fmort = o%fmort
     n%hmort = o%hmort
     n%frmort = o%frmort
 

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
@@ -139,7 +140,6 @@ subroutine disturbance_rates( site_in, bc_in)
           currentCohort%bmort = bmort
           currentCohort%hmort = hmort
           currentCohort%frmort = frmort
-          currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed
 
           call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, &
                 lmort_direct,lmort_collateral,lmort_infra )
@@ -216,7 +216,6 @@ subroutine disturbance_rates( site_in, bc_in)
           currentCohort => currentPatch%shortest
           do while(associated(currentCohort))
              if(currentCohort%canopy_layer == 1)then
-                currentCohort%fmort = 0.0_r8
                 currentCohort%cmort = currentCohort%cmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction)
@@ -267,7 +266,6 @@ subroutine disturbance_rates( site_in, bc_in)
                 currentCohort%lmort_direct    = 0.0_r8
                 currentCohort%lmort_collateral = 0.0_r8
                 currentCohort%lmort_infra      = 0.0_r8
-                currentCohort%fmort            = 0.0_r8
              end if
              currentCohort => currentCohort%taller
           enddo !currentCohort
@@ -323,13 +321,15 @@ subroutine spawn_patches( currentSite, bc_in)
     real(r8) :: leaf_litter_local(maxpft)    ! initial value of leaf litter. KgC/m2
     real(r8) :: cwd_ag_local(ncwd)           ! initial value of above ground coarse woody debris. KgC/m2
     real(r8) :: cwd_bg_local(ncwd)           ! initial value of below ground coarse woody debris. KgC/m2
+    integer  :: levcan                       ! canopy level
     real(r8) :: leaf_c                       ! leaf carbon [kg]
-    real(r8) :: fnrt_c               ! fineroot carbon [kg]
-    real(r8) :: sapw_c               ! sapwood carbon [kg]
-    real(r8) :: store_c              ! storage carbon [kg]
-    real(r8) :: struct_c             ! structure carbon [kg]
-    real(r8) :: total_c              ! total carbon of plant [kg]
+    real(r8) :: fnrt_c                       ! fineroot carbon [kg]
+    real(r8) :: sapw_c                       ! sapwood carbon [kg]
+    real(r8) :: store_c                      ! storage carbon [kg]
+    real(r8) :: struct_c                     ! structure carbon [kg]
+    real(r8) :: total_c                      ! total carbon of plant [kg]
 
+    !---------------------------------------------------------------------
 
     storesmallcohort => null() ! storage of the smallest cohort for insertion routine
     storebigcohort   => null() ! storage of the largest cohort for insertion routine 
@@ -449,7 +449,6 @@ subroutine spawn_patches( currentSite, bc_in)
                    nc%cmort = nan     ! The mortality diagnostics are set to nan because the cohort should dissappear
                    nc%hmort = nan
                    nc%bmort = nan
-                   nc%fmort = nan
                    nc%frmort = nan
                    nc%lmort_direct     = nan
                    nc%lmort_collateral = nan
@@ -492,7 +491,6 @@ subroutine spawn_patches( currentSite, bc_in)
                       ! number density in EDCLMLink, and the number density of this new patch is donated
                       ! so with the number density must come the effective mortality rates.
 
-                      nc%fmort            = 0.0_r8               ! Should had also been zero in the donor
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -517,7 +515,6 @@ subroutine spawn_patches( currentSite, bc_in)
                       ! Those remaining in the existing
                       currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)
 
-                      nc%fmort            = 0.0_r8
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -540,11 +537,40 @@ subroutine spawn_patches( currentSite, bc_in)
                 ! loss of individuals from source patch due to area shrinking
                 currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) 
 
+                levcan = currentCohort%canopy_layer 
+
+                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
+               
                 ! loss of individual from fire in new patch.
                 nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) 
 
-                nc%fmort            = currentCohort%fire_mort/hlm_freq_day
-                
                 nc%cmort            = currentCohort%cmort
                 nc%hmort            = currentCohort%hmort
                 nc%bmort            = currentCohort%bmort
@@ -553,7 +579,8 @@ subroutine spawn_patches( currentSite, bc_in)
                 nc%lmort_direct     = currentCohort%lmort_direct
                 nc%lmort_collateral = currentCohort%lmort_collateral
                 nc%lmort_infra      = currentCohort%lmort_infra
-                
+
+
              ! Logging is the dominant disturbance  
              elseif (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. &
                      currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifire)) then  ! Logging 
@@ -574,7 +601,6 @@ subroutine spawn_patches( currentSite, bc_in)
                    nc%hmort            = nan
                    nc%bmort            = nan
                    nc%frmort           = nan
-                   nc%fmort            = nan
                    nc%lmort_direct     = nan
                    nc%lmort_collateral = nan
                    nc%lmort_infra      = nan
@@ -616,7 +642,6 @@ subroutine spawn_patches( currentSite, bc_in)
                       currentCohort%n = currentCohort%n * (1._r8 -  patch_site_areadis/currentPatch%area)
 
 
-                      nc%fmort            = 0.0_r8
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -636,7 +661,6 @@ subroutine spawn_patches( currentSite, bc_in)
                       currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)
 
                       ! No grass impact mortality imposed on the newly created patch
-                      nc%fmort            = 0.0_r8
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort

fire/SFMainMod.F90

@@ -158,14 +158,15 @@ end subroutine fire_danger_index
   subroutine  charecteristics_of_fuel ( currentSite )
   !*****************************************************************
 
-    use SFParamsMod, only  : SF_val_alpha_FMC, SF_val_SAV, SF_val_FBD
+    use SFParamsMod, only  : SF_val_drying_ratio, SF_val_SAV, SF_val_FBD
 
     type(ed_site_type), intent(in), target :: currentSite
 
     type(ed_patch_type),  pointer :: currentPatch
     type(ed_cohort_type), pointer :: currentCohort
 
-    real(r8) timeav_swc 
+    real(r8) timeav_swc
+    real(r8) alpha_FMC(nfsc)     ! Relative fuel moisture adjusted per drying ratio
     real(r8) fuel_moisture(nfsc) ! Scaled moisture content of small litter fuels. 
     real(r8) MEF(nfsc)           ! Moisture extinction factor of fuels     integer n 
 
@@ -235,13 +236,15 @@ subroutine  charecteristics_of_fuel ( currentSite )
           ! Equation 6 in Thonicke et al. 2010. across leaves,twig, small branch, and large branch
           ! dead leaves and twigs included in 1hr pool per Thonicke (2010) 
           ! Calculate fuel moisture for trunks to hold value for fuel consumption
-          fuel_moisture(dl_sf:tr_sf)  = exp(-1.0_r8 * SF_val_alpha_FMC(dl_sf:tr_sf) * currentSite%acc_NI) 
+          alpha_FMC(dl_sf:tr_sf) = SF_val_SAV(dl_sf:tr_sf)/SF_val_drying_ratio
+          
+          fuel_moisture(dl_sf:tr_sf)  = exp(-1.0_r8 * alpha_FMC(dl_sf:tr_sf) * currentSite%acc_NI) 
 
           if(write_SF == itrue)then
              if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac
              if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fm ',fuel_moisture
              if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%acc_NI
-             if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',SF_val_alpha_FMC
+             if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC
           endif
           ! FIX(RF,032414): needs refactoring. 
           ! average water content !is this the correct metric?         

fire/SFParamsMod.F90

@@ -23,7 +23,7 @@ module SFParamsMod
    real(r8),protected :: SF_val_miner_damp
    real(r8),protected :: SF_val_max_durat
    real(r8),protected :: SF_val_durat_slope
-   real(r8),protected :: SF_val_alpha_FMC(NFSC)
+   real(r8),protected :: SF_val_drying_ratio
    real(r8),protected :: SF_val_CWD_frac(NCWD)
    real(r8),protected :: SF_val_max_decomp(NFSC)
    real(r8),protected :: SF_val_SAV(NFSC)
@@ -44,7 +44,7 @@ module SFParamsMod
    character(len=param_string_length),parameter :: SF_name_miner_damp = "fates_miner_damp"
    character(len=param_string_length),parameter :: SF_name_max_durat = "fates_max_durat"
    character(len=param_string_length),parameter :: SF_name_durat_slope = "fates_durat_slope"
-   character(len=param_string_length),parameter :: SF_name_alpha_FMC = "fates_alpha_FMC"
+   character(len=param_string_length),parameter :: SF_name_drying_ratio = "fates_drying_ratio"
    character(len=param_string_length),parameter :: SF_name_CWD_frac = "fates_CWD_frac"
    character(len=param_string_length),parameter :: SF_name_max_decomp = "fates_max_decomp"
    character(len=param_string_length),parameter :: SF_name_SAV = "fates_SAV"
@@ -88,10 +88,10 @@ subroutine SpitFireParamsInit()
     SF_val_miner_damp = nan
     SF_val_max_durat = nan
     SF_val_durat_slope = nan
+    SF_val_drying_ratio = nan
 
     SF_val_CWD_frac(:) = nan
 
-    SF_val_alpha_FMC(:) = nan
     SF_val_max_decomp(:) = nan
 
     SF_val_SAV(:) = nan
@@ -174,6 +174,9 @@ subroutine SpitFireRegisterScalars(fates_params)
     call fates_params%RegisterParameter(name=SF_name_durat_slope, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
+    call fates_params%RegisterParameter(name=SF_name_drying_ratio, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
   end subroutine SpitFireRegisterScalars
 
  !-----------------------------------------------------------------------
@@ -212,6 +215,9 @@ subroutine SpitFireReceiveScalars(fates_params)
     call fates_params%RetreiveParameter(name=SF_name_durat_slope, &
          data=SF_val_durat_slope)
 
+    call fates_params%RetreiveParameter(name=SF_name_drying_ratio, &
+         data=SF_val_drying_ratio)
+
   end subroutine SpitFireReceiveScalars
 
   !-----------------------------------------------------------------------
@@ -279,9 +285,6 @@ subroutine SpitFireRegisterNFSC(fates_params)
     call fates_params%RegisterParameter(name=SF_name_mid_moisture_Slope, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
-    call fates_params%RegisterParameter(name=SF_name_alpha_FMC, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names)
-
     call fates_params%RegisterParameter(name=SF_name_max_decomp, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
@@ -321,9 +324,6 @@ subroutine SpitFireReceiveNFSC(fates_params)
     call fates_params%RetreiveParameter(name=SF_name_mid_moisture_Slope, &
          data=SF_val_mid_moisture_Slope)
 
-    call fates_params%RetreiveParameter(name=SF_name_alpha_FMC, &
-         data=SF_val_alpha_FMC)
-
     call fates_params%RetreiveParameter(name=SF_name_max_decomp, &
          data=SF_val_max_decomp)
 

main/EDInitMod.F90

@@ -70,11 +70,17 @@ subroutine init_site_vars( site_in )
     ! !LOCAL VARIABLES:
     !----------------------------------------------------------------------
     !
-    allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,2))
+    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_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))
+
     !
     end subroutine init_site_vars
 
@@ -123,11 +129,19 @@ 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_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
 
     ! fusoin-induced growth flux of individuals
     site_in%growthflux_fusion(:,:) = 0._r8

main/EDMainMod.F90

@@ -686,7 +686,6 @@ subroutine bypass_dynamics(currentSite)
           currentCohort%bmort = 0.0_r8
           currentCohort%hmort = 0.0_r8
           currentCohort%cmort = 0.0_r8
-          currentCohort%fmort = 0.0_r8
           currentCohort%frmort = 0.0_r8
 
           currentCohort%dndt      = 0.0_r8