Add NBP, TOTECOSYSC, and related summary variables for FATES configurations

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

@@ -24,7 +24,7 @@ hist_fincl1       = 'FATES_NCOHORTS', 'FATES_TRIMMING', 'FATES_AREA_PLANTS',
 'FATES_DEMOTION_CARBONFLUX', 'FATES_PROMOTION_CARBONFLUX',
 'FATES_MORTALITY_CFLUX_CANOPY', 'FATES_MORTALITY_CFLUX_USTORY',
 'FATES_NEP', 'FATES_HET_RESP', 'FATES_FIRE_CLOSS', 'FATES_FIRE_FLUX_EL',
-'FATES_CBALANCE_ERROR', 'FATES_ERROR_EL', 'FATES_LEAF_ALLOC',
+'FATES_CBALANCE_ERROR', 'FATES_LEAF_ALLOC',
 'FATES_SEED_ALLOC', 'FATES_STEM_ALLOC', 'FATES_FROOT_ALLOC',
 'FATES_CROOT_ALLOC', 'FATES_STORE_ALLOC',
 'FATES_PATCHAREA_LU', 'FATES_DISTURBANCE_RATE_MATRIX_LULU'

components/elm/cime_config/testdefs/testmods_dirs/elm/fates_cold_twostream/user_nl_elm

@@ -1 +1 @@
-fates_paramfile = '$CASEROOT/fates_params_twostream.nc'
+fates_radiation_model = 'twostream'

components/elm/src/biogeochem/EcosystemBalanceCheckMod.F90

@@ -219,8 +219,9 @@ subroutine ColCBalanceCheck(bounds, &
          col_endcb                 =>    col_cs%endcb                    , & ! Output: [real(r8) (:) ]  carbon mass, end of time step (gC/m**2)
          col_errcb                 =>    col_cs%errcb                    ,  & ! Output: [real(r8) (:) ]  carbon balance error for the timestep (gC/m**2)
          hr                        =>    col_cf%hr                       , &  ! Input: heterotrophic respiration flux (gC/m2/s)
-         litfall                   =>    col_cf%litfall )                     ! Input: carbon flux from litterfall (particularly fates) ( gc/m2/s)
-
+         litfall                   =>    col_cf%litfall                  , &  ! Input: carbon flux from litterfall (particularly fates) ( gc/m2/s)
+         hrv_deadstemc_to_prod10c  =>    col_cf%hrv_deadstemc_to_prod10c , &  ! Input: harvest of wood to 10-year product pool (gC/m2/s)
+         hrv_deadstemc_to_prod100c =>    col_cf%hrv_deadstemc_to_prod100c )   ! Input: harvest of wood to 100-year product pool (gC/m2/s)
 
       ! set time steps
       dt = real( get_step_size(), r8 )
