Merge pull request #825 from sshu88/fatesluc

Pass Harvested C from FATES to ELM

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
@@ -2060,6 +2061,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_out(s))
+
     end do
 
     ! This call to RecruitWaterStorage() makes an accounting of

biogeochem/EDLoggingMortalityMod.F90

@@ -54,6 +54,7 @@ module EDLoggingMortalityMod
    use FatesConstantsMod , only : itrue,ifalse
    use FatesGlobals      , only : endrun => fates_endrun 
    use FatesGlobals      , only : fates_log
+   use FatesGlobals      , only : fates_global_verbose
    use shr_log_mod       , only : errMsg => shr_log_errMsg
    use FatesPlantHydraulicsMod, only : AccumulateMortalityWaterStorage
    use PRTGenericMod     , only : all_carbon_elements,carbon12_element
@@ -62,7 +63,7 @@ module EDLoggingMortalityMod
    use FatesAllometryMod , only : set_root_fraction
    use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold
    use FatesConstantsMod , only : fates_tiny
-   use FatesConstantsMod , only : months_per_year
+   use FatesConstantsMod , only : months_per_year, days_per_sec, years_per_day, g_per_kg
    use FatesConstantsMod , only : hlm_harvest_area_fraction
    use FatesConstantsMod , only : hlm_harvest_carbon
    use FatesConstantsMod, only : fates_check_param_set
@@ -85,6 +86,11 @@ module EDLoggingMortalityMod
 
    real(r8), parameter :: harvest_litter_localization = 0.0_r8
 
+   ! ! transfer factor from kg biomass (dry matter) to kg carbon
+   ! ! now we applied a simple fraction of 50% based on the IPCC
+   ! ! guideline
+   ! real(r8), parameter :: carbon_per_kg_biomass = 0.5_r8
+
    character(len=*), parameter, private :: sourcefile = &
          __FILE__
 
@@ -93,6 +99,7 @@ module EDLoggingMortalityMod
    public :: logging_time
    public :: IsItLoggingTime
    public :: get_harvest_rate_area
+   public :: UpdateHarvestC
 
 contains
 
@@ -250,11 +257,26 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
             call get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, &
                  hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate)
 
+            if (fates_global_verbose()) then
+               write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.'
+            end if
+
          else if (hlm_use_lu_harvest == itrue .and. hlm_harvest_units == hlm_harvest_carbon) then
             ! 2=use carbon from hlm
-            ! not implemented yet
+            ! Shijie: Shall call another function, which transfer biomass/carbon into fraction?
+            ! Is it the correct place to call the function?
+            ! Inputs: patch_area, patch_biomass, what else?
+
+            ! call get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_catnames, &
+            !       hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate)
+
+            ! if (fates_global_verbose()) then
+            !    write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.', hlm_harvest_rates(:), harvest_rate 
+            ! end if
+
             write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.'
             call endrun(msg=errMsg(sourcefile, __LINE__))
+
          endif
 
          ! transfer of area to secondary land is based on overall area affected, not just logged crown area
@@ -358,7 +380,7 @@ subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_ca
 
      !  Normalize by site-level primary or secondary forest fraction
      !  since harvest_rate is specified as a fraction of the gridcell
-     ! also need to put a cap so as not to harvest more primary or secondary area than there is in a gridcell
+     !  also need to put a cap so as not to harvest more primary or secondary area than there is in a gridcell
      if (patch_anthro_disturbance_label .eq. primaryforest) then
         if (frac_site_primary .gt. fates_tiny) then
            harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary)
@@ -691,7 +713,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site
             ! This is for checking the total mass balance [kg/site/day]
             site_mass%wood_product = site_mass%wood_product + &
                   ag_wood * logging_export_frac
-            
+
             new_litt%ag_cwd(ncwd) = new_litt%ag_cwd(ncwd) + ag_wood * &
                   (1._r8-logging_export_frac)*donate_frac/newPatch%area
 
@@ -797,4 +819,45 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site
       return
    end subroutine logging_litter_fluxes
 
