diff --git a/models/clm/DART_SourceMods/cesm2_2_0/SourceMods/src.clm/biogeochem/CNVegetationFacade.F90 b/models/clm/DART_SourceMods/cesm2_2_0/SourceMods/src.clm/biogeochem/CNVegetationFacade.F90
new file mode 100644
index 000000000..691ba5023
--- /dev/null
+++ b/models/clm/DART_SourceMods/cesm2_2_0/SourceMods/src.clm/biogeochem/CNVegetationFacade.F90
@@ -0,0 +1,1552 @@
+module CNVegetationFacade
+
+  !-----------------------------------------------------------------------
+  ! !DESCRIPTION:
+  ! Facade for the CN Vegetation subsystem.
+  !
+  ! (A "facade", in software engineering terms, is a unified interface to a set of
+  ! interfaces in a subsystem. The facade defines a higher-level interface that makes the
+  ! subsystem easier to use.)
+  !
+  ! NOTE(wjs, 2016-02-19) I envision that we will introduce an abstract base class
+  ! (VegBase). Then both CNVeg and EDVeg will extend VegBase. The rest of the CLM code can
+  ! then have an instance of VegBase, which depending on the run, can be either a CNVeg or
+  ! EDVeg instance.
+  !
+  ! In addition, we probably want an implementation when running without CN or fates - i.e.,
+  ! an SPVeg inst. This would provide implementations for get_leafn_patch,
+  ! get_downreg_patch, etc., so that we don't need to handle the non-cn case here (note
+  ! that, currently, we return NaN for most of these getters, because these arrays are
+  ! invalid and shouldn't be used when running in SP mode). Also, in its EcosystemDynamics
+  ! routine, it would call SatellitePhenology (but note that the desired interface for
+  ! EcosystemDynamics would be quite different... could just pass everything needed by any
+  ! model, and ignore unneeded arguments). Then we can get rid of comments in this module
+  ! like, "only call if use_cn is true", as well as use_cn conditionals in this module.
+  !
+  ! NOTE(wjs, 2016-02-23) Currently, SatellitePhenology is called even when running with
+  ! CN, for the sake of dry deposition. This seems weird to me, and my gut feeling -
+  ! without understanding it well - is that this should be rewritten to depend on LAI from
+  ! CN rather than from satellite phenology. Until that is done, the separation between SP
+  ! and other Veg modes will be messier.
+  !
+  ! NOTE(wjs, 2016-02-23) Currently, this class coordinates calls to soil BGC routines as
+  ! well as veg BGC routines (even though it doesn't contain any soil BGC types). This is
+  ! because CNDriver coordinates both the veg & soil BGC. We should probably split up
+  ! CNDriver so that there is a cleaner separation between veg BGC and soil BGC, to allow
+  ! easier swapping of (for example) CN and ED. At that point, this class could
+  ! coordinate just the calls to veg BGC routines, with a similar facade class
+  ! coordinating the calls to soil BGC routines.
+  !
+  ! !USES:
+#include "shr_assert.h"
+  use shr_kind_mod                    , only : r8 => shr_kind_r8
+  use shr_infnan_mod                  , only : nan => shr_infnan_nan, assignment(=)
+  use shr_log_mod                     , only : errMsg => shr_log_errMsg
+  use perf_mod                        , only : t_startf, t_stopf
+  use decompMod                       , only : bounds_type
+  use clm_varctl                      , only : iulog, use_cn, use_cndv, use_c13, use_c14
+  use abortutils                      , only : endrun
+  use spmdMod                         , only : masterproc
+  use clm_time_manager                , only : get_curr_date, get_ref_date
+  use clm_time_manager                , only : get_nstep, is_end_curr_year, is_first_step
+  use CNBalanceCheckMod               , only : cn_balance_type
+  use CNVegStateType                  , only : cnveg_state_type
+  use CNVegCarbonFluxType             , only : cnveg_carbonflux_type
+  use CNVegCarbonStateType            , only : cnveg_carbonstate_type
+  use CNVegNitrogenFluxType           , only : cnveg_nitrogenflux_type
+  use CNVegNitrogenStateType          , only : cnveg_nitrogenstate_type
+  use FireMethodType                  , only : fire_method_type
+  use CNProductsMod                   , only : cn_products_type
+  use NutrientCompetitionMethodMod    , only : nutrient_competition_method_type
+  use SpeciesIsotopeType              , only : species_isotope_type
+  use SpeciesNonIsotopeType           , only : species_non_isotope_type
+  use CanopyStateType                 , only : canopystate_type
+  use PhotosynthesisMod               , only : photosyns_type
+  use atm2lndType                     , only : atm2lnd_type
+  use WaterStateBulkType                  , only : waterstatebulk_type
+  use WaterDiagnosticBulkType                  , only : waterdiagnosticbulk_type
+  use WaterFluxBulkType                   , only : waterfluxbulk_type
+  use Wateratm2lndBulkType                   , only : wateratm2lndbulk_type
+  use SoilStateType                   , only : soilstate_type
+  use TemperatureType                 , only : temperature_type 
+  use CropType                        , only : crop_type
+  use ch4Mod                          , only : ch4_type
+  use CNDVType                        , only : dgvs_type
+  use CNDVDriverMod                   , only : CNDVDriver, CNDVHIST
+  use EnergyFluxType                  , only : energyflux_type
+  use SaturatedExcessRunoffMod        , only : saturated_excess_runoff_type
+  use FrictionVelocityMod             , only : frictionvel_type
+  use ActiveLayerMod                  , only : active_layer_type
+  use SoilBiogeochemStateType         , only : soilBiogeochem_state_type
+  use SoilBiogeochemCarbonStateType   , only : soilbiogeochem_carbonstate_type
+  use SoilBiogeochemCarbonFluxType    , only : soilBiogeochem_carbonflux_type
+  use SoilBiogeochemNitrogenStateType , only : soilbiogeochem_nitrogenstate_type
+  use SoilBiogeochemNitrogenFluxType  , only : soilbiogeochem_nitrogenflux_type
+  use CNFireEmissionsMod              , only : fireemis_type, CNFireEmisUpdate
+  use CNDriverMod                     , only : CNDriverInit
+  use CNDriverMod                     , only : CNDriverSummarizeStates, CNDriverSummarizeFluxes
+  use CNDriverMod                     , only : CNDriverNoLeaching, CNDriverLeaching
+  use CNCStateUpdate1Mod              , only : CStateUpdateDynPatch
+  use CNNStateUpdate1Mod              , only : NStateUpdateDynPatch
+  use CNVegStructUpdateMod            , only : CNVegStructUpdate
+  use CNAnnualUpdateMod               , only : CNAnnualUpdate
+  use dynConsBiogeochemMod            , only : dyn_cnbal_patch, dyn_cnbal_col
+  use dynCNDVMod                      , only : dynCNDV_init, dynCNDV_interp
+  use CNPrecisionControlMod           , only: CNPrecisionControl
+  use SoilBiogeochemPrecisionControlMod , only: SoilBiogeochemPrecisionControl
+  !
+  implicit none
+  private
+
+  ! !PUBLIC TYPES:
+
+  type, public :: cn_vegetation_type
+     ! FIXME(bja, 2016-06) These need to be public for use when fates is
+     ! turned on. Should either be moved out of here or create some ED
+     ! version of the facade....
+     type(cnveg_state_type)         :: cnveg_state_inst
+     type(cnveg_carbonstate_type)   :: cnveg_carbonstate_inst
+     type(cnveg_carbonflux_type)    :: cnveg_carbonflux_inst
+
+     !X!private
+
+     type(cnveg_carbonstate_type)   :: c13_cnveg_carbonstate_inst
+     type(cnveg_carbonstate_type)   :: c14_cnveg_carbonstate_inst
+     type(cnveg_carbonflux_type)    :: c13_cnveg_carbonflux_inst
+     type(cnveg_carbonflux_type)    :: c14_cnveg_carbonflux_inst
+     type(cnveg_nitrogenstate_type) :: cnveg_nitrogenstate_inst
+     type(cnveg_nitrogenflux_type)  :: cnveg_nitrogenflux_inst
+
+     type(cn_products_type)         :: c_products_inst
+     type(cn_products_type)         :: c13_products_inst
+     type(cn_products_type)         :: c14_products_inst
+     type(cn_products_type)         :: n_products_inst
+
+     type(cn_balance_type)          :: cn_balance_inst
+     class(fire_method_type), allocatable :: cnfire_method
+     type(dgvs_type)                :: dgvs_inst
+
+     ! Control variables
+     logical, private :: reseed_dead_plants              ! Flag to indicate if should reseed dead plants when starting up the model
+     logical, private :: dribble_crophrv_xsmrpool_2atm = .False. ! Flag to indicate if should harvest xsmrpool to the atmosphere
+
+     ! TODO(wjs, 2016-02-19) Evaluate whether some other variables should be moved in
+     ! here. Whether they should be moved in depends on how tightly they are tied in with
+     ! the other CN Vegetation stuff. A question to ask is: Is this module used when
+     ! running with SP or ED? If so, then it should probably remain outside of CNVeg.
+     !
+     ! From the clm_instMod section on "CN vegetation types":
+     ! - nutrient_competition_method
+     !   - I'm pretty sure this should be moved into here; it's just a little messy to do
+     !     so, because of how it's initialized (specifically, the call to readParameters
+     !     in clm_initializeMod).
+     !
+     ! From the clm_instMod section on "general biogeochem types":
+     ! - ch4_inst
+     !   - probably not: really seems to belong in soilbiogeochem
+     ! - crop_inst
+     ! - dust_inst
+     ! - vocemis_inst
+     ! - fireemis_inst
+     ! - drydepvel_inst
+     
+   contains
+     procedure, public :: Init
+     procedure, public :: InitAccBuffer
+     procedure, public :: InitAccVars
+     procedure, public :: UpdateAccVars
+     procedure, public :: Restart
+
+     procedure, public :: Init2                         ! Do initialization in initialize phase, after subgrid weights are determined
+     procedure, public :: InitEachTimeStep              ! Do initializations at the start of each time step
+     procedure, public :: InterpFileInputs              ! Interpolate inputs from files
+     procedure, public :: UpdateSubgridWeights          ! Update subgrid weights if running with prognostic patch weights
+     procedure, public :: DynamicAreaConservation       ! Conserve C & N with updates in subgrid weights
+     procedure, public :: InitColumnBalance             ! Set the starting point for col-level balance checks
+     procedure, public :: InitGridcellBalance           ! Set the starting point for gridcell-level balance checks
+     procedure, public :: EcosystemDynamicsPreDrainage  ! Do the main science that needs to be done before hydrology-drainage
+     procedure, public :: EcosystemDynamicsPostDrainage ! Do the main science that needs to be done after hydrology-drainage
+     procedure, public :: BalanceCheck                  ! Check the carbon and nitrogen balance
+     procedure, public :: EndOfTimeStepVegDynamics      ! Do vegetation dynamics that should be done at the end of each time step
+     procedure, public :: WriteHistory                  ! Do any history writes that are specific to veg dynamics
+
+     procedure, public :: get_net_carbon_exchange_grc   ! Get gridcell-level net carbon exchange array
+     procedure, public :: get_leafn_patch               ! Get patch-level leaf nitrogen array
+     procedure, public :: get_downreg_patch             ! Get patch-level downregulation array
+     procedure, public :: get_root_respiration_patch    ! Get patch-level root respiration array
+     procedure, public :: get_annsum_npp_patch          ! Get patch-level annual sum NPP array
+     procedure, public :: get_agnpp_patch               ! Get patch-level aboveground NPP array
+     procedure, public :: get_bgnpp_patch               ! Get patch-level belowground NPP array
+     procedure, public :: get_froot_carbon_patch        ! Get patch-level fine root carbon array
+     procedure, public :: get_croot_carbon_patch        ! Get patch-level coarse root carbon array
+     procedure, public :: get_totvegc_col               ! Get column-level total vegetation carbon array
+
+     procedure, private :: CNReadNML                    ! Read in the CN general namelist
+  end type cn_vegetation_type
+
+  ! !PRIVATE DATA MEMBERS:
+
+  integer, private :: skip_steps    ! Number of steps to skip at startup
+  character(len=*), parameter, private :: sourcefile = &
+       __FILE__
+
+contains
+
+  !-----------------------------------------------------------------------
+  subroutine Init(this, bounds, NLFilename, nskip_steps, params_ncid)
+    !
+    ! !DESCRIPTION:
+    ! Initialize a CNVeg object.
+    !
+    ! Should be called regardless of whether use_cn is true
+    !
+    ! !USES:
+    use CNFireFactoryMod , only : create_cnfire_method
+    use clm_varcon       , only : c13ratio, c14ratio
+    use ncdio_pio        , only : file_desc_t
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type), intent(in)    :: bounds
+    character(len=*) , intent(in)    :: NLFilename  ! namelist filename
+    integer          , intent(in)    :: nskip_steps ! Number of steps to skip at startup
+    type(file_desc_t), intent(inout) :: params_ncid ! NetCDF handle to parameter file
+    !
+    ! !LOCAL VARIABLES:
+    integer :: begp, endp
+
+    character(len=*), parameter :: subname = 'Init'
+    !-----------------------------------------------------------------------
+
+    begp = bounds%begp
+    endp = bounds%endp
+
+    ! Note - always initialize the memory for cnveg_state_inst (used in biogeophys/)
+    call this%cnveg_state_inst%Init(bounds)
+    
+    skip_steps = nskip_steps
+
+    if (use_cn) then
+
+       ! Read in the general CN namelist
+       call this%CNReadNML( NLFilename )    ! MUST be called first as passes down control information to others
+
+       call this%cnveg_carbonstate_inst%Init(bounds, carbon_type='c12', ratio=1._r8, &
+                                             NLFilename=NLFilename,  dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm )
+       if (use_c13) then
+          call this%c13_cnveg_carbonstate_inst%Init(bounds, carbon_type='c13', ratio=c13ratio, &
+               NLFilename=NLFilename, dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm,        &
+               c12_cnveg_carbonstate_inst=this%cnveg_carbonstate_inst)
+       end if
+       if (use_c14) then
+          call this%c14_cnveg_carbonstate_inst%Init(bounds, carbon_type='c14', ratio=c14ratio, &
+               NLFilename=NLFilename, dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm,        &
+               c12_cnveg_carbonstate_inst=this%cnveg_carbonstate_inst)
+       end if
+       call this%cnveg_carbonflux_inst%Init(bounds, carbon_type='c12', dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm )
+       if (use_c13) then
+          call this%c13_cnveg_carbonflux_inst%Init(bounds, carbon_type='c13', dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm)
+       end if
+       if (use_c14) then
+          call this%c14_cnveg_carbonflux_inst%Init(bounds, carbon_type='c14', dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm)
+       end if
+       call this%cnveg_nitrogenstate_inst%Init(bounds,                   &
+            this%cnveg_carbonstate_inst%leafc_patch(begp:endp),          &
+            this%cnveg_carbonstate_inst%leafc_storage_patch(begp:endp),  &
+            this%cnveg_carbonstate_inst%frootc_patch(begp:endp),         &
+            this%cnveg_carbonstate_inst%frootc_storage_patch(begp:endp), &
+            this%cnveg_carbonstate_inst%deadstemc_patch(begp:endp) )
+       call this%cnveg_nitrogenflux_inst%Init(bounds) 
+
+       call this%c_products_inst%Init(bounds, species_non_isotope_type('C'))
+       if (use_c13) then
+          call this%c13_products_inst%Init(bounds, species_isotope_type('C', '13'))
+       end if
+       if (use_c14) then
+          call this%c14_products_inst%Init(bounds, species_isotope_type('C', '14'))
+       end if
+       call this%n_products_inst%Init(bounds, species_non_isotope_type('N'))
+
+       call this%cn_balance_inst%Init(bounds)
+
+       ! Initialize the memory for the dgvs_inst data structure regardless of whether
+       ! use_cndv is true so that it can be used in associate statements (nag compiler
+       ! complains otherwise)
+       call this%dgvs_inst%Init(bounds)
+    end if
+
+    call create_cnfire_method(NLFilename, this%cnfire_method)
+    call this%cnfire_method%CNFireReadParams( params_ncid )
+
+  end subroutine Init
+
+  !-----------------------------------------------------------------------
+  subroutine CNReadNML( this, NLFilename )
+    !
+    ! !DESCRIPTION:
+    ! Read in the general CN control namelist
+    !
+    ! !USES:
+    use fileutils      , only : getavu, relavu, opnfil
+    use shr_nl_mod     , only : shr_nl_find_group_name
+    use spmdMod        , only : masterproc, mpicom
+    use shr_mpi_mod    , only : shr_mpi_bcast
+    use clm_varctl     , only : iulog
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    character(len=*)         , intent(in)    :: NLFilename                 ! Namelist filename
+    !
+    ! !LOCAL VARIABLES:
+    integer :: ierr                 ! error code
+    integer :: unitn                ! unit for namelist file
+
+    character(len=*), parameter :: subname = 'CNReadNML'
+    character(len=*), parameter :: nmlname = 'cn_general'   ! MUST match what is in namelist below
+    !-----------------------------------------------------------------------
+    logical :: reseed_dead_plants
+    logical :: dribble_crophrv_xsmrpool_2atm
+    namelist /cn_general/ reseed_dead_plants, dribble_crophrv_xsmrpool_2atm
+
+    reseed_dead_plants    = this%reseed_dead_plants
+    dribble_crophrv_xsmrpool_2atm = this%dribble_crophrv_xsmrpool_2atm
+
+    if (masterproc) then
+       unitn = getavu()
+       write(iulog,*) 'Read in '//nmlname//'  namelist'
+       call opnfil (NLFilename, unitn, 'F')
+       call shr_nl_find_group_name(unitn, nmlname, status=ierr)
+       if (ierr == 0) then
+          read(unitn, nml=cn_general, iostat=ierr)   ! Namelist name here MUST be the same as in nmlname above!
+          if (ierr /= 0) then
+             call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__))
+          end if
+       else
+          call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__))
+       end if
+       call relavu( unitn )
+    end if
+
+    call shr_mpi_bcast (reseed_dead_plants     , mpicom)
+    call shr_mpi_bcast (dribble_crophrv_xsmrpool_2atm  , mpicom)
+
+    this%reseed_dead_plants = reseed_dead_plants
+    this%dribble_crophrv_xsmrpool_2atm = dribble_crophrv_xsmrpool_2atm
+
+    if (masterproc) then
+       write(iulog,*) ' '
+       write(iulog,*) nmlname//' settings:'
+       write(iulog,nml=cn_general)    ! Name here MUST be the same as in nmlname above!
+       write(iulog,*) ' '
+    end if
+
+    !-----------------------------------------------------------------------
+
+  end subroutine CNReadNML
+
+
+  !-----------------------------------------------------------------------
+  subroutine InitAccBuffer(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Initialize accumulation buffer for types contained here 
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type), intent(in)    :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitAccBuffer'
+    !-----------------------------------------------------------------------
+
+    if (use_cndv) then
+       call this%dgvs_inst%InitAccBuffer(bounds)
+    end if
+
+  end subroutine InitAccBuffer
+
+  !-----------------------------------------------------------------------
+  subroutine InitAccVars(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Initialize variables that are associated with accumulated fields
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type), intent(in)    :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitAccVars'
+    !-----------------------------------------------------------------------
+
+    if (use_cndv) then
+       call this%dgvs_inst%initAccVars(bounds)
+    end if
+
+  end subroutine InitAccVars
+
+  !-----------------------------------------------------------------------
+  subroutine UpdateAccVars(this, bounds, t_a10_patch, t_ref2m_patch)
+    !
+    ! !DESCRIPTION:
+    ! Update accumulated variables
+    !
+    ! Should be called every time step
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type), intent(in)    :: bounds
+    ! NOTE(wjs, 2016-02-23) These need to be pointers to agree with the interface of
+    ! UpdateAccVars in CNDVType (they are pointers there as a workaround for a compiler
+    ! bug).
+    real(r8), pointer , intent(in)   :: t_a10_patch(:)      ! 10-day running mean of the 2 m temperature (K)
+    real(r8), pointer , intent(in)   :: t_ref2m_patch(:)    ! 2 m height surface air temperature (K)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'UpdateAccVars'
+    !-----------------------------------------------------------------------
+
+    SHR_ASSERT_ALL_FL((ubound(t_a10_patch) == (/bounds%endp/)), sourcefile, __LINE__)
+    SHR_ASSERT_ALL_FL((ubound(t_ref2m_patch) == (/bounds%endp/)), sourcefile, __LINE__)
+
+    if (use_cndv) then
+       call this%dgvs_inst%UpdateAccVars(bounds, &
+            t_a10_patch = t_a10_patch, &
+            t_ref2m_patch = t_ref2m_patch)
+    end if
+
+  end subroutine UpdateAccVars
+
+
+  !-----------------------------------------------------------------------
+  subroutine Restart(this, bounds, ncid, flag)
+    !
+    ! !DESCRIPTION:
+    ! Handle restart (read / write) for CNVeg
+    !
+    ! Should be called regardless of whether use_cn is true
+    !
+    ! !USES:
+    use ncdio_pio, only : file_desc_t
+    use clm_varcon, only : c3_r2, c14ratio
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type), intent(in)    :: bounds 
+    type(file_desc_t), intent(inout) :: ncid   
+    character(len=*) , intent(in)    :: flag   
+    integer  :: reseed_patch(bounds%endp-bounds%begp+1)
+    integer  :: num_reseed_patch
+    !
+    ! !LOCAL VARIABLES:
+
+    integer :: begp, endp
+
+    character(len=*), parameter :: subname = 'Restart'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       begp = bounds%begp
+       endp = bounds%endp
+       call this%cnveg_carbonstate_inst%restart(bounds, ncid, flag=flag, carbon_type='c12', &
+               reseed_dead_plants=this%reseed_dead_plants, filter_reseed_patch=reseed_patch, &
+               num_reseed_patch=num_reseed_patch )
+       if ( flag /= 'read' .and. num_reseed_patch /= 0 )then
+          call endrun(msg="ERROR num_reseed should be zero and is not"//errmsg(sourcefile, __LINE__))
+       end if
+       if (use_c13) then
+          call this%c13_cnveg_carbonstate_inst%restart(bounds, ncid, flag=flag, carbon_type='c13', &
+               reseed_dead_plants=this%reseed_dead_plants, c12_cnveg_carbonstate_inst=this%cnveg_carbonstate_inst)
+       end if
+       if (use_c14) then
+          call this%c14_cnveg_carbonstate_inst%restart(bounds, ncid, flag=flag, carbon_type='c14', &
+               reseed_dead_plants=this%reseed_dead_plants, c12_cnveg_carbonstate_inst=this%cnveg_carbonstate_inst)
+       end if
+
+       call this%cnveg_carbonflux_inst%restart(bounds, ncid, flag=flag, carbon_type='c12')
+       if (use_c13) then
+          call this%c13_cnveg_carbonflux_inst%restart(bounds, ncid, flag=flag, carbon_type='c13')
+       end if
+       if (use_c14) then
+          call this%c14_cnveg_carbonflux_inst%restart(bounds, ncid, flag=flag, carbon_type='c14')
+       end if
+
+       call this%cnveg_nitrogenstate_inst%restart(bounds, ncid, flag=flag,  &
+            leafc_patch=this%cnveg_carbonstate_inst%leafc_patch(begp:endp),         &
+            leafc_storage_patch=this%cnveg_carbonstate_inst%leafc_storage_patch(begp:endp), &
+            frootc_patch=this%cnveg_carbonstate_inst%frootc_patch(begp:endp), &
+            frootc_storage_patch=this%cnveg_carbonstate_inst%frootc_storage_patch(begp:endp), &
+            deadstemc_patch=this%cnveg_carbonstate_inst%deadstemc_patch(begp:endp), &
+            filter_reseed_patch=reseed_patch, num_reseed_patch=num_reseed_patch)
+       call this%cnveg_nitrogenflux_inst%restart(bounds, ncid, flag=flag)
+       call this%cnveg_state_inst%restart(bounds, ncid, flag=flag, &
+            cnveg_carbonstate=this%cnveg_carbonstate_inst, &
+            cnveg_nitrogenstate=this%cnveg_nitrogenstate_inst, &
+            filter_reseed_patch=reseed_patch, num_reseed_patch=num_reseed_patch)
+
+       call this%c_products_inst%restart(bounds, ncid, flag)
+       if (use_c13) then
+          call this%c13_products_inst%restart(bounds, ncid, flag, &
+               template_for_missing_fields = this%c_products_inst, &
+               template_multiplier = c3_r2)
+       end if
+       if (use_c14) then
+          call this%c14_products_inst%restart(bounds, ncid, flag, &
+               template_for_missing_fields = this%c_products_inst, &
+               template_multiplier = c14ratio)
+       end if
+       call this%n_products_inst%restart(bounds, ncid, flag)
+
+    end if
+
+    if (use_cndv) then
+       call this%dgvs_inst%Restart(bounds, ncid, flag=flag)
+    end if
+
+  end subroutine Restart
+
+  !-----------------------------------------------------------------------
+  subroutine Init2(this, bounds, NLFilename)
+    !
+    ! !DESCRIPTION:
+    ! Do initialization that is needed in the initialize phase, after subgrid weights are
+    ! determined
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type) , intent(inout) :: this
+    type(bounds_type) , intent(in)    :: bounds
+    character(len=*)  , intent(in)    :: NLFilename ! namelist filename
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'Init2'
+    !-----------------------------------------------------------------------
+
+    call CNDriverInit(bounds, NLFilename, this%cnfire_method)
+
+    if (use_cndv) then
+       call dynCNDV_init(bounds, this%dgvs_inst)
+    end if
+
+  end subroutine Init2
+
+
+  !-----------------------------------------------------------------------
+  subroutine InitEachTimeStep(this, bounds, num_soilc, filter_soilc)
+    !
+    ! !DESCRIPTION:
+    ! Do initializations that need to be done at the start of every time step
+    !
+    ! This includes zeroing fluxes
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type) , intent(inout) :: this
+    type(bounds_type) , intent(in)    :: bounds
+    integer           , intent(in)    :: num_soilc       ! number of soil columns filter
+    integer           , intent(in)    :: filter_soilc(:) ! filter for soil columns
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitEachTimeStep'
+    !-----------------------------------------------------------------------
+
+    call this%cnveg_carbonflux_inst%ZeroDWT(bounds)
+    if (use_c13) then
+       call this%c13_cnveg_carbonflux_inst%ZeroDWT(bounds)
+    end if
+    if (use_c14) then
+       call this%c14_cnveg_carbonflux_inst%ZeroDWT(bounds)
+    end if
+    call this%cnveg_nitrogenflux_inst%ZeroDWT(bounds)
+    call this%cnveg_carbonstate_inst%ZeroDWT(bounds)
+    call this%cnveg_nitrogenstate_inst%ZeroDWT(bounds)
+
+  end subroutine InitEachTimeStep
+
+  !-----------------------------------------------------------------------
+  subroutine InterpFileInputs(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Interpolate inputs from files
+    !
+    ! NOTE(wjs, 2016-02-23) Stuff done here could probably be done at the end of
+    ! InitEachTimeStep, rather than in this separate routine, except for the fact that
+    ! (currently) this Interp stuff is done with proc bounds rather thna clump bounds. I
+    ! think that is needed so that you don't update a given stream multiple times. If we
+    ! rework the handling of threading / clumps so that there is a separate object for
+    ! each clump, then I think this problem would disappear - at which point we could
+    ! remove this Interp routine, moving its body to the end of InitEachTimeStep.
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type) , intent(inout) :: this
+    type(bounds_type) , intent(in)    :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InterpFileInputs'
+    !-----------------------------------------------------------------------
+
+    call this%cnfire_method%FireInterp(bounds)
+
+  end subroutine InterpFileInputs
+
+
+  !-----------------------------------------------------------------------
+  subroutine UpdateSubgridWeights(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Update subgrid weights if running with prognostic patch weights
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type) , intent(inout) :: this
+    type(bounds_type) , intent(in)    :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'UpdateSubgridWeights'
+    !-----------------------------------------------------------------------
+
+    if (use_cndv) then
+       call dynCNDV_interp(bounds, this%dgvs_inst)
+    end if
+
+  end subroutine UpdateSubgridWeights
+
+
+  !-----------------------------------------------------------------------
+  subroutine DynamicAreaConservation(this, bounds, clump_index, &
+       num_soilp_with_inactive, filter_soilp_with_inactive, &
+       num_soilc_with_inactive, filter_soilc_with_inactive, &
+       prior_weights, patch_state_updater, column_state_updater, &
+       canopystate_inst, photosyns_inst, &
+       soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, &
+       c13_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonstate_inst, &
+       soilbiogeochem_nitrogenstate_inst, ch4_inst, soilbiogeochem_state_inst)
+    !
+    ! !DESCRIPTION:
+    ! Conserve C & N with updates in subgrid weights
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    use dynPriorWeightsMod      , only : prior_weights_type
+    use dynPatchStateUpdaterMod, only : patch_state_updater_type
+    use dynColumnStateUpdaterMod, only : column_state_updater_type
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type)                       , intent(in)    :: bounds        
+
+    ! Index of clump on which we're currently operating. Note that this implies that this
+    ! routine must be called from within a clump loop.
+    integer                                 , intent(in)    :: clump_index
+
+    integer                                 , intent(in)    :: num_soilp_with_inactive ! number of points in filter_soilp_with_inactive
+    integer                                 , intent(in)    :: filter_soilp_with_inactive(:) ! soil patch filter that includes inactive points
+    integer                                 , intent(in)    :: num_soilc_with_inactive ! number of points in filter_soilc_with_inactive
+    integer                                 , intent(in)    :: filter_soilc_with_inactive(:) ! soil column filter that includes inactive points
+    type(prior_weights_type)                , intent(in)    :: prior_weights         ! weights prior to the subgrid weight updates
+    type(patch_state_updater_type)          , intent(in)    :: patch_state_updater
+    type(column_state_updater_type)         , intent(in)    :: column_state_updater
+    type(canopystate_type)                  , intent(inout) :: canopystate_inst
+    type(photosyns_type)                    , intent(inout) :: photosyns_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c13_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c14_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+    type(ch4_type)                          , intent(inout) :: ch4_inst
+    type(soilbiogeochem_state_type)         , intent(in)    :: soilbiogeochem_state_inst
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'DynamicAreaConservation'
+    !-----------------------------------------------------------------------
+
+    call t_startf('dyn_cnbal_patch')
+    call dyn_cnbal_patch(bounds, &
+         num_soilp_with_inactive, filter_soilp_with_inactive, &
+         prior_weights, patch_state_updater, &
+         canopystate_inst, photosyns_inst, &
+         this%cnveg_state_inst, &
+         this%cnveg_carbonstate_inst, this%c13_cnveg_carbonstate_inst, this%c14_cnveg_carbonstate_inst, &
+         this%cnveg_carbonflux_inst, this%c13_cnveg_carbonflux_inst, this%c14_cnveg_carbonflux_inst, &
+         this%cnveg_nitrogenstate_inst, this%cnveg_nitrogenflux_inst, &
+         soilbiogeochem_carbonflux_inst, soilbiogeochem_state_inst)
+    call t_stopf('dyn_cnbal_patch')
+
+    ! It is important to update column-level state variables based on the fluxes
+    ! generated by dyn_cnbal_patch (which handles the change in aboveground / patch-level
+    ! C/N due to shrinking patches), before calling dyn_cnbal_col (which handles the
+    ! change in belowground / column-level C/N due to changing column areas). This way,
+    ! any aboveground biomass which is sent to litter or soil due to shrinking patch
+    ! areas is accounted for by the column-level conservation. This is important if
+    ! column weights on the grid cell are changing at the same time as patch weights on
+    ! the grid cell (which will typically be the case when columns change in area).
+    !
+    ! The filters here need to include inactive points as well as active points so that
+    ! we correctly update column states in columns that have just shrunk to 0 area -
+    ! since those column states are still important in the following dyn_cnbal_col.
+    call t_startf('CNUpdateDynPatch')
+    call CStateUpdateDynPatch(bounds, num_soilc_with_inactive, filter_soilc_with_inactive, &
+         this%cnveg_carbonflux_inst, this%cnveg_carbonstate_inst, &
+         soilbiogeochem_carbonstate_inst)
+    if (use_c13) then
+       call CStateUpdateDynPatch(bounds, num_soilc_with_inactive, filter_soilc_with_inactive, &
+            this%c13_cnveg_carbonflux_inst, this%c13_cnveg_carbonstate_inst, &
+            c13_soilbiogeochem_carbonstate_inst)
+    end if
+    if (use_c14) then
+       call CStateUpdateDynPatch(bounds, num_soilc_with_inactive, filter_soilc_with_inactive, &
+            this%c14_cnveg_carbonflux_inst, this%c14_cnveg_carbonstate_inst, &
+            c14_soilbiogeochem_carbonstate_inst)
+    end if
+    call NStateUpdateDynPatch(bounds, num_soilc_with_inactive, filter_soilc_with_inactive, &
+         this%cnveg_nitrogenflux_inst, this%cnveg_nitrogenstate_inst, &
+         soilbiogeochem_nitrogenstate_inst)
+    call t_stopf('CNUpdateDynPatch')
+
+    ! This call fixes issue #741 by performing precision control on decomp_cpools_vr_col
+    call t_startf('SoilBiogeochemPrecisionControl')
+    call SoilBiogeochemPrecisionControl(num_soilc_with_inactive, filter_soilc_with_inactive,  &
+         soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonstate_inst,soilbiogeochem_nitrogenstate_inst)
+    call t_stopf('SoilBiogeochemPrecisionControl')
+
+    call t_startf('dyn_cnbal_col')
+    call dyn_cnbal_col(bounds, clump_index, column_state_updater, &
+         soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonstate_inst, soilbiogeochem_nitrogenstate_inst, &
+         ch4_inst)
+    call t_stopf('dyn_cnbal_col')
+
+  end subroutine DynamicAreaConservation
+
+  !-----------------------------------------------------------------------
+  subroutine InitColumnBalance(this, bounds, num_allc, filter_allc, &
+       num_soilc, filter_soilc, num_soilp, filter_soilp, &
+       soilbiogeochem_carbonstate_inst, &
+       c13_soilbiogeochem_carbonstate_inst, &
+       c14_soilbiogeochem_carbonstate_inst, &
+       soilbiogeochem_nitrogenstate_inst)
+    !
+    ! !DESCRIPTION:
+    ! Set the starting point for column-level balance checks.
+    !
+    ! This should be called after DynamicAreaConservation, since the changes made by
+    ! DynamicAreaConservation can break column-level conservation checks.
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type)               , intent(inout) :: this
+    type(bounds_type)                       , intent(in)    :: bounds  
+    integer                                 , intent(in)    :: num_allc          ! number of columns in allc filter
+    integer                                 , intent(in)    :: filter_allc(:)    ! filter for all active columns
+    integer                                 , intent(in)    :: num_soilc         ! number of soil columns in filter
+    integer                                 , intent(in)    :: filter_soilc(:)   ! filter for soil columns
+    integer                                 , intent(in)    :: num_soilp         ! number of soil patches in filter
+    integer                                 , intent(in)    :: filter_soilp(:)   ! filter for soil patches
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c13_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c14_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitColumnBalance'
+    !-----------------------------------------------------------------------
+
+    call CNDriverSummarizeStates(bounds, &
+         num_allc, filter_allc, &
+         num_soilc, filter_soilc, &
+         num_soilp, filter_soilp, &
+         this%cnveg_carbonstate_inst, &
+         this%c13_cnveg_carbonstate_inst, &
+         this%c14_cnveg_carbonstate_inst, &
+         this%cnveg_nitrogenstate_inst, &
+         soilbiogeochem_carbonstate_inst, &
+         c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonstate_inst, &
+         soilbiogeochem_nitrogenstate_inst)
+
+    call this%cn_balance_inst%BeginCNColumnBalance( &
+         bounds, num_soilc, filter_soilc, &
+         this%cnveg_carbonstate_inst, this%cnveg_nitrogenstate_inst)
+
+  end subroutine InitColumnBalance
+
+
+  !-----------------------------------------------------------------------
+  subroutine InitGridcellBalance(this, bounds, num_allc, filter_allc, &
+       num_soilc, filter_soilc, num_soilp, filter_soilp, &
+       soilbiogeochem_carbonstate_inst, &
+       c13_soilbiogeochem_carbonstate_inst, &
+       c14_soilbiogeochem_carbonstate_inst, &
+       soilbiogeochem_nitrogenstate_inst)
+    !
+    ! !DESCRIPTION:
+    ! Set the starting point for gridcell-level balance checks.
+    !
+    ! Gridcell level:
+    ! Called before DynamicAreaConservation.
+    !
+    ! !USES:
+    use subgridAveMod, only : c2g
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type)               , intent(inout) :: this
+    type(bounds_type)                       , intent(in)    :: bounds
+    integer                                 , intent(in)    :: num_allc          ! number of columns in allc filter
+    integer                                 , intent(in)    :: filter_allc(:)    ! filter for all active columns
+    integer                                 , intent(in)    :: num_soilc         ! number of soil columns in filter
+    integer                                 , intent(in)    :: filter_soilc(:)   ! filter for soil columns
+    integer                                 , intent(in)    :: num_soilp         ! number of soil patches in filter
+    integer                                 , intent(in)    :: filter_soilp(:)   ! filter for soil patches
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c13_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c14_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitGridcellBalance'
+    !-----------------------------------------------------------------------
+
+    call CNDriverSummarizeStates(bounds, &
+         num_allc, filter_allc, &
+         num_soilc, filter_soilc, &
+         num_soilp, filter_soilp, &
+         this%cnveg_carbonstate_inst, &
+         this%c13_cnveg_carbonstate_inst, &
+         this%c14_cnveg_carbonstate_inst, &
+         this%cnveg_nitrogenstate_inst, &
+         soilbiogeochem_carbonstate_inst, &
+         c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonstate_inst, &
+         soilbiogeochem_nitrogenstate_inst)
+
+    ! total gridcell carbon (TOTGRIDCELLC)
+    call c2g( bounds = bounds, &
+         carr = this%cnveg_carbonstate_inst%totc_col(bounds%begc:bounds%endc), &
+         garr = this%cnveg_carbonstate_inst%totc_grc(bounds%begg:bounds%endg), &
+         c2l_scale_type = 'unity', &
+         l2g_scale_type = 'unity')
+    ! total gridcell nitrogen (TOTGRIDCELLN)
+    call c2g( bounds = bounds, &
+         carr = this%cnveg_nitrogenstate_inst%totn_col(bounds%begc:bounds%endc), &
+         garr = this%cnveg_nitrogenstate_inst%totn_grc(bounds%begg:bounds%endg), &
+         c2l_scale_type = 'unity', &
+         l2g_scale_type = 'unity')
+
+    call this%cn_balance_inst%BeginCNGridcellBalance( bounds, &
+         this%cnveg_carbonflux_inst, &
+         this%cnveg_carbonstate_inst, this%cnveg_nitrogenstate_inst, &
+         this%c_products_inst, this%n_products_inst)
+
+  end subroutine InitGridcellBalance
+
+
+  !-----------------------------------------------------------------------
+  subroutine EcosystemDynamicsPreDrainage(this, bounds, &
+       num_soilc, filter_soilc, &
+       num_soilp, filter_soilp, &
+       num_pcropp, filter_pcropp, &
+       doalb, &
+       soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst,         &
+       c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, &
+       c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, &
+       soilbiogeochem_state_inst,                                               &
+       soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst,     &
+       active_layer_inst, &
+       atm2lnd_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst,                           &
+       wateratm2lndbulk_inst, canopystate_inst, soilstate_inst, temperature_inst, crop_inst, ch4_inst, &
+       photosyns_inst, saturated_excess_runoff_inst, energyflux_inst,          &
+       nutrient_competition_method, fireemis_inst)
+    !
+    ! !DESCRIPTION:
+    ! Do the main science for CN vegetation that needs to be done before hydrology-drainage
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type)               , intent(inout) :: this
+    type(bounds_type)                       , intent(in)    :: bounds  
+    integer                                 , intent(in)    :: num_soilc         ! number of soil columns in filter
+    integer                                 , intent(in)    :: filter_soilc(:)   ! filter for soil columns
+    integer                                 , intent(in)    :: num_soilp         ! number of soil patches in filter
+    integer                                 , intent(in)    :: filter_soilp(:)   ! filter for soil patches
+    integer                                 , intent(in)    :: num_pcropp        ! number of prog. crop patches in filter
+    integer                                 , intent(in)    :: filter_pcropp(:)  ! filter for prognostic crop patches
+    logical                                 , intent(in)    :: doalb             ! true = surface albedo calculation time step
+    type(soilbiogeochem_state_type)         , intent(inout) :: soilbiogeochem_state_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: c13_soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c13_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: c14_soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c14_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_nitrogenflux_type)  , intent(inout) :: soilbiogeochem_nitrogenflux_inst
+    type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+    type(active_layer_type)                 , intent(in)    :: active_layer_inst
+    type(atm2lnd_type)                      , intent(in)    :: atm2lnd_inst
+    type(waterstatebulk_type)                   , intent(in)    :: waterstatebulk_inst
+    type(waterdiagnosticbulk_type)                   , intent(in)    :: waterdiagnosticbulk_inst
+    type(waterfluxbulk_type)                    , intent(inout) :: waterfluxbulk_inst
+    type(wateratm2lndbulk_type)                    , intent(inout) :: wateratm2lndbulk_inst
+    type(canopystate_type)                  , intent(inout) :: canopystate_inst
+    type(soilstate_type)                    , intent(inout) :: soilstate_inst
+    type(temperature_type)                  , intent(inout) :: temperature_inst
+    type(crop_type)                         , intent(inout) :: crop_inst
+    type(ch4_type)                          , intent(in)    :: ch4_inst
+    type(photosyns_type)                    , intent(in)    :: photosyns_inst
+    type(saturated_excess_runoff_type)      , intent(in)    :: saturated_excess_runoff_inst
+    type(energyflux_type)                   , intent(in)    :: energyflux_inst
+    class(nutrient_competition_method_type) , intent(inout) :: nutrient_competition_method
+    type(fireemis_type)                     , intent(inout) :: fireemis_inst
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'EcosystemDynamicsPreDrainage'
+    !-----------------------------------------------------------------------
+
+    call crop_inst%CropIncrementYear(num_pcropp, filter_pcropp)
+
+    call CNDriverNoLeaching(bounds,                                         &
+         num_soilc, filter_soilc,                       &
+         num_soilp, filter_soilp,                       &
+         num_pcropp, filter_pcropp, doalb,              &
+         this%cnveg_state_inst,                                                        &
+         this%cnveg_carbonflux_inst, this%cnveg_carbonstate_inst,                           &
+         this%c13_cnveg_carbonflux_inst, this%c13_cnveg_carbonstate_inst,                   &
+         this%c14_cnveg_carbonflux_inst, this%c14_cnveg_carbonstate_inst,                   &
+         this%cnveg_nitrogenflux_inst, this%cnveg_nitrogenstate_inst,                       &
+         this%c_products_inst, this%c13_products_inst, this%c14_products_inst,    &
+         this%n_products_inst,                                                    &
+         soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst,         &
+         c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, &
+         soilbiogeochem_state_inst,                                               &
+         soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst,     &
+         active_layer_inst, &
+         atm2lnd_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst,                           &
+         wateratm2lndbulk_inst, canopystate_inst, soilstate_inst, temperature_inst, crop_inst, ch4_inst, &
+         this%dgvs_inst, photosyns_inst, saturated_excess_runoff_inst, energyflux_inst,          &
+         nutrient_competition_method, this%cnfire_method, this%dribble_crophrv_xsmrpool_2atm)
+
+    ! fire carbon emissions 
+    call CNFireEmisUpdate(bounds, num_soilp, filter_soilp, &
+         this%cnveg_carbonflux_inst, this%cnveg_carbonstate_inst, fireemis_inst )
+
+    call CNAnnualUpdate(bounds,            &
+         num_soilc, filter_soilc, &
+         num_soilp, filter_soilp, &
+         this%cnveg_state_inst, this%cnveg_carbonflux_inst)
+
+  end subroutine EcosystemDynamicsPreDrainage
+
+  !-----------------------------------------------------------------------
+  subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, &
+       num_soilc, filter_soilc, num_soilp, filter_soilp, doalb, crop_inst, &
+       waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, frictionvel_inst, canopystate_inst, &
+       soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, &
+       c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, &
+       c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, &
+       soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst)
+    !
+    ! !DESCRIPTION:
+    ! Do the main science for CN vegetation that needs to be done after hydrology-drainage
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type)               , intent(inout) :: this
+    type(bounds_type)                       , intent(in)    :: bounds  
+    integer                                 , intent(in)    :: num_allc          ! number of columns in allc filter
+    integer                                 , intent(in)    :: filter_allc(:)    ! filter for all active columns
+    integer                                 , intent(in)    :: num_soilc         ! number of soil columns in filter
+    integer                                 , intent(in)    :: filter_soilc(:)   ! filter for soil columns
+    integer                                 , intent(in)    :: num_soilp         ! number of soil patches in filter
+    integer                                 , intent(in)    :: filter_soilp(:)   ! filter for soil patches
+    logical                                 , intent(in)    :: doalb             ! true = surface albedo calculation time step
+    type(crop_type)                         , intent(in)    :: crop_inst
+    type(waterstatebulk_type)                   , intent(in)    :: waterstatebulk_inst
+    type(waterdiagnosticbulk_type)                   , intent(in)    :: waterdiagnosticbulk_inst
+    type(waterfluxbulk_type)                    , intent(inout) :: waterfluxbulk_inst
+    type(frictionvel_type)                  , intent(in)    :: frictionvel_inst
+    type(canopystate_type)                  , intent(inout) :: canopystate_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: c13_soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c13_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: c14_soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_carbonstate_type)   , intent(inout) :: c14_soilbiogeochem_carbonstate_inst
+    type(soilbiogeochem_nitrogenflux_type)  , intent(inout) :: soilbiogeochem_nitrogenflux_inst
+    type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'EcosystemDynamicsPostDrainage'
+    !-----------------------------------------------------------------------
+
+    ! Update the nitrogen leaching rate as a function of soluble mineral N 
+    ! and total soil water outflow.
+    
+    call CNDriverLeaching(bounds, &
+         num_soilc, filter_soilc, &
+         num_soilp, filter_soilp, &
+         waterstatebulk_inst, waterfluxbulk_inst, &
+         this%cnveg_nitrogenflux_inst, this%cnveg_nitrogenstate_inst, &
+         soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst)
+
+    ! Set controls on very low values in critical state variables 
+
+    call t_startf('CNPrecisionControl')
+    call CNPrecisionControl(bounds, num_soilp, filter_soilp, &
+         this%cnveg_carbonstate_inst, this%c13_cnveg_carbonstate_inst, &
+         this%c14_cnveg_carbonstate_inst, this%cnveg_nitrogenstate_inst)
+    call t_stopf('CNPrecisionControl')
+
+    call t_startf('SoilBiogeochemPrecisionControl')
+    call SoilBiogeochemPrecisionControl(num_soilc, filter_soilc,  &
+         soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonstate_inst,soilbiogeochem_nitrogenstate_inst)
+    call t_stopf('SoilBiogeochemPrecisionControl')
+
+    ! Call to all CN summary routines
+
+    call  CNDriverSummarizeStates(bounds, &
+         num_allc, filter_allc, &
+         num_soilc, filter_soilc, &
+         num_soilp, filter_soilp, &
+         this%cnveg_carbonstate_inst, &
+         this%c13_cnveg_carbonstate_inst, &
+         this%c14_cnveg_carbonstate_inst, &
+         this%cnveg_nitrogenstate_inst, &
+         soilbiogeochem_carbonstate_inst, &
+         c13_soilbiogeochem_carbonstate_inst, &
+         c14_soilbiogeochem_carbonstate_inst, &
+         soilbiogeochem_nitrogenstate_inst)
+
+    call  CNDriverSummarizeFluxes(bounds, &
+         num_soilc, filter_soilc, &
+         num_soilp, filter_soilp, &
+         this%cnveg_carbonflux_inst, &
+         this%c13_cnveg_carbonflux_inst, &
+         this%c14_cnveg_carbonflux_inst, &
+         this%cnveg_nitrogenflux_inst, &
+         this%c_products_inst, this%c13_products_inst, this%c14_products_inst, &
+         soilbiogeochem_carbonflux_inst, &
+         c13_soilbiogeochem_carbonflux_inst, &
+         c14_soilbiogeochem_carbonflux_inst, &
+         soilbiogeochem_nitrogenflux_inst)
+
+    ! On the radiation time step, use C state variables to calculate
+    ! vegetation structure (LAI, SAI, height)
+
+    if (doalb) then   
+       call CNVegStructUpdate(num_soilp, filter_soilp, &
+            waterdiagnosticbulk_inst, frictionvel_inst, this%dgvs_inst, this%cnveg_state_inst, &
+            crop_inst, this%cnveg_carbonstate_inst, canopystate_inst)
+    end if
+
+  end subroutine EcosystemDynamicsPostDrainage
+
+  !-----------------------------------------------------------------------
+  subroutine BalanceCheck(this, bounds, num_soilc, filter_soilc, &
+       soilbiogeochem_carbonflux_inst, soilbiogeochem_nitrogenflux_inst, &
+       atm2lnd_inst)
+    !
+    ! !DESCRIPTION:
+    ! Check the carbon and nitrogen balance
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    use clm_time_manager   , only : get_nstep_since_startup_or_lastDA_restart_or_pause
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type)               , intent(inout) :: this
+    type(bounds_type)                       , intent(in)    :: bounds  
+    integer                                 , intent(in)    :: num_soilc         ! number of soil columns in filter
+    integer                                 , intent(in)    :: filter_soilc(:)   ! filter for soil columns
+    type(soilbiogeochem_carbonflux_type)    , intent(inout) :: soilbiogeochem_carbonflux_inst
+    type(soilbiogeochem_nitrogenflux_type)  , intent(inout) :: soilbiogeochem_nitrogenflux_inst
+    type(atm2lnd_type)                      , intent(in)    :: atm2lnd_inst
+    !
+    ! !LOCAL VARIABLES:
+    integer              :: DA_nstep                   ! time step number
+
+    character(len=*), parameter :: subname = 'BalanceCheck'
+    !-----------------------------------------------------------------------
+
+    DA_nstep = get_nstep_since_startup_or_lastDA_restart_or_pause()
+
+    ! DART SourceMod edit
+    ! Always call C and N balances, allow CNBalanceCheckMod subroutines to skip/perform balance checks
+
+    !if (DA_nstep <= skip_steps )then
+    !   if (masterproc) then
+    !      write(iulog,*) '--WARNING-- skipping CN balance check for first timesteps after startup or data assimilation'
+    !   end if
+    !else
+
+       call this%cn_balance_inst%CBalanceCheck( &
+            bounds, num_soilc, filter_soilc, &
+            soilbiogeochem_carbonflux_inst, &
+            this%cnveg_carbonflux_inst, &
+            this%cnveg_carbonstate_inst, &
+            this%c_products_inst)
+
+       call this%cn_balance_inst%NBalanceCheck( &
+            bounds, num_soilc, filter_soilc, &
+            soilbiogeochem_nitrogenflux_inst, &
+            this%cnveg_nitrogenflux_inst, &
+            this%cnveg_nitrogenstate_inst, &
+            this%n_products_inst, &
+            atm2lnd_inst)
+
+    !end if
+
+  end subroutine BalanceCheck
+
+  !-----------------------------------------------------------------------
+  subroutine EndOfTimeStepVegDynamics(this, bounds, num_natvegp, filter_natvegp, &
+       atm2lnd_inst, wateratm2lndbulk_inst)
+    !
+    ! !DESCRIPTION:
+    ! Do vegetation dynamics that should be done at the end of each time step
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(inout) :: this
+    type(bounds_type)  , intent(in)    :: bounds                  
+    integer            , intent(inout) :: num_natvegp       ! number of naturally-vegetated patches in filter
+    integer            , intent(inout) :: filter_natvegp(:) ! filter for naturally-vegetated patches
+    type(atm2lnd_type) , intent(inout) :: atm2lnd_inst
+    type(wateratm2lndbulk_type) , intent(inout) :: wateratm2lndbulk_inst
+    !
+    ! !LOCAL VARIABLES:
+    integer  :: nstep  ! time step number
+    integer  :: yr     ! year (0, ...)
+    integer  :: mon    ! month (1, ..., 12)
+    integer  :: day    ! day of month (1, ..., 31)
+    integer  :: sec    ! seconds of the day
+    integer  :: ncdate ! current date
+    integer  :: nbdate ! base date (reference date)
+    integer  :: kyr    ! thousand years, equals 2 at end of first year
+
+    character(len=*), parameter :: subname = 'EndOfTimeStepVegDynamics'
+    !-----------------------------------------------------------------------
+
+    if (use_cndv) then
+       ! Call dv (dynamic vegetation) at last time step of year
+
+       call t_startf('d2dgvm')
+       if (is_end_curr_year() .and. .not. is_first_step())  then
+
+          ! Get date info.  kyr is used in lpj().  At end of first year, kyr = 2.
+          call get_curr_date(yr, mon, day, sec)
+          ncdate = yr*10000 + mon*100 + day
+          call get_ref_date(yr, mon, day, sec)
+          nbdate = yr*10000 + mon*100 + day
+          kyr = ncdate/10000 - nbdate/10000 + 1
+
+          if (masterproc) then
+             nstep = get_nstep()
+             write(iulog,*) 'End of year. CNDV called now: ncdate=', &
+                  ncdate,' nbdate=',nbdate,' kyr=',kyr,' nstep=', nstep
+          end if
+
+          call CNDVDriver(bounds, &
+               num_natvegp, filter_natvegp, kyr,  &
+               atm2lnd_inst, wateratm2lndbulk_inst, &
+               this%cnveg_carbonflux_inst, this%cnveg_carbonstate_inst, this%dgvs_inst)
+       end if
+       call t_stopf('d2dgvm')
+    end if
+
+  end subroutine EndOfTimeStepVegDynamics
+
+  !-----------------------------------------------------------------------
+  subroutine WriteHistory(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Do any history writes that are specific to vegetation dynamics
+    !
+    ! NOTE(wjs, 2016-02-23) This could probably be combined with
+    ! EndOfTimeStepVegDynamics, except for the fact that (currently) history writes are
+    ! done with proc bounds rather than clump bounds. If that were changed, then the body
+    ! of this could be moved into EndOfTimeStepVegDynamics, inside a "if (.not.
+    ! use_noio)" conditional.
+    !
+    ! Should only be called if use_cn is true
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type)  , intent(in) :: bounds                  
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'WriteHistory'
+    !-----------------------------------------------------------------------
+
+    ! Write to CNDV history buffer if appropriate
+    if (use_cndv) then
+       if (is_end_curr_year() .and. .not. is_first_step())  then
+          call t_startf('clm_drv_io_hdgvm')
+          call CNDVHist( bounds, this%dgvs_inst )
+          if (masterproc) write(iulog,*) 'Annual CNDV calculations are complete'
+          call t_stopf('clm_drv_io_hdgvm')
+       end if
+    end if
+
+  end subroutine WriteHistory
+
+
+  !-----------------------------------------------------------------------
+  function get_net_carbon_exchange_grc(this, bounds) result(net_carbon_exchange_grc)
+    !
+    ! !DESCRIPTION:
+    ! Get gridcell-level net carbon exchange array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: net_carbon_exchange_grc(bounds%begg:bounds%endg)  ! function result: net carbon exchange between land and atmosphere, includes fire, landuse, harvest and hrv_xsmrpool flux, positive for source (gC/m2/s)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_net_carbon_exchange_grc'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       net_carbon_exchange_grc(bounds%begg:bounds%endg) = &
+            -this%cnveg_carbonflux_inst%nbp_grc(bounds%begg:bounds%endg)
+    else
+       net_carbon_exchange_grc(bounds%begg:bounds%endg) = 0._r8
+    end if
+
+  end function get_net_carbon_exchange_grc
+
+
+  !-----------------------------------------------------------------------
+  function get_leafn_patch(this, bounds) result(leafn_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level leaf nitrogen array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: leafn_patch(bounds%begp:bounds%endp)  ! function result: leaf N (gN/m2)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_leafn_patch'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       leafn_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_nitrogenstate_inst%leafn_patch(bounds%begp:bounds%endp)
+    else
+       leafn_patch(bounds%begp:bounds%endp) = nan
+    end if
+
+  end function get_leafn_patch
+
+  !-----------------------------------------------------------------------
+  function get_downreg_patch(this, bounds) result(downreg_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level downregulation array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: downreg_patch(bounds%begp:bounds%endp)  ! function result: fractional reduction in GPP due to N limitation (dimensionless)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_downreg_patch'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       downreg_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_state_inst%downreg_patch(bounds%begp:bounds%endp)
+    else
+       downreg_patch(bounds%begp:bounds%endp) = nan
+    end if
+
+  end function get_downreg_patch
+
+  !-----------------------------------------------------------------------
+  function get_root_respiration_patch(this, bounds) result(root_respiration_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level root respiration array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: root_respiration_patch(bounds%begp:bounds%endp)  ! function result: root respiration (fine root MR + total root GR) (gC/m2/s)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_root_respiration_patch'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       root_respiration_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_carbonflux_inst%rr_patch(bounds%begp:bounds%endp)
+    else
+       root_respiration_patch(bounds%begp:bounds%endp) = nan
+    end if
+
+  end function get_root_respiration_patch
+
+  ! TODO(wjs, 2016-02-19) annsum_npp, agnpp and bgnpp are all needed for the estimation
+  ! of tillers in ch4Mod. Rather than providing getters for these three things so that
+  ! ch4Mod can estimate tillers, it would probably be better if the tiller estimation
+  ! algorithm was moved into some CNVeg-specific module, and then tillers could be
+  ! queried directly.
+
+  !-----------------------------------------------------------------------
+  function get_annsum_npp_patch(this, bounds) result(annsum_npp_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level annual sum NPP array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: annsum_npp_patch(bounds%begp:bounds%endp)  ! function result: annual sum NPP (gC/m2/yr)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_annsum_npp_patch'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       annsum_npp_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_carbonflux_inst%annsum_npp_patch(bounds%begp:bounds%endp)
+    else
+       annsum_npp_patch(bounds%begp:bounds%endp) = nan
+    end if
+
+  end function get_annsum_npp_patch
+
+  !-----------------------------------------------------------------------
+  function get_agnpp_patch(this, bounds) result(agnpp_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level aboveground NPP array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: agnpp_patch(bounds%begp:bounds%endp)  ! function result: aboveground NPP (gC/m2/s)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_agnpp_patch'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       agnpp_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_carbonflux_inst%agnpp_patch(bounds%begp:bounds%endp)
+    else
+       agnpp_patch(bounds%begp:bounds%endp) = nan
+    end if
+
+  end function get_agnpp_patch
+
+  !-----------------------------------------------------------------------
+  function get_bgnpp_patch(this, bounds) result(bgnpp_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level belowground NPP array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: bgnpp_patch(bounds%begp:bounds%endp)  ! function result: belowground NPP (gC/m2/s)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_bgnpp_patch'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       bgnpp_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_carbonflux_inst%bgnpp_patch(bounds%begp:bounds%endp)
+    else
+       bgnpp_patch(bounds%begp:bounds%endp) = nan
+    end if
+
+  end function get_bgnpp_patch
+
+  !-----------------------------------------------------------------------
+  function get_froot_carbon_patch(this, bounds, tlai) result(froot_carbon_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level fine root carbon array
+    !
+    ! !USES:
+    use pftconMod           , only : pftcon
+    use PatchType           , only : patch
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8)         , intent(in) :: tlai( bounds%begp: )
+    real(r8) :: froot_carbon_patch(bounds%begp:bounds%endp)  ! function result: (gC/m2)
+    !
+    ! !LOCAL VARIABLES:
+    character(len=*), parameter :: subname = 'get_froot_carbon_patch'
+    integer :: p
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       froot_carbon_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_carbonstate_inst%frootc_patch(bounds%begp:bounds%endp)
+    else
+! To get leaf biomass:
+! bleaf = LAI / slatop
+! g/m2 =  m2/m2  / m2/g 
+! To get root biomass: 
+! broot = bleaf * froot_leaf(ivt(p))
+! g/m2 = g/m2 * g/g
+       do p=bounds%begp, bounds%endp
+          if (pftcon%slatop(patch%itype(p)) > 0._r8) then 
+             froot_carbon_patch(p) = tlai(p) &
+                  / pftcon%slatop(patch%itype(p)) &
+                  *pftcon%froot_leaf(patch%itype(p))
+          else
+             froot_carbon_patch(p) = 0._r8
+          endif
+       enddo
+    end if
+
+  end function get_froot_carbon_patch
+
+  !-----------------------------------------------------------------------
+  function get_croot_carbon_patch(this, bounds, tlai) result(croot_carbon_patch)
+    !
+    ! !DESCRIPTION:
+    ! Get patch-level live coarse root carbon array
+    !
+    ! !USES:
+    use pftconMod           , only : pftcon
+    use PatchType           , only : patch
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8)         , intent(in) :: tlai( bounds%begp: )
+    real(r8) :: croot_carbon_patch(bounds%begp:bounds%endp)  ! function result: (gC/m2)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_croot_carbon_patch'
+    integer :: p
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       croot_carbon_patch(bounds%begp:bounds%endp) = &
+            this%cnveg_carbonstate_inst%livecrootc_patch(bounds%begp:bounds%endp)
+    else
+! To get leaf biomass:
+! bleaf = LAI / slatop
+! g/m2 =  m2/m2  / m2/g 
+! To get root biomass: 
+! broot = bleaf * froot_leaf(ivt(p))
+! g/m2 = g/m2 * g/g
+       do p=bounds%begp, bounds%endp
+          if (pftcon%slatop(patch%itype(p)) > 0._r8) then 
+             croot_carbon_patch(p) = tlai(p) &
+                  / pftcon%slatop(patch%itype(p)) &
+                  *pftcon%stem_leaf(patch%itype(p)) &
+                  *pftcon%croot_stem(patch%itype(p))
+          else
+             croot_carbon_patch(p) = 0._r8
+          endif
+       enddo
+    end if
+
+  end function get_croot_carbon_patch
+
+  !-----------------------------------------------------------------------
+  function get_totvegc_col(this, bounds) result(totvegc_col)
+    !
+    ! !DESCRIPTION:
+    ! Get column-level total vegetation carbon array
+    !
+    ! !USES:
+    !
+    ! !ARGUMENTS:
+    class(cn_vegetation_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    real(r8) :: totvegc_col(bounds%begc:bounds%endc)  ! function result: (gC/m2)
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'get_totvegc_col'
+    !-----------------------------------------------------------------------
+
+    if (use_cn) then
+       totvegc_col(bounds%begc:bounds%endc) = &
+            this%cnveg_carbonstate_inst%totvegc_col(bounds%begc:bounds%endc)
+    else
+       totvegc_col(bounds%begc:bounds%endc) = nan
+    end if
+
+  end function get_totvegc_col
+
+
+end module CNVegetationFacade
diff --git a/models/clm/DART_SourceMods/cesm2_2_0/SourceMods/src.clm/biogeophys/BalanceCheckMod.F90 b/models/clm/DART_SourceMods/cesm2_2_0/SourceMods/src.clm/biogeophys/BalanceCheckMod.F90
new file mode 100644
index 000000000..da8644716
--- /dev/null
+++ b/models/clm/DART_SourceMods/cesm2_2_0/SourceMods/src.clm/biogeophys/BalanceCheckMod.F90
@@ -0,0 +1,743 @@
+module BalanceCheckMod
+
+  !-----------------------------------------------------------------------
+  ! !DESCRIPTION:
+  ! Water and energy balance check.
+  !
+  ! !USES:
+#include "shr_assert.h"
+  use shr_kind_mod       , only : r8 => shr_kind_r8
+  use shr_log_mod        , only : errMsg => shr_log_errMsg
+  use decompMod          , only : bounds_type
+  use abortutils         , only : endrun
+  use clm_varctl         , only : iulog
+  use clm_varcon         , only : namep, namec
+  use clm_varpar         , only : nlevsoi
+  use GetGlobalValuesMod , only : GetGlobalIndex
+  use atm2lndType        , only : atm2lnd_type
+  use EnergyFluxType     , only : energyflux_type
+  use SolarAbsorbedType  , only : solarabs_type
+  use SoilHydrologyType  , only : soilhydrology_type
+  use SurfaceAlbedoType  , only : surfalb_type
+  use WaterStateType     , only : waterstate_type
+  use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type
+  use WaterDiagnosticType, only : waterdiagnostic_type
+  use Wateratm2lndType   , only : wateratm2lnd_type
+  use WaterBalanceType   , only : waterbalance_type
+  use WaterFluxType      , only : waterflux_type
+  use WaterType          , only : water_type
+  use TotalWaterAndHeatMod, only : ComputeWaterMassNonLake, ComputeWaterMassLake
+  use GridcellType       , only : grc                
+  use LandunitType       , only : lun                
+  use ColumnType         , only : col                
+  use PatchType          , only : patch                
+  use landunit_varcon    , only : istdlak, istsoil,istcrop,istwet,istice_mec
+  use column_varcon      , only : icol_roof, icol_sunwall, icol_shadewall
+  use column_varcon      , only : icol_road_perv, icol_road_imperv
+  use clm_time_manager   , only : is_first_restart_step
+  !
+  ! !PUBLIC TYPES:
+  implicit none
+  save
+  private
+  !
+  ! !PUBLIC MEMBER FUNCTIONS:
+
+  public :: BalanceCheckInit         ! Initialization of Water and energy balance check
+  public :: BeginWaterBalance        ! Initialize water balance check
+  public :: BalanceCheck             ! Water and energy balance check
+  public :: GetBalanceCheckSkipSteps ! Get the number of steps to skip for the balance check
+  public :: BalanceCheckClean        ! Clean up for BalanceCheck
+
+  ! !PRIVATE MEMBER DATA:
+  real(r8), private, parameter :: skip_size = 3600.0_r8   ! Time steps to skip the balance check at startup (sec)
+  integer,  private            :: skip_steps = -999       ! Number of time steps to skip the balance check for
+
+  !
+  ! !PRIVATE MEMBER FUNCTIONS:
+  private :: BeginWaterBalanceSingle ! Initialize water balance check for bulk or a single tracer
+
+  character(len=*), parameter, private :: sourcefile = &
+       __FILE__
+  !-----------------------------------------------------------------------
+
+contains
+
+  !-----------------------------------------------------------------------
+  subroutine BalanceCheckInit( )
+  !-----------------------------------------------------------------------
+    !
+    ! !DESCRIPTION:
+    ! Initialize balance check
+    !
+    ! !USES:
+    use spmdMod           , only : masterproc
+    use clm_time_manager  , only : get_step_size_real
+    ! !ARGUMENTS:
+    !
+    ! !LOCAL VARIABLES:
+    real(r8) :: dtime                    ! land model time step (sec)
+    !-----------------------------------------------------------------------
+    dtime = get_step_size_real()
+    ! Skip a minimum of two time steps, but otherwise skip the number of time-steps in the skip_size rounded to the nearest integer
+    skip_steps = max(2, nint( (skip_size / dtime) ) )
+
+    if ( masterproc ) write(iulog,*) ' Skip balance checking for the first ', skip_steps, ' time steps'
+
+  end subroutine BalanceCheckInit
+
+  !-----------------------------------------------------------------------
+  subroutine BalanceCheckClean( )
+  !-----------------------------------------------------------------------
+    !
+    ! !DESCRIPTION:
+    ! Clean up BalanceCheck
+    !
+    ! !USES:
+    ! !ARGUMENTS:
+    !
+    ! !LOCAL VARIABLES:
+    !-----------------------------------------------------------------------
+    skip_steps = -999
+  end subroutine BalanceCheckClean
+
+  !-----------------------------------------------------------------------
+  function GetBalanceCheckSkipSteps( ) result( get_skip )
+  !-----------------------------------------------------------------------
+    !
+    ! !DESCRIPTION:
+    ! Get the number of steps to skip for the balance check
+    !
+    ! !ARGUMENTS:
+    integer :: get_skip    ! function result
+    ! !LOCAL VARIABLES:
+    if ( skip_steps > 0 )then
+       get_skip = skip_steps
+    else
+       get_skip = -999
+       call endrun('ERROR:: GetBalanceCheckSkipSteps called before BalanceCheckInit')
+    end if
+  end function GetBalanceCheckSkipSteps
+  !-----------------------------------------------------------------------
+
+  !-----------------------------------------------------------------------
+  subroutine BeginWaterBalance(bounds, &
+       num_nolakec, filter_nolakec, num_lakec, filter_lakec, &
+       water_inst, soilhydrology_inst, &
+       use_aquifer_layer)
+    !
+    ! !DESCRIPTION:
+    ! Initialize column-level water balance at beginning of time step, for bulk water and
+    ! each water tracer
+    !
+    ! !ARGUMENTS:
+    type(bounds_type)         , intent(in)    :: bounds
+    integer                   , intent(in)    :: num_nolakec          ! number of column non-lake points in column filter
+    integer                   , intent(in)    :: filter_nolakec(:)    ! column filter for non-lake points
+    integer                   , intent(in)    :: num_lakec            ! number of column lake points in column filter
+    integer                   , intent(in)    :: filter_lakec(:)      ! column filter for lake points
+    type(water_type)          , intent(inout) :: water_inst
+    type(soilhydrology_type)  , intent(in)    :: soilhydrology_inst
+    logical                   , intent(in)    :: use_aquifer_layer    ! whether an aquifer layer is used in this run
+    !
+    ! !LOCAL VARIABLES:
+    integer :: i
+
+    character(len=*), parameter :: subname = 'BeginWaterBalance'
+    !-----------------------------------------------------------------------
+
+    do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end
+       call BeginWaterBalanceSingle(bounds, &
+            num_nolakec, filter_nolakec, &
+            num_lakec, filter_lakec, &
+            soilhydrology_inst, &
+            water_inst%bulk_and_tracers(i)%waterstate_inst, &
+            water_inst%bulk_and_tracers(i)%waterdiagnostic_inst, &
+            water_inst%bulk_and_tracers(i)%waterbalance_inst, &
+            use_aquifer_layer = use_aquifer_layer)
+    end do
+
+  end subroutine BeginWaterBalance
+
+  !-----------------------------------------------------------------------
+  subroutine BeginWaterBalanceSingle(bounds, &
+       num_nolakec, filter_nolakec, num_lakec, filter_lakec, &
+       soilhydrology_inst, waterstate_inst, waterdiagnostic_inst, waterbalance_inst, &
+       use_aquifer_layer)
+    !
+    ! !DESCRIPTION:
+    ! Initialize column-level water balance at beginning of time step, for bulk or a
+    ! single tracer
+    !
+    ! !ARGUMENTS:
+    type(bounds_type)         , intent(in)    :: bounds     
+    integer                   , intent(in)    :: num_nolakec          ! number of column non-lake points in column filter
+    integer                   , intent(in)    :: filter_nolakec(:)    ! column filter for non-lake points
+    integer                   , intent(in)    :: num_lakec            ! number of column lake points in column filter
+    integer                   , intent(in)    :: filter_lakec(:)      ! column filter for lake points
+    type(soilhydrology_type)  , intent(in)    :: soilhydrology_inst
+    class(waterstate_type)    , intent(inout) :: waterstate_inst
+    class(waterdiagnostic_type), intent(in)   :: waterdiagnostic_inst
+    class(waterbalance_type)  , intent(inout) :: waterbalance_inst
+    logical                   , intent(in)    :: use_aquifer_layer    ! whether an aquifer layer is used in this run
+    !
+    ! !LOCAL VARIABLES:
+    integer :: c, j, fc                  ! indices
+    !-----------------------------------------------------------------------
+
+    associate(                                               &
+         zi         => col%zi                           , & ! Input:  [real(r8) (:,:) ]  interface level below a "z" level (m)
+         zwt        => soilhydrology_inst%zwt_col       , & ! Input:  [real(r8) (:)   ]  water table depth (m)
+         aquifer_water_baseline => waterstate_inst%aquifer_water_baseline, & ! Input: [real(r8)] baseline value for water in the unconfined aquifer (wa_col) for this bulk / tracer (mm)
+         wa         => waterstate_inst%wa_col           , & ! Output: [real(r8) (:)   ]  water in the unconfined aquifer (mm)
+         begwb      => waterbalance_inst%begwb_col      , & ! Output: [real(r8) (:)   ]  water mass begining of the time step
+         h2osno_old => waterbalance_inst%h2osno_old_col   & ! Output: [real(r8) (:)   ]  snow water (mm H2O) at previous time step
+         )
+
+      if(use_aquifer_layer) then
+         do fc = 1, num_nolakec
+            c = filter_nolakec(fc)
+            if (col%hydrologically_active(c)) then
+               if(zwt(c) <= zi(c,nlevsoi)) then
+                  wa(c) = aquifer_water_baseline
+               end if
+            end if
+         end do
+      endif
+
+    call ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, &
+         waterstate_inst, waterdiagnostic_inst, &
+         subtract_dynbal_baselines = .false., &
+         water_mass = begwb(bounds%begc:bounds%endc))
+
+    call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
+         waterstate_inst, &
+         subtract_dynbal_baselines = .false., &
+         water_mass = begwb(bounds%begc:bounds%endc))
+
+    call waterstate_inst%CalculateTotalH2osno(bounds, num_nolakec, filter_nolakec, &
+         caller = 'BeginWaterBalanceSingle-nolake', &
+         h2osno_total = h2osno_old(bounds%begc:bounds%endc))
+    call waterstate_inst%CalculateTotalH2osno(bounds, num_lakec, filter_lakec, &
+         caller = 'BeginWaterBalanceSingle-lake', &
+         h2osno_total = h2osno_old(bounds%begc:bounds%endc))
+
+    end associate 
+
+  end subroutine BeginWaterBalanceSingle
+
+   !-----------------------------------------------------------------------
+   subroutine BalanceCheck( bounds, &
+        num_allc, filter_allc, &
+        atm2lnd_inst, solarabs_inst, waterflux_inst, waterstate_inst, &
+        waterdiagnosticbulk_inst, waterbalance_inst, wateratm2lnd_inst, &
+        surfalb_inst, energyflux_inst, canopystate_inst)
+     !
+     ! !DESCRIPTION:
+     ! This subroutine accumulates the numerical truncation errors of the water
+     ! and energy balance calculation. It is helpful to see the performance of
+     ! the process of integration.
+     !
+     ! The error for energy balance:
+     !
+     ! error = abs(Net radiation - change of internal energy - Sensible heat
+     !             - Latent heat)
+     !
+     ! The error for water balance:
+     !
+     ! error = abs(precipitation - change of water storage - evaporation - runoff)
+     !
+     ! !USES:
+     use clm_varcon        , only : spval
+     use clm_time_manager  , only : get_step_size_real, get_nstep
+     use clm_time_manager  , only : get_nstep_since_startup_or_lastDA_restart_or_pause
+     use CanopyStateType   , only : canopystate_type
+     use subgridAveMod
+     !
+     ! !ARGUMENTS:
+     type(bounds_type)     , intent(in)    :: bounds  
+     integer               , intent(in)    :: num_allc        ! number of columns in allc filter
+     integer               , intent(in)    :: filter_allc(:)  ! filter for all columns
+     type(atm2lnd_type)    , intent(in)    :: atm2lnd_inst
+     type(solarabs_type)   , intent(in)    :: solarabs_inst
+     class(waterflux_type) , intent(in)    :: waterflux_inst
+     class(waterstate_type), intent(in)    :: waterstate_inst
+     type(waterdiagnosticbulk_type), intent(in) :: waterdiagnosticbulk_inst
+     class(waterbalance_type), intent(inout) :: waterbalance_inst
+     class(wateratm2lnd_type) , intent(in) :: wateratm2lnd_inst
+     type(surfalb_type)    , intent(in)    :: surfalb_inst
+     type(energyflux_type) , intent(inout) :: energyflux_inst
+     type(canopystate_type), intent(inout) :: canopystate_inst
+     !
+     ! !LOCAL VARIABLES:
+     integer  :: p,c,l,g,fc                             ! indices
+     real(r8) :: dtime                                  ! land model time step (sec)
+     integer  :: nstep                                  ! time step number
+     integer  :: DAnstep                                ! time step number since last Data Assimilation (DA)
+     integer  :: indexp,indexc,indexl,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) :: h2osno_total(bounds%begc:bounds%endc)  ! total snow water [mm H2O]
+
+     real(r8) :: errh2o_max_val                         ! Maximum value of error in water conservation error  over all columns [mm H2O]
+     real(r8) :: errh2osno_max_val                      ! Maximum value of error in h2osno conservation error over all columns [kg m-2]
+     real(r8) :: errsol_max_val                         ! Maximum value of error in solar radiation conservation error over all columns [W m-2]
+     real(r8) :: errlon_max_val                         ! Maximum value of error in longwave radiation conservation error over all columns [W m-2]
+     real(r8) :: errseb_max_val                         ! Maximum value of error in surface energy conservation error over all columns [W m-2]
+     real(r8) :: errsoi_col_max_val                     ! Maximum value of column-level soil/lake energy conservation error over all columns [W m-2]
+
+     real(r8), parameter :: h2o_warning_thresh       = 1.e-9_r8                       ! Warning threshhold for error in errh2o and errh2osnow 
+     real(r8), parameter :: energy_warning_thresh    = 1.e-7_r8                       ! Warning threshhold for error in errsol, errsol, errseb, errlonv
+     real(r8), parameter :: error_thresh             = 1.e-5_r8                       ! Error threshhold for conservation error
+
+     !-----------------------------------------------------------------------
+
+     associate(                                                                   & 
+          forc_solad              =>    atm2lnd_inst%forc_solad_grc             , & ! Input:  [real(r8) (:,:) ]  direct beam radiation (vis=forc_sols , nir=forc_soll )
+          forc_solai              =>    atm2lnd_inst%forc_solai_grc             , & ! Input:  [real(r8) (:,:) ]  diffuse radiation     (vis=forc_solsd, nir=forc_solld)
+          forc_rain               =>    wateratm2lnd_inst%forc_rain_downscaled_col   , & ! Input:  [real(r8) (:)   ]  rain rate [mm/s]
+          forc_snow               =>    wateratm2lnd_inst%forc_snow_downscaled_col   , & ! Input:  [real(r8) (:)   ]  snow rate [mm/s]
+          forc_lwrad              =>    atm2lnd_inst%forc_lwrad_downscaled_col  , & ! Input:  [real(r8) (:)   ]  downward infrared (longwave) radiation (W/m**2)
+
+          h2osno_old              =>    waterbalance_inst%h2osno_old_col          , & ! Input:  [real(r8) (:)   ]  snow water (mm H2O) at previous time step
+          frac_sno_eff            =>    waterdiagnosticbulk_inst%frac_sno_eff_col        , & ! Input:  [real(r8) (:)   ]  effective snow fraction                 
+          frac_sno                =>    waterdiagnosticbulk_inst%frac_sno_col            , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by snow (0 to 1)
+          snow_depth              =>    waterdiagnosticbulk_inst%snow_depth_col          , & ! Input:  [real(r8) (:)   ]  snow height (m)                         
+          begwb                   =>    waterbalance_inst%begwb_col               , & ! Input:  [real(r8) (:)   ]  water mass begining of the time step    
+          errh2o                  =>    waterbalance_inst%errh2o_col              , & ! Output: [real(r8) (:)   ]  water conservation error (mm H2O)       
+          errh2osno               =>    waterbalance_inst%errh2osno_col           , & ! Output: [real(r8) (:)   ]  error in h2osno (kg m-2)                
+          endwb                   =>    waterbalance_inst%endwb_col               , & ! Output: [real(r8) (:)   ]  water mass end of the time step         
+          snow_sources            =>    waterbalance_inst%snow_sources_col         , & ! Output: [real(r8) (:)   ]  snow sources (mm H2O /s)
+          snow_sinks              =>    waterbalance_inst%snow_sinks_col           , & ! Output: [real(r8) (:)   ]  snow sinks (mm H2O /s)
+          qflx_liq_grnd_col       =>    waterflux_inst%qflx_liq_grnd_col        , & ! Input:  [real(r8) (:)   ]  liquid on ground after interception (mm H2O/s) [+]
+          qflx_snow_grnd_col      =>    waterflux_inst%qflx_snow_grnd_col       , & ! Input:  [real(r8) (:)   ]  snow on ground after interception (mm H2O/s) [+]
+          qflx_snwcp_liq          =>    waterflux_inst%qflx_snwcp_liq_col       , & ! Input:  [real(r8) (:)   ]  excess liquid h2o due to snow capping (outgoing) (mm H2O /s) [+]`
+          qflx_snwcp_ice          =>    waterflux_inst%qflx_snwcp_ice_col       , & ! Input:  [real(r8) (:)   ]  excess solid h2o due to snow capping (outgoing) (mm H2O /s) [+]`
+          qflx_snwcp_discarded_liq =>   waterflux_inst%qflx_snwcp_discarded_liq_col, & ! Input:  [real(r8) (:)   ]  excess liquid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) [+]`
+          qflx_snwcp_discarded_ice =>   waterflux_inst%qflx_snwcp_discarded_ice_col, & ! Input:  [real(r8) (:)   ]  excess solid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) [+]`
+          qflx_evap_tot           =>    waterflux_inst%qflx_evap_tot_col        , & ! Input:  [real(r8) (:)   ]  qflx_evap_soi + qflx_evap_can + qflx_tran_veg
+          qflx_soliddew_to_top_layer    => waterflux_inst%qflx_soliddew_to_top_layer_col   , & ! Input:  [real(r8) (:)   ]  rate of solid water deposited on top soil or snow layer (frost) (mm H2O /s) [+]
+          qflx_solidevap_from_top_layer => waterflux_inst%qflx_solidevap_from_top_layer_col, & ! Input:  [real(r8) (:)   ]  rate of ice evaporated from top soil or snow layer (sublimation) (mm H2O /s) [+]
+          qflx_liqevap_from_top_layer   => waterflux_inst%qflx_liqevap_from_top_layer_col  , & ! Input:  [real(r8) (:)   ]  rate of liquid water evaporated from top soil or snow layer (mm H2O/s) [+]
+          qflx_liqdew_to_top_layer      => waterflux_inst%qflx_liqdew_to_top_layer_col     , & ! Input:  [real(r8) (:)   ]  rate of liquid water deposited on top soil or snow layer (dew) (mm H2O /s) [+]
+          qflx_prec_grnd          =>    waterdiagnosticbulk_inst%qflx_prec_grnd_col, & ! Input:  [real(r8) (:)   ]  water onto ground including canopy runoff [kg/(m2 s)]
+          qflx_snow_h2osfc        =>    waterflux_inst%qflx_snow_h2osfc_col     , & ! Input:  [real(r8) (:)   ]  snow falling on surface water (mm/s)
+          qflx_h2osfc_to_ice      =>    waterflux_inst%qflx_h2osfc_to_ice_col   , & ! Input:  [real(r8) (:)   ]  conversion of h2osfc to ice             
+          qflx_drain_perched      =>    waterflux_inst%qflx_drain_perched_col   , & ! Input:  [real(r8) (:)   ]  sub-surface runoff (mm H2O /s)          
+          qflx_floodc             =>    waterflux_inst%qflx_floodc_col          , & ! Input:  [real(r8) (:)   ]  total runoff due to flooding            
+          qflx_snow_drain         =>    waterflux_inst%qflx_snow_drain_col      , & ! Input:  [real(r8) (:)   ]  drainage from snow pack                         
+          qflx_surf               =>    waterflux_inst%qflx_surf_col            , & ! Input:  [real(r8) (:)   ]  surface runoff (mm H2O /s)              
+          qflx_qrgwl              =>    waterflux_inst%qflx_qrgwl_col           , & ! Input:  [real(r8) (:)   ]  qflx_surf at glaciers, wetlands, lakes  
+          qflx_drain              =>    waterflux_inst%qflx_drain_col           , & ! Input:  [real(r8) (:)   ]  sub-surface runoff (mm H2O /s)          
+          qflx_ice_runoff_snwcp   =>    waterflux_inst%qflx_ice_runoff_snwcp_col, & ! Input:  [real(r8) (:)   ] solid runoff from snow capping (mm H2O /s)
+          qflx_ice_runoff_xs      =>    waterflux_inst%qflx_ice_runoff_xs_col   , & ! Input:  [real(r8) (:)   ] solid runoff from excess ice in soil (mm H2O /s)
+          qflx_sl_top_soil        =>    waterflux_inst%qflx_sl_top_soil_col     , & ! Input:  [real(r8) (:)   ]  liquid water + ice from layer above soil to top soil layer or sent to qflx_qrgwl (mm H2O/s)
+
+          qflx_sfc_irrig          =>    waterflux_inst%qflx_sfc_irrig_col       , & ! Input:  [real(r8) (:)   ]  irrigation flux (mm H2O /s)             
+          qflx_glcice_dyn_water_flux => waterflux_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)]  water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system)
+
+          eflx_lwrad_out          =>    energyflux_inst%eflx_lwrad_out_patch    , & ! Input:  [real(r8) (:)   ]  emitted infrared (longwave) radiation (W/m**2)
+          eflx_lwrad_net          =>    energyflux_inst%eflx_lwrad_net_patch    , & ! Input:  [real(r8) (:)   ]  net infrared (longwave) rad (W/m**2) [+ = to atm]
+          eflx_sh_tot             =>    energyflux_inst%eflx_sh_tot_patch       , & ! Input:  [real(r8) (:)   ]  total sensible heat flux (W/m**2) [+ to atm]
+          eflx_lh_tot             =>    energyflux_inst%eflx_lh_tot_patch       , & ! Input:  [real(r8) (:)   ]  total latent heat flux (W/m**2)  [+ to atm]
+          eflx_soil_grnd          =>    energyflux_inst%eflx_soil_grnd_patch    , & ! Input:  [real(r8) (:)   ]  soil heat flux (W/m**2) [+ = into soil] 
+          eflx_wasteheat_patch    =>    energyflux_inst%eflx_wasteheat_patch    , & ! Input:  [real(r8) (:)   ]  sensible heat flux from urban heating/cooling sources of waste heat (W/m**2)
+          eflx_heat_from_ac_patch =>    energyflux_inst%eflx_heat_from_ac_patch , & ! Input:  [real(r8) (:)   ]  sensible heat flux put back into canyon due to removal by AC (W/m**2)
+          eflx_traffic_patch      =>    energyflux_inst%eflx_traffic_patch      , & ! Input:  [real(r8) (:)   ]  traffic sensible heat flux (W/m**2)     
+          eflx_dynbal             =>    energyflux_inst%eflx_dynbal_grc         , & ! Input:  [real(r8) (:)   ]  energy conversion flux due to dynamic land cover change(W/m**2) [+ to atm]
+          errsoi_col              =>    energyflux_inst%errsoi_col              , & ! Output: [real(r8) (:)   ]  column-level soil/lake energy conservation error (W/m**2)
+          errsol                  =>    energyflux_inst%errsol_patch            , & ! Output: [real(r8) (:)   ]  solar radiation conservation error (W/m**2)
+          errseb                  =>    energyflux_inst%errseb_patch            , & ! Output: [real(r8) (:)   ]  surface energy conservation error (W/m**2)
+          errlon                  =>    energyflux_inst%errlon_patch            , & ! Output: [real(r8) (:)   ]  longwave radiation conservation error (W/m**2)
+
+          sabg_soil               =>    solarabs_inst%sabg_soil_patch           , & ! Input:  [real(r8) (:)   ]  solar radiation absorbed by soil (W/m**2)
+          sabg_snow               =>    solarabs_inst%sabg_snow_patch           , & ! Input:  [real(r8) (:)   ]  solar radiation absorbed by snow (W/m**2)
+          sabg_chk                =>    solarabs_inst%sabg_chk_patch            , & ! Input:  [real(r8) (:)   ]  sum of soil/snow using current fsno, for balance check
+          fsa                     =>    solarabs_inst%fsa_patch                 , & ! Input:  [real(r8) (:)   ]  solar radiation absorbed (total) (W/m**2)
+          fsr                     =>    solarabs_inst%fsr_patch                 , & ! Input:  [real(r8) (:)   ]  solar radiation reflected (W/m**2)      
+          sabv                    =>    solarabs_inst%sabv_patch                , & ! Input:  [real(r8) (:)   ]  solar radiation absorbed by vegetation (W/m**2)
+          sabg                    =>    solarabs_inst%sabg_patch                , & ! Input:  [real(r8) (:)   ]  solar radiation absorbed by ground (W/m**2)
+          
+          elai                    =>    canopystate_inst%elai_patch             , & ! Input:  [real(r8) (:,:)]  
+          esai                    =>    canopystate_inst%esai_patch             , & ! Input:  [real(r8) (:,:)]  
+
+          fabd                    =>    surfalb_inst%fabd_patch                 , & ! Input:  [real(r8) (:,:)]  flux absorbed by canopy per unit direct flux
+          fabi                    =>    surfalb_inst%fabi_patch                 , & ! Input:  [real(r8) (:,:)]  flux absorbed by canopy per unit indirect flux
+          albd                    =>    surfalb_inst%albd_patch                 , & ! Input:  [real(r8) (:,:)]  surface albedo (direct)
+          albi                    =>    surfalb_inst%albi_patch                 , & ! Input:  [real(r8) (:,:)]  surface albedo (diffuse)
+          ftdd                    =>    surfalb_inst%ftdd_patch                 , & ! Input:  [real(r8) (:,:)]  down direct flux below canopy per unit direct flux
+          ftid                    =>    surfalb_inst%ftid_patch                 , & ! Input:  [real(r8) (:,:)]  down diffuse flux below canopy per unit direct flux
+          ftii                    =>    surfalb_inst%ftii_patch                 , & ! Input:  [real(r8) (:,:)]  down diffuse flux below canopy per unit diffuse flux
+
+          netrad                  =>    energyflux_inst%netrad_patch              & ! Output: [real(r8) (:)   ]  net radiation (positive downward) (W/m**2)
+          )
+
+       ! Get step size and time step
+
+       nstep = get_nstep()
+       DAnstep = get_nstep_since_startup_or_lastDA_restart_or_pause()
+       dtime = get_step_size_real()
+
+       ! Determine column level incoming snow and rain
+       ! Assume no incident precipitation on urban wall columns (as in CanopyHydrologyMod.F90).
+
+       do c = bounds%begc,bounds%endc
+          g = col%gridcell(c)
+          l = col%landunit(c)       
+
+          if (col%itype(c) == icol_sunwall .or.  col%itype(c) == icol_shadewall) then
+             forc_rain_col(c) = 0.
+             forc_snow_col(c) = 0.
+          else
+             forc_rain_col(c) = forc_rain(c)
+             forc_snow_col(c) = forc_snow(c)
+          end if
+       end do
+
+       ! Water balance check
+
+       do c = bounds%begc, bounds%endc
+
+          ! add qflx_drain_perched and qflx_flood
+          if (col%active(c)) then
+
+             errh2o(c) = endwb(c) - begwb(c) &
+                  - (forc_rain_col(c)        &
+                  + forc_snow_col(c)         &
+                  + qflx_floodc(c)           &
+                  + qflx_sfc_irrig(c)        &
+                  + qflx_glcice_dyn_water_flux(c) &
+                  - qflx_evap_tot(c)         &
+                  - qflx_surf(c)             &
+                  - qflx_qrgwl(c)            &
+                  - qflx_drain(c)            &
+                  - qflx_drain_perched(c)    &
+                  - qflx_ice_runoff_snwcp(c) &
+                  - qflx_ice_runoff_xs(c)    &
+                  - qflx_snwcp_discarded_liq(c) &
+                  - qflx_snwcp_discarded_ice(c)) * dtime
+
+          else
+
+             errh2o(c) = 0.0_r8
+
+          end if
+
+       end do
+       
+       errh2o_max_val = maxval(abs(errh2o(bounds%begc:bounds%endc)))
+
+       if (errh2o_max_val > h2o_warning_thresh) then
+
+           indexc = maxloc( abs(errh2o(bounds%begc:bounds%endc)), 1 ) + bounds%begc -1
+           write(iulog,*)'WARNING:  water balance error ',&
+             ' nstep= ',nstep, &
+             ' local indexc= ',indexc,&
+           ! ' global indexc= ',GetGlobalIndex(decomp_index=indexc, clmlevel=namec), &
+             ' errh2o= ',errh2o(indexc)
+
+         ! DART SourceMod edits    
+         if ((errh2o_max_val > error_thresh) .and. (.not. is_first_restart_step()) ) then    
+              write(iulog,*)'clm urban model is stopping - error is greater than 1e-5 (mm)'
+              write(iulog,*)'nstep                     = ',nstep
+              write(iulog,*)'errh2o                    = ',errh2o(indexc)
+              write(iulog,*)'forc_rain                 = ',forc_rain_col(indexc)*dtime
+              write(iulog,*)'forc_snow                 = ',forc_snow_col(indexc)*dtime
+              write(iulog,*)'endwb                     = ',endwb(indexc)
+              write(iulog,*)'begwb                     = ',begwb(indexc)
+
+              write(iulog,*)'qflx_evap_tot             = ',qflx_evap_tot(indexc)*dtime
+              write(iulog,*)'qflx_sfc_irrig            = ',qflx_sfc_irrig(indexc)*dtime
+              write(iulog,*)'qflx_surf                 = ',qflx_surf(indexc)*dtime
+              write(iulog,*)'qflx_qrgwl                = ',qflx_qrgwl(indexc)*dtime
+              write(iulog,*)'qflx_drain                = ',qflx_drain(indexc)*dtime
+
+              write(iulog,*)'qflx_ice_runoff_snwcp     = ',qflx_ice_runoff_snwcp(indexc)*dtime
+              write(iulog,*)'qflx_ice_runoff_xs        = ',qflx_ice_runoff_xs(indexc)*dtime
+
+              write(iulog,*)'qflx_snwcp_discarded_ice  = ',qflx_snwcp_discarded_ice(indexc)*dtime
+              write(iulog,*)'qflx_snwcp_discarded_liq  = ',qflx_snwcp_discarded_liq(indexc)*dtime
+              write(iulog,*)'deltawb                   = ',endwb(indexc)-begwb(indexc)
+              write(iulog,*)'deltawb/dtime             = ',(endwb(indexc)-begwb(indexc))/dtime
+
+              if (.not.(col%itype(indexc) == icol_roof .or. &
+                   col%itype(indexc) == icol_road_imperv .or. &
+                   col%itype(indexc) == icol_road_perv)) then
+                   write(iulog,*)'qflx_drain_perched         = ',qflx_drain_perched(indexc)*dtime
+                   write(iulog,*)'qflx_flood                 = ',qflx_floodc(indexc)*dtime
+                   write(iulog,*)'qflx_glcice_dyn_water_flux = ', qflx_glcice_dyn_water_flux(indexc)*dtime
+              end if
+              
+              write(iulog,*)'clm model is stopping'
+              call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
+         end if
+       
+       end if
+       
+
+       ! Snow balance check
+
+       call waterstate_inst%CalculateTotalH2osno(bounds, num_allc, filter_allc, &
+            caller = 'BalanceCheck', &
+            h2osno_total = h2osno_total(bounds%begc:bounds%endc))
+
+       do c = bounds%begc,bounds%endc
+          if (col%active(c)) then
+             g = col%gridcell(c)
+             l = col%landunit(c)
+
+             ! As defined here, snow_sources - snow_sinks will equal the change in h2osno at 
+             ! any given time step but only if there is at least one snow layer.  h2osno 
+             ! also includes snow that is part of the soil column (an initial snow layer is 
+             ! only created if h2osno > 10mm).
+
+             if (col%snl(c) < 0) then
+                snow_sources(c) = qflx_prec_grnd(c) + qflx_soliddew_to_top_layer(c) &
+                     + qflx_liqdew_to_top_layer(c)
+                snow_sinks(c)  = qflx_solidevap_from_top_layer(c) + qflx_liqevap_from_top_layer(c) &
+                     + qflx_snow_drain(c) + qflx_snwcp_ice(c) + qflx_snwcp_liq(c) &
+                     + qflx_snwcp_discarded_ice(c) + qflx_snwcp_discarded_liq(c) &
+                     + qflx_sl_top_soil(c)
+
+                if (lun%itype(l) == istdlak) then 
+                   snow_sources(c) = qflx_snow_grnd_col(c) &
+                        + frac_sno_eff(c) * (qflx_liq_grnd_col(c) &
+                        +  qflx_soliddew_to_top_layer(c) + qflx_liqdew_to_top_layer(c) ) 
+                   snow_sinks(c)   = frac_sno_eff(c) * (qflx_solidevap_from_top_layer(c) &
+                        + qflx_liqevap_from_top_layer(c) ) + qflx_snwcp_ice(c) + qflx_snwcp_liq(c)  &
+                        + qflx_snwcp_discarded_ice(c) + qflx_snwcp_discarded_liq(c)  &
+                        + qflx_snow_drain(c)  + qflx_sl_top_soil(c)
+                endif
+
+                 if (col%itype(c) == icol_road_perv .or. lun%itype(l) == istsoil .or. &
+                      lun%itype(l) == istcrop .or. lun%itype(l) == istwet .or. &
+                      lun%itype(l) == istice_mec) then
+                   snow_sources(c) = (qflx_snow_grnd_col(c) - qflx_snow_h2osfc(c) ) &
+                          + frac_sno_eff(c) * (qflx_liq_grnd_col(c) &
+                          + qflx_soliddew_to_top_layer(c) + qflx_liqdew_to_top_layer(c) ) &
+                          + qflx_h2osfc_to_ice(c)
+                   snow_sinks(c) = frac_sno_eff(c) * (qflx_solidevap_from_top_layer(c) &
+                          + qflx_liqevap_from_top_layer(c)) + qflx_snwcp_ice(c) + qflx_snwcp_liq(c) &
+                          + qflx_snwcp_discarded_ice(c) + qflx_snwcp_discarded_liq(c) &
+                          + qflx_snow_drain(c) + qflx_sl_top_soil(c)
+                endif
+
+                errh2osno(c) = (h2osno_total(c) - h2osno_old(c)) - (snow_sources(c) - snow_sinks(c)) * dtime
+             else
+                snow_sources(c) = 0._r8
+                snow_sinks(c) = 0._r8
+                errh2osno(c) = 0._r8
+             end if
+          else
+             errh2osno(c) = 0._r8
+          end if
+       end do
+
+
+       errh2osno_max_val = maxval( abs(errh2osno(bounds%begc:bounds%endc)))
+       
+       if (errh2osno_max_val > h2o_warning_thresh) then
+            indexc = maxloc( abs(errh2osno(bounds%begc:bounds%endc)), 1) + bounds%begc -1
+            write(iulog,*)'WARNING:  snow balance error '
+            write(iulog,*)'nstep= ',nstep, &
+                 ' local indexc= ',indexc, &
+                 ! ' global indexc= ',GetGlobalIndex(decomp_index=indexc, clmlevel=namec), &
+                 ' col%itype= ',col%itype(indexc), &
+                 ' lun%itype= ',lun%itype(col%landunit(indexc)), &
+                 ' errh2osno= ',errh2osno(indexc)
+
+            ! DART SourceMod edits
+            if ((errh2osno_max_val > error_thresh) .and. (.not. is_first_restart_step()) ) then
+                 write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
+                 write(iulog,*)'nstep              = ',nstep
+                 write(iulog,*)'errh2osno          = ',errh2osno(indexc)
+                 write(iulog,*)'snl                = ',col%snl(indexc)
+                 write(iulog,*)'snow_depth         = ',snow_depth(indexc)
+                 write(iulog,*)'frac_sno_eff       = ',frac_sno_eff(indexc)
+                 write(iulog,*)'h2osno             = ',h2osno_total(indexc)
+                 write(iulog,*)'h2osno_old         = ',h2osno_old(indexc)
+                 write(iulog,*)'snow_sources       = ',snow_sources(indexc)*dtime
+                 write(iulog,*)'snow_sinks         = ',snow_sinks(indexc)*dtime
+                 write(iulog,*)'qflx_prec_grnd     = ',qflx_prec_grnd(indexc)*dtime
+                 write(iulog,*)'qflx_snow_grnd_col = ',qflx_snow_grnd_col(indexc)*dtime
+                 write(iulog,*)'qflx_liq_grnd_col  = ',qflx_liq_grnd_col(indexc)*dtime
+                 write(iulog,*)'qflx_solidevap_from_top_layer = ',qflx_solidevap_from_top_layer(indexc)*dtime
+                 write(iulog,*)'qflx_snow_drain    = ',qflx_snow_drain(indexc)*dtime
+                 write(iulog,*)'qflx_liqevap_from_top_layer = ',qflx_liqevap_from_top_layer(indexc)*dtime
+                 write(iulog,*)'qflx_soliddew_to_top_layer  = ',qflx_soliddew_to_top_layer(indexc)*dtime
+                 write(iulog,*)'qflx_liqdew_to_top_layer    = ',qflx_liqdew_to_top_layer(indexc)*dtime
+                 write(iulog,*)'qflx_snwcp_ice     = ',qflx_snwcp_ice(indexc)*dtime
+                 write(iulog,*)'qflx_snwcp_liq     = ',qflx_snwcp_liq(indexc)*dtime
+                 write(iulog,*)'qflx_snwcp_discarded_ice = ',qflx_snwcp_discarded_ice(indexc)*dtime
+                 write(iulog,*)'qflx_snwcp_discarded_liq = ',qflx_snwcp_discarded_liq(indexc)*dtime
+                 write(iulog,*)'qflx_sl_top_soil   = ',qflx_sl_top_soil(indexc)*dtime
+                 write(iulog,*)'clm model is stopping'
+                 call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
+            end if
+
+       end if
+
+       ! Energy balance checks
+
+       do p = bounds%begp, bounds%endp
+          if (patch%active(p)) then
+             c = patch%column(p)
+             l = patch%landunit(p)
+             g = patch%gridcell(p)
+
+             ! Solar radiation energy balance
+             ! Do not do this check for an urban patch since it will not balance on a per-column
+             ! level because of interactions between columns and since a separate check is done
+             ! in the urban radiation module
+             if (.not. lun%urbpoi(l)) then
+                errsol(p) = fsa(p) + fsr(p) &
+                     - (forc_solad(g,1) + forc_solad(g,2) + forc_solai(g,1) + forc_solai(g,2))
+             else
+                errsol(p) = spval
+             end if
+
+             ! Longwave radiation energy balance
+             ! Do not do this check for an urban patch since it will not balance on a per-column
+             ! level because of interactions between columns and since a separate check is done
+             ! in the urban radiation module
+             if (.not. lun%urbpoi(l)) then
+                errlon(p) = eflx_lwrad_out(p) - eflx_lwrad_net(p) - forc_lwrad(c)
+             else
+                errlon(p) = spval
+             end if
+
+             ! Surface energy balance
+             ! Changed to using (eflx_lwrad_net) here instead of (forc_lwrad - eflx_lwrad_out) because
+             ! there are longwave interactions between urban columns (and therefore patches). 
+             ! For surfaces other than urban, (eflx_lwrad_net) equals (forc_lwrad - eflx_lwrad_out),
+             ! and a separate check is done above for these terms.
+
+             if (.not. lun%urbpoi(l)) then
+                errseb(p) = sabv(p) + sabg_chk(p) + forc_lwrad(c) - eflx_lwrad_out(p) &
+                     - eflx_sh_tot(p) - eflx_lh_tot(p) - eflx_soil_grnd(p)
+             else
+                errseb(p) = sabv(p) + sabg(p) &
+                     - eflx_lwrad_net(p) &
+                     - eflx_sh_tot(p) - eflx_lh_tot(p) - eflx_soil_grnd(p) &
+                     + eflx_wasteheat_patch(p) + eflx_heat_from_ac_patch(p) + eflx_traffic_patch(p)
+             end if
+             !TODO MV - move this calculation to a better place - does not belong in BalanceCheck 
+             netrad(p) = fsa(p) - eflx_lwrad_net(p)
+          else
+             errsol(p) = 0._r8
+             errlon(p) = 0._r8
+             errseb(p) = 0._r8
+          end if
+       end do
+
+       ! Solar radiation energy balance check
+
+       errsol_max_val = maxval( abs(errsol(bounds%begp:bounds%endp)), mask = (errsol(bounds%begp:bounds%endp) /= spval) ) 
+
+       ! DART SourceMods
+
+       if  ((errsol_max_val > energy_warning_thresh) .and. (.not. is_first_restart_step()) ) then
+
+           indexp = maxloc( abs(errsol(bounds%begp:bounds%endp)), 1 , mask = (errsol(bounds%begp:bounds%endp) /= spval) ) + bounds%begp -1
+           indexg = patch%gridcell(indexp)
+           write(iulog,*)'WARNING:: BalanceCheck, solar radiation balance error (W/m2)'
+           write(iulog,*)'nstep         = ',nstep
+           write(iulog,*)'errsol        = ',errsol(indexp)
+          
+
+           if (errsol_max_val > error_thresh) then
+               write(iulog,*)'clm model is stopping - error is greater than 1e-5 (W/m2)'
+               write(iulog,*)'fsa           = ',fsa(indexp)
+               write(iulog,*)'fsr           = ',fsr(indexp)
+               write(iulog,*)'forc_solad(1) = ',forc_solad(indexg,1)
+               write(iulog,*)'forc_solad(2) = ',forc_solad(indexg,2)
+               write(iulog,*)'forc_solai(1) = ',forc_solai(indexg,1)
+               write(iulog,*)'forc_solai(2) = ',forc_solai(indexg,2)
+               write(iulog,*)'forc_tot      = ',forc_solad(indexg,1)+forc_solad(indexg,2) &
+                  +forc_solai(indexg,1)+forc_solai(indexg,2)
+               write(iulog,*)'clm model is stopping'
+               call endrun(decomp_index=indexp, clmlevel=namep, msg=errmsg(sourcefile, __LINE__))
+           end if
+
+       end if
+       
+       ! Longwave radiation energy balance check
+
+       errlon_max_val = maxval( abs(errlon(bounds%begp:bounds%endp)), mask = (errlon(bounds%begp:bounds%endp) /= spval) )
+
+       ! DART SourceMods
+
+       if ((errlon_max_val > energy_warning_thresh) .and. (.not. is_first_restart_step()) ) then        
+            indexp = maxloc( abs(errlon(bounds%begp:bounds%endp)), 1 , mask = (errlon(bounds%begp:bounds%endp) /= spval) ) + bounds%begp -1
+            write(iulog,*)'indexp         = ',indexp
+            write(iulog,*)'WARNING: BalanceCheck: longwave energy balance error (W/m2)'
+            write(iulog,*)'nstep        = ',nstep
+            write(iulog,*)'errlon       = ',errlon(indexp)
+            if (errlon_max_val > error_thresh ) then
+                 write(iulog,*)'clm model is stopping - error is greater than 1e-5 (W/m2)'
+                 call endrun(decomp_index=indexp, clmlevel=namep, msg=errmsg(sourcefile, __LINE__))
+            end if
+       end if
+
+       ! Surface energy balance check
+
+       errseb_max_val = maxval( abs(errseb(bounds%begp:bounds%endp)))
+ 
+       ! DART SourceMods
+
+       if ((errseb_max_val > energy_warning_thresh) .and. (.not. is_first_restart_step()) ) then
+
+           indexp = maxloc( abs(errseb(bounds%begp:bounds%endp)), 1 ) + bounds%begp -1
+           indexc = patch%column(indexp)
+           indexg = patch%gridcell(indexp)
+
+           write(iulog,*)'WARNING: BalanceCheck: surface flux energy balance error (W/m2)'
+           write(iulog,*)'nstep          = ' ,nstep
+           write(iulog,*)'errseb         = ' ,errseb(indexp)
+
+           if ( errseb_max_val > error_thresh ) then
+              write(iulog,*)'clm model is stopping - error is greater than 1e-5 (W/m2)'
+              write(iulog,*)'sabv           = ' ,sabv(indexp)
+              write(iulog,*)'sabg           = ' ,sabg(indexp), ((1._r8- frac_sno(indexc))*sabg_soil(indexp) + &
+                   frac_sno(indexc)*sabg_snow(indexp)),sabg_chk(indexp)
+              write(iulog,*)'forc_tot      = '  ,forc_solad(indexg,1) + forc_solad(indexg,2) + &
+                   forc_solai(indexg,1) + forc_solai(indexg,2)
+
+              write(iulog,*)'eflx_lwrad_net = ' ,eflx_lwrad_net(indexp)
+              write(iulog,*)'eflx_sh_tot    = ' ,eflx_sh_tot(indexp)
+              write(iulog,*)'eflx_lh_tot    = ' ,eflx_lh_tot(indexp)
+              write(iulog,*)'eflx_soil_grnd = ' ,eflx_soil_grnd(indexp)
+              write(iulog,*)'fsa fsr = '        ,fsa(indexp),    fsr(indexp)
+              write(iulog,*)'fabd fabi = '      ,fabd(indexp,:), fabi(indexp,:)
+              write(iulog,*)'albd albi = '      ,albd(indexp,:), albi(indexp,:)
+              write(iulog,*)'ftii ftdd ftid = ' ,ftii(indexp,:), ftdd(indexp,:),ftid(indexp,:)
+              write(iulog,*)'elai esai = '      ,elai(indexp),   esai(indexp)
+              write(iulog,*)'clm model is stopping'
+              call endrun(decomp_index=indexp, clmlevel=namep, msg=errmsg(sourcefile, __LINE__))
+           end if
+
+       end if
+
+       ! Soil energy balance check
+
+       errsoi_col_max_val  =  maxval( abs(errsoi_col(bounds%begc:bounds%endc)) , mask = col%active(bounds%begc:bounds%endc))
+
+       if (errsoi_col_max_val > 1.0e-5_r8 ) then
+           indexc = maxloc( abs(errsoi_col(bounds%begc:bounds%endc)), 1 , mask = col%active(bounds%begc:bounds%endc) ) + bounds%begc -1
+           write(iulog,*)'WARNING: BalanceCheck: soil balance error (W/m2)'
+           write(iulog,*)'nstep         = ',nstep
+           write(iulog,*)'errsoi_col    = ',errsoi_col(indexc)
+
+       ! DART SourceMods
+
+           if ((errsoi_col_max_val > 1.e-4_r8) .and.  (.not. is_first_restart_step()) ) then
+              write(iulog,*)'clm model is stopping'
+              call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
+           end if
+       end if 
+       
+     end associate
+
+   end subroutine BalanceCheck
+
+end module BalanceCheckMod
diff --git a/models/clm/shell_scripts/cesm2_2/CLM5_hybrid_freerun b/models/clm/shell_scripts/cesm2_2/CLM5_hybrid_freerun
index 465cd5142..25934e731 100755
--- a/models/clm/shell_scripts/cesm2_2/CLM5_hybrid_freerun
+++ b/models/clm/shell_scripts/cesm2_2/CLM5_hybrid_freerun
@@ -253,12 +253,20 @@ endif
 ./xmlchange INFO_DBUG=0
 
 # ==============================================================================
