Merge pull request #599 from JessicaNeedham/JessicaNeedham-IO-sizeage…

…mortality Size- and age- dependent mortality
NGEET · Mar 26, 2020 · 3b6ac22 · 3b6ac22
2 parents ab1647a + 10809b9
commit 3b6ac22
Show file tree

Hide file tree

Showing 23 changed files with 1,254 additions and 533 deletions.
diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90
@@ -26,10 +26,11 @@ module EDCanopyStructureMod
   use FatesGlobals          , only : endrun => fates_endrun
   use FatesInterfaceMod     , only : hlm_days_per_year
   use FatesInterfaceMod     , only : hlm_use_planthydro
+  use FatesInterfaceMod     , only : hlm_use_cohort_age_tracking
   use FatesInterfaceMod     , only : numpft
   use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
   use EDTypesMod            , only : maxCohortsPerPatch
-
+  
   use PRTGenericMod,          only : leaf_organ
   use PRTGenericMod,          only : all_carbon_elements
   use PRTGenericMod,          only : leaf_organ
@@ -1256,11 +1257,14 @@ subroutine canopy_summarization( nsites, sites, bc_in )
      ! ---------------------------------------------------------------------------------
 
     use FatesInterfaceMod    , only : bc_in_type
+    use FatesInterfaceMod    , only : hlm_use_cohort_age_tracking
     use EDPatchDynamicsMod   , only : set_patchno
     use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index
+    use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index
     use EDtypesMod           , only : area
     use EDPftvarcon          , only : EDPftvarcon_inst
-
+    use FatesConstantsMod    , only : itrue
+
     ! !ARGUMENTS    
     integer                 , intent(in)            :: nsites
     type(ed_site_type)      , intent(inout), target :: sites(nsites)
@@ -1279,8 +1283,9 @@ subroutine canopy_summarization( nsites, sites, bc_in )
     real(r8) :: sapw_c           ! sapwood carbon [kg]
     real(r8) :: store_c          ! storage carbon [kg]
     real(r8) :: struct_c         ! structure carbon [kg]
-    !----------------------------------------------------------------------
 
+    !----------------------------------------------------------------------
+
     if ( debug ) then
        write(fates_log(),*) 'in canopy_summarization'
     endif
@@ -1319,8 +1324,13 @@ subroutine canopy_summarization( nsites, sites, bc_in )
              ! Update the cohort's index within the size bin classes
              ! Update the cohort's index within the SCPF classification system
              call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
-                                       currentCohort%size_class,currentCohort%size_by_pft_class)
+                  currentCohort%size_class,currentCohort%size_by_pft_class)
 
+             if (hlm_use_cohort_age_tracking .eq. itrue) then
+             call coagetype_class_index(currentCohort%coage,currentCohort%pft, &
+                  currentCohort%coage_class,currentCohort%coage_by_pft_class)
+          end if
+
              call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,&
                   currentCohort%pft,currentCohort%c_area)
 

diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90
diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90
@@ -42,29 +42,36 @@ module EDMortalityFunctionsMod
 
 
 
-  subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
+  subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmort )
 
     ! ============================================================================
     !  Calculate mortality rates from carbon storage, hydraulic cavitation, 
-    !  background and freezing
+    !  background and freezing and size and age dependent senescence
     ! ============================================================================
 
     use FatesConstantsMod,  only : tfrz => t_water_freeze_k_1atm
-
-
+    use FatesInterfaceMod        , only : hlm_hio_ignore_val   
+    use FatesConstantsMod,  only : fates_check_param_set
+
     type (ed_cohort_type), intent(in) :: cohort_in 
     type (bc_in_type), intent(in) :: bc_in
     real(r8),intent(out) :: bmort ! background mortality : Fraction per year
     real(r8),intent(out) :: cmort  ! carbon starvation mortality
     real(r8),intent(out) :: hmort  ! hydraulic failure mortality
     real(r8),intent(out) :: frmort ! freezing stress mortality
+    real(r8),intent(out) :: smort  ! size dependent senescence term
+    real(r8),intent(out) :: asmort ! age dependent senescence term 
 
 
     real(r8) :: frac  ! relativised stored carbohydrate
     real(r8) :: leaf_c_target      ! target leaf biomass kgC
     real(r8) :: store_c
     real(r8) :: hf_sm_threshold    ! hydraulic failure soil moisture threshold 
     real(r8) :: hf_flc_threshold   ! hydraulic failure fractional loss of conductivity threshold