+  ! =====================================================================================
+
+   subroutine UpdateHarvestC(currentSite,bc_out)
+
+      ! ----------------------------------------------------------------------------------
+      ! Added by Shijie Shu.
+      ! This subroutine is called when logging is completed and need to update 
+      ! Harvested C flux in HLM.
+      ! ----------------------------------------------------------------------------------
+      use EDtypesMod             , only : ed_site_type
+      use EDTypesMod             , only : AREA_INV
+      use PRTGenericMod          , only : element_pos
+      use PRTGenericMod          , only : carbon12_element
+      use FatesInterfaceTypesMod , only : bc_out_type
+      use EDParamsMod            , only : pprodharv10_forest_mean
+
+      ! Arguments
+      type(ed_site_type), intent(inout), target :: currentSite     ! site structure
+      type(bc_out_type), intent(inout)          :: bc_out
+
+      integer :: icode
+      real(r8) :: unit_trans_factor
+
+
+      ! Flush the older value before update
+      bc_out%hrv_deadstemc_to_prod10c = 0._r8
+      bc_out%hrv_deadstemc_to_prod100c = 0._r8
+
+      ! Calculate the unit transfer factor (from kgC m-2 day-1 to gC m-2 s-1)
+      unit_trans_factor = g_per_kg * days_per_sec
+
+      bc_out%hrv_deadstemc_to_prod10c = bc_out%hrv_deadstemc_to_prod10c + &
+          currentSite%mass_balance(element_pos(carbon12_element))%wood_product * &
+          AREA_INV * pprodharv10_forest_mean * unit_trans_factor
+      bc_out%hrv_deadstemc_to_prod100c = bc_out%hrv_deadstemc_to_prod100c + &
+          currentSite%mass_balance(element_pos(carbon12_element))%wood_product * &
+          AREA_INV * (1._r8 - pprodharv10_forest_mean) * unit_trans_factor  
+
+      return
+   end subroutine UpdateHarvestC
+
 end module EDLoggingMortalityMod

biogeochem/EDPatchDynamicsMod.F90

@@ -197,8 +197,6 @@ 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)
 
-    site_in%harvest_carbon_flux = 0._r8
-
     currentPatch => site_in%oldest_patch
     do while (associated(currentPatch))   
 
@@ -234,16 +232,6 @@ subroutine disturbance_rates( site_in, bc_in)
           currentCohort%lmort_infra      = lmort_infra
           currentCohort%l_degrad         = l_degrad
 
-          ! estimate the wood product (trunk_product_site)
-          if (currentCohort%canopy_layer>=1) then
-             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
-          endif
-
           currentCohort => currentCohort%taller
        end do
        currentPatch%disturbance_mode = fates_unset_int
@@ -307,7 +295,7 @@ 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. 

main/EDInitMod.F90

@@ -209,7 +209,6 @@ subroutine zero_site( site_in )
 
     ! Disturbance rates tracking
     site_in%primary_land_patchfusion_error = 0.0_r8
-    site_in%harvest_carbon_flux = 0.0_r8
     site_in%potential_disturbance_rates(:) = 0.0_r8
     site_in%disturbance_rates_secondary_to_secondary(:) = 0.0_r8
     site_in%disturbance_rates_primary_to_secondary(:) = 0.0_r8

main/EDParamsMod.F90

@@ -233,6 +233,10 @@ module EDParamsMod
                                                     ! leftovers will be left onsite as large CWD
    character(len=param_string_length),parameter,public :: logging_name_export_frac ="fates_landuse_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_landuse_pprodharv10_forest_mean"
+
    real(r8),protected,public :: eca_plant_escalar  ! scaling factor for plant fine root biomass to 
                                                ! calculate nutrient carrier enzyme abundance (ECA)
 
@@ -303,6 +307,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
@@ -481,6 +486,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)
 
@@ -685,6 +693,9 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetrieveParameter(name=logging_name_export_frac, &
           data=logging_export_frac)
 
+    call fates_params%RetrieveParameter(name=logging_name_pprodharv10, &
+         data=pprodharv10_forest_mean)
+
     call fates_params%RetrieveParameter(name=eca_name_plant_escalar, &
           data=eca_plant_escalar)
 

main/EDTypesMod.F90

@@ -802,8 +802,6 @@ module EDTypesMod
 
      real(r8) :: demotion_carbonflux                     ! biomass of demoted individuals from canopy to understory [kgC/ha/day]
      real(r8) :: promotion_carbonflux                    ! biomass of promoted individuals from understory to canopy [kgC/ha/day]
