Merge pull request #1 from billsacks/dynlakes

Merge dynbal baseline changes into dynlake branch

bld/CLMBuildNamelist.pm

@@ -2503,6 +2503,11 @@ sub setup_logic_dynamic_subgrid {
    setup_logic_do_transient_crops($opts, $nl_flags, $definition, $defaults, $nl, $physv);
    setup_logic_do_harvest($opts, $nl_flags, $definition, $defaults, $nl, $physv);
 
+   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'reset_dynbal_baselines');
+   if ( &value_is_true($nl->get_value('reset_dynbal_baselines')) &&
+        &remove_leading_and_trailing_quotes($nl_flags->{'clm_start_type'}) eq "branch") {
+      $log->fatal_error("reset_dynbal_baselines has no effect in a branch run");
+   }
 }
 
 sub setup_logic_do_transient_pfts {

bld/namelist_files/namelist_defaults_clm4_5.xml

@@ -2619,6 +2619,13 @@ lnd/clm2/surfdata_map/landuse.timeseries_ne30np4_hist_16pfts_Irrig_CMIP6_simyr18
 <use_fates_inventory_init      use_fates=".true.">.false.</use_fates_inventory_init>
 <fates_parteh_mode             use_fates=".true.">1</fates_parteh_mode>
 
+<!-- =========================================  -->
+<!-- Defaults for dynamic subgrid               -->
+<!-- (some other defaults for dynamic subgrid   -->
+<!-- are set in the BuildNamelist code)         -->
+<!-- =========================================  -->
+<reset_dynbal_baselines>.false.</reset_dynbal_baselines>
+
 <!-- =========================================  -->
 <!-- Defaults for init_interp                   -->
 <!-- =========================================  -->

bld/namelist_files/namelist_definition_clm4_5.xml

@@ -2281,6 +2281,52 @@ If TRUE, apply harvest from flanduse_timeseries file.
 (Also, only valid for use_cn = true.)
 </entry>
 
+<entry id="reset_dynbal_baselines" type="logical" category="physics"
+       group="dynamic_subgrid" valid_values="">
+If TRUE, reset baseline values of total column water and energy in the
+first step of the run. This should typically be set when transitioning
+from an offline spinup to a coupled run with transient glaciers, and
+*possibly* when transitioning from the spinup to the transient portion
+of a coupled run with transient glaciers; it should typically remain
+unset at other times.
+
+These baseline values are computed only for particular columns
+(currently, only for glacier columns). They provide values that are
+subtracted from the current state when counting total column water and
+energy. For glacier columns, these discount the water and energy
+contents in the "virtual" glacier ice, while adding representative
+amounts of water and energy in the non-explicitly-modeled soil beneath
+the ice. Subtracting these baselines reduces the fictitious dynbal
+fluxes generated when total grid cell water and energy changes as a
+result of dynamic column/landunit areas.
+
+These baseline values are initially computed based on cold start
+states. If this flag remains unset (.false.), these baseline values will
+remain fixed at their cold start values. This will conserve mass and
+energy, but may result in larger-than-desired dynbal energy fluxes (and,
+in principle, also larger dynbal water fluxes; but currently, the mass
+of glacier ice remains fixed over time, so dynbal water fluxes are
+fairly small regardless of whether this flag is ever set). To further
+reduce these dynbal fluxes, you can set this flag to .true. when
+starting a startup or hybrid run from a partially or entirely spun-up
+state. This will reset the baseline values based on the state at the
+start of this run.
+
+Note that setting this flag can break water and energy conservation!
+Specifically, any water and energy that has previously been added to or
+removed from states that contribute to these baselines (currently, (a)
+glacier ice and (b) soil water and energy in the vegetated landunit in
+the same grid cell as glaciers) will effectively be ignored when
+computing conservation corrections due to land cover change. Instead,
+only the change in states from this point forward will be
+considered. So, for example, this flag should NOT be set when
+transitioning from a historical run to a future scenario.
+
+This setting only impacts startup and hybrid runs; it has no effect in a
+continue run; it is an error for this to be set in a branch
+run. Furthermore, this setting has no effect in a cold start run.
+</entry>
+
 <entry id="for_testing_allow_non_annual_changes" type="logical" category="physics"
        group="dynamic_subgrid" valid_values="" >
 If TRUE, allow area changes at times other than the year boundary. This should

cime_config/testdefs/testlist_clm.xml

@@ -785,13 +785,13 @@
       <option name="wallclock">00:20:00</option>
     </options>
   </test>
-  <test name="ERS_D_Ld10" grid="f10_f10_musgs" compset="IHistClm50SpG" testmods="clm/glcMEC_decrease">
+  <test name="ERP_P36x2_D_Ld10" grid="f10_f10_musgs" compset="IHistClm50SpG" testmods="clm/glcMEC_decrease">
     <machines>
       <machine name="cheyenne" compiler="intel" category="aux_clm"/>
     </machines>
     <options>
       <option name="wallclock">00:20:00</option>
-      <option name="comment"  >test transient PFTs (via HIST) in conjunction with changing glacier area</option>
+      <option name="comment"  >Test transient PFTs (via HIST) in conjunction with changing glacier area. This test also covers the reset_dynbal_baselines option. CISM is not answer preserving across processor changes, but short test length should be OK.</option>
     </options>
   </test>
   <test name="ERS_D_Ld12" grid="f10_f10_musgs" compset="I1850Clm50BgcCropG" testmods="clm/glcMEC_spunup_inc_dec_bgc">

cime_config/testdefs/testmods_dirs/clm/glcMEC_decrease/user_nl_clm

@@ -1 +1,7 @@
 for_testing_allow_non_annual_changes = .true.
+
+! Also test the reset_dynbal_baselines flag. We want at least one restart test that uses
+! this flag; for this test to have power, it should have some change in glacier area
+! after the restart. (So, in this decrease test, the glacier decrease should be slow
+! enough that there is some decrease after the restart.)
+reset_dynbal_baselines = .true.

cime_config/testdefs/testmods_dirs/clm/glcMEC_increase/README

@@ -1,2 +1,2 @@
 This testmods directory is like the standard glcMEC testmods directory, except
-it forces a quick increas in glacier area.
+it forces a quick increase in glacier area.

src/biogeophys/BalanceCheckMod.F90

@@ -150,10 +150,14 @@ subroutine BeginWaterBalance(bounds, &
     end do
 
     call ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, &
-         soilhydrology_inst, waterstate_inst, begwb(bounds%begc:bounds%endc))
+         soilhydrology_inst, waterstate_inst, &
+         subtract_dynbal_baselines = .false., &
+         water_mass = begwb(bounds%begc:bounds%endc))
 
     call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
