Squashed commit of the following:

commit 194f073
Author: Jessica Needham <[email protected]>
Date:   Tue Dec 17 14:47:53 2019 -0800

    [ Fix bug in error message]

    [Error message if user tries to have age dependent mortality
    and inventory init on. ]

    Fixes: [NGT-ED Github issue #]

    User interface changes?: [ No]

    Code review: [Names]

    Test suite: [suite name, machine, compilers]
    Test baseline:
    Test namelist changes:
    Test answer changes: [bit for bit, roundoff, climate changing]
    Test summary:No testing

commit 34019b7
Author: Jessica Needham <[email protected]>
Date:   Tue Dec 17 12:09:06 2019 -0800

    [ Update FatesPFTIndexswapper to be compatible with python3]

    [Change iteritems to items ]

    Fixes: [NGT-ED Github issue #]

    User interface changes?: [Yes (describe what changes), No]

    Code review: [Names]

    Test suite: [suite name, machine, compilers]
    Test baseline:
    Test namelist changes:
    Test answer changes: [bit for bit, roundoff, climate changing]
    Test summary:Doesn't change fates. No testing.

commit ed84697
Author: Jessica Needham <[email protected]>
Date:   Tue Dec 17 11:44:47 2019 -0800

    [Size- and age-dependent mortality functionality ]

    [These changes introduce two new mortality terms: M9 -
    size dependent and M10 age dependent mortality. They are
    logistic functions that gradually increase the rate of
    mortality with size or age. Each is controlled by two
    parameters a rate and an inflection point. To turn off
    size dependent mortality the inflection point can be set
    to a very large number (>10000). To turn off age mortality
    set the fates_cohort_age_fusion_tol parameter to 0. Cohort
    age will simply stay at 0. Otherwise, cohort age is tracked.
    Cohorts cannot be fused if they are too different in age. ]

    Fixes: [NGT-ED Github issue #]

    User interface changes?: [Yes - five new parameters. User has to
    specify age and size mortality as either on or off by setting
    fates_mort_ip_size_senescence and fates_cohort_age_fusion_tol. If
    they decide to have size or age on, then sensible parameter values
    need to be selected for inflection point and rate.]

    Code review: [Names]

    Test suite: [suite name, machine, compilers]
    Test baseline:
    Test namelist changes:
    Test answer changes: [bit for bit, roundoff, climate changing]
    Test summary: No testing yet.

Conflicts:
	main/EDParamsMod.F90
	parameter_files/fates_params_default.cdl

biogeochem/EDCanopyStructureMod.F90

@@ -29,7 +29,8 @@ module EDCanopyStructureMod
   use FatesInterfaceMod     , only : numpft
   use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
   use EDTypesMod            , only : maxCohortsPerPatch
-
+  use EDParamsMod           , only : ED_val_cohort_age_fusion_tol
+
   use PRTGenericMod,          only : leaf_organ
   use PRTGenericMod,          only : all_carbon_elements
   use PRTGenericMod,          only : leaf_organ
@@ -1258,9 +1259,11 @@ subroutine canopy_summarization( nsites, sites, bc_in )
     use FatesInterfaceMod    , only : bc_in_type
     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 EDParamsMod          , only : ED_val_cohort_age_fusion_tol
+
     ! !ARGUMENTS    
     integer                 , intent(in)            :: nsites
     type(ed_site_type)      , intent(inout), target :: sites(nsites)
@@ -1279,8 +1282,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 +1323,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 (ED_val_cohort_age_fusion_tol > 0.0_r8) 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)
 

biogeochem/EDMortalityFunctionsMod.F90

@@ -19,6 +19,7 @@ module EDMortalityFunctionsMod
    use FatesInterfaceMod     , only : hlm_use_planthydro
    use EDLoggingMortalityMod , only : LoggingMortality_frac
    use EDParamsMod           , only : fates_mortality_disturbance_fraction
+   use EDParamsMod           , only : ED_val_cohort_age_fusion_tol
    use FatesInterfaceMod     , only : bc_in_type
 
    use PRTGenericMod,          only : all_carbon_elements
@@ -42,11 +43,11 @@ 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 dependent senescence
     ! ============================================================================
 
     use FatesConstantsMod,  only : tfrz => t_water_freeze_k_1atm
@@ -58,12 +59,18 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     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
@@ -75,12 +82,42 @@ 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 the user has set the inflection point to be a huge number then
+    ! we zero smort (even though it would essentially be zero anyway)
+    if (mort_ip_size_senescence < 10000 ) 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 user has set cohort age fusion param > 0 we calculate age
+    ! dependent mortality
+    if (ED_val_cohort_age_fusion_tol > 0.0_r8) then
+       ! Age Dependent Senescence
+       ! rate and inflection point define the change in mortality with age
+       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)
+       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
+
+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 +178,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 +195,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 +214,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 +226,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 +236,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 +249,36 @@ 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
+
+       ! this  caps bmort so that daily mortality cannot exceed 0.995
+       if ((cmort+hmort+bmort+frmort+smort+asmort + dndt_logging)*hlm_freq_day &
+            > 0.995_r8)then
+          bmort = (0.995_r8 - ((cmort+hmort+asmort+frmort+smort+dndt_logging)*hlm_freq_day)) &
+               /hlm_freq_day
+       endif
+
+       currentCohort%dndt = -1.0_r8 * &
+            (cmort+hmort+bmort+frmort+smort+asmort + dndt_logging) &
+            * currentCohort%n
     else
+
+        ! cap on bmort for canopy layers - bmort is adjusted so that total mortality
+       ! including disturbance fraction does not exceed 0.995
+       if(((cmort+hmort+bmort+frmort+smort+asmort)*hlm_freq_day) > &
+            0.995_r8 - fates_mortality_disturbance_fraction)then
+          bmort = (0.995_r8 - fates_mortality_disturbance_fraction - &
+               ((cmort+hmort+asmort+frmort+smort)*hlm_freq_day))/hlm_freq_day
+       endif
+
+
        ! 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+dndt_logging) * &
+            currentCohort%n
+
     endif
 
     return

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

biogeochem/EDPhysiologyMod.F90

@@ -1361,6 +1361,7 @@ 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
        call h2d_allom(temp_cohort%hite,ft,temp_cohort%dbh)
 
        ! Initialize live pools
@@ -1460,7 +1461,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)
 
@@ -1560,7 +1560,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, cohortstatus, recruitstatus, &
                 temp_cohort%canopy_trim, currentPatch%NCL_p, currentSite%spread, bc_in)