-# Use the stream template that has all CAM reanalysis years available.
-
+# If the experiment only spans one year, copy a stream template for a single year
+# otherwise, use 'all' the years.
+
+if (${stream_year_first} == ${stream_year_last}) then
+   set STREAMFILE_SOLAR        = datm.streams.txt.CPLHISTForcing.Solar_single_year
+   set STREAMFILE_STATE1HR     = datm.streams.txt.CPLHISTForcing.State1hr_single_year
+   set STREAMFILE_STATE3HR     = datm.streams.txt.CPLHISTForcing.State3hr_single_year
+   set STREAMFILE_NONSOLARFLUX = datm.streams.txt.CPLHISTForcing.nonSolarFlux_single_year
+else
    set STREAMFILE_SOLAR        = datm.streams.txt.CPLHISTForcing.Solar_complete
    set STREAMFILE_STATE1HR     = datm.streams.txt.CPLHISTForcing.State1hr_complete
    set STREAMFILE_STATE3HR     = datm.streams.txt.CPLHISTForcing.State3hr_complete
    set STREAMFILE_NONSOLARFLUX = datm.streams.txt.CPLHISTForcing.nonSolarFlux_complete
+endif
 
 # ==============================================================================
 # Modify namelist templates for each instance.
@@ -553,7 +561,7 @@ if ( ${use_SourceMods} == TRUE ) then
    endif
 endif
 