+    real(r8) :: mort_ip_size_senescence ! inflection point for increase in mortality with dbh 
+    real(r8) :: mort_r_size_senescence  ! rate of mortality increase with dbh in senesence term
+    real(r8) :: mort_ip_age_senescence ! inflection point for increase in mortality with age
+    real(r8) :: mort_r_age_senescence ! rate of mortality increase with age in senescence term 
     real(r8) :: temp_dep_fraction  ! Temp. function (freezing mortality)
     real(r8) :: temp_in_C          ! Daily averaged temperature in Celcius
     real(r8) :: min_fmc_ag         ! minimum fraction of maximum conductivity for aboveground
@@ -76,12 +83,45 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     logical, parameter :: test_zero_mortality = .false. ! Developer test which
                                                         ! may help to debug carbon imbalances
                                                         ! and the like
+
+   ! Size Dependent Senescence
+    ! rate (r) and inflection point (ip) define the increase in mortality rate with dbh
+    mort_r_size_senescence = EDPftvarcon_inst%mort_r_size_senescence(cohort_in%pft)
+    mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft)
+
+    ! if param values have been set then calculate smort
+    if ( mort_ip_size_senescence < fates_check_param_set ) then 
+       smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * &
+            (cohort_in%dbh - mort_ip_size_senescence) ) ) 
+    else
+       smort = 0.0_r8
+    end if
 
-    if (hlm_use_ed_prescribed_phys .eq. ifalse) then
+    ! if param values have been set then calculate asmort
+
+
+
+    mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft)
+    mort_ip_age_senescence = EDPftvarcon_inst%mort_ip_age_senescence(cohort_in%pft)
+
+    if ( mort_ip_age_senescence < fates_check_param_set ) then
+       ! Age Dependent Senescence
+       ! rate and inflection point define the change in mortality with age
+
+       asmort = 1.0_r8 / (1.0_r8 + exp(-1.0_r8 * mort_r_age_senescence * &
+            (cohort_in%coage - mort_ip_age_senescence ) ) )
+    else
+       asmort = 0.0_r8
+    end if
+
+
+
+if (hlm_use_ed_prescribed_phys .eq. ifalse) then
 
     ! 'Background' mortality (can vary as a function of 
     !  density as in ED1.0 and ED2.0, but doesn't here for tractability) 
-    bmort = EDPftvarcon_inst%bmort(cohort_in%pft) 
+
+       bmort = EDPftvarcon_inst%bmort(cohort_in%pft)
 
     ! Proxy for hydraulic failure induced mortality. 
     hf_sm_threshold = EDPftvarcon_inst%hf_sm_threshold(cohort_in%pft)
@@ -141,9 +181,10 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
 
     !mortality_rates = bmort + hmort + cmort
 
-    else ! i.e. hlm_use_ed_prescribed_phys is true
+ else ! i.e. hlm_use_ed_prescribed_phys is true 
+
        if ( cohort_in%canopy_layer .eq. 1) then
-          bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft)
+          bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft) 
        else
           bmort = EDPftvarcon_inst%prescribed_mortality_understory(cohort_in%pft)
        endif
@@ -157,6 +198,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
        hmort = 0.0_r8
        frmort = 0.0_r8
        bmort = 0.0_r8
+       smort = 0.0_r8
+       asmort = 0.0_r8
     end if
 
     return
@@ -174,6 +217,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in)
     !
     ! !USES:
 
+    use FatesInterfaceMod, only : hlm_freq_day
     !
     ! !ARGUMENTS    
     type(ed_site_type), intent(inout), target  :: currentSite
@@ -185,6 +229,8 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in)
     real(r8) :: bmort    ! background mortality rate (fraction per year)
     real(r8) :: hmort    ! hydraulic failure mortality rate (fraction per year)
     real(r8) :: frmort   ! freezing mortality rate (fraction per year)
+    real(r8) :: smort    ! size dependent senescence mortality rate (fraction per year)
+    real(r8) :: asmort   ! age dependent senescence mortality rate (fraction per year)
     real(r8) :: dndt_logging      ! Mortality rate (per day) associated with the a logging event
     integer  :: ipft              ! local copy of the pft index
     !----------------------------------------------------------------------
@@ -193,7 +239,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in)
 
     ! Mortality for trees in the understorey. 
     !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology
