diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90
index 65f102a548..69f7cdb167 100644
--- a/src/biogeophys/BalanceCheckMod.F90
+++ b/src/biogeophys/BalanceCheckMod.F90
@@ -216,7 +216,7 @@ subroutine BeginWaterBalanceSingle(bounds, &
 
     call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
          waterstate_inst, lakestate_inst, &
-         subtract_dynbal_baselines = .false., &
+         add_lake_water_and_subtract_dynbal_baselines = .false., &
          water_mass = begwb(bounds%begc:bounds%endc))
 
     call waterstate_inst%CalculateTotalH2osno(bounds, num_nolakec, filter_nolakec, &
diff --git a/src/biogeophys/LakeHydrologyMod.F90 b/src/biogeophys/LakeHydrologyMod.F90
index c078e958e1..f6f83d8956 100644
--- a/src/biogeophys/LakeHydrologyMod.F90
+++ b/src/biogeophys/LakeHydrologyMod.F90
@@ -650,7 +650,7 @@ subroutine LakeHydrology(bounds, &
 
     call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
          b_waterstate_inst, lakestate_inst, &
-         subtract_dynbal_baselines = .false., &
+         add_lake_water_and_subtract_dynbal_baselines = .false., &
          water_mass = endwb(bounds%begc:bounds%endc))
 
     do j = 1, nlevgrnd
diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90
index 3bec595bf0..cd7316af81 100644
--- a/src/biogeophys/TotalWaterAndHeatMod.F90
+++ b/src/biogeophys/TotalWaterAndHeatMod.F90
@@ -143,7 +143,7 @@ end subroutine ComputeWaterMassNonLake
   !-----------------------------------------------------------------------
   subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
        waterstate_inst, lakestate_inst, &
-       subtract_dynbal_baselines, &
+       add_lake_water_and_subtract_dynbal_baselines, &
        water_mass)
     !
     ! !DESCRIPTION:
@@ -159,9 +159,12 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
     class(waterstate_type)   , intent(in)    :: waterstate_inst
     type(lakestate_type)     , intent(in)    :: lakestate_inst
 
+    ! Whether to (1) add lake water/ice to total accounting, and (2) subtract
+    ! dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
+    !
     ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658
     ! is resolved, remove this argument, always assuming it's true.
-    logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
+    logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines
 
     real(r8)                 , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2)
     !
@@ -181,7 +184,7 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
          filter_lakec = filter_lakec, &
          waterstate_inst = waterstate_inst, &
          lakestate_inst = lakestate_inst, &         
-         subtract_dynbal_baselines = subtract_dynbal_baselines, &
+         add_lake_water_and_subtract_dynbal_baselines = add_lake_water_and_subtract_dynbal_baselines, &
          liquid_mass = liquid_mass(bounds%begc:bounds%endc), &
          ice_mass = ice_mass(bounds%begc:bounds%endc))
 
@@ -387,8 +390,8 @@ end subroutine AccumulateSoilLiqIceMassNonLake
 
   !-----------------------------------------------------------------------
   subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
-       waterstate_inst, lakestate_inst,&
-       subtract_dynbal_baselines, &
+       waterstate_inst, lakestate_inst, &
+       add_lake_water_and_subtract_dynbal_baselines, &
        liquid_mass, ice_mass)
     !
     ! !DESCRIPTION:
@@ -407,10 +410,12 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
     class(waterstate_type), intent(in)    :: waterstate_inst
     type(lakestate_type)  , intent(in)    :: lakestate_inst
 
-
+    ! Whether to (1) add lake water/ice to total accounting, and (2) subtract
+    ! dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
+    !
     ! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658
     ! is resolved, remove this argument, always assuming it's true.
-    logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
+    logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines
 
     real(r8)              , intent(inout) :: liquid_mass( bounds%begc: ) ! computed liquid water mass (kg m-2)
     real(r8)              , intent(inout) :: ice_mass( bounds%begc: )    ! computed ice mass (kg m-2)
@@ -439,6 +444,34 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
        ice_mass(c) = 0._r8
     end do
 
+    ! ------------------------------------------------------------------------
+    ! Start with some large terms that will often cancel (the negative baselines and the
+    ! lake water content): In order to maintain precision in the other terms, it can help if
+    ! we deal with these large, often-canceling terms first. (If we accumulated some
+    ! small terms, then added a big term and then subtracted a big term, we would have
+    ! lost some precision in the small terms.)
+    ! ------------------------------------------------------------------------
+
+    if (add_lake_water_and_subtract_dynbal_baselines) then
+       ! Subtract baselines set in initialization
+       do fc = 1, num_lakec
+          c = filter_lakec(fc)
+          liquid_mass(c) = liquid_mass(c) - dynbal_baseline_liq(c)
+          ice_mass(c) = ice_mass(c) - dynbal_baseline_ice(c)
+       end do
+
+       ! Lake water content
+       call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, &
+            lakestate_inst, &
+            tracer_ratio = waterstate_inst%info%get_ratio(), &
+            liquid_mass = liquid_mass(bounds%begc:bounds%endc), &
+            ice_mass = ice_mass(bounds%begc:bounds%endc))
+    end if
+
+    ! ------------------------------------------------------------------------
+    ! Now add some other terms
+    ! ------------------------------------------------------------------------
+
     ! Snow water content
     do fc = 1, num_lakec
        c = filter_lakec(fc)