@@ -238,10 +239,10 @@ subroutine ColCBalanceCheck(bounds, &
          ! a change in the total stock. So hwere we assume that
          ! the fates stocks are 0 for simplicity, and are only
          ! concerned that the changes in the soil stocks (including
-         ! fragmented litter, match the fluxes in and out of fates
+         ! fragmented litter) and wood product pools, match the fluxes in and out of fates
          if(use_fates) then
 
-            col_cinputs  = litfall(c)
+            col_cinputs  = litfall(c) + hrv_deadstemc_to_prod10c(c) + hrv_deadstemc_to_prod100c(c)
 
             col_coutputs = er(c)
 
@@ -899,13 +900,15 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
          grc_er                    => grc_cf%er                    , & ! Output: [real(r8) (:) ] (gC/m2/s) total ecosystem respiration, autotrophic + heterotrophic
          grc_fire_closs            => grc_cf%fire_closs            , & ! Output: [real(r8) (:) ] (gC/m2/s) total column-level fire C loss
          grc_prod1c_loss           => grc_cf%prod1_loss            , & ! Output: [real(r8) (:) ] (gC/m2/s) crop leafc harvested
-         grc_prod10c_loss          => grc_cf%prod10_loss           , & ! Output: [real(r8) (:) ] (gC/m2/s) 10-year wood C harvested
-         grc_prod100c_loss         => grc_cf%prod100_loss          , & ! Output: [real(r8) (:) ] (gC/m2/s) 100-year wood C harvested 
+         grc_prod10c_loss          => grc_cf%prod10_loss           , & ! Output: [real(r8) (:) ] (gC/m2/s) 10-year wood C lost
+         grc_prod100c_loss         => grc_cf%prod100_loss          , & ! Output: [real(r8) (:) ] (gC/m2/s) 100-year wood C lost
          grc_hrv_xsmrpool_to_atm   => grc_cf%hrv_xsmrpool_to_atm   , & ! Output: [real(r8) (:) ] (gC/m2/s) excess MR pool harvest mortality
          grc_som_c_leached         => grc_cf%som_c_leached         , & ! Output: [real(r8) (:) ] (gC/m^2/s)total SOM C loss from vertical transport
          grc_som_c_yield           => grc_cf%somc_yield            , & ! Output: [real(r8) (:) ] (gC/m^2/s)total SOM C loss by erosion
          grc_cinputs               => grc_cf%cinputs               , & ! Output: [real(r8) (:) ] (gC/m2/s) column-level C inputs
          grc_coutputs              => grc_cf%coutputs              , & ! Output: [real(r8) (:) ] (gC/m2/s) column-level C outputs
+         grc_hrv_deadstemc_to_prod10c  => grc_cf%hrv_deadstemc_to_prod10c , & ! Output: [real(r8) (:) ] (gC/m2/s) 10-year wood C harvested
+         grc_hrv_deadstemc_to_prod100c => grc_cf%hrv_deadstemc_to_prod100c, & ! Output: [real(r8) (:) ] (gC/m2/s) 100-year wood C harvested
          beg_totpftc               => grc_cs%beg_totpftc          , & ! Input: [real(r8) (:)] (gC/m2) patch-level carbon aggregated to column level, incl veg and cpool
          beg_cwdc                  => grc_cs%beg_cwdc             , & ! Input: [real(r8) (:)] (gC/m2) total column coarse woody debris carbon
          beg_totsomc               => grc_cs%beg_totsomc          , & ! Input: [real(r8) (:)] (gC/m2) total column soil organic matter carbon
@@ -955,7 +958,13 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
 
       if (use_fates) then 
         call c2g(bounds, col_cf%litfall(bounds%begc:bounds%endc), grc_cinputs(bounds%begg:bounds%endg), &
-               c2l_scale_type = 'unity', l2g_scale_type = 'unity')
+             c2l_scale_type = 'unity', l2g_scale_type = 'unity')
+
+        call c2g(bounds, col_cf%hrv_deadstemc_to_prod10c(bounds%begc:bounds%endc), grc_hrv_deadstemc_to_prod10c(bounds%begg:bounds%endg), &
+             c2l_scale_type = 'unity', l2g_scale_type = 'unity')
+
+        call c2g(bounds, col_cf%hrv_deadstemc_to_prod100c(bounds%begc:bounds%endc), grc_hrv_deadstemc_to_prod100c(bounds%begg:bounds%endg), &
+             c2l_scale_type = 'unity', l2g_scale_type = 'unity')
       end if 
 
       dt = real( get_step_size(), r8 )
@@ -965,7 +974,9 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
 
          if (.not. use_fates) then 
            grc_cinputs(g) = grc_gpp(g) + grc_dwt_seedc_to_leaf(g) + grc_dwt_seedc_to_deadstem(g)
-         end if 
+         else
+           grc_cinputs(g) = grc_cinputs(g) + grc_hrv_deadstemc_to_prod10c(g) + grc_hrv_deadstemc_to_prod100c(g)
+         end if
 
          grc_coutputs(g) = grc_er(g) + grc_fire_closs(g) + grc_hrv_xsmrpool_to_atm(g) + &
               grc_prod1c_loss(g) + grc_prod10c_loss(g) + grc_prod100c_loss(g) - grc_som_c_leached(g) + &

components/elm/src/biogeochem/EcosystemDynMod.F90

@@ -274,6 +274,11 @@ subroutine EcosystemDynLeaching(bounds, num_soilc, filter_soilc, &
 
     call t_stop_lnd(event)
 
+    if (use_fates) then
+       call alm_fates%wrap_FatesAtmosphericCarbonFluxes(bounds, num_soilc, filter_soilc)
+       call alm_fates%wrap_FatesCarbonStocks(bounds, num_soilc, filter_soilc)
+    endif
+
   end subroutine EcosystemDynLeaching
 
 

components/elm/src/data_types/ColumnDataType.F90

@@ -6318,6 +6318,11 @@ subroutine col_cf_init(this, begc, endc, carbon_type)
                //' harvest, and hrv_xsmrpool flux, positive for source', &
                 ptr_col=this%nee)
 