-    call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort)
+    call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort,smort, asmort)
     call LoggingMortality_frac(ipft, currentCohort%dbh, currentCohort%canopy_layer, &
                                currentCohort%lmort_direct,                       &
                                currentCohort%lmort_collateral,                    &
@@ -206,14 +252,23 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in)
     if (currentCohort%canopy_layer > 1)then 
        ! Include understory logging mortality rates not associated with disturbance
        dndt_logging = (currentCohort%lmort_direct     + &
-                       currentCohort%lmort_collateral + &
-                       currentCohort%lmort_infra)/hlm_freq_day
-       currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort+frmort+dndt_logging) * currentCohort%n
+            currentCohort%lmort_collateral + &
+            currentCohort%lmort_infra)/hlm_freq_day
+
+
+       currentCohort%dndt = -1.0_r8 * &
+            (cmort+hmort+bmort+frmort+smort+asmort + dndt_logging) &
+            * currentCohort%n
     else
+
+
+
        ! Mortality from logging in the canopy is ONLY disturbance generating, don't
        ! update number densities via non-disturbance inducing death
-       currentCohort%dndt = -(1.0_r8 - fates_mortality_disturbance_fraction) &
-            * (cmort+hmort+bmort+frmort) * currentCohort%n
+       currentCohort%dndt= -(1.0_r8-fates_mortality_disturbance_fraction) &
+            * (cmort+hmort+bmort+frmort+smort+asmort) * &
+            currentCohort%n
+
     endif
 
     return

diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90
@@ -158,6 +158,8 @@ subroutine disturbance_rates( site_in, bc_in)
     real(r8) :: bmort
     real(r8) :: hmort
     real(r8) :: frmort
+    real(r8) :: smort
+    real(r8) :: asmort
 
     real(r8) :: lmort_direct
     real(r8) :: lmort_collateral
@@ -181,8 +183,8 @@ subroutine disturbance_rates( site_in, bc_in)
           ! Mortality for trees in the understorey.
           currentCohort%patchptr => currentPatch
 
-          call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort)
-          currentCohort%dmort  = cmort+hmort+bmort+frmort
+          call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort,smort,asmort)
+          currentCohort%dmort  = cmort+hmort+bmort+frmort+smort+asmort
           call carea_allom(currentCohort%dbh,currentCohort%n,site_in%spread,currentCohort%pft, &
                currentCohort%c_area)
 
@@ -191,6 +193,8 @@ subroutine disturbance_rates( site_in, bc_in)
           currentCohort%bmort = bmort
           currentCohort%hmort = hmort
           currentCohort%frmort = frmort
+          currentCohort%smort = smort
+          currentCohort%asmort = asmort
 
           call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_layer, &
                 lmort_direct,lmort_collateral,lmort_infra,l_degrad )
@@ -287,6 +291,8 @@ subroutine disturbance_rates( site_in, bc_in)
                 currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%frmort = currentCohort%frmort*(1.0_r8 - fates_mortality_disturbance_fraction)
+                currentCohort%smort = currentCohort%smort*(1.0_r8 - fates_mortality_disturbance_fraction)
+                currentCohort%asmort = currentCohort%asmort*(1.0_r8 - fates_mortality_disturbance_fraction)
              end if
              currentCohort => currentCohort%taller
           enddo !currentCohort
@@ -307,6 +313,8 @@ subroutine disturbance_rates( site_in, bc_in)
                 currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%frmort = currentCohort%frmort*(1.0_r8 - fates_mortality_disturbance_fraction)
+                currentCohort%smort = currentCohort%smort*(1.0_r8 - fates_mortality_disturbance_fraction)
+                currentCohort%asmort = currentCohort%asmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%lmort_direct    = 0.0_r8
                 currentCohort%lmort_collateral = 0.0_r8
                 currentCohort%lmort_infra      = 0.0_r8
@@ -621,6 +629,8 @@ subroutine spawn_patches( currentSite, bc_in)
                       nc%hmort = nan
                       nc%bmort = nan
                       nc%frmort = nan
+                      nc%smort = nan
+                      nc%asmort = nan
                       nc%lmort_direct     = nan
                       nc%lmort_collateral = nan
                       nc%lmort_infra      = nan
@@ -673,6 +683,8 @@ subroutine spawn_patches( currentSite, bc_in)
                          nc%hmort            = currentCohort%hmort
                          nc%bmort            = currentCohort%bmort
                          nc%frmort           = currentCohort%frmort