@@ -449,11 +482,7 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
           ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j)
        end do
     end do
-    
-    ! Lake water content
-    call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, &
-       lakestate_inst, liquid_mass, ice_mass)
-    
+
     ! Soil water content of the soil under the lake
     do j = 1, nlevgrnd
        do fc = 1, num_lakec
@@ -463,22 +492,13 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
        end do
     end do
 
-    if (subtract_dynbal_baselines) then
-       ! Subtract baselines set in initialization
-       do fc = 1, num_lakec
-          c = filter_lakec(fc)
-          liquid_mass(c) = liquid_mass(c) - dynbal_baseline_liq(c)
-          ice_mass(c) = ice_mass(c) - dynbal_baseline_ice(c)
-       end do
-    end if
-
     end associate
 
   end subroutine ComputeLiqIceMassLake
 
     !-----------------------------------------------------------------------
   subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, &
-       lakestate_inst, liquid_mass, ice_mass)
+       lakestate_inst, tracer_ratio, liquid_mass, ice_mass)
     !
     ! !DESCRIPTION:
     ! Accumulate lake water mass of lake columns, separated into liquid and ice.
@@ -493,6 +513,7 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, &
     integer                , intent(in)    :: num_c                       ! number of column points in column filter (should not include lake points; can be a subset of the no-lake filter)
     integer                , intent(in)    :: filter_c(:)                 ! column filter (should not include lake points; can be a subset of the no-lake filter)
     type(lakestate_type)   , intent(in)    :: lakestate_inst
+    real(r8)               , intent(in)    :: tracer_ratio                ! for water tracers, standard ratio of this tracer to bulk water (1 for bulk water)
     real(r8)               , intent(inout) :: liquid_mass( bounds%begc: ) ! accumulated liquid mass (kg m-2)
     real(r8)               , intent(inout) :: ice_mass( bounds%begc: )    ! accumulated ice mass (kg m-2)
     !
@@ -516,12 +537,14 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, &
     do j = 1, nlevlak
        do fc = 1, num_c
           c = filter_c(fc)
-          ! calculate lake liq and ice content per lake layer first
-          h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j))
+          ! Lake water volume isn't tracked explicitly, so we calculate it from lake
+          ! depth. Because it isn't tracked explicitly, we also don't have any water
+          ! tracer information, so we assume a fixed, standard tracer ratio.
+          h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) * tracer_ratio
 
           ! use water density rather than ice density because lake layer depths are not
           ! adjusted when the layer freezes
-          h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j)
+          h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) * tracer_ratio
           
           liquid_mass(c) = liquid_mass(c) + h2olak_liq
           ice_mass(c) = ice_mass(c) + h2olak_ice
@@ -936,6 +959,42 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, &
        latent_heat_liquid(c) = 0._r8
     end do
 
+    ! ------------------------------------------------------------------------
+    ! Start with some large terms that will often cancel (the negative baselines and the
+    ! lake water content): In order to maintain precision in the other terms, it can help if
+    ! we deal with these large, often-canceling terms first. (If we accumulated some
+    ! small terms, then added a big term and then subtracted a big term, we would have
+    ! lost some precision in the small terms.)
+    ! ------------------------------------------------------------------------
+
+    ! Subtract baselines set in initialization
+    !
+    ! NOTE(wjs, 2019-03-01) I haven't given enough thought to how (if at all) we should
+    ! correct for heat_liquid and cv_liquid, which are used to determine the weighted
+    ! average liquid water temperature. For example, if we're subtracting out a baseline
+    ! water amount because a particular water state is fictitious, we probably shouldn't
+    ! include that particular state when determining the weighted average temperature of
+    ! liquid water. And conversely, if we're adding a state via these baselines, should
+    ! we also add some water temperature of that state? The tricky thing here is what to
+    ! do when we end up subtracting water due to the baselines, so for now I'm simply not
+    ! trying to account for the temperature of these baselines at all. This amounts to
+    ! assuming that the baselines that we add / subtract are at the average temperature
+    ! of the real liquid water in the column.
+    do fc = 1, num_lakec
+       c = filter_lakec(fc)
+       heat(c) = -dynbal_baseline_heat(c)
+    end do
+
+    ! Lake water heat content
+    ! Note that we do NOT accumulate heat_liquid and cv_liquid in this call. See the
+    ! comments at the top of AccumulateHeatLake for rationale.
+    call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, &
+       heat)
+
+    ! ------------------------------------------------------------------------
+    ! Now add some other terms
+    ! ------------------------------------------------------------------------
+
     ! Snow heat content
     do fc = 1, num_lakec
        c = filter_lakec(fc)
