Merge pull request #331 from rgknox/rgknox-allocation-integrator

Now that this has passed all tests and we understand all the reasons why it gives answer changes, pulling.  Thanks @rgknox!

Author: ckoven

biogeochem/EDCanopyStructureMod.F90

@@ -4,27 +4,37 @@ module EDMortalityFunctionsMod
   ! Functions that control mortality.
   ! ============================================================================
 
-  use FatesConstantsMod, only : r8 => fates_r8
-  use FatesGlobals     , only : fates_log
-  use EDPftvarcon      , only : EDPftvarcon_inst
-  use EDTypesMod       , only : ed_cohort_type
-  use FatesConstantsMod, only : itrue,ifalse
-  use FatesAllometryMod, only : bleaf
-  use EDParamsMod      , only : ED_val_stress_mort
-  use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys
-  use FatesInterfaceMod, only : bc_in_type
-
-  implicit none
-  private
-
-
-  public ::  mortality_rates
+   use FatesConstantsMod     , only : r8 => fates_r8
+   use FatesGlobals          , only : fates_log
+   use EDPftvarcon           , only : EDPftvarcon_inst
+   use EDTypesMod            , only : ed_cohort_type
+   use EDTypesMod            , only : ed_site_type
+   use EDTypesMod            , only : ed_patch_type
+   use FatesConstantsMod     , only : itrue,ifalse
+   use FatesAllometryMod     , only : bleaf
+   use EDParamsMod           , only : ED_val_stress_mort
+   use FatesInterfaceMod     , only : bc_in_type
+   use FatesInterfaceMod     , only : hlm_use_ed_prescribed_phys
+   use FatesInterfaceMod     , only : hlm_freq_day
+   use EDLoggingMortalityMod , only : LoggingMortality_frac
+   use EDParamsMod           , only : fates_mortality_disturbance_fraction
+   use FatesInterfaceMod     , only : bc_in_type
+
+
+   implicit none
+   private
+   
+   
+   public :: mortality_rates
+   public :: Mortality_Derivative
+   
+   logical :: DEBUG_growth = .false.
+   
+   ! ============================================================================
+   ! 10/30/09: Created by Rosie Fisher
+   ! 02/20/18: Refactored Ryan Knox
+   ! ============================================================================
 
-  logical :: DEBUG_growth = .false.
-
-  ! ============================================================================
-  ! 10/30/09: Created by Rosie Fisher
-  ! ============================================================================
 
 contains
 
@@ -55,6 +65,10 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     real(r8),parameter :: frost_mort_scaler = 3.0_r8  ! Scaling factor for freezing mortality
     real(r8),parameter :: frost_mort_buffer = 5.0_r8  ! 5deg buffer for freezing mortality
 
+    logical, parameter :: test_zero_mortality = .false. ! Developer test which
+                                                        ! may help to debug carbon imbalances
+                                                        ! and the like
+
     if (hlm_use_ed_prescribed_phys .eq. ifalse) then
 
     ! 'Background' mortality (can vary as a function of 
@@ -72,7 +86,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
 
     ! Carbon Starvation induced mortality.
     if ( cohort_in%dbh  >  0._r8 ) then
-       call bleaf(cohort_in%dbh,cohort_in%hite,cohort_in%pft,cohort_in%canopy_trim,b_leaf)
+       call bleaf(cohort_in%dbh,cohort_in%pft,cohort_in%canopy_trim,b_leaf)
        if( b_leaf > 0._r8 .and. cohort_in%bstore <= b_leaf )then
           frac = cohort_in%bstore/ b_leaf
           cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac))
@@ -111,8 +125,68 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
        frmort = 0._r8
     endif
 
+    if (test_zero_mortality) then
+       cmort = 0.0_r8
+       hmort = 0.0_r8
+       frmort = 0.0_r8
+       bmort = 0.0_r8
+    end if
+       
+    return
  end subroutine mortality_rates
 