+                         nc%smort            = currentCohort%smort
+                         nc%asmort           = currentCohort%asmort
                          nc%dmort            = currentCohort%dmort
                          nc%lmort_direct     = currentCohort%lmort_direct
                          nc%lmort_collateral = currentCohort%lmort_collateral
@@ -697,6 +709,8 @@ subroutine spawn_patches( currentSite, bc_in)
                          nc%hmort            = currentCohort%hmort
                          nc%bmort            = currentCohort%bmort
                          nc%frmort           = currentCohort%frmort
+                         nc%smort            = currentCohort%smort
+                         nc%asmort           = currentCohort%asmort
                          nc%dmort            = currentCohort%dmort
                          nc%lmort_direct    = currentCohort%lmort_direct
                          nc%lmort_collateral = currentCohort%lmort_collateral
@@ -752,6 +766,8 @@ subroutine spawn_patches( currentSite, bc_in)
                    nc%hmort            = currentCohort%hmort
                    nc%bmort            = currentCohort%bmort
                    nc%frmort           = currentCohort%frmort
+                   nc%smort            = currentCohort%smort
+                   nc%asmort           = currentCohort%asmort
                    nc%dmort            = currentCohort%dmort
                    nc%lmort_direct     = currentCohort%lmort_direct
                    nc%lmort_collateral = currentCohort%lmort_collateral
@@ -825,6 +841,8 @@ subroutine spawn_patches( currentSite, bc_in)
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
                       nc%frmort           = currentCohort%frmort
+                      nc%smort            = currentCohort%smort
+                      nc%asmort           = currentCohort%asmort
                       nc%dmort            = currentCohort%dmort
 
                       ! since these are the ones that weren't logged, 
@@ -884,6 +902,8 @@ subroutine spawn_patches( currentSite, bc_in)
                          nc%hmort            = currentCohort%hmort
                          nc%bmort            = currentCohort%bmort
                          nc%frmort           = currentCohort%frmort
+                         nc%smort            = currentCohort%smort
+                         nc%asmort           = currentCohort%asmort
                          nc%dmort            = currentCohort%dmort
                          nc%lmort_direct     = currentCohort%lmort_direct
                          nc%lmort_collateral = currentCohort%lmort_collateral
@@ -904,6 +924,8 @@ subroutine spawn_patches( currentSite, bc_in)
                          nc%hmort            = currentCohort%hmort
                          nc%bmort            = currentCohort%bmort
                          nc%frmort           = currentCohort%frmort
+                         nc%smort            = currentCohort%smort
+                         nc%asmort           = currentCohort%asmort
                          nc%dmort            = currentCohort%dmort
                          nc%lmort_direct     = currentCohort%lmort_direct
                          nc%lmort_collateral = currentCohort%lmort_collateral
@@ -1959,7 +1981,7 @@ subroutine zero_patch(cp_p)
 
 
     ! FIRE
-    currentPatch%litter_moisture(:)         = 0.0_r8 ! litter moisture
+   currentPatch%litter_moisture(:)         = 0.0_r8 ! litter moisture
     currentPatch%fuel_eff_moist             = 0.0_r8 ! average fuel moisture content of the ground fuel 
     ! (incl. live grasses. omits 1000hr fuels)
     currentPatch%livegrass                  = 0.0_r8 ! total ag grass biomass in patch. 1=c3 grass, 2=c4 grass. gc/m2

diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90
@@ -1471,7 +1471,9 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
        temp_cohort%canopy_trim = 0.8_r8  !starting with the canopy not fully expanded 
        temp_cohort%pft         = ft
        temp_cohort%hite        = EDPftvarcon_inst%hgt_min(ft)
+       temp_cohort%coage       = 0.0_r8
        stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(ft)
+
        call h2d_allom(temp_cohort%hite,ft,temp_cohort%dbh)
 
        ! Initialize live pools
@@ -1589,7 +1591,6 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
           ! Initialize the PARTEH object, and determine the initial masses of all
           ! organs and elements.
           ! -----------------------------------------------------------------------------
-
           prt => null()
           call InitPRTObject(prt)
 
@@ -1689,7 +1690,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
 
            ! This initializes the cohort
            call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, & 
-                temp_cohort%hite, temp_cohort%dbh, prt, & 
+                temp_cohort%hite, temp_cohort%coage, temp_cohort%dbh, prt, & 
                 temp_cohort%laimemory, temp_cohort%sapwmemory, temp_cohort%structmemory, &
                 cohortstatus, recruitstatus, &
                 temp_cohort%canopy_trim, currentPatch%NCL_p, currentSite%spread, bc_in)