[Land Use] C-based wood harvest

biogeochem/EDCanopyStructureMod.F90

@@ -1784,7 +1784,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
     use EDTypesMod        , only : ed_patch_type, ed_cohort_type, &
          ed_site_type, AREA
-    use FatesInterfaceTypesMod , only : bc_out_type
+    use FatesInterfaceTypesMod , only : bc_in_type, bc_out_type
 
     !
     ! !ARGUMENTS
@@ -1983,7 +1983,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
        ! Pass FATES Harvested C to bc_out.
        call UpdateHarvestC(sites(s),bc_out(s))
-       
+
     end do
 
     ! This call to RecruitWaterStorage() makes an accounting of
@@ -1997,7 +1997,6 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        call RecruitWaterStorage(nsites,sites,bc_out)
     end if
 
-
   end subroutine update_hlm_dynamics
 
   ! =====================================================================================

biogeochem/EDMortalityFunctionsMod.F90

@@ -225,7 +225,8 @@ end subroutine mortality_rates
 
  ! ============================================================================
 
- subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary)
+ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary, &
+         harvestable_forest_c, harvest_tag)
 
     !
     ! !DESCRIPTION:
@@ -234,14 +235,21 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr
     ! elsewhere).
     !
     ! !USES:
-
     use FatesInterfaceTypesMod, only : hlm_freq_day
     !
     ! !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
     real(r8), intent(in)                       :: frac_site_primary
+
+    real(r8), intent(in) :: harvestable_forest_c(:)   ! total carbon available for logging, kgC site-1
+    integer, intent(out) :: harvest_tag(:)    ! tag to record the harvest status
+                                              ! for the calculation of harvest debt in C-based
+                                              ! harvest mode
+                                              ! 0 - successful;
+                                              ! 1 - unsuccessful since not enough carbon
+                                              ! 2 - not applicable
     !
     ! !LOCAL VARIABLES:
     real(r8) :: cmort    ! starvation mortality rate (fraction per year)
@@ -271,10 +279,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr
                                bc_in%hlm_harvest_units, &
                                currentCohort%patchptr%anthro_disturbance_label, &
                                currentCohort%patchptr%age_since_anthro_disturbance, &
-                               frac_site_primary)
-
-
-
+                               frac_site_primary, harvestable_forest_c, harvest_tag)
 
     if (currentCohort%canopy_layer > 1)then 
        ! Include understory logging mortality rates not associated with disturbance

biogeochem/EDPatchDynamicsMod.F90

@@ -47,6 +47,7 @@ module EDPatchDynamicsMod
   use FatesInterfaceTypesMod    , only : hlm_use_sp
   use FatesInterfaceTypesMod    , only : hlm_use_nocomp
   use FatesInterfaceTypesMod    , only : hlm_use_fixed_biogeog
+  use FatesInterfaceTypesMod    , only : hlm_num_lu_harvest_cats
   use FatesGlobals         , only : endrun => fates_endrun
   use FatesConstantsMod    , only : r8 => fates_r8
   use FatesConstantsMod    , only : itrue, ifalse
@@ -57,6 +58,8 @@ module EDPatchDynamicsMod
   use EDLoggingMortalityMod, only : logging_litter_fluxes 
   use EDLoggingMortalityMod, only : logging_time
   use EDLoggingMortalityMod, only : get_harvest_rate_area
+  use EDLoggingMortalityMod, only : get_harvest_rate_carbon
+  use EDLoggingMortalityMod, only : get_harvestable_carbon
   use EDParamsMod          , only : fates_mortality_disturbance_fraction
   use FatesAllometryMod    , only : carea_allom
   use FatesAllometryMod    , only : set_root_fraction
@@ -69,6 +72,7 @@ module EDPatchDynamicsMod
   use FatesConstantsMod    , only : n_anthro_disturbance_categories
   use FatesConstantsMod    , only : fates_unset_r8
   use FatesConstantsMod    , only : fates_unset_int