@@ -976,41 +1035,16 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, &
     
     do fc = 1, num_lakec
        c = filter_lakec(fc)
-       heat(c) = heat_dry_mass(c) + heat_ice(c) + heat_liquid(c) + latent_heat_liquid(c)
+       heat(c) = heat(c) + heat_dry_mass(c) + heat_ice(c) + heat_liquid(c) + latent_heat_liquid(c)
     end do
     
-    ! Lake water heat content
-    ! Note that we do NOT accumulate heat_liquid and cv_liquid in this call. See the
-    ! comments at the top of AccumulateHeatLake for rationale.
-    call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, &
-       heat)
-
-    ! Subtract baselines set in initialization
-    !
-    !
-    ! NOTE(wjs, 2019-03-01) I haven't given enough thought to how (if at all) we should
-    ! correct for heat_liquid and cv_liquid, which are used to determine the weighted
-    ! average liquid water temperature. For example, if we're subtracting out a baseline
-    ! water amount because a particular water state is fictitious, we probably shouldn't
-    ! include that particular state when determining the weighted average temperature of
-    ! liquid water. And conversely, if we're adding a state via these baselines, should
-    ! we also add some water temperature of that state? The tricky thing here is what to
-    ! do when we end up subtracting water due to the baselines, so for now I'm simply not
-    ! trying to account for the temperature of these baselines at all. This amounts to
-    ! assuming that the baselines that we add / subtract are at the average temperature
-    ! of the real liquid water in the column.
-    do fc = 1, num_lakec
-       c = filter_lakec(fc)
-       heat(c) = heat(c) - dynbal_baseline_heat(c)
-    end do
-
     end associate
 
   end subroutine ComputeHeatLake
 
   !-----------------------------------------------------------------------
   subroutine AccumulateHeatLake(bounds, num_c, filter_c, &
-        temperature_inst, lakestate_inst, &
+       temperature_inst, lakestate_inst, &
        heat)
     !
     ! !DESCRIPTION:
@@ -1097,11 +1131,11 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, &
         end do
     end do 
 
-   ! add ice heat and liquid heat together
+    ! add ice heat and liquid heat together
     do fc = 1, num_c
        c = filter_c(fc)
-       heat(c) = heat(c) + lake_heat_ice(c) + &
-            lake_heat_liquid(c) + lake_latent_heat_liquid(c)
+       heat(c) = heat(c) + (lake_heat_ice(c) + &
+            lake_heat_liquid(c) + lake_latent_heat_liquid(c))
     end do
 
     end associate
diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90
index d8bb809bba..089715233d 100644
--- a/src/dyn_subgrid/dynConsBiogeophysMod.F90
+++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90
@@ -66,7 +66,8 @@ module dynConsBiogeophysMod
   !-----------------------------------------------------------------------
   subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
        num_lakec, filter_lakec, &
-       urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst)
+       urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, &
+       reset_all_baselines, reset_lake_baselines)
     !
     ! !DESCRIPTION:
     ! Set start-of-run baseline values for heat and water content in some columns.
@@ -79,13 +80,18 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
     ! represented). These corrections will typically reduce the fictitious dynbal
     ! conservation fluxes.
     !
-    ! At a minimum, this should be called during cold start initialization, to initialize
-    ! the baseline values based on cold start states. In addition, this can be called when
-    ! starting a new run from existing initial conditions, if the user desires. This
-    ! optional resetting can further reduce the dynbal fluxes; however, it can break
-    ! conservation. (So, for example, it can be done when transitioning from an offline
-    ! spinup to a coupled run, but it should not be done when transitioning from a
-    ! coupled historical run to a future scenario.)
+    ! At a minimum, this should be called during cold start initialization with
+    ! reset_all_baselines set to true, to initialize the baseline values based on cold
+    ! start states. This should also be called after reading initial conditions, but the
+    ! various reset_* flags should be set appropriately; setting all of these flags to
+    ! false will result in no baselines being reset in this call.
+
+    ! Setting reset_all_baselines can be done when starting a new run from existing
+    ! initial conditions if the user desires. This optional resetting can further reduce
+    ! the dynbal fluxes; however, it can break conservation. (So, for example, it can be
+    ! done when transitioning from an offline spinup to a coupled run, but it should not
+    ! be done when transitioning from a coupled historical run to a future scenario.)
+    ! Other reset_* flags are described below.
     !
     ! !ARGUMENTS:
     type(bounds_type), intent(in) :: bounds
@@ -98,12 +104,25 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
     integer, intent(in) :: num_lakec  ! number of points in filter_lakec
     integer, intent(in) :: filter_lakec(:) ! filter for lake columns
 
-    
     type(urbanparams_type), intent(in) :: urbanparams_inst
     type(soilstate_type), intent(in) :: soilstate_inst
     type(lakestate_type), intent(in) :: lakestate_inst
     type(water_type), intent(inout) :: water_inst
     type(temperature_type), intent(inout) :: temperature_inst
