C-based harvest in ELM-FATES

components/elm/bld/namelist_files/namelist_defaults.xml

@@ -125,7 +125,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <!-- ================================================================== -->
 <!-- FATES default parameter file                                       -->
 <!-- ================================================================== -->
-<fates_paramfile >lnd/clm2/paramdata/fates_params_api.24.0.0_12pft_c220719.nc</fates_paramfile>
+<fates_paramfile >lnd/clm2/paramdata/fates_params_api.25.0.0_12pft_c221212.nc</fates_paramfile>
 
 <!-- soil order related parameters (relative to {csmdata}) -->
 <fsoilordercon             >lnd/clm2/paramdata/CNP_parameters_c131108.nc</fsoilordercon>

components/elm/cime_config/testdefs/testmods_dirs/elm/fates/user_nl_elm

@@ -16,7 +16,7 @@ hist_fincl1       = 'FATES_NCOHORTS', 'FATES_TRIMMING', 'FATES_AREA_PLANTS',
 'FATES_FUEL_AMOUNT', 'FATES_LITTER_IN', 'FATES_LITTER_OUT',
 'FATES_SEED_BANK', 'FATES_SEEDS_IN', 'FATES_STOREC', 'FATES_VEGC',
 'FATES_SAPWOODC', 'FATES_LEAFC', 'FATES_FROOTC', 'FATES_REPROC',
-'FATES_CEFFLUX', 'FATES_STRUCTC', 'FATES_NONSTRUCTC', 'FATES_VEGC_ABOVEGROUND',
+'FATES_STRUCTC', 'FATES_NONSTRUCTC', 'FATES_VEGC_ABOVEGROUND',
 'FATES_CANOPY_VEGC', 'FATES_USTORY_VEGC', 'FATES_PRIMARY_PATCHFUSION_ERR',
 'FATES_DISTURBANCE_RATE_P2P', 'FATES_DISTURBANCE_RATE_P2S',
 'FATES_DISTURBANCE_RATE_S2S', 'FATES_DISTURBANCE_RATE_FIRE',

components/elm/cime_config/testdefs/testmods_dirs/elm/fates_eca/user_nl_elm

@@ -18,7 +18,7 @@ hist_fincl1       = 'FATES_NCOHORTS', 'FATES_TRIMMING', 'FATES_AREA_PLANTS',
 'FATES_FUEL_AMOUNT', 'FATES_LITTER_IN', 'FATES_LITTER_OUT',
 'FATES_SEED_BANK', 'FATES_SEEDS_IN', 'FATES_STOREC', 'FATES_VEGC',
 'FATES_SAPWOODC', 'FATES_LEAFC', 'FATES_FROOTC', 'FATES_REPROC',
-'FATES_CEFFLUX', 'FATES_STRUCTC', 'FATES_NONSTRUCTC', 'FATES_VEGC_ABOVEGROUND',
+'FATES_STRUCTC', 'FATES_NONSTRUCTC', 'FATES_VEGC_ABOVEGROUND',
 'FATES_CANOPY_VEGC', 'FATES_USTORY_VEGC', 'FATES_PRIMARY_PATCHFUSION_ERR',
 'FATES_DISTURBANCE_RATE_P2P', 'FATES_DISTURBANCE_RATE_P2S',
 'FATES_DISTURBANCE_RATE_S2S', 'FATES_DISTURBANCE_RATE_FIRE',
@@ -33,7 +33,16 @@ hist_fincl1       = 'FATES_NCOHORTS', 'FATES_TRIMMING', 'FATES_AREA_PLANTS',
 'FATES_CBALANCE_ERROR', 'FATES_ERROR_EL', 'FATES_LEAF_ALLOC',
 'FATES_SEED_ALLOC', 'FATES_STEM_ALLOC', 'FATES_FROOT_ALLOC',
 'FATES_CROOT_ALLOC', 'FATES_STORE_ALLOC','FATES_LITTER_IN_EL','FATES_LITTER_OUT_EL',