-! ============================================================================
+ ! ============================================================================
+
+ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in)
+
+    !
+    ! !DESCRIPTION:
+    ! Calculate the change in number density per unit time from the contributing
+    ! rates.  These rates are not disturbance-inducing rates (that is handled
+    ! elsewhere).
+    !
+    ! !USES:
+
+    !
+    ! !ARGUMENTS    
+    type(ed_site_type), intent(inout), target  :: currentSite
+    type(ed_cohort_type),intent(inout), target :: currentCohort
+    type(bc_in_type), intent(in)               :: bc_in
+    !
+    ! !LOCAL VARIABLES:
+    real(r8) :: cmort    ! starvation mortality rate (fraction per year)
+    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) :: dndt_logging      ! Mortality rate (per day) associated with the a logging event
+    integer  :: ipft              ! local copy of the pft index
+    !----------------------------------------------------------------------
+
+    ipft = currentCohort%pft
+    
+    ! 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 LoggingMortality_frac(ipft, currentCohort%dbh, &
+                               currentCohort%lmort_direct,                       &
+                               currentCohort%lmort_collateral,                    &
+                               currentCohort%lmort_infra )
+
+    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
+    else
+       currentCohort%dndt = -(1.0_r8 - fates_mortality_disturbance_fraction) &
+            * (cmort+hmort+bmort+frmort) * currentCohort%n
+    endif
+
+    return
+
+ end subroutine Mortality_Derivative
 
 end module EDMortalityFunctionsMod

biogeochem/EDCohortDynamicsMod.F90

@@ -432,8 +432,8 @@ subroutine spawn_patches( currentSite, bc_in)
                            nc%n * ED_val_understorey_death / hlm_freq_day
                       currentSite%imort_carbonflux = currentSite%imort_carbonflux + &
                            (nc%n * ED_val_understorey_death / hlm_freq_day ) * &
-                           currentCohort%b * g_per_kg * days_per_sec * years_per_day * ha_per_m2
-
+                           currentCohort%b_total() * 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...  
                       nc%n = nc%n * (1.0_r8 - ED_val_understorey_death)
@@ -556,7 +556,8 @@ subroutine spawn_patches( currentSite, bc_in)
                            nc%n * logging_coll_under_frac / hlm_freq_day
                       currentSite%imort_carbonflux = currentSite%imort_carbonflux + &
                            (nc%n * logging_coll_under_frac/ hlm_freq_day ) * &
-                           currentCohort%b * g_per_kg * days_per_sec * years_per_day * ha_per_m2
+                           currentCohort%b_total() * g_per_kg * days_per_sec * years_per_day * ha_per_m2
+
 
                       ! Step 2:  Apply survivor ship function based on the understory death fraction
 
@@ -956,14 +957,14 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si
 
           call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area)
           if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
-             burned_leaves = (currentCohort%bl+currentCohort%bsw) * currentCohort%cfa
+             burned_leaves = min(currentCohort%bl, (currentCohort%bl+currentCohort%bsw) * currentCohort%cfa)
           else
-             burned_leaves = (currentCohort%bl+currentCohort%bsw) * currentPatch%burnt_frac_litter(6)
+             burned_leaves = min(currentCohort%bl, (currentCohort%bl+currentCohort%bsw) * currentPatch%burnt_frac_litter(6))
           endif
           if (burned_leaves > 0.0_r8) then
 
-             currentCohort%balive = max(currentCohort%br,currentCohort%balive - burned_leaves)
-             currentCohort%bl     = max(0.00001_r8,   currentCohort%bl - burned_leaves)
+             currentCohort%bl     = currentCohort%bl - burned_leaves
+
              !KgC/gridcell/day
              currentSite%flux_out = currentSite%flux_out + burned_leaves * currentCohort%n * &
                   patch_site_areadis/currentPatch%area * AREA 
@@ -1337,7 +1338,6 @@ subroutine fuse_patches( csite, bc_in )
     integer  :: fuse_flag   !do patches get fused (1) or not (0). 
     !---------------------------------------------------------------------
 
-    !maxpatch = 4  
     maxpatch = maxPatchesPerSite
 
     currentSite => csite