+
+    ! Whether to reset baselines for all columns
+    logical, intent(in) :: reset_all_baselines
+
+    ! Whether to reset baselines for lake columns. Ignored if reset_all_baselines is
+    ! true.
+    !
+    ! BACKWARDS_COMPATIBILITY(wjs, 2020-09-02) This is needed when reading old initial
+    ! conditions files created before https://github.com/ESCOMP/CTSM/issues/1140 was
+    ! resolved via https://github.com/ESCOMP/CTSM/pull/1109: The definition of total
+    ! column water content has been changed for lakes, so we need to reset baseline values
+    ! for lakes on older initial conditions files.
+    logical, intent(in) :: reset_lake_baselines
+
     !
     ! !LOCAL VARIABLES:
     integer :: i
@@ -125,22 +144,26 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
 
        call dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
             num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
-            bulk_or_tracer%waterstate_inst, lakestate_inst)
+            bulk_or_tracer%waterstate_inst, lakestate_inst, &
+            reset_all_baselines = reset_all_baselines, &
+            reset_lake_baselines = reset_lake_baselines)
        end associate
     end do
        
     call dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
          num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
          urbanparams_inst, soilstate_inst, lakestate_inst, water_inst%waterstatebulk_inst, &
-         temperature_inst)
-         
+         temperature_inst, &
+         reset_all_baselines = reset_all_baselines, &
+         reset_lake_baselines = reset_lake_baselines)
 
   end subroutine dyn_hwcontent_set_baselines
 
   !-----------------------------------------------------------------------
   subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
        num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
-       waterstate_inst, lakestate_inst)
+       waterstate_inst, lakestate_inst, &
+       reset_all_baselines, reset_lake_baselines)
     !
     ! !DESCRIPTION:
     ! Set start-of-run baseline values for water content, for a single water tracer or
@@ -155,9 +178,13 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
     integer, intent(in) :: filter_icemecc(:) ! filter for icemec (i.e., glacier) columns
     integer, intent(in) :: num_lakec  ! number of points in filter_lakec
     integer, intent(in) :: filter_lakec(:) ! filter for lake columns
-    type(lakestate_type), intent(in) :: lakestate_inst
 
+    type(lakestate_type), intent(in) :: lakestate_inst
     class(waterstate_type), intent(inout) :: waterstate_inst
+
+    ! See dyn_hwcontent_set_baselines for documentation of these arguments
+    logical, intent(in) :: reset_all_baselines
+    logical, intent(in) :: reset_lake_baselines
     !
     ! !LOCAL VARIABLES:
     integer  :: c, fc  ! indices
@@ -174,46 +201,50 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
          dynbal_baseline_ice => waterstate_inst%dynbal_baseline_ice_col  & ! Output: [real(r8) (:)   ]  baseline ice content subtracted from each column's total ice calculation (mm H2O)
          )
 
-    soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
-    soil_ice_mass_col   (bounds%begc:bounds%endc) = 0._r8
-    lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
-    lake_ice_mass_col   (bounds%begc:bounds%endc) = 0._r8
-    
-    
-    call AccumulateSoilLiqIceMassNonLake(bounds, &
-         natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
-         waterstate_inst, &
-         liquid_mass = soil_liquid_mass_col(bounds%begc:bounds%endc), &
-         ice_mass = soil_ice_mass_col(bounds%begc:bounds%endc))
-
-    ! For glacier columns, the liquid and ice in the "soil" (i.e., in the glacial ice) is
-    ! a virtual pool. So we'll subtract this amount when determining the dynbal
-    ! fluxes. (Note that a positive value in these baseline variables indicates an extra
-    ! stock that will be subtracted later.) But glacier columns do not represent the soil
-    ! under the glacial ice. Let's assume that the soil state under each glacier column is
-    ! the same as the soil state in the natural vegetation landunit on that grid cell. We
-    ! subtract this from the dynbal baseline variables to indicate a missing stock.
-    call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
-         vals_col = soil_liquid_mass_col(bounds%begc:bounds%endc), &
-         baselines_col = dynbal_baseline_liq(bounds%begc:bounds%endc))
-    call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
-         vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), &
-         baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc))
-         
-         
-    ! set baselines for lake columns
-    
-    ! Calculate the total water volume of the lake column
-    call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, &
-        liquid_mass = lake_liquid_mass_col(bounds%begc:bounds%endc), &
-        ice_mass = lake_ice_mass_col(bounds%begc:bounds%endc))
-     
-    do fc = 1, num_lakec
-       c = filter_lakec(fc)
-       dynbal_baseline_liq(c) = lake_liquid_mass_col(c)       
-       dynbal_baseline_ice(c) = lake_ice_mass_col(c)
-    end do
+    if (reset_all_baselines) then
+       soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
+       soil_ice_mass_col   (bounds%begc:bounds%endc) = 0._r8
+
+       call AccumulateSoilLiqIceMassNonLake(bounds, &
+            natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
+            waterstate_inst, &
+            liquid_mass = soil_liquid_mass_col(bounds%begc:bounds%endc), &
+            ice_mass = soil_ice_mass_col(bounds%begc:bounds%endc))
+
+       ! For glacier columns, the liquid and ice in the "soil" (i.e., in the glacial ice) is
+       ! a virtual pool. So we'll subtract this amount when determining the dynbal
+       ! fluxes. (Note that a positive value in these baseline variables indicates an extra
+       ! stock that will be subtracted later.) But glacier columns do not represent the soil
+       ! under the glacial ice. Let's assume that the soil state under each glacier column is
+       ! the same as the soil state in the natural vegetation landunit on that grid cell. We
+       ! subtract this from the dynbal baseline variables to indicate a missing stock.
+       call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
+            vals_col = soil_liquid_mass_col(bounds%begc:bounds%endc), &
+            baselines_col = dynbal_baseline_liq(bounds%begc:bounds%endc))
+       call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
+            vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), &
+            baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc))
+    end if
+
+    if (reset_all_baselines .or. reset_lake_baselines) then
+       ! set baselines for lake columns
        