-'FATES_LEAFN','FATES_LEAFP','FATES_FROOTN','FATES_FROOTP','FATES_NH4UPTAKE',
-'FATES_NO3UPTAKE','FATES_NEFFLUX','FATES_PUPTAKE','FATES_PEFFLUX',
 'NET_NMIN','NET_PMIN','LITR1P_vr','LITR1C_vr','LITR1N_vr','SOLUTIONP',
-'SMIN_NO3','SMIN_NH4','TOTSOMN','TOTSOMP'
+'SMIN_NO3','SMIN_NH4','TOTSOMN','TOTSOMP',
+'FATES_L2FR','FATES_L2FR_CANOPY_REC_PF','FATES_L2FR_USTORY_REC_PF',
+'FATES_NH4UPTAKE_SZPF','FATES_NO3UPTAKE_SZPF','FATES_NEFFLUX_SZPF'	,
+'FATES_NDEMAND_SZPF','FATES_NFIX_SYM_SZPF','FATES_NH4UPTAKE','FATES_NO3UPTAKE',
+'FATES_NEFFLUX','FATES_NDEMAND','FATES_NFIX_SYM','FATES_STOREN','FATES_STOREN_TF',
+'FATES_VEGN','FATES_SAPWOODN','FATES_LEAFN','FATES_FROOTN','FATES_REPRON','FATES_VEGN_SZPF',
+'FATES_LEAFN_SZPF','FATES_FROOTN_SZPF','FATES_SAPWOODN_SZPF','FATES_STOREN_SZPF','FATES_STOREN_TF_CANOPY_SZPF',
+'FATES_STOREN_TF_USTORY_SZPF','FATES_REPRON_SZPF','FATES_STOREP','FATES_STOREP_TF','FATES_VEGP','FATES_SAPWOODP',
+'FATES_LEAFP','FATES_FROOTP','FATES_REPROP','FATES_PUPTAKE','FATES_PEFFLUX','FATES_PDEMAND',
+'FATES_VEGP_SZPF','FATES_LEAFP_SZPF','FATES_FROOTP_SZPF','FATES_SAPWOODP_SZPF','FATES_STOREP_SZPF',
+'FATES_STOREP_TF_CANOPY_SZPF','FATES_STOREP_TF_USTORY_SZPF','FATES_REPROP_SZPF','FATES_PUPTAKE_SZPF',
+'FATES_PEFFLUX_SZPF','FATES_PDEMAND_SZPF'

components/elm/cime_config/testdefs/testmods_dirs/elm/fates_rd/user_nl_elm

@@ -18,7 +18,7 @@ hist_fincl1       = 'FATES_NCOHORTS', 'FATES_TRIMMING', 'FATES_AREA_PLANTS',
 'FATES_FUEL_AMOUNT', 'FATES_LITTER_IN', 'FATES_LITTER_OUT',
 'FATES_SEED_BANK', 'FATES_SEEDS_IN', 'FATES_STOREC', 'FATES_VEGC',
 'FATES_SAPWOODC', 'FATES_LEAFC', 'FATES_FROOTC', 'FATES_REPROC',
-'FATES_CEFFLUX', 'FATES_STRUCTC', 'FATES_NONSTRUCTC', 'FATES_VEGC_ABOVEGROUND',
+'FATES_STRUCTC', 'FATES_NONSTRUCTC', 'FATES_VEGC_ABOVEGROUND',
 'FATES_CANOPY_VEGC', 'FATES_USTORY_VEGC', 'FATES_PRIMARY_PATCHFUSION_ERR',
 'FATES_DISTURBANCE_RATE_P2P', 'FATES_DISTURBANCE_RATE_P2S',
 'FATES_DISTURBANCE_RATE_S2S', 'FATES_DISTURBANCE_RATE_FIRE',
@@ -33,7 +33,16 @@ hist_fincl1       = 'FATES_NCOHORTS', 'FATES_TRIMMING', 'FATES_AREA_PLANTS',
 'FATES_CBALANCE_ERROR', 'FATES_ERROR_EL', 'FATES_LEAF_ALLOC',
 'FATES_SEED_ALLOC', 'FATES_STEM_ALLOC', 'FATES_FROOT_ALLOC',
 'FATES_CROOT_ALLOC', 'FATES_STORE_ALLOC','FATES_LITTER_IN_EL','FATES_LITTER_OUT_EL',
