C-based harvest in ELM-FATES

components/elm/src/biogeochem/EcosystemDynMod.F90

@@ -805,6 +805,13 @@ subroutine EcosystemDynNoLeaching2(bounds,                                  &
 
        call t_stop_lnd(event)
 
+   else
+       call alm_fates%wrap_WoodProducts(bounds, num_soilc, filter_soilc)
+
+       call WoodProducts(num_soilc, filter_soilc )
+
+       call CropHarvestPools(num_soilc, filter_soilc, dt)
+
    end if
 
    if ( use_erosion ) then

components/elm/src/biogeophys/BalanceCheckMod.F90

@@ -616,7 +616,10 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
        found = .false.
        do p = bounds%begp, bounds%endp
           if (veg_pp%active(p)) then
-             if ( (errsol(p) /= spval) .and. (abs(errsol(p)) > 1.e-7_r8) ) then
+             !if ( (errsol(p) /= spval) .and. (abs(errsol(p)) > 1.e-7_r8) ) then
+             ! solar radiation balance error is high when running FATES
+             ! adjust the threshold to 5.e-7
+             if ( (errsol(p) /= spval) .and. (abs(errsol(p)) > 5.e-7_r8) ) then
                 found = .true.
                 indexp = p
                 indext = veg_pp%topounit(indexp)

components/elm/src/data_types/ColumnDataType.F90

@@ -2054,6 +2054,16 @@ subroutine col_cs_init(this, begc, endc, carbon_type, ratio, c12_carbonstate_var
             avgflag='A', long_name='total column carbon, incl veg and cpool but excl product pools', &
             ptr_col=this%totcolc)
 
+       this%prod10c(begc:endc) = spval
+       call hist_addfld1d (fname='PROD10C', units='gC/m^2', &
+            avgflag='A', long_name='10-yr wood product C', &
+            ptr_col=this%prod10c, default='inactive')
+
+       this%prod100c(begc:endc) = spval
+       call hist_addfld1d (fname='PROD100C', units='gC/m^2', &
+            avgflag='A', long_name='100-yr wood product C', &
+            ptr_col=this%prod100c, default='inactive')
+
        if(.not.use_fates)then
 
           this%seedc(begc:endc) = spval
@@ -5717,6 +5727,34 @@ subroutine col_cf_init(this, begc, endc, carbon_type)
            call hist_addfld1d (fname='HR', units='gC/m^2/s', &
                 avgflag='A', long_name='total heterotrophic respiration', &
                  ptr_col=this%hr)
+! Shijie
+          this%hrv_deadstemc_to_prod10c(begc:endc) = spval
+          call hist_addfld1d (fname='HRV_DEADSTEMC_TO_PROD10C', units='gC/m^2/s', &
+               avgflag='A', long_name='flux into 10-yr wood product C', &
+               ptr_col=this%hrv_deadstemc_to_prod10c, default='inactive')
+
+          this%hrv_deadstemc_to_prod100c(begc:endc) = spval
+          call hist_addfld1d (fname='HRV_DEADSTEMC_TO_PROD100C', units='gC/m^2/s', &
+               avgflag='A', long_name='flux into 100-yr wood product C', &
+               ptr_col=this%hrv_deadstemc_to_prod100c, default='inactive')
+
+          this%nep(begc:endc) = spval
+          call hist_addfld1d (fname='NEP', units='gC/m^2/s', &
+               avgflag='A', long_name='net ecosystem production, excludes fire, landuse, and harvest flux, positive for sink', &
+                ptr_col=this%nep)
+
+          this%nbp(begc:endc) = spval
+          call hist_addfld1d (fname='NBP', units='gC/m^2/s', &
+               avgflag='A', long_name='net biome production, includes fire, landuse, and harvest flux, positive for sink', &
+                ptr_col=this%nbp)
+
+          this%nee(begc:endc) = spval
+          call hist_addfld1d (fname='NEE', units='gC/m^2/s', &
+               avgflag='A', long_name='net ecosystem exchange of carbon, includes fire, landuse,'&
+               //' harvest, and hrv_xsmrpool flux, positive for source', &
+                ptr_col=this%nee)
+! Shijie
+
        end if
        ! end of use_fates (C12) block
 
@@ -6902,9 +6940,13 @@ subroutine col_cf_summary(this, bounds, num_soilc, filter_soilc, isotope)
 
        ! total product loss
        this%product_closs(c) = &
-            this%prod10c_loss(c)  + &
+            this%prod10c_loss(c) + &
             this%prod100c_loss(c) + &
             this%prod1c_loss(c)
+       ! if (masterproc) then
+       !    write(iulog,*) "this%prod10c_loss?", this%prod10c_loss(c), "this%prod100c_loss?", &
+       !        this%prod100c_loss(c), "this%prod1c_loss?", this%prod1c_loss(c)
+       ! endif
 
        ! soil organic matter fire losses (SOMFIRE)
        this%somfire(c) = 0._r8
@@ -7405,18 +7447,18 @@ subroutine col_cf_zero_forfates_veg(this, bounds, num_soilc, filter_soilc)
 
     if(.not.use_fates) return
 
-    do fc = 1,num_soilc
-       c = filter_soilc(fc)
-       this%gpp(c) = 0._r8
-       this%ar(c) = 0._r8
-       this%npp(c) = 0._r8
-       this%vegfire(c) = 0._r8
-       this%wood_harvestc(c) = 0._r8
-       this%fire_closs_p2c(c) = 0._r8
-       !this%litfall(c) = 0._r8 (overwritten)
-       this%hrv_xsmrpool_to_atm(c) = 0._r8
-
-    end do
+     do fc = 1,num_soilc
+        c = filter_soilc(fc)
+        this%gpp(c) = 0._r8
+        this%ar(c) = 0._r8
+        this%npp(c) = 0._r8
+        this%vegfire(c) = 0._r8
+        this%wood_harvestc(c) = 0._r8
+        this%fire_closs_p2c(c) = 0._r8
+        !this%litfall(c) = 0._r8 (overwritten)
+        this%hrv_xsmrpool_to_atm(c) = 0._r8
+
+     end do
 
 
   end subroutine col_cf_zero_forfates_veg

components/elm/src/dyn_subgrid/dynHarvestMod.F90

@@ -25,8 +25,7 @@ module dynHarvestMod
   use ColumnDataType        , only : col_cf, col_nf, col_pf  
   use VegetationType        , only : veg_pp                
   use VegetationDataType    , only : veg_cs, veg_cf, veg_ns, veg_nf  
-  use topounit_varcon      , only : max_topounits
-  use VegetationDataType    , only : veg_ps, veg_pf  
+  use topounit_varcon       , only : max_topounits
   use VegetationDataType    , only : veg_ps, veg_pf
   use elm_varctl            , only : use_cn, use_fates, iulog
   use FatesConstantsMod      , only : hlm_harvest_area_fraction
@@ -66,22 +65,23 @@ module dynHarvestMod
   ! for FATES: capacity for passing harvest data in units of carbon harvested per year (per grid cell) has been added
   ! but these data are not yet included in the input file
   ! the code here can be changed to wood_harvest_units = harvest_carbon to pass carbon data to FATES if:
-  !  the carbon data are in the same input variables as listed below
-  !  and the carbon units in the input file match that expected by FATES
+  ! the carbon data are in the same input variables as listed below
+  ! and the carbon units in the input file match that expected by FATES
 
   integer, public, parameter :: num_harvest_vars = 5
   character(len=64), public, parameter :: harvest_varnames(num_harvest_vars) = &
-       [character(len=64)  :: 'HARVEST_VH1', 'HARVEST_VH2', 'HARVEST_SH1', 'HARVEST_SH2', 'HARVEST_SH3']
+       [character(len=64) :: 'HARVEST_VH1', 'HARVEST_VH2', 'HARVEST_SH1', 'HARVEST_SH2', 'HARVEST_SH3']
+
+  type(dyn_var_time_uninterp_type) :: harvest_vars(num_harvest_vars)   ! value of each harvest variable
+
   ! the units flag must match the units of harvest_varnames
   ! set this here because dynHarvest_init is called after alm_fates%init
   ! this flag is accessed only if namelist do_harvest is TRUE
-  integer, public          :: wood_harvest_units = hlm_harvest_area_fraction
-
-  type(dyn_var_time_uninterp_type) :: harvest_vars(num_harvest_vars)   ! value of each harvest variable
 
-  real(r8) , allocatable, public   :: harvest_rates(:,:) ! category harvest rates (d1) in each gridcell (d2)
+  integer, public, parameter    :: wood_harvest_units = 2    ! 1 = area fraction, 2 = carbon
+  real(r8), allocatable, public :: harvest_rates(:,:) ! harvest rates
+  logical, private              :: do_harvest ! whether we're in a period when we should do harvest
 
-  logical, private         :: do_harvest ! whether we're in a period when we should do harvest
   !---------------------------------------------------------------------------
 
 contains
@@ -92,8 +92,9 @@ subroutine dynHarvest_init(bounds, harvest_filename)
     ! !DESCRIPTION:
     ! Initialize data structures for harvest information.
     ! This should be called once, during model initialization.
-    !
-    use elm_varctl            , only : use_cn
+
+    ! !USES:
+    use elm_varctl            , only : use_cn, use_fates
     use dynVarTimeUninterpMod , only : dyn_var_time_uninterp_type
     use dynTimeInfoMod        , only : YEAR_POSITION_START_OF_TIMESTEP
     use dynTimeInfoMod        , only : YEAR_POSITION_END_OF_TIMESTEP
@@ -156,7 +157,7 @@ subroutine dynHarvest_interp_harvest_types(bounds)
     !
     ! !USES:
     use dynTimeInfoMod , only : time_info_type
-    use elm_varctl            , only : use_cn, use_fates
+    use elm_varctl     , only : use_cn, use_fates
     !
     ! !ARGUMENTS:
     type(bounds_type), intent(in) :: bounds  ! proc-level bounds
@@ -181,7 +182,8 @@ subroutine dynHarvest_interp_harvest_types(bounds)
           ! means that harvest rates will be maintained at the rate given in the last
           ! year of the file for all years past the end of this specified time series.
           do_harvest = .true.
-          allocate(this_data(bounds%begg:bounds%endg))   
+          ! Right now we don't account for the topounit in plant harvest
+          allocate(this_data(bounds%begg:bounds%endg))
           do varnum = 1, num_harvest_vars
              call harvest_vars(varnum)%get_current_data(this_data)
              harvest_rates(varnum,bounds%begg:bounds%endg) = this_data(bounds%begg:bounds%endg)

components/elm/src/main/elmfates_interfaceMod.F90

@@ -152,7 +152,7 @@ module ELMFatesInterfaceMod
    use FatesPlantHydraulicsMod, only : InitHydrSites
    use FatesPlantHydraulicsMod, only : RestartHydrStates
 
-   use dynHarvestMod          , only : num_harvest_vars, harvest_varnames
+   use dynHarvestMod          , only : num_harvest_vars, harvest_varnames, wood_harvest_units
    use dynHarvestMod          , only : harvest_rates ! these are dynamic in space and time
    use dynHarvestMod          , only : num_harvest_vars, harvest_varnames, wood_harvest_units
 
@@ -211,6 +211,7 @@ module ELMFatesInterfaceMod
       procedure, public :: wrap_accumulatefluxes
       procedure, public :: prep_canopyfluxes
       procedure, public :: wrap_canopy_radiation
+      procedure, public :: wrap_WoodProducts
       procedure, public :: wrap_update_hifrq_hist
       procedure, public :: TransferZ0mDisp
       procedure, public :: UpdateLitterFluxes
@@ -271,6 +272,7 @@ subroutine ELMFatesGlobals()
      integer                                        :: pass_biogeog
      integer                                        :: pass_num_lu_harvest_types
      integer                                        :: pass_lu_harvest
+     integer                                        :: pass_harvest_bypass_criteria
      integer                                        :: pass_nocomp
      integer                                        :: pass_sp
 
@@ -294,11 +296,17 @@ subroutine ELMFatesGlobals()
 
      ! We will use this switch temporarily, until  we complete
      ! the ELM-FATES harvest integration
-     logical, parameter :: do_elm_fates_harvest = .false.
-
+     logical, parameter :: do_elm_fates_harvest = .true.
+     ! Another switch to determine if we want to harvest forest 
+     ! by ignoring the harvest criteria
+     ! This is only useful when performing carbon-based harvest
+     ! Need to update as a namelist option in the future
+     logical, parameter :: allow_harvest_bypass_criteria = .false.
+
      if (use_fates) then
 
-verbose_output = .false.
+        ! Shijie: Turn on verbose_output for test
+        verbose_output = .false.
         call FatesInterfaceInit(iulog, verbose_output)
 
         ! Force FATES parameters that are recieve type, to the unset value
@@ -313,7 +321,7 @@ subroutine ELMFatesGlobals()
         call set_fates_ctrlparms('hlm_name',cval='ELM')
         call set_fates_ctrlparms('hio_ignore_val',rval=spval)
         call set_fates_ctrlparms('soilwater_ipedof',ival=get_ipedof(0))
-        call set_fates_ctrlparms('max_patch_per_site',ival=(natpft_size-1))
+        call set_fates_ctrlparms('max_patch_per_site',ival=natpft_size-1)
 
 
         call set_fates_ctrlparms('parteh_mode',ival=fates_parteh_mode)
@@ -411,12 +419,19 @@ subroutine ELMFatesGlobals()
            pass_logging = 1
            pass_num_lu_harvest_types = num_harvest_vars
            pass_lu_harvest = 1
+           if(allow_harvest_bypass_criteria) then
+              pass_harvest_bypass_criteria = 1
+           else
+              pass_harvest_bypass_criteria = 0
+           end if
         else
            pass_lu_harvest = 0
            pass_num_lu_harvest_types = 0
+           pass_harvest_bypass_criteria = 0
         end if
 
         call set_fates_ctrlparms('use_lu_harvest',ival=pass_lu_harvest)
+        call set_fates_ctrlparms('use_harvest_bypass_criteria',ival=pass_harvest_bypass_criteria)
         call set_fates_ctrlparms('num_lu_harvest_cats',ival=pass_num_lu_harvest_types)
         call set_fates_ctrlparms('use_logging',ival=pass_logging)
 
@@ -692,7 +707,6 @@ subroutine init(this, bounds_proc )
 
          end do
 
-
          ! Initialize site-level static quantities dictated by the HLM
          ! currently ground layering depth
 
@@ -886,6 +900,7 @@ subroutine dynamics_driv(this, bounds_clump, top_as_inst,          &
             this%fates(nc)%bc_in(s)%hlm_harvest_catnames = harvest_varnames
             this%fates(nc)%bc_in(s)%hlm_harvest_units = wood_harvest_units
          end if
+         this%fates(nc)%bc_in(s)%site_area=col_pp%wtgcell(c)*grc_pp%area(g)*1e6_r8
 
       end do
 
@@ -1108,10 +1123,11 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
             this%fates(nc)%sites,  &
             this%fates(nc)%bc_in)
 
-       ! Canopy diagnostic outputs for HLM
+       ! Canopy diagnostic outputs for HLM, including LUC
        call update_hlm_dynamics(this%fates(nc)%nsites, &
             this%fates(nc)%sites,  &
             this%f2hmap(nc)%fcolumn, &
+            this%fates(nc)%bc_in, &
             this%fates(nc)%bc_out )
 
        !---------------------------------------------------------------------------------
@@ -1180,7 +1196,6 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
 
           ! initialize SP mode pft order index to 0.  Below ground is the 0th patch
           veg_pp%sp_pftorder_index(col_pp%pfti(c)) = 0
-
           areacheck = veg_pp%wt_ed(col_pp%pfti(c))
 
           do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno
@@ -2213,6 +2228,46 @@ end subroutine wrap_accumulatefluxes
 
  ! ======================================================================================
 
+ subroutine wrap_WoodProducts(this, bounds_clump, fc, filterc)
+
+   ! !ARGUMENTS:
+   class(hlm_fates_interface_type), intent(inout) :: this
+   type(bounds_type)              , intent(in)    :: bounds_clump
+   integer                        , intent(in)    :: fc                   ! size of column filter
+   integer                        , intent(in)    :: filterc(fc)          ! column filter
+
+   ! Locacs
+   integer                                        :: s,c,icc
+   integer                                        :: nc
+
+   associate(&
+         gpp     => col_cf%gpp    , &
+         ar     => col_cf%ar    , &
+         hrv_deadstemc_to_prod10c     => col_cf%hrv_deadstemc_to_prod10c    , &
+         hrv_deadstemc_to_prod100c    => col_cf%hrv_deadstemc_to_prod100c)
+
+    nc = bounds_clump%clump_index
+    ! Loop over columns
+    do icc = 1,fc
+       c = filterc(icc)
+       s = this%f2hmap(nc)%hsites(c)
+
+       ! Shijie: Pass harvested wood products to ELM variable
+       hrv_deadstemc_to_prod10c(c)  = this%fates(nc)%bc_out(s)%hrv_deadstemc_to_prod10c
+       hrv_deadstemc_to_prod100c(c) = this%fates(nc)%bc_out(s)%hrv_deadstemc_to_prod100c
+
+       ! Pass LUC related C fluxes which are calculated in FATES [gC m-2 s-1]
+       gpp(c) = this%fates(nc)%bc_out(s)%gpp_site*1e3
+       ar(c) = this%fates(nc)%bc_out(s)%ar_site*1e3
+
+    end do
+
+    end associate
+    return
+ end subroutine wrap_WoodProducts
+
+ ! ======================================================================================
+
  subroutine wrap_canopy_radiation(this, bounds_clump, &
          num_vegsol, filter_vegsol, coszen, surfalb_inst)