Dynamical lakes and new lake cover map

bld/CLMBuildNamelist.pm

@@ -2524,6 +2524,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

@@ -2647,6 +2647,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

@@ -2337,13 +2337,65 @@ If TRUE, apply transient crops from flanduse_timeseries file.
 (Only valid for transient runs, where there is a flanduse_timeseries file.)
 </entry>
 
+<entry id="do_transient_lakes" type="logical" category="physics"
+       group="dynamic_subgrid" valid_values="" >
+If TRUE, apply transient lakes from flanduse_timeseries file.
+(Only valid for transient runs, where there is a flanduse_timeseries file.)
+</entry>
+
 <entry id="do_harvest" type="logical" category="physics"
        group="dynamic_subgrid" valid_values="" >
 If TRUE, apply harvest from flanduse_timeseries file.
 (Only valid for transient runs, where there is a 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

@@ -769,13 +769,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

@@ -19,6 +19,7 @@ module BalanceCheckMod
   use SolarAbsorbedType  , only : solarabs_type
   use SoilHydrologyType  , only : soilhydrology_type  
   use WaterstateType     , only : waterstate_type
+  use LakestateType      , only : lakestate_type
   use WaterfluxType      , only : waterflux_type
   use IrrigationMod      , only : irrigation_type
   use GlacierSurfaceMassBalanceMod, only : glacier_smb_type
@@ -115,7 +116,7 @@ end function GetBalanceCheckSkipSteps
   !-----------------------------------------------------------------------
   subroutine BeginWaterBalance(bounds, &
        num_nolakec, filter_nolakec, num_lakec, filter_lakec, &
-       soilhydrology_inst, waterstate_inst)
+       soilhydrology_inst, waterstate_inst, lakestate_inst)
     !
     ! !DESCRIPTION:
     ! Initialize column-level water balance at beginning of time step
@@ -126,6 +127,7 @@ subroutine BeginWaterBalance(bounds, &
     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(lakestate_type)      , intent(in)    :: lakestate_inst
     type(soilhydrology_type)  , intent(inout) :: soilhydrology_inst
     type(waterstate_type)     , intent(inout) :: waterstate_inst
     !
@@ -150,10 +152,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, lakestate_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

@@ -142,7 +142,7 @@ subroutine LakeHydrology(bounds, &
          t_soisno             =>  temperature_inst%t_soisno_col         , & ! Output: [real(r8) (:,:) ]  snow temperature (Kelvin)             
          dTdz_top             =>  temperature_inst%dTdz_top_col         , & ! Output: [real(r8) (:)   ]  temperature gradient in top layer K m-1] !TOD 
          snot_top             =>  temperature_inst%snot_top_col         , & ! Output: [real(r8) (:)   ]  snow temperature in top layer [K]  !TODO
-         t_sno_mul_mss        => temperature_inst%t_sno_mul_mss_col     , & ! Output: [real(r8) (:)   ]  col snow temperature multiplied by layer mass, layer sum (K * kg/m2) 
+         t_sno_mul_mss        =>  temperature_inst%t_sno_mul_mss_col     , & ! Output: [real(r8) (:)   ]  col snow temperature multiplied by layer mass, layer sum (K * kg/m2) 
 
          begwb                =>  waterstate_inst%begwb_col             , & ! Input:  [real(r8) (:)   ]  water mass begining of the time step    
          endwb                =>  waterstate_inst%endwb_col             , & ! Output: [real(r8) (:)   ]  water mass end of the time step         
@@ -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, lakestate_inst, &
+           subtract_dynbal_baselines = .false., &
+           water_mass = endwb(bounds%begc:bounds%endc))
 
       do j = 1, nlevgrnd
          do fc = 1, num_lakec

src/biogeophys/LakeTemperatureMod.F90

@@ -133,6 +133,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
     type(energyflux_type)  , intent(inout) :: energyflux_inst
     type(temperature_type) , intent(inout) :: temperature_inst
     type(lakestate_type)   , intent(inout) :: lakestate_inst
+
     !
     ! !LOCAL VARIABLES:
     real(r8), parameter :: p0 = 1._r8                                      ! neutral value of turbulent prandtl number
@@ -246,7 +247,8 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
          t_grnd          =>   temperature_inst%t_grnd_col          , & ! Input:  [real(r8) (:)   ]  ground temperature (Kelvin)             
          t_soisno        =>   temperature_inst%t_soisno_col        , & ! Output: [real(r8) (:,:) ]  soil (or snow) temperature (Kelvin)   
          t_lake          =>   temperature_inst%t_lake_col          , & ! Output: [real(r8) (:,:) ]  col lake temperature (Kelvin)             
-
+         lake_heat       =>   temperature_inst%lake_heat           , & ! Output: [real(r8) (:) ]    col lake heat (J/m²) 
+
          beta            =>   lakestate_inst%betaprime_col         , & ! Output: [real(r8) (:)   ]  col effective beta: sabg_lyr(p,jtop) for snow layers, beta otherwise
          lake_icefrac    =>   lakestate_inst%lake_icefrac_col      , & ! Output: [real(r8) (:,:) ]  col mass fraction of lake layer that is frozen
          lake_icefracsurf => lakestate_inst%lake_icefracsurf_col   , & ! Output: [real(r8) (:,:) ]  col mass fraction of surface lake layer that is frozen
@@ -995,8 +997,13 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
           ncvts(c) = ncvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
                    + cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) 
           fin(c) = fin(c) + phi(c,j)
+
        end do
     end do
+	! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts
+
+	! IV: currently commented out: caused crash. To do: look at this part of the code!!!
+    ! lake_heat(c) = ncvts(c) 
 
     do j = -nlevsno + 1, nlevgrnd
        do fc = 1, num_lakec
@@ -1015,6 +1022,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
     end do
 
 
+
     ! Check energy conservation.
     do fp = 1, num_lakep
        p = filter_lakep(fp)

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)
@@ -113,6 +117,9 @@ module TemperatureType
      real(r8), pointer    :: fact_col              (:,:) ! used in computing tridiagonal matrix
      real(r8), pointer    :: c_h2osfc_col          (:)   ! heat capacity of surface water
 
+     ! lake heat
+     real(r8), pointer    :: lake_heat              (:)  ! total heat of lake water (J/m²)
+
    contains
 
      procedure, public  :: Init
@@ -257,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
@@ -274,6 +282,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%fact_col                 (begc:endc, -nlevsno+1:nlevgrnd)) ; this%fact_col                 (:,:) = nan
     allocate(this%c_h2osfc_col             (begc:endc))                      ; this%c_h2osfc_col             (:)   = nan
 
+    allocate(this%lake_heat                (begc:endc))                      ; this%lake_heat                (:)   = nan
   end subroutine InitAllocate
 
   !------------------------------------------------------------------------
@@ -613,6 +622,11 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp )
             ptr_patch=this%t_veg10_night_patch, default='inactive')
     endif
 
+    ! add lake heat history field here
+    this%lake_heat(begc:endc) = spval
+    call hist_addfld1d (fname='LAKE_HEAT', units='J/m^2', &
+            avgflag='A', long_name='Heat content of gridcell lake water', &
+            ptr_col=this%lake_heat, default='active')
 
   end subroutine InitHistory
 
@@ -840,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
 
   !------------------------------------------------------------------------
@@ -982,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', &