+          this%product_closs(begc:endc) = spval
+          call hist_addfld1d (fname='PRODUCT_CLOSS', units='gC/m^2/s', &
+               avgflag='A', long_name='total carbon loss from wood product pools', &
+               ptr_col=this%product_closs, default='inactive')
+
        end if
        ! end of use_fates (C12) block
 

components/elm/src/main/elmfates_interfaceMod.F90

@@ -270,6 +270,8 @@ module ELMFatesInterfaceMod
       procedure, public :: prep_canopyfluxes
       procedure, public :: wrap_canopy_radiation
       procedure, public :: wrap_WoodProducts
+      procedure, public :: wrap_FatesAtmosphericCarbonFluxes
+      procedure, public :: wrap_FatesCarbonStocks
       procedure, public :: wrap_update_hifrq_hist
       procedure, public :: TransferZ0mDisp
       procedure, public :: InterpFileInputs  ! Interpolate inputs from files
@@ -1434,6 +1436,10 @@ subroutine UpdateLitterFluxes(this,bounds_clump)
       do s = 1, this%fates(nc)%nsites
          c = this%f2hmap(nc)%fcolumn(s)
 
+         call FluxIntoLitterPools(this%fates(nc)%sites(s), &
+                                  this%fates(nc)%bc_in(s), &
+                                  this%fates(nc)%bc_out(s))
+
          col_cf%decomp_cpools_sourcesink(c,1:nlevdecomp,i_met_lit) = &
               col_cf%decomp_cpools_sourcesink(c,1:nlevdecomp,i_met_lit) + &
               this%fates(nc)%bc_out(s)%litt_flux_lab_c_si(1:nlevdecomp) * dtime
@@ -1967,13 +1973,6 @@ subroutine restart( this, bounds_proc, ncid, flag, &
                         this%fates(nc)%bc_out(s), &
                         is_restarting = .true.)
 
-                  ! This call sends internal fates variables into the
-                  ! output boundary condition structures. Note: this is called
-                  ! internally in fates dynamics as well.
-
-                  call FluxIntoLitterPools(this%fates(nc)%sites(s), &
-                       this%fates(nc)%bc_in(s), &
-                       this%fates(nc)%bc_out(s))
                end do
 
                if(use_fates_sp)then
@@ -2221,13 +2220,6 @@ subroutine init_coldstart(this, canopystate_inst, soilstate_inst, frictionvel_in
                    this%fates(nc)%bc_out(s), &
                    is_restarting = .false.)
 
-              ! This call sends internal fates variables into the
-              ! output boundary condition structures. Note: this is called
-              ! internally in fates dynamics as well.
-
-              call FluxIntoLitterPools(this%fates(nc)%sites(s), &
-                   this%fates(nc)%bc_in(s), &
-                   this%fates(nc)%bc_out(s))
            end do
 
            ! ------------------------------------------------------------------------
@@ -2752,8 +2744,6 @@ subroutine wrap_WoodProducts(this, bounds_clump, fc, filterc)
    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)
 
@@ -2767,17 +2757,104 @@ subroutine wrap_WoodProducts(this, bounds_clump, fc, filterc)
        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*g_per_kg
-       ar(c) = this%fates(nc)%bc_out(s)%ar_site*g_per_kg
-
     end do
 
     end associate
     return
  end subroutine wrap_WoodProducts
 
  ! ======================================================================================