+       lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
+       lake_ice_mass_col   (bounds%begc:bounds%endc) = 0._r8
+
+       ! Calculate the total water volume of the lake column
+       call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, &
+            tracer_ratio = waterstate_inst%info%get_ratio(), &
+            liquid_mass = lake_liquid_mass_col(bounds%begc:bounds%endc), &
+            ice_mass = lake_ice_mass_col(bounds%begc:bounds%endc))
+
+       do fc = 1, num_lakec
+          c = filter_lakec(fc)
+          dynbal_baseline_liq(c) = lake_liquid_mass_col(c)       
+          dynbal_baseline_ice(c) = lake_ice_mass_col(c)
+       end do
+    end if
+
     end associate
 
   end subroutine dyn_water_content_set_baselines
@@ -222,7 +253,8 @@ end subroutine dyn_water_content_set_baselines
   subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
        num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
        urbanparams_inst, soilstate_inst, lakestate_inst, waterstatebulk_inst, &
-       temperature_inst)
+       temperature_inst, &
+       reset_all_baselines, reset_lake_baselines)
     !
     ! !DESCRIPTION:
     ! Set start-of-run baseline values for heat content.
@@ -242,6 +274,10 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
     type(lakestate_type)  , intent(in)    :: lakestate_inst
     type(waterstatebulk_type), intent(in) :: waterstatebulk_inst
     type(temperature_type), intent(inout) :: temperature_inst
+
+    ! See dyn_hwcontent_set_baselines for documentation of these arguments
+    logical, intent(in) :: reset_all_baselines
+    logical, intent(in) :: reset_lake_baselines
     !
     ! !LOCAL VARIABLES:
     integer  :: c, fc  ! indices
@@ -257,51 +293,52 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
          dynbal_baseline_heat => temperature_inst%dynbal_baseline_heat_col & ! Output: [real(r8) (:)   ]  baseline heat content subtracted from each column's total heat calculation (J/m2)
          )
 
-    soil_heat_col(bounds%begc:bounds%endc) = 0._r8
-    soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8
-    soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8
-
-    lake_heat_col(bounds%begc:bounds%endc) = 0._r8
-      
-    call AccumulateSoilHeatNonLake(bounds, &
-         natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
-         urbanparams_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, &
-         heat = soil_heat_col(bounds%begc:bounds%endc), &
-         heat_liquid = soil_heat_liquid_col(bounds%begc:bounds%endc), &
-         cv_liquid = soil_cv_liquid_col(bounds%begc:bounds%endc))
-        
-
-    ! See comments in dyn_water_content_set_baselines for rationale for these glacier
-    ! baselines. Even though the heat in glacier ice can interact with the rest of the
-    ! system (e.g., giving off heat to the atmosphere or receiving heat from the
-    ! atmosphere), it is still a virtual state in the sense that the glacier ice column
-    ! is virtual. And, as for water, we subtract the heat of the soil in the associated
-    ! natural vegetated landunit to account for the fact that we don't explicitly model
-    ! the soil under glacial ice.
-    !
-    ! Aside from these considerations of what is virtual vs. real, another rationale
-    ! driving the use of these baselines is the desire to minimize the dynbal fluxes. For
-    ! the sake of conservation, it seems okay to pick any fixed baseline when summing the
-    ! energy (or water) content of a given column (as long as that baseline doesn't
-    ! change over time). By using the baselines computed here, we reduce the dynbal
-    ! fluxes to more reasonable values.
-    call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
-         vals_col = soil_heat_col(bounds%begc:bounds%endc), &
-         baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc))
-         
-    
-    ! Set baselines for lake columns
-    call AccumulateHeatLake(bounds, num_lakec, filter_lakec, &
-        temperature_inst, lakestate_inst, &
-        heat = lake_heat_col)
-    
-    do fc = 1, num_lakec
-       c = filter_lakec(fc)
+    if (reset_all_baselines) then
+       soil_heat_col(bounds%begc:bounds%endc) = 0._r8
+       soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8
+       soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8
+
+       call AccumulateSoilHeatNonLake(bounds, &
+            natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
+            urbanparams_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, &
+            heat = soil_heat_col(bounds%begc:bounds%endc), &
+            heat_liquid = soil_heat_liquid_col(bounds%begc:bounds%endc), &
+            cv_liquid = soil_cv_liquid_col(bounds%begc:bounds%endc))
+
+       ! See comments in dyn_water_content_set_baselines for rationale for these glacier
+       ! baselines. Even though the heat in glacier ice can interact with the rest of the
+       ! system (e.g., giving off heat to the atmosphere or receiving heat from the
+       ! atmosphere), it is still a virtual state in the sense that the glacier ice column
+       ! is virtual. And, as for water, we subtract the heat of the soil in the associated
+       ! natural vegetated landunit to account for the fact that we don't explicitly model
+       ! the soil under glacial ice.
+       !
+       ! Aside from these considerations of what is virtual vs. real, another rationale
+       ! driving the use of these baselines is the desire to minimize the dynbal fluxes. For
+       ! the sake of conservation, it seems okay to pick any fixed baseline when summing the
+       ! energy (or water) content of a given column (as long as that baseline doesn't
+       ! change over time). By using the baselines computed here, we reduce the dynbal
+       ! fluxes to more reasonable values.
+       call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
+            vals_col = soil_heat_col(bounds%begc:bounds%endc), &
+            baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc))
+    end if
+
+    if (reset_all_baselines .or. reset_lake_baselines) then
+       ! Set baselines for lake columns
+
+       lake_heat_col(bounds%begc:bounds%endc) = 0._r8
+
+       call AccumulateHeatLake(bounds, num_lakec, filter_lakec, &
+            temperature_inst, lakestate_inst, &
+            heat = lake_heat_col)
+
+       do fc = 1, num_lakec
+          c = filter_lakec(fc)
+          dynbal_baseline_heat(c) = lake_heat_col(c)
+       end do
+    end if
 