-     real(r8) :: harvest_carbon_flux                     ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day]
-
      real(r8) :: recruitment_rate(1:maxpft)              ! number of individuals that were recruited into new cohorts
      real(r8), allocatable :: demotion_rate(:)           ! rate of individuals demoted from canopy to understory per FATES timestep
      real(r8), allocatable :: promotion_rate(:)          ! rate of individuals promoted from understory to canopy per FATES timestep

main/FatesHistoryInterfaceMod.F90

@@ -2146,10 +2146,8 @@ subroutine update_history_dyn(this,nc,nsites,sites)
 
       hio_potential_disturbance_rate_si(io_si) = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) * days_per_year
 
-      ! harvest carbon flux in [kgC/m2/d] -> [kgC/m2/yr]
-      hio_harvest_carbonflux_si(io_si) = sites(s)%harvest_carbon_flux *      &
-         days_per_year
-
+      hio_harvest_carbonflux_si(io_si) = sites(s)%mass_balance(element_pos(carbon12_element))%wood_product * AREA_INV
+
       ! Loop through patches to sum up diagonistics
       ipa = 0
       cpatch => sites(s)%oldest_patch

main/FatesInterfaceMod.F90

@@ -374,7 +374,6 @@ subroutine zero_bcs(fates,s)
     end if
     fates%bc_out(s)%plant_stored_h2o_si = 0.0_r8
 
-
     fates%bc_out(s)%hrv_deadstemc_to_prod10c = 0.0_r8
     fates%bc_out(s)%hrv_deadstemc_to_prod100c = 0.0_r8
 

main/FatesInterfaceTypesMod.F90

@@ -531,7 +531,7 @@ module FatesInterfaceTypesMod
       character(len=64), allocatable :: hlm_harvest_catnames(:)  ! names of hlm_harvest d1
 
       integer :: hlm_harvest_units  ! what units are the harvest rates specified in? [area vs carbon]
-
+    
       ! Fixed biogeography mode 
       real(r8), allocatable :: pft_areafrac(:)     ! Fractional area of the FATES column occupied by each PFT  
 
@@ -722,7 +722,6 @@ module FatesInterfaceTypesMod
                                                        ! small fluxes for various reasons
                                                        ! [mm H2O/s]
 
-
       ! FATES LULCC
       real(r8) :: hrv_deadstemc_to_prod10c   ! Harvested C flux to 10-yr wood product pool [Site-Level, gC m-2 s-1]
       real(r8) :: hrv_deadstemc_to_prod100c  ! Harvested C flux to 100-yr wood product pool [Site-Level, gC m-2 s-1]

main/FatesRestartInterfaceMod.F90

@@ -219,6 +219,7 @@ module FatesRestartInterfaceMod
   integer :: ir_uptake_flxdg
   integer :: ir_oldstock_mbal
   integer :: ir_errfates_mbal
+  integer :: ir_woodprod_mbal
   integer :: ir_prt_base     ! Base index for all PRT variables
 
   ! Hydraulic indices
@@ -1056,7 +1057,11 @@ subroutine define_restart_vars(this, initialize_variables)
          units='kg/ha', veclength=num_elements, flushval = flushzero, &
          hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_errfates_mbal)
 
-
+    call this%RegisterCohortVector(symbol_base='fates_woodproduct', vtype=site_r8, &
+         long_name_base='Current wood product flux', &
+         units='kg/m2/day', veclength=num_elements, flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_woodprod_mbal)
+
     ! Only register satellite phenology related restart variables if it is turned on!
 
     if(hlm_use_sp .eq. itrue) then
@@ -1910,7 +1915,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
 
              this%rvars(ir_oldstock_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%old_stock
              this%rvars(ir_errfates_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%err_fates
-
+             this%rvars(ir_woodprod_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%wood_product
+
           end do
 
 
@@ -2738,7 +2744,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
 
              sites(s)%mass_balance(el)%old_stock = this%rvars(ir_oldstock_mbal+el-1)%r81d(io_idx_si)
              sites(s)%mass_balance(el)%err_fates = this%rvars(ir_errfates_mbal+el-1)%r81d(io_idx_si)
-
+             sites(s)%mass_balance(el)%wood_product = this%rvars(ir_woodprod_mbal+el-1)%r81d(io_idx_si)
           end do
 
           sites(s)%spread = rio_spread_si(io_idx_si)