-${BUILD_WRAPPER} ./case.build || exit 9
+./case.build || exit 9
 
 # ==============================================================================
 # What to do next
diff --git a/models/clm/shell_scripts/cesm2_2/CLM5_setup_pmo b/models/clm/shell_scripts/cesm2_2/CLM5_setup_pmo
index f0156f93c..55027b5ab 100755
--- a/models/clm/shell_scripts/cesm2_2/CLM5_setup_pmo
+++ b/models/clm/shell_scripts/cesm2_2/CLM5_setup_pmo
@@ -224,12 +224,20 @@ endif
 ./xmlchange INFO_DBUG=0
 
 # ==============================================================================
-# Use the stream template that has all CAM reanalysis years available.
-
+# If the experiment only spans one year, copy a stream template for a single year
+# otherwise, use 'all' the years.
+
+if (${stream_year_first} == ${stream_year_last}) then
+   set STREAMFILE_SOLAR        = datm.streams.txt.CPLHISTForcing.Solar_single_year
+   set STREAMFILE_STATE1HR     = datm.streams.txt.CPLHISTForcing.State1hr_single_year
+   set STREAMFILE_STATE3HR     = datm.streams.txt.CPLHISTForcing.State3hr_single_year
+   set STREAMFILE_NONSOLARFLUX = datm.streams.txt.CPLHISTForcing.nonSolarFlux_single_year
+else
    set STREAMFILE_SOLAR        = datm.streams.txt.CPLHISTForcing.Solar_complete
    set STREAMFILE_STATE1HR     = datm.streams.txt.CPLHISTForcing.State1hr_complete
    set STREAMFILE_STATE3HR     = datm.streams.txt.CPLHISTForcing.State3hr_complete
    set STREAMFILE_NONSOLARFLUX = datm.streams.txt.CPLHISTForcing.nonSolarFlux_complete
+endif
 
 # ==============================================================================
 # Modify namelist templates for each instance.
@@ -528,7 +536,7 @@ if ( ${use_SourceMods} == TRUE ) then
    endif
 endif
 
-${BUILD_WRAPPER} ./case.build || exit 7
+./case.build || exit 7
 
 # ==============================================================================
 # What to do next
diff --git a/models/clm/shell_scripts/cesm2_2/CLM5_startup_freerun b/models/clm/shell_scripts/cesm2_2/CLM5_startup_freerun
index 5502e80ad..5b1343274 100755
--- a/models/clm/shell_scripts/cesm2_2/CLM5_startup_freerun
+++ b/models/clm/shell_scripts/cesm2_2/CLM5_startup_freerun
@@ -364,7 +364,7 @@ if ( ${use_SourceMods} == TRUE ) then
    endif
 endif
 
-${BUILD_WRAPPER} ./case.build || exit 7
+./case.build || exit 7
 
 # ==============================================================================
 # What to do next