+  use FatesConstantsMod    , only : hlm_harvest_carbon
   use EDCohortDynamicsMod  , only : InitPRTObject
   use EDCohortDynamicsMod  , only : InitPRTBoundaryConditions
   use ChecksBalancesMod,      only : SiteMassStock
@@ -187,9 +191,16 @@ subroutine disturbance_rates( site_in, bc_in)
     real(r8) :: dist_rate_ldist_notharvested
     integer  :: threshold_sizeclass
     integer  :: i_dist
+    integer  :: h_index
     real(r8) :: frac_site_primary
     real(r8) :: harvest_rate
     real(r8) :: tempsum
+    real(r8) :: harvest_debt_pri
+    real(r8) :: harvest_debt_sec_mature
+    real(r8) :: harvest_debt_sec_young
+    real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats)
+    integer  :: harvest_tag(hlm_num_lu_harvest_cats)
+    integer  :: patch_no_secondary
 
     !----------------------------------------------------------------------------------------------
     ! Calculate Mortality Rates (these were previously calculated during growth derivatives)
@@ -198,7 +209,16 @@ subroutine disturbance_rates( site_in, bc_in)
 
     ! first calculate the fractino of the site that is primary land
     call get_frac_site_primary(site_in, frac_site_primary)
+
+    ! get available biomass for harvest for all patches
+    call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c)
 
+    ! Initialize local variables
+    patch_no_secondary = 0
+    harvest_debt_pri = 0._r8
+    harvest_debt_sec_mature = 0._r8
+    harvest_debt_sec_young = 0._r8
+
     currentPatch => site_in%oldest_patch
     do while (associated(currentPatch))   
 
@@ -228,7 +248,9 @@ subroutine disturbance_rates( site_in, bc_in)
                 bc_in%hlm_harvest_units, &
                 currentPatch%anthro_disturbance_label, &
                 currentPatch%age_since_anthro_disturbance, &
-                frac_site_primary)
+                frac_site_primary, &
+                harvestable_forest_c, &
+                harvest_tag)
 
           currentCohort%lmort_direct     = lmort_direct
           currentCohort%lmort_collateral = lmort_collateral
@@ -237,9 +259,60 @@ subroutine disturbance_rates( site_in, bc_in)
 
           currentCohort => currentCohort%taller
        end do
+
+       ! ! Counter of secondary patch used in logging debt calculation
+       ! if ( currentPatch%anthro_disturbance_label .eq. secondaryforest ) then
+       !    patch_no_secondary = patch_no_secondary + 1
+       ! end if
+
        currentPatch => currentPatch%younger
     end do
 
+    ! Calculate if we have harvest debt for primary and secondary land
+    ! Harvest debt is the accumulated total carbon 
+    ! deficiency once the carbon amount available for harvest 
+    ! is smaller than the harvest rate of forcing data.
+    ! Harvest debt is calculated on site level
+    ! TODO: we can define harvest debt as a fraction of the 
+    ! harvest rate in the future
+    ! Note: Non-forest harvest is accounted for under forest
+    ! harvest, thus the harvest tag for non-forest is not applicable (= 2)
+
+    if(logging_time) then
+       ! First we need to get harvest rate for all three categories
+       do h_index = 1, hlm_num_lu_harvest_cats
+          ! Primary forest harvest rate
+          if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. &
+              bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2" ) then
+                harvest_debt_pri = harvest_debt_pri + bc_in%hlm_harvest_rates(h_index)
+          else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then
+              harvest_debt_sec_mature = harvest_debt_sec_mature + bc_in%hlm_harvest_rates(h_index)
+          else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. &
+                   bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then
+              harvest_debt_sec_mature = harvest_debt_sec_mature + bc_in%hlm_harvest_rates(h_index)
+          end if
+       end do
+       ! Next we get the harvest debt through the harvest tag 
+       do h_index = 1, hlm_num_lu_harvest_cats
+          if (harvest_tag(h_index) .eq. 1) then
+             if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1") then
+                site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + &
+                    harvest_debt_pri
+             else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then
+                site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + &
+                    harvest_debt_sec_mature
+                site_in%resources_management%harvest_debt_sec = site_in%resources_management%harvest_debt_sec + &
+                    harvest_debt_sec_mature
+             else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2") then
+                site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + &
+                    harvest_debt_sec_young
+                site_in%resources_management%harvest_debt_sec = site_in%resources_management%harvest_debt_sec + &
+                    harvest_debt_sec_young
+             end if
+          end if
+       end do
+    end if
+
     ! ---------------------------------------------------------------------------------------------
     ! Calculate Disturbance Rates based on the mortality rates just calculated
     ! ---------------------------------------------------------------------------------------------
