[WIP] Improved forest management / wood harvest strategies for FATES

biogeochem/EDCohortDynamicsMod.F90

@@ -1124,6 +1124,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                                            nextc%n*nextc%lmort_infra)/newn
                                       currentCohort%l_degrad = (currentCohort%n*currentCohort%l_degrad + &
                                            nextc%n*nextc%l_degrad)/newn
+                                      currentCohort%harv_c = (currentCohort%n*currentCohort%harv_c + &
+                                           nextc%n*nextc%harv_c)/newn
 
                                       ! biomass and dbh tendencies
                                       currentCohort%ddbhdt     = (currentCohort%n*currentCohort%ddbhdt  + &

biogeochem/EDMortalityFunctionsMod.F90

@@ -22,6 +22,7 @@ module EDMortalityFunctionsMod
    use FatesInterfaceTypesMod     , only : hlm_use_tree_damage
    use EDLoggingMortalityMod , only : LoggingMortality_frac
    use EDParamsMod           , only : fates_mortality_disturbance_fraction
+   use EDParamsMod           , only : logging_preference_options
 
    use PRTGenericMod,          only : carbon12_element
    use PRTGenericMod,          only : store_organ
@@ -235,7 +236,8 @@ end subroutine mortality_rates
 
  subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
       mean_temp, land_use_label, age_since_anthro_disturbance,       &
-      frac_site_primary, harvestable_forest_c, harvest_tag)
+      frac_site_primary, frac_site_secondary_mature, harvestable_forest_c, &
+      harvest_rate_scale, harvest_tag)
 
     !
     ! !DESCRIPTION:
@@ -255,8 +257,10 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
     integer,          intent(in)               :: land_use_label
     real(r8),         intent(in)               :: age_since_anthro_disturbance
     real(r8),         intent(in)               :: frac_site_primary
+    real(r8),         intent(in)               :: frac_site_secondary_mature
 
     real(r8), intent(in) :: harvestable_forest_c(:)   ! total carbon available for logging, kgC site-1
+    real(r8), intent(in) :: harvest_rate_scale   ! scale factor for patch-specific harvest rate 
     integer, intent(out) :: harvest_tag(:)    ! tag to record the harvest status
                                               ! for the calculation of harvest debt in C-based
                                               ! harvest mode
@@ -293,7 +297,10 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
                                bc_in%hlm_harvest_units, &
                                land_use_label, &
                                age_since_anthro_disturbance, &
-                               frac_site_primary, harvestable_forest_c, harvest_tag)
+                               frac_site_primary, frac_site_secondary_mature, &
+                               harvestable_forest_c, &
+                               harvest_rate_scale, harvest_tag, &
+                               currentSite%resources_management%harvest_wp_scale, logging_preference_options)
 
     if (currentCohort%canopy_layer > 1)then 
        ! Include understory logging mortality rates not associated with disturbance

biogeochem/EDPhysiologyMod.F90

@@ -488,6 +488,10 @@ subroutine PreDisturbanceLitterFluxes( currentSite, currentPatch, bc_in )
 
        site_mass => currentSite%mass_balance(el)
 