-       dynbal_baseline_heat(c) = lake_heat_col(c)
-       
-    end do
-    
     end associate
 
   end subroutine dyn_heat_content_set_baselines
@@ -645,7 +682,7 @@ subroutine dyn_water_content(bounds, &
     call ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
          waterstate_inst, &
          lakestate_inst, &
-         subtract_dynbal_baselines = .true., &
+         add_lake_water_and_subtract_dynbal_baselines = .true., &
          liquid_mass = liquid_mass_col(bounds%begc:bounds%endc), &
          ice_mass = ice_mass_col(bounds%begc:bounds%endc))
 
diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90
index c07a47c7ed..912bbb293b 100644
--- a/src/main/clm_initializeMod.F90
+++ b/src/main/clm_initializeMod.F90
@@ -339,6 +339,8 @@ subroutine initialize2( )
     logical               :: lexist
     integer               :: closelatidx,closelonidx
     real(r8)              :: closelat,closelon
+    logical               :: reset_dynbal_baselines_all_columns
+    logical               :: reset_dynbal_baselines_lake_columns
     integer               :: begp, endp
     integer               :: begc, endc
     integer               :: begl, endl
@@ -490,8 +492,12 @@ subroutine initialize2( )
             filter_inactive_and_active(nc)%icemecc, &
             filter_inactive_and_active(nc)%num_lakec, &
             filter_inactive_and_active(nc)%lakec, &
-            urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst)
+            urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, &
+            reset_all_baselines = .true., &
+            ! reset_lake_baselines is irrelevant since reset_all_baselines is true
+            reset_lake_baselines = .false.)
     end do
+    !$OMP END PARALLEL DO
 
     ! ------------------------------------------------------------------------
     ! Initialize modules (after time-manager initialization in most cases)
@@ -543,6 +549,7 @@ subroutine initialize2( )
 
     is_cold_start = .false.
     is_interpolated_start = .false.
+    reset_dynbal_baselines_lake_columns = .false.
 
     if (nsrest == nsrStartup) then
 
@@ -563,7 +570,8 @@ subroutine initialize2( )
              write(iulog,*)'Reading initial conditions from ',trim(finidat)
           end if
           call getfil( finidat, fnamer, 0 )
-          call restFile_read(bounds_proc, fnamer, glc_behavior)
+          call restFile_read(bounds_proc, fnamer, glc_behavior, &
+               reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)
        end if
 
     else if ((nsrest == nsrContinue) .or. (nsrest == nsrBranch)) then
@@ -571,8 +579,8 @@ subroutine initialize2( )
        if (masterproc) then
           write(iulog,*)'Reading restart file ',trim(fnamer)
        end if
-       call restFile_read(bounds_proc, fnamer, glc_behavior)
-
+       call restFile_read(bounds_proc, fnamer, glc_behavior, &
+            reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)
     end if
 
     ! ------------------------------------------------------------------------
@@ -598,7 +606,8 @@ subroutine initialize2( )
             glc_behavior=glc_behavior)
 
        ! Read new interpolated conditions file back in