-'FATES_LEAFN','FATES_LEAFP','FATES_FROOTN','FATES_FROOTP','FATES_NH4UPTAKE',
-'FATES_NO3UPTAKE','FATES_NEFFLUX','FATES_PUPTAKE','FATES_PEFFLUX',
 'NET_NMIN','NET_PMIN','LITR1P_vr','LITR1C_vr','LITR1N_vr','SOLUTIONP',
-'SMIN_NO3','SMIN_NH4','TOTSOMN','TOTSOMP'
+'SMIN_NO3','SMIN_NH4','TOTSOMN','TOTSOMP',
+'FATES_L2FR','FATES_L2FR_CANOPY_REC_PF','FATES_L2FR_USTORY_REC_PF',
+'FATES_NH4UPTAKE_SZPF','FATES_NO3UPTAKE_SZPF','FATES_NEFFLUX_SZPF'	,
+'FATES_NDEMAND_SZPF','FATES_NFIX_SYM_SZPF','FATES_NH4UPTAKE','FATES_NO3UPTAKE',
+'FATES_NEFFLUX','FATES_NDEMAND','FATES_NFIX_SYM','FATES_STOREN','FATES_STOREN_TF',
+'FATES_VEGN','FATES_SAPWOODN','FATES_LEAFN','FATES_FROOTN','FATES_REPRON','FATES_VEGN_SZPF',
+'FATES_LEAFN_SZPF','FATES_FROOTN_SZPF','FATES_SAPWOODN_SZPF','FATES_STOREN_SZPF','FATES_STOREN_TF_CANOPY_SZPF',
+'FATES_STOREN_TF_USTORY_SZPF','FATES_REPRON_SZPF','FATES_STOREP','FATES_STOREP_TF','FATES_VEGP','FATES_SAPWOODP',
+'FATES_LEAFP','FATES_FROOTP','FATES_REPROP','FATES_PUPTAKE','FATES_PEFFLUX','FATES_PDEMAND',
+'FATES_VEGP_SZPF','FATES_LEAFP_SZPF','FATES_FROOTP_SZPF','FATES_SAPWOODP_SZPF','FATES_STOREP_SZPF',
+'FATES_STOREP_TF_CANOPY_SZPF','FATES_STOREP_TF_USTORY_SZPF','FATES_REPROP_SZPF','FATES_PUPTAKE_SZPF',
+'FATES_PEFFLUX_SZPF','FATES_PDEMAND_SZPF'

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 ( ero_ccycle ) then

components/elm/src/biogeophys/BalanceCheckMod.F90

@@ -191,6 +191,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
      integer  :: indexp,indexc,indexl,indext,indexg     ! index of first found in search loop
      real(r8) :: forc_rain_col(bounds%begc:bounds%endc) ! column level rain rate [mm/s]
      real(r8) :: forc_snow_col(bounds%begc:bounds%endc) ! column level snow rate [mm/s]
+     real(r8) :: sol_err_th                             ! solar radiation imbalance threshold
      !-----------------------------------------------------------------------
 
      associate(                                                                         &
@@ -616,7 +617,14 @@ 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
+             ! solar radiation balance error is high when running FATES
+             ! adjust the threshold to 5.e-7
+             if (veg_pp%is_fates(p)) then
+                sol_err_th = 5.e-7_r8
+             else
+                sol_err_th = 1.e-7_r8
+             end if
+             if ( (errsol(p) /= spval) .and. (abs(errsol(p)) > sol_err_th) ) then
                 found = .true.
                 indexp = p
                 indext = veg_pp%topounit(indexp)

components/elm/src/data_types/ColumnDataType.F90

@@ -2090,6 +2090,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
@@ -5763,6 +5773,33 @@ 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)
+
+          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)
+
        end if
        ! end of use_fates (C12) block
 
@@ -6948,9 +6985,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
@@ -7451,18 +7492,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)