Merge pull request #888 from sshu88/fatesluc_cbasedhrv_v2

[Land Use] C-based wood harvest

biogeochem/EDCanopyStructureMod.F90

@@ -1995,7 +1995,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
@@ -2009,7 +2009,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

@@ -48,6 +48,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
@@ -58,6 +59,9 @@ 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 EDLoggingMortalityMod, only : get_harvest_debt
   use EDParamsMod          , only : fates_mortality_disturbance_fraction
   use FatesAllometryMod    , only : carea_allom
   use FatesAllometryMod    , only : set_root_fraction
@@ -70,6 +74,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 +192,12 @@ 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) :: harvestable_forest_c(hlm_num_lu_harvest_cats)
+    integer  :: harvest_tag(hlm_num_lu_harvest_cats)
 
     !----------------------------------------------------------------------------------------------
     ! Calculate Mortality Rates (these were previously calculated during growth derivatives)
@@ -198,6 +206,9 @@ 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)
 
     currentPatch => site_in%oldest_patch
     do while (associated(currentPatch))   
@@ -228,7 +239,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 +250,12 @@ subroutine disturbance_rates( site_in, bc_in)
 
           currentCohort => currentCohort%taller
        end do
+
        currentPatch => currentPatch%younger
     end do
 
+    call get_harvest_debt(site_in, bc_in, harvest_tag)
+
     ! ---------------------------------------------------------------------------------------------
     ! Calculate Disturbance Rates based on the mortality rates just calculated
     ! ---------------------------------------------------------------------------------------------
@@ -289,6 +305,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 +320,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 +342,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

@@ -1753,7 +1753,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

@@ -2520,7 +2520,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

@@ -303,6 +303,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

main/EDMainMod.F90

@@ -87,6 +87,7 @@ module EDMainMod
   use FatesPlantHydraulicsMod  , only : AccumulateMortalityWaterStorage
   use FatesAllometryMod        , only : h_allom,tree_sai,tree_lai
   use EDLoggingMortalityMod    , only : IsItLoggingTime
+  use EDLoggingMortalityMod    , only : get_harvestable_carbon
   use DamageMainMod            , only : IsItDamageTime
   use EDPatchDynamicsMod       , only : get_frac_site_primary
   use FatesGlobals             , only : endrun => fates_endrun
@@ -322,6 +323,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
     ! FIX(SPM,032414) refactor so everything goes through interface
     !
     ! !USES:
+    use FatesInterfaceTypesMod, only : hlm_num_lu_harvest_cats
     use FatesInterfaceTypesMod, only : nlevdamage
     use FatesAllometryMod    , only : bleaf
     use FatesAllometryMod    , only : carea_allom
@@ -381,6 +383,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
     real(r8) :: target_leaf_c
     real(r8) :: frac_site_primary
 
+    real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats)
+    integer  :: harvest_tag(hlm_num_lu_harvest_cats)
+
     real(r8) :: n_old
     real(r8) :: n_recover
     real(r8) :: sapw_c
@@ -414,6 +419,13 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
 
     call get_frac_site_primary(currentSite, frac_site_primary)
 
+    ! Clear site GPP and AR passing to HLM
+    bc_out%gpp_site = 0._r8
+    bc_out%ar_site = 0._r8
+
+    ! Patch level biomass are required for C-based harvest
+    call get_harvestable_carbon(currentSite, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c)
+
     ! Set a pointer to this sites carbon12 mass balance
     site_cmass => currentSite%mass_balance(element_pos(carbon12_element))
 
@@ -456,7 +468,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
        do while(associated(currentCohort))
 
           ft = currentCohort%pft
-
           ! Some cohorts are created and inserted to the list while
           ! the loop is going. These are pointed to the "taller" position
           ! of current, and then inherit properties of their donor (current)
@@ -465,8 +476,10 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
 
           if_not_newlyrecovered: if(.not.newly_recovered) then
 
+
              ! Calculate the mortality derivatives
-             call Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary )
+             call Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary, &
+                 harvestable_forest_c, harvest_tag )
 
              ! -----------------------------------------------------------------------------
              ! Apply Plant Allocation and Reactive Transport
@@ -477,6 +490,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
              !    decrement the available carbon pool to zero.
              ! -----------------------------------------------------------------------------
 
+
              if (hlm_use_ed_prescribed_phys .eq. itrue) then
                 if (currentCohort%canopy_layer .eq. 1) then
                    currentCohort%npp_acc = EDPftvarcon_inst%prescribed_npp_canopy(ft) &
@@ -485,14 +499,14 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
                    currentCohort%npp_acc = EDPftvarcon_inst%prescribed_npp_understory(ft) &
                         * currentCohort%c_area / currentCohort%n / hlm_days_per_year
                 endif