-       call restFile_read(bounds_proc, finidat_interp_dest, glc_behavior)
+       call restFile_read(bounds_proc, finidat_interp_dest, glc_behavior, &
+            reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)
 
        ! Reset finidat to now be finidat_interp_dest
        ! (to be compatible with routines still using finidat)
@@ -613,33 +622,49 @@ subroutine initialize2( )
     ! interpolated restart file, if applicable).
     ! ------------------------------------------------------------------------
 
-    if (get_reset_dynbal_baselines()) then
-       if (nsrest == nsrStartup) then
-          if (masterproc) then
-             write(iulog,*) ' '
-             write(iulog,*) 'Resetting dynbal baselines'
-             write(iulog,*) ' '
-          end if
-
-          !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
-          do nc = 1,nclumps
-             call get_clump_bounds(nc, bounds_clump)
-
-             call dyn_hwcontent_set_baselines(bounds_clump, &
-                  filter_inactive_and_active(nc)%num_icemecc, &
-                  filter_inactive_and_active(nc)%icemecc, &
-                  filter_inactive_and_active(nc)%num_lakec, &
-                  filter_inactive_and_active(nc)%lakec, &
-                  urbanparams_inst, soilstate_inst, lakestate_inst, &
-                  water_inst, temperature_inst)
-          end do
-       else if (nsrest == nsrBranch) then
+    reset_dynbal_baselines_all_columns = get_reset_dynbal_baselines()
+    if (nsrest == nsrBranch) then
+       if (reset_dynbal_baselines_all_columns) then
           call endrun(msg='ERROR clm_initializeMod: '//&
                'Cannot set reset_dynbal_baselines in a branch run')
        end if
-       ! nsrContinue not explicitly handled: it's okay for reset_dynbal_baselines to
-       ! remain set in a continue run, but it has no effect
+    else if (nsrest == nsrContinue) then
+       ! It's okay for the reset_dynbal_baselines flag to remain set in a continue
+       ! run, but we'll ignore it. (This way, the user doesn't have to change their
+       ! namelist file for the continue run.)
+       reset_dynbal_baselines_all_columns = .false.
     end if
+    ! Note that we will still honor reset_dynbal_baselines_lake_columns even in a branch
+    ! or continue run: even in these runs, we want to reset those baselines if they are
+    ! wrong on the restart file.
+
+    if (masterproc) then
+       if (reset_dynbal_baselines_all_columns) then
+          write(iulog,*) ' '
+          write(iulog,*) 'Resetting dynbal baselines for all columns'
+          write(iulog,*) ' '
+       else if (reset_dynbal_baselines_lake_columns) then
+          write(iulog,*) ' '
+          write(iulog,*) 'Resetting dynbal baselines for lake columns'
+          write(iulog,*) ' '
+       end if
+    end if
+
+    !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
+    do nc = 1,nclumps
+       call get_clump_bounds(nc, bounds_clump)
+
+       call dyn_hwcontent_set_baselines(bounds_clump, &
+            filter_inactive_and_active(nc)%num_icemecc, &
+            filter_inactive_and_active(nc)%icemecc, &
+            filter_inactive_and_active(nc)%num_lakec, &
+            filter_inactive_and_active(nc)%lakec, &
+            urbanparams_inst, soilstate_inst, lakestate_inst, &
+            water_inst, temperature_inst, &
+            reset_all_baselines = reset_dynbal_baselines_all_columns, &
+            reset_lake_baselines = reset_dynbal_baselines_lake_columns)
+    end do
+    !$OMP END PARALLEL DO
 
     ! ------------------------------------------------------------------------
     ! Initialize nitrogen deposition
diff --git a/src/main/restFileMod.F90 b/src/main/restFileMod.F90
index 83be13835b..4302cbc1c0 100644
--- a/src/main/restFileMod.F90
+++ b/src/main/restFileMod.F90
@@ -25,6 +25,7 @@ module restFileMod
   use ncdio_pio        , only : check_att, ncd_getatt
   use glcBehaviorMod   , only : glc_behavior_type
   use reweightMod      , only : reweight_wrapup
+  use IssueFixedMetadataHandler, only : write_issue_fixed_metadata, read_issue_fixed_metadata
   !
   ! !PUBLIC TYPES:
   implicit none
@@ -44,6 +45,8 @@ module restFileMod
   private :: restFile_write_pfile       ! Writes restart pointer file
   private :: restFile_closeRestart      ! Close restart file and write restart pointer file
   private :: restFile_dimset
+  private :: restFile_write_issues_fixed ! Write metadata for issues fixed
+  private :: restFile_read_issues_fixed ! Read and process metadata for issues fixed
   private :: restFile_add_flag_metadata ! Add global metadata for some logical flag
   private :: restFile_add_ilun_metadata ! Add global metadata defining landunit types
   private :: restFile_add_icol_metadata ! Add global metadata defining column types
@@ -56,6 +59,9 @@ module restFileMod
   !
   ! !PRIVATE TYPES:
 
