Merge branch 'fatesluc_cbasedhrv' into HEAD

biogeochem/EDCanopyStructureMod.F90

@@ -23,6 +23,7 @@ module EDCanopyStructureMod
   use EDTypesMod            , only : nclmax
   use EDTypesMod            , only : nlevleaf
   use EDtypesMod            , only : AREA
+  use EDLoggingMortalityMod , only : UpdateHarvestC
   use FatesGlobals          , only : endrun => fates_endrun
   use FatesInterfaceTypesMod     , only : hlm_days_per_year
   use FatesInterfaceTypesMod     , only : hlm_use_planthydro
@@ -1891,7 +1892,7 @@ end subroutine leaf_area_profile
 
   ! ======================================================================================
 
-  subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
+  subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_in,bc_out)
 
     ! ----------------------------------------------------------------------------------
     ! The purpose of this routine is to package output boundary conditions related
@@ -2110,6 +2111,9 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
           call UpdateH2OVeg(sites(s),bc_out(s),bc_out(s)%plant_stored_h2o_si,1)
        end if
 
+       ! Pass FATES Harvested C to bc_out.
+       call UpdateHarvestC(sites(s),bc_in(s),bc_out(s))
+
     end do
 
     ! This call to RecruitWaterStorage() makes an accounting of
@@ -2123,7 +2127,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

@@ -206,7 +206,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, available_forest_c, harvest_tag)
 
     !
     ! !DESCRIPTION:
@@ -215,14 +216,17 @@ 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(:)
+    real(r8), intent(in) :: available_forest_c(:)
+    integer,  intent(inout) :: harvest_tag(:)
+
     !
     ! !LOCAL VARIABLES:
     real(r8) :: cmort    ! starvation mortality rate (fraction per year)
@@ -233,6 +237,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr
     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
+
     !----------------------------------------------------------------------
 
     ipft = currentCohort%pft
@@ -250,10 +255,11 @@ 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, &
+                               available_forest_c, &
+                               harvest_tag)
 
-
-
 
     if (currentCohort%canopy_layer > 1)then 
        ! Include understory logging mortality rates not associated with disturbance

biogeochem/EDPatchDynamicsMod.F90

@@ -49,6 +49,8 @@ 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 FatesInterfaceTypesMod    , only : hlm_harvest_bypass_criteria
   use FatesGlobals         , only : endrun => fates_endrun
   use FatesConstantsMod    , only : r8 => fates_r8
   use FatesConstantsMod    , only : itrue, ifalse
@@ -59,6 +61,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
@@ -71,6 +75,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
@@ -157,7 +162,7 @@ subroutine disturbance_rates( site_in, bc_in)
     use EDMortalityFunctionsMod , only : mortality_rates
     ! loging flux
     use EDLoggingMortalityMod , only : LoggingMortality_frac
-
+    use EDTypesMod , only : ed_resources_management_type
 
     ! !ARGUMENTS:
     type(ed_site_type) , intent(inout), target :: site_in
@@ -183,8 +188,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) :: harvestable_forest_c(hlm_num_lu_harvest_cats)
+    real(r8) :: available_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)
@@ -193,6 +202,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, available_forest_c)
 
     site_in%harvest_carbon_flux = 0._r8
 
@@ -224,7 +236,10 @@ 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, &
+                available_forest_c, &
+                harvest_tag)
 
           currentCohort%lmort_direct     = lmort_direct
           currentCohort%lmort_collateral = lmort_collateral
@@ -233,12 +248,13 @@ subroutine disturbance_rates( site_in, bc_in)
 
           ! estimate the wood product (trunk_product_site)
           if (currentCohort%canopy_layer>=1) then
+             ! kgC m-2 day-1
              site_in%harvest_carbon_flux = site_in%harvest_carbon_flux + &
                   currentCohort%lmort_direct * currentCohort%n * &
                   ( currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + &
                   currentCohort%prt%GetState(struct_organ, all_carbon_elements)) * &
                   prt_params%allom_agb_frac(currentCohort%pft) * &
-                  SF_val_CWD_frac(ncwd) * logging_export_frac
+                  SF_val_CWD_frac(ncwd) * logging_export_frac * AREA_INV
           endif
 
           currentCohort => currentCohort%taller
@@ -247,6 +263,17 @@ subroutine disturbance_rates( site_in, bc_in)
        currentPatch => currentPatch%younger
     end do
 