-                
+
                 ! We don't explicitly define a respiration rate for prescribe phys
                 ! but we do need to pass mass balance. So we say it is zero respiration
                 currentCohort%gpp_acc  = currentCohort%npp_acc
                 currentCohort%resp_acc = 0._r8
-                
+
              end if
-             
+
              ! -----------------------------------------------------------------------------
              ! Save NPP/GPP/R in these "hold" style variables. These variables
              ! persist after this routine is complete, and used in I/O diagnostics.
@@ -504,21 +518,26 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
              ! <x>_acc will be reset soon and will be accumulated on the next leaf
              !         photosynthesis step
              ! -----------------------------------------------------------------------------
-             
+
              currentCohort%npp_acc_hold  = currentCohort%npp_acc  * real(hlm_days_per_year,r8)
              currentCohort%gpp_acc_hold  = currentCohort%gpp_acc  * real(hlm_days_per_year,r8)
              currentCohort%resp_acc_hold = currentCohort%resp_acc * real(hlm_days_per_year,r8)
-
+
+             ! Passing gpp_acc_hold to HLM 
+             bc_out%gpp_site = bc_out%gpp_site + currentCohort%gpp_acc_hold * &
+                  AREA_INV * currentCohort%n / hlm_days_per_year / sec_per_day
+             bc_out%ar_site = bc_out%ar_site + currentCohort%resp_acc_hold * & 
+                  AREA_INV * currentCohort%n / hlm_days_per_year / sec_per_day
+
              ! Conduct Maintenance Turnover (parteh)
              if(debug) call currentCohort%prt%CheckMassConservation(ft,3)
              if(any(currentSite%dstatus == [phen_dstat_moiston,phen_dstat_timeon])) then
                 is_drought = .false.
              else
                 is_drought = .true.
              end if
-
              call PRTMaintTurnover(currentCohort%prt,ft,is_drought)
-
+             
              ! -----------------------------------------------------------------------------------
              ! Call the routine that advances leaves in age.
              ! This will move a portion of the leaf mass in each
@@ -538,7 +557,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
              currentCohort%resp_excess = 0._r8
 
           end if if_not_newlyrecovered
-             
+
           ! If the current diameter of a plant is somehow less than what is consistent
           ! with what is allometrically consistent with the stuctural biomass, then
           ! correct the dbh to match.

main/EDTypesMod.F90

@@ -613,6 +613,9 @@ module EDTypesMod
   type, public :: ed_resources_management_type
 
      real(r8) ::  trunk_product_site                       ! Actual  trunk product at site level KgC/site
+     real(r8) ::  harvest_debt                             ! the amount of kgC per site that did not successfully harvested 
+     real(r8) ::  harvest_debt_sec                         ! the amount of kgC per site from secondary patches that did
+                                                           ! not successfully harvested
 
      !debug variables
      real(r8) ::  delta_litter_stock                       ! kgC/site = kgC/ha
@@ -781,7 +784,6 @@ module EDTypesMod
                                                                ! in runs that are restarted, regardless of
                                                                ! the conditions of restart
 
-
      real(r8) ::  water_memory(numWaterMem)                             ! last 10 days of soil moisture memory...
 
 
@@ -1166,6 +1168,10 @@ subroutine dump_cohort(ccohort)
      write(fates_log(),*) 'co%dgmort                 = ', ccohort%dgmort
      write(fates_log(),*) 'co%hmort                  = ', ccohort%hmort
      write(fates_log(),*) 'co%frmort                 = ', ccohort%frmort
+     write(fates_log(),*) 'co%asmort                 = ', ccohort%asmort
+     write(fates_log(),*) 'co%lmort_direct           = ', ccohort%lmort_direct
+     write(fates_log(),*) 'co%lmort_collateral       = ', ccohort%lmort_collateral
+     write(fates_log(),*) 'co%lmort_infra            = ', ccohort%lmort_infra
      write(fates_log(),*) 'co%isnew                  = ', ccohort%isnew
      write(fates_log(),*) 'co%dndt                   = ', ccohort%dndt
      write(fates_log(),*) 'co%dhdt                   = ', ccohort%dhdt
@@ -1177,6 +1183,7 @@ subroutine dump_cohort(ccohort)
      write(fates_log(),*) 'co%cambial_mort           = ', ccohort%cambial_mort
      write(fates_log(),*) 'co%size_class             = ', ccohort%size_class
      write(fates_log(),*) 'co%size_by_pft_class      = ', ccohort%size_by_pft_class
+
      if (associated(ccohort%co_hydr) ) then
         call dump_cohort_hydr(ccohort)
      endif