+
+ subroutine wrap_FatesAtmosphericCarbonFluxes(this, bounds_clump, fc, filterc)
+
+   ! summarize the high-level fluxes that integrate information from both
+   ! FATES and outside-of-FATES decomposition and product decay code.
+
+   use FatesConstantsMod     , only : g_per_kg
+
+   ! !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(&
+        nep     => col_cf%nep    , &
+        nee     => col_cf%nee    , &
+        nbp     => col_cf%nbp    , &
+        product_closs => col_cf%product_closs ,  &
+        hr     => col_cf%hr)
+
+    nc = bounds_clump%clump_index
+    ! Loop over columns
+    do icc = 1,fc
+       c = filterc(icc)
+       s = this%f2hmap(nc)%hsites(c)
+
+       nep(c) = this%fates(nc)%bc_out(s)%gpp_site*g_per_kg &
+            - this%fates(nc)%bc_out(s)%ar_site*g_per_kg &
+            - hr(c)
+
+       nbp(c) = nep(c) &
+            - this%fates(nc)%bc_out(s)%grazing_closs_to_atm_si*g_per_kg &
+            - this%fates(nc)%bc_out(s)%fire_closs_to_atm_si*g_per_kg &
+            - product_closs(c)
+
+       nee(c) = -nbp(c)
+
+    end do
+
+    end associate
+    return
+ end subroutine wrap_FatesAtmosphericCarbonFluxes
+
+ ! ======================================================================================
+
+ subroutine wrap_FatesCarbonStocks(this, bounds_clump, fc, filterc)
+
+   ! summarize the high-level fluxes that integrate information from both
+   ! FATES and outside-of-FATES decomposition and product decay code.
+
+   use FatesConstantsMod     , only : g_per_kg
+
+   ! !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(&
+        totecosysc     => col_cs%totecosysc, &
+        totlitc        => col_cs%totlitc, &
+        totsomc        => col_cs%totsomc, &
+        totprodc       => col_cs%totprodc)
+
+    nc = bounds_clump%clump_index
+    ! Loop over columns
+    do icc = 1,fc
+       c = filterc(icc)
+       s = this%f2hmap(nc)%hsites(c)
+
+       totecosysc(c) = totsomc(c) + totlitc(c) + totprodc(c) + &
+            this%fates(nc)%bc_out(s)%veg_c_si + &
+            this%fates(nc)%bc_out(s)%litter_cwd_c_si + &
+            this%fates(nc)%bc_out(s)%seed_c_si
+
+    end do
+
+    end associate
+    return
+  end subroutine wrap_FatesCarbonStocks
+
+ ! ======================================================================================
 
  subroutine wrap_canopy_radiation(this, bounds_clump, surfalb_inst,nextsw_cday,declinp1)
 
@@ -3260,7 +3337,7 @@ subroutine init_history_io(this,bounds_proc)
    use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8
    use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8
    use FatesIOVariableKindMod, only : site_cdpf_r8, site_cdsc_r8, site_cdam_r8
-   use FatesIOVariableKindMod, only : site_landuse_r8, site_lulu_r8
+   use FatesIOVariableKindMod, only : site_landuse_r8, site_lulu_r8, site_lupft_r8
    use FatesIODimensionsMod, only : fates_bounds_type
 
 
@@ -3362,7 +3439,7 @@ subroutine init_history_io(this,bounds_proc)
              site_scagpft_r8, site_agepft_r8, site_elem_r8, site_elpft_r8, &
              site_elcwd_r8, site_elage_r8, site_coage_r8, site_coage_pft_r8, &
              site_agefuel_r8,site_cdsc_r8, site_cdpf_r8, site_cdam_r8, &
-             site_landuse_r8, site_lulu_r8)
+             site_landuse_r8, site_lulu_r8, site_lupft_r8)
 
            d_index = fates_hist%dim_kinds(dk_index)%dim2_index
            dim2name = fates_hist%dim_bounds(d_index)%name
@@ -3720,6 +3797,9 @@ subroutine hlm_bounds_to_fates_bounds(hlm, fates)
    fates%lulu_begin = 1
    fates%lulu_end   = n_landuse_cats * n_landuse_cats
 
+   fates%lupft_begin = 1
+   fates%lupft_end   = n_landuse_cats * numpft_fates
+
  end subroutine hlm_bounds_to_fates_bounds
 
 ! ======================================================================================

components/elm/src/main/histFileMod.F90

@@ -1958,6 +1958,7 @@ subroutine htape_create (t, histrest)
        call ncd_defdim(lnfid, 'fates_levagefuel', nlevage_fates * nfc_fates, dimid)
        call ncd_defdim(lnfid, 'fates_levlanduse', n_landuse_cats, dimid)
        call ncd_defdim(lnfid, 'fates_levlulu', n_landuse_cats * n_landuse_cats, dimid)
+       call ncd_defdim(lnfid, 'fates_levlupft', n_landuse_cats * numpft_fates, dimid)
     end if
 
     if ( .not. lhistrest )then
@@ -4856,6 +4857,8 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out,
        num2d = n_landuse_cats
     case ('fates_levlulu')
        num2d = n_landuse_cats * n_landuse_cats
+    case ('fates_levlupft')
+       num2d = n_landuse_cats * numpft_fates
     case ('fates_levage')
        num2d = nlevage_fates
     case ('fates_levfuel')