@@ -289,6 +362,11 @@ subroutine disturbance_rates( site_in, bc_in)
                                currentCohort%lmort_infra +                           &
                                currentCohort%l_degrad ) *                            &
                                currentCohort%c_area/currentPatch%area
+
+             if(currentPatch%disturbance_rates(dtype_ilog)>1.0) then
+                 write(fates_log(),*) 'See luc mortalities:', currentCohort%lmort_direct, &
+                     currentCohort%lmort_collateral, currentCohort%lmort_infra, currentCohort%l_degrad
+             end if
 
              ! Non-harvested part of the logging disturbance rate
              dist_rate_ldist_notharvested = dist_rate_ldist_notharvested + currentCohort%l_degrad * &
@@ -299,13 +377,19 @@ subroutine disturbance_rates( site_in, bc_in)
        enddo !currentCohort
 
        ! for non-closed-canopy areas subject to logging, add an additional increment of area disturbed
-       ! equivalent to the fradction logged to account for transfer of interstitial ground area to new secondary lands
+       ! equivalent to the fraction logged to account for transfer of interstitial ground area to new secondary lands
        if ( logging_time .and. &
             (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then
           ! The canopy is NOT closed. 
 
-          call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, &
-               bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate)
+          if(bc_in%hlm_harvest_units == hlm_harvest_carbon) then
+             call get_harvest_rate_carbon (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, &
+                   bc_in%hlm_harvest_rates, currentPatch%age_since_anthro_disturbance, harvestable_forest_c, &
+                   harvest_rate, harvest_tag)
+          else
+             call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, &
+                  bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate)
+          end if
 
           currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + &
                (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area
@@ -315,6 +399,12 @@ subroutine disturbance_rates( site_in, bc_in)
                (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area
        endif
 
+       ! For nocomp mode, we need to prevent producing too small patches, which may produce small patches
+       if ((hlm_use_nocomp .eq. itrue) .and. &
+           (currentPatch%disturbance_rates(dtype_ilog)*currentPatch%area .lt. min_patch_area_forced)) then
+          currentPatch%disturbance_rates(dtype_ilog) = 0._r8
+       end if
+
        ! fraction of the logging disturbance rate that is non-harvested
        if (currentPatch%disturbance_rates(dtype_ilog) .gt. nearzero) then
           currentPatch%fract_ldist_not_harvested = dist_rate_ldist_notharvested / &

biogeochem/EDPhysiologyMod.F90

@@ -1749,7 +1749,6 @@ subroutine SeedIn( currentSite, bc_in )
     ! !USES:
     use EDTypesMod, only : area
     use EDTypesMod, only : homogenize_seed_pfts
-
     !
     ! !ARGUMENTS
     type(ed_site_type), intent(inout), target  :: currentSite

biogeochem/FatesAllometryMod.F90

@@ -2534,7 +2534,6 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, d, h, bdead, bl )
      end if
 
      call h_allom(d,ipft,h)
-
      if(counter>20)then
         write(fates_log(),*) 'dbh counter: ',counter,' is woody: ',&
              (prt_params%woody(ipft) == itrue)

main/EDInitMod.F90

@@ -302,6 +302,8 @@ subroutine zero_site( site_in )
     site_in%fmort_cflux_ustory_damage(:,:) = 0._r8
 
     ! Resources management (logging/harvesting, etc)
+    site_in%resources_management%harvest_debt = 0.0_r8
+    site_in%resources_management%harvest_debt_sec = 0.0_r8
     site_in%resources_management%trunk_product_site  = 0.0_r8
 
     ! canopy spread