-         waterstate_inst, begwb(bounds%begc:bounds%endc))
+         waterstate_inst, &
+         subtract_dynbal_baselines = .false., &
+         water_mass = begwb(bounds%begc:bounds%endc))
 
     end associate 
 

src/biogeophys/HydrologyDrainageMod.F90

@@ -158,9 +158,9 @@ subroutine HydrologyDrainage(bounds,               &
       end do
 
       call ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, &
-           soilhydrology_inst, waterstate_inst, endwb(bounds%begc:bounds%endc))
-
-
+           soilhydrology_inst, waterstate_inst, &
+           subtract_dynbal_baselines = .false., &
+           water_mass = endwb(bounds%begc:bounds%endc))
 
 
       ! Determine wetland and land ice hydrology (must be placed here

src/biogeophys/LakeHydrologyMod.F90

@@ -624,7 +624,9 @@ subroutine LakeHydrology(bounds, &
       ! Determine ending water balance and volumetric soil water
 
       call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
-           waterstate_inst, endwb(bounds%begc:bounds%endc))
+           waterstate_inst, &
+           subtract_dynbal_baselines = .false., &
+           water_mass = endwb(bounds%begc:bounds%endc))
 
       do j = 1, nlevgrnd
          do fc = 1, num_lakec

src/biogeophys/TemperatureType.F90

@@ -95,6 +95,10 @@ module TemperatureType
 
      ! Heat content
      real(r8), pointer :: beta_col                 (:)   ! coefficient of convective velocity [-]
+     ! For the following dynbal baseline variable: positive values are subtracted to avoid
+     ! counting liquid water content of "virtual" states; negative values are added to
+     ! account for missing states in the model.
+     real(r8), pointer :: dynbal_baseline_heat_col (:)   ! baseline heat content subtracted from each column's total heat calculation [J/m^2]
      real(r8), pointer :: heat1_grc                (:)   ! grc initial gridcell total heat content
      real(r8), pointer :: heat2_grc                (:)   ! grc post land cover change total heat content
      real(r8), pointer :: liquid_water_temp1_grc   (:)   ! grc initial weighted average liquid water temperature (K)
@@ -260,6 +264,7 @@ subroutine InitAllocate(this, bounds)
 
     ! Heat content
     allocate(this%beta_col                 (begc:endc))                      ; this%beta_col                 (:)   = nan
+    allocate(this%dynbal_baseline_heat_col (begc:endc))                      ; this%dynbal_baseline_heat_col (:)   = nan
     allocate(this%heat1_grc                (begg:endg))                      ; this%heat1_grc                (:)   = nan
     allocate(this%heat2_grc                (begg:endg))                      ; this%heat2_grc                (:)   = nan
     allocate(this%liquid_water_temp1_grc   (begg:endg))                      ; this%liquid_water_temp1_grc   (:)   = nan
@@ -849,6 +854,11 @@ subroutine InitCold(this, bounds, &
        if (col%itype(c) == icol_road_perv  ) this%emg_col(c) = em_perroad_lun(l)
     end do
 
+    ! Initialize dynbal_baseline_heat_col: for some columns, this is set elsewhere in
+    ! initialization, but we need it to be 0 for columns for which it is not explicitly
+    ! set.
+    this%dynbal_baseline_heat_col(bounds%begc:bounds%endc) = 0._r8
+
   end subroutine InitCold
 
   !------------------------------------------------------------------------
@@ -991,6 +1001,12 @@ subroutine Restart(this, bounds, ncid, flag, is_simple_buildtemp, is_prog_buildt
          long_name='urban canopy air temperature', units='K',                                                         &
          interpinic_flag='interp', readvar=readvar, data=this%taf_lun)
 
+    call restartvar(ncid=ncid, flag=flag, varname='DYNBAL_BASELINE_HEAT', xtype=ncd_double, &
+         dim1name='column', &
+         long_name="baseline heat content subtracted from each column's total heat calculation", &
+         units='J/m2', &
+         interpinic_flag='interp', readvar=readvar, data=this%dynbal_baseline_heat_col)
+
     if (use_crop) then
        call restartvar(ncid=ncid, flag=flag,  varname='gdd1020', xtype=ncd_double,  &
             dim1name='pft', long_name='20 year average of growing degree-days base 10C from planting', units='ddays', &