+  ! Issue numbers for issue-fixed metadata
+  integer, parameter :: lake_dynbal_baseline_issue = 1140
+
   character(len=*), parameter, private :: sourcefile = &
        __FILE__
   !-----------------------------------------------------------------------
@@ -108,6 +114,9 @@ subroutine restFile_write( bounds, file, writing_finidat_interp_dest_file, rdate
        call hist_restart_ncd (bounds, ncid, flag='define', rdate=rdate )
     end if
 
+    call restFile_write_issues_fixed(ncid, &
+         writing_finidat_interp_dest_file = writing_finidat_interp_dest_file)
+
     call restFile_enddef( ncid )
 
     ! Write variables
@@ -142,7 +151,7 @@ subroutine restFile_write( bounds, file, writing_finidat_interp_dest_file, rdate
   end subroutine restFile_write
 
   !-----------------------------------------------------------------------
-  subroutine restFile_read( bounds_proc, file, glc_behavior )
+  subroutine restFile_read( bounds_proc, file, glc_behavior, reset_dynbal_baselines_lake_columns )
     !
     ! !DESCRIPTION:
     ! Read a CLM restart file.
@@ -151,6 +160,13 @@ subroutine restFile_read( bounds_proc, file, glc_behavior )
     type(bounds_type) , intent(in) :: bounds_proc      ! processor-level bounds
     character(len=*)  , intent(in) :: file             ! output netcdf restart file
     type(glc_behavior_type), intent(in) :: glc_behavior
+
+    ! BACKWARDS_COMPATIBILITY(wjs, 2020-09-02) This is needed when reading old initial
+    ! conditions files created before https://github.com/ESCOMP/CTSM/issues/1140 was
+    ! resolved via https://github.com/ESCOMP/CTSM/pull/1109: The definition of total
+    ! column water content has been changed for lakes, so we need to reset baseline values
+    ! for lakes on older initial conditions files.
+    logical, intent(out) :: reset_dynbal_baselines_lake_columns
     !
     ! !LOCAL VARIABLES:
     type(file_desc_t) :: ncid    ! netcdf id
@@ -201,6 +217,9 @@ subroutine restFile_read( bounds_proc, file, glc_behavior )
 
     call hist_restart_ncd (bounds_proc, ncid, flag='read' )
 
+    call restFile_read_issues_fixed(ncid, &
+         reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)
+
     ! Do error checking on file
     
     call restFile_check_consistency(bounds_proc, ncid)
@@ -563,6 +582,63 @@ subroutine restFile_dimset( ncid )
 
   end subroutine restFile_dimset
 
+  !-----------------------------------------------------------------------
+  subroutine restFile_write_issues_fixed(ncid, writing_finidat_interp_dest_file)
+    !
+    ! !DESCRIPTION:
+    ! Write metadata for issues fixed
+    !
+    ! !ARGUMENTS:
+    type(file_desc_t), intent(inout) :: ncid ! local file id
+    logical          , intent(in)    :: writing_finidat_interp_dest_file ! true if we are writing a finidat_interp_dest file
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'restFile_write_issues_fixed'
+    !-----------------------------------------------------------------------
+
+    ! See comment associated with reset_dynbal_baselines_lake_columns in restFile_read
+    call write_issue_fixed_metadata( &
+         ncid = ncid, &
+         writing_finidat_interp_dest_file = writing_finidat_interp_dest_file, &
+         issue_num = lake_dynbal_baseline_issue)
+
+  end subroutine restFile_write_issues_fixed
+
+  !-----------------------------------------------------------------------
+  subroutine restFile_read_issues_fixed(ncid, reset_dynbal_baselines_lake_columns)
+    !
+    ! !DESCRIPTION:
+    ! Read and process metadata for issues fixed
+    !
+    ! !ARGUMENTS:
+    type(file_desc_t), intent(inout) :: ncid ! local file id
+    logical, intent(out) :: reset_dynbal_baselines_lake_columns ! see comment in restFile_read for details
+    !
+    ! !LOCAL VARIABLES:
+    integer :: attribute_value
+
+    character(len=*), parameter :: subname = 'restFile_read_issues_fixed'
+    !-----------------------------------------------------------------------
+
+    ! See comment associated with reset_dynbal_baselines_lake_columns in restFile_read
+    call read_issue_fixed_metadata( &
+         ncid = ncid, &
+         issue_num = lake_dynbal_baseline_issue, &
+         attribute_value = attribute_value)
+    if (attribute_value == 0) then
+       ! Old restart file, from before lake_dynbal_baseline_issue was resolved, so we
+       ! need to reset dynbal baselines for lake columns later in initialization.
+       reset_dynbal_baselines_lake_columns = .true.
+    else
+       ! Recent restart file, so no need to reset dynbal baselines for lake columns,
+       ! because they are already up to date.
+       reset_dynbal_baselines_lake_columns = .false.
+    end if
+
+  end subroutine restFile_read_issues_fixed
+
+
   !-----------------------------------------------------------------------
   subroutine restFile_add_flag_metadata(ncid, flag, flag_name)
     !