+       site_mass%frag_in = site_mass%frag_in + currentPatch%area * &
+            ( sum(litt%ag_cwd_in) + sum(litt%bg_cwd_in) + &
+            sum(litt%leaf_fines_in) + sum(litt%root_fines_in))
+
        ! Fragmentation flux to soil decomposition model [kg/site/day]
        site_mass%frag_out = site_mass%frag_out + currentPatch%area * &
             ( sum(litt%ag_cwd_frag) + sum(litt%bg_cwd_frag) + &
@@ -2877,7 +2881,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
 
        flux_diags%leaf_litter_input(pft) = &
             flux_diags%leaf_litter_input(pft) +  &
-            leaf_m_turnover * currentCohort%n
+            (leaf_m_turnover+repro_m_turnover) * currentCohort%n
 
        root_fines_tot = (fnrt_m_turnover + store_m_turnover ) * &
             plant_dens
@@ -2898,7 +2902,6 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
             flux_diags%root_litter_input(pft) +  &
             (fnrt_m_turnover + store_m_turnover ) * currentCohort%n
 
-
        ! Assumption: turnover from deadwood and sapwood are lumped together in CWD pool
 
        !update partitioning of stem wood (struct + sapw) to cwd based on cohort dbh
@@ -2926,10 +2929,8 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
 
           flux_diags%cwd_bg_input(c)  = flux_diags%cwd_bg_input(c) + &
                bg_cwd_tot*currentPatch%area
-
        enddo
 
-
        ! ---------------------------------------------------------------------------------
        ! PART 2 Litter fluxes from non-disturbance inducing  mortality. Kg/m2/day
        ! ---------------------------------------------------------------------------------
@@ -2963,8 +2964,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
 
        flux_diags%leaf_litter_input(pft) = &
             flux_diags%leaf_litter_input(pft) +  &
-            leaf_m * dead_n*currentPatch%area
-
+            (leaf_m+repro_m) * dead_n*currentPatch%area
 
        ! %n has not been updated due to mortality yet, thus
        ! the litter flux has already been counted since it captured
@@ -3022,6 +3022,13 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
              site_mass%wood_product = site_mass%wood_product + &
                   trunk_wood * currentPatch%area * logging_export_frac
 
+             site_mass%wood_product_sz(currentCohort%size_class) = &
+                 site_mass%wood_product_sz(currentCohort%size_class) + &
+                  trunk_wood * currentPatch%area * logging_export_frac
+
+             currentCohort%harv_c = currentCohort%harv_c + &
+                   (trunk_wood * currentPatch%area * logging_export_frac) / currentCohort%n
+
              ! Add AG wood to litter from the non-exported fraction of wood
              ! from direct anthro sources
 
@@ -3037,7 +3044,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
                   SF_val_CWD_frac_adj(c) * (dead_n_natural+dead_n_ilogging)  * &
                   prt_params%allom_agb_frac(pft)
 
-             flux_diags%cwd_ag_input(c)  = flux_diags%cwd_ag_input(c) + &
+             flux_diags%cwd_ag_input(c)  = flux_diags%cwd_ag_input(c) + (struct_m + sapw_m) * &
                   SF_val_CWD_frac_adj(c) * (dead_n_natural+dead_n_ilogging) * &
                   currentPatch%area * prt_params%allom_agb_frac(pft)
 
@@ -3050,12 +3057,10 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
              flux_diags%cwd_ag_input(c)  = flux_diags%cwd_ag_input(c) + &
                   SF_val_CWD_frac_adj(c) * dead_n * (struct_m + sapw_m) * &
                   currentPatch%area * prt_params%allom_agb_frac(pft)
-
           end if
 
        end do
 
-
        ! Update diagnostics that track resource management
 
        if( element_id .eq. carbon12_element ) then
@@ -3097,7 +3102,6 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
        currentCohort => currentCohort%taller
     enddo  ! end loop over cohorts
 
-
     return
   end subroutine CWDInput
 

biogeochem/FatesCohortMod.F90

@@ -239,6 +239,7 @@ module FatesCohortMod
     real(r8) :: l_degrad         ! rate of trees that are not killed but suffer from forest degradation
                                  !  (i.e. they are moved to newly-anthro-disturbed secondary 
                                  !  forest patch)  [fraction/logging activity]
+    real(r8) :: harv_c           ! harvested carbon [kgC/indiv]
 
     !---------------------------------------------------------------------------
 
@@ -434,6 +435,7 @@ subroutine NanValues(this)
       this%lmort_collateral        = nan
       this%lmort_infra             = nan
       this%l_degrad                = nan
+      this%harv_c                  = nan
 
       ! GROWTH DERIVATIVES
       this%dndt                    = nan 
@@ -530,6 +532,7 @@ subroutine ZeroValues(this)
       this%lmort_collateral        = 0._r8
       this%lmort_infra             = 0._r8
       this%l_degrad                = 0._r8
+      this%harv_c                  = 0._r8
       this%fraction_crown_burned   = 0._r8
       this%cambial_mort            = 0._r8
       this%crownfire_mort          = 0._r8
@@ -769,6 +772,7 @@ subroutine Copy(this, copyCohort)
       copyCohort%lmort_collateral        = this%lmort_collateral
       copyCohort%lmort_infra             = this%lmort_infra
       copyCohort%l_degrad                = this%l_degrad
+      copyCohort%harv_c                  = this%harv_c
 
       ! GROWTH DERIVATIVES
       copyCohort%dndt                    = this%dndt
@@ -1074,6 +1078,7 @@ subroutine Dump(this)
       write(fates_log(),*) 'cohort%lmort_direct           = ', this%lmort_direct
       write(fates_log(),*) 'cohort%lmort_collateral       = ', this%lmort_collateral
       write(fates_log(),*) 'cohort%lmort_infra            = ', this%lmort_infra
+      write(fates_log(),*) 'cohort%harv_c                 = ', this%harv_c
       write(fates_log(),*) 'cohort%isnew                  = ', this%isnew
       write(fates_log(),*) 'cohort%dndt                   = ', this%dndt
       write(fates_log(),*) 'cohort%dhdt                   = ', this%dhdt