+    ! Determine harvest debt from all three categories
+    do h_index = 1, hlm_num_lu_harvest_cats
+       if (harvest_tag(h_index) == 2 .or. &
+           (harvest_tag(h_index) == 1 .and. .not. (hlm_harvest_bypass_criteria))) then
+           if(logging_time) then
+              site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + &
+                  bc_in%hlm_harvest_rates(h_index)
+           end if
+       end if
+    end do
+
     ! ---------------------------------------------------------------------------------------------
     ! Calculate Disturbance Rates based on the mortality rates just calculated
     ! ---------------------------------------------------------------------------------------------
@@ -294,6 +321,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 * &
@@ -304,13 +336,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 loged to account for transfer of interstitial ground area to new secondary lands
+       ! equivalent to the fradction 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, &
+                   available_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

biogeochem/FatesAllometryMod.F90

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

main/EDInitMod.F90

@@ -248,6 +248,7 @@ subroutine zero_site( site_in )
     site_in%promotion_carbonflux = 0._r8
 
     ! Resources management (logging/harvesting, etc)
+    site_in%resources_management%harvest_debt = 0.0_r8
     site_in%resources_management%trunk_product_site  = 0.0_r8
 
     ! canopy spread

main/EDMainMod.F90

@@ -80,6 +80,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 EDPatchDynamicsMod       , only : get_frac_site_primary
   use FatesGlobals             , only : endrun => fates_endrun
   use ChecksBalancesMod        , only : SiteMassStock
@@ -303,6 +304,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 : hlm_use_cohort_age_tracking
     use FatesConstantsMod, only : itrue
     ! !ARGUMENTS:
@@ -333,9 +335,16 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
     !-----------------------------------------------------------------------
     real(r8) :: frac_site_primary
 
+    real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats)
+    real(r8) :: available_forest_c(hlm_num_lu_harvest_cats)
+    integer  :: harvest_tag(hlm_num_lu_harvest_cats)
+
 
     call get_frac_site_primary(currentSite, frac_site_primary)
 
+    ! 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, available_forest_c)
+
     ! Set a pointer to this sites carbon12 mass balance
     site_cmass => currentSite%mass_balance(element_pos(carbon12_element))
 
@@ -368,7 +377,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
           ft = currentCohort%pft
 
           ! 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, available_forest_c, harvest_tag)
 
           ! -----------------------------------------------------------------------------
           ! Apply Plant Allocation and Reactive Transport

main/EDParamsMod.F90

@@ -199,6 +199,9 @@ module EDParamsMod
    real(r8),protected,public :: logging_export_frac        ! "fraction of trunk product being shipped offsite, the 
                                                     ! leftovers will be left onsite as large CWD
    character(len=param_string_length),parameter,public :: logging_name_export_frac ="fates_logging_export_frac"   
+   real(r8),protected,public :: pprodharv10_forest_mean ! "mean harvest mortality proportion of deadstem to 10-yr
+                                                        ! product pool (pprodharv10) of all woody PFT types
+   character(len=param_string_length),parameter,public :: logging_name_pprodharv10="fates_pprodharv10_forest_mean"
 
    real(r8),protected,public :: eca_plant_escalar  ! scaling factor for plant fine root biomass to 
                                                ! calculate nutrient carrier enzyme abundance (ECA)
@@ -263,6 +266,7 @@ subroutine FatesParamsInit()
     logging_event_code                    = nan
     logging_dbhmax_infra                  = nan
     logging_export_frac                   = nan
+    pprodharv10_forest_mean               = nan
     eca_plant_escalar                     = nan
     q10_mr                                = nan
     q10_froz                              = nan
@@ -426,6 +430,9 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=logging_name_export_frac, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
+    call fates_params%RegisterParameter(name=logging_name_pprodharv10, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
     call fates_params%RegisterParameter(name=eca_name_plant_escalar, dimension_shape=dimension_shape_scalar, & 
          dimension_names=dim_names_scalar)
 
@@ -603,6 +610,9 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=logging_name_export_frac, &
           data=logging_export_frac)
 
+    call fates_params%RetreiveParameter(name=logging_name_pprodharv10, &
+          data=pprodharv10_forest_mean)
+
     call fates_params%RetreiveParameter(name=eca_name_plant_escalar, &
           data=eca_plant_escalar)
 

main/EDTypesMod.F90

@@ -593,6 +593,7 @@ 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 
 
      !debug variables
      real(r8) ::  delta_litter_stock                       ! kgC/site = kgC/ha
@@ -1102,6 +1103,9 @@ subroutine dump_cohort(ccohort)
      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
@@ -1113,6 +1117,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