ESCOMP · billsacks · Aug 6, 2019 · Jul 10, 2019 · Jul 10, 2019 · Jul 22, 2019
diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm
@@ -1682,7 +1682,12 @@ sub process_namelist_inline_logic {
   #######################################################################
   # namelist groups: clm_hydrology1_inparm and clm_soilhydrology_inparm #
   #######################################################################
-  setup_logic_hydrology_switches($nl);
+  setup_logic_hydrology_switches($opts, $nl_flags, $definition, $defaults, $nl);
+
+  #######################################################################
+  # namelist group: scf_clm5_inparm                                     #
+  #######################################################################
+  setup_logic_scf_clm5($opts, $nl_flags, $definition, $defaults, $nl);
 
   #########################################
   # namelist group: clm_initinterp_inparm #
@@ -2720,16 +2725,22 @@ sub setup_logic_hydrology_switches {
   #
   # Check on Switches for hydrology
   #
-  my ($nl) = @_;
+  my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
 
-  my $subgrid    = $nl->get_value('subgridflag' );
+  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_subgrid_fluxes');
+  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_cover_fraction_method');
+  my $subgrid    = $nl->get_value('use_subgrid_fluxes' );
   my $origflag   = $nl->get_value('origflag'    );
   my $h2osfcflag = $nl->get_value('h2osfcflag'  );
-  if ( $origflag == 1 && $subgrid == 1 ) {
-    $log->fatal_error("if origflag is ON, subgridflag can NOT also be on!");
+  my $scf_method = $nl->get_value('snow_cover_fraction_method');
+  if ( $origflag == 1 && &value_is_true($subgrid) ) {
+    $log->fatal_error("if origflag is ON, use_subgrid_fluxes can NOT also be on!");
+  }
+  if ( $h2osfcflag == 1 && ! &value_is_true($subgrid) ) {
+    $log->fatal_error("if h2osfcflag is ON, use_subgrid_fluxes can NOT be off!");
   }
-  if ( $h2osfcflag == 1 && $subgrid != 1 ) {
-    $log->fatal_error("if h2osfcflag is ON, subgridflag can NOT be off!");
+  if ( remove_leading_and_trailing_quotes($scf_method) eq 'ny07' && &value_is_true($subgrid) ) {
+     $log->fatal_error("snow_cover_fraction_method ny07 is incompatible with use_subgrid_fluxes");
   }
   # Test bad configurations
   my $lower   = $nl->get_value( 'lower_boundary_condition'  );
@@ -3498,8 +3509,6 @@ sub setup_logic_snowpack {
               'structure'=>$nl_flags->{'structure'});
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'h2osno_max',
               'structure'=>$nl_flags->{'structure'});
-  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'int_snow_max');
-  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'n_melt_glcmec');
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'wind_dependent_snow_density');
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_overburden_compaction_method');
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'lotmp_snowdensity_method');
@@ -3523,6 +3532,26 @@ sub setup_logic_snowpack {
 
 #-------------------------------------------------------------------------------
 
+sub setup_logic_scf_clm5 {
+   # Options related to the clm5 snow cover fraction method
+  my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
+
+  if (remove_leading_and_trailing_quotes($nl->get_value('snow_cover_fraction_method')) eq 'clm5') {
+     add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'int_snow_max');
+     add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'n_melt_glcmec');
+  }
+  else {
+     if (defined($nl->get_value('int_snow_max'))) {
+        $log->fatal_error('int_snow_max is set, but only applies for snow_cover_fraction_method=clm5');
+     }
+     if (defined($nl->get_value('n_melt_glcmec'))) {
+        $log->fatal_error('n_melt_glcmec is set, but only applies for snow_cover_fraction_method=clm5');
+     }
+  }
+}
+
+#-------------------------------------------------------------------------------
+
 sub setup_logic_atm_forcing {
    #
    # Options related to atmospheric forcings
@@ -3667,6 +3696,9 @@ sub write_output_files {
   push @groups, "lifire_inparm";
   push @groups, "ch4finundated";
   push @groups, "clm_canopy_inparm";
+  if (remove_leading_and_trailing_quotes($nl->get_value('snow_cover_fraction_method')) eq 'clm5') {
+     push @groups, "scf_clm5_inparm";
+  }
 
   my $outfile;
   $outfile = "$opts->{'dir'}/lnd_in";

diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml
@@ -329,6 +329,8 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
      to .true. -->
 <reset_snow_glc_ela>1.e9</reset_snow_glc_ela>
 
+<snow_cover_fraction_method>clm5</snow_cover_fraction_method>
+
 <!-- Default glacier behavior is:
      Mountain glaciers: single_at_atm_topo
      Greenland - inside CISM grid but outside Greenland itself: virtual
@@ -436,6 +438,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <constrain_stress_deciduous_onset phys="clm4_5">.false.</constrain_stress_deciduous_onset>
 <constrain_stress_deciduous_onset phys="clm5_0">.true.</constrain_stress_deciduous_onset>
 
+<!-- Whether to use subgrid fluxes for snow -->
+<use_subgrid_fluxes>.true.</use_subgrid_fluxes>
+
 <!-- dynamic roots -->
 <use_dynroot phys="clm4_5">.false.</use_dynroot>
 <use_dynroot phys="clm5_0" bgc_mode="sp">.false.</use_dynroot>

diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml
@@ -910,9 +910,9 @@ Human heat stress indices:
 If TRUE, write diagnostic of global radiative temperature written to CLM log file.
 </entry>
 
-<entry id="subgridflag" type="integer" category="clm_physics"
-       group="clm_inparm" valid_values="0,1" >
-Subgrid fluxes for snow
+<entry id="use_subgrid_fluxes" type="logical" category="clm_physics"
+       group="clm_inparm" >
+Whether to use subgrid fluxes for snow
 </entry>
 
 <entry id="snowveg_affects_radiation" type="logical" category="clm_physics"
@@ -2443,18 +2443,20 @@ Changes in this value should possibly be accompanied by changes in:
 </entry>
 
 <entry id="int_snow_max" type="real" category="clm_physics"
-       group="clm_inparm" valid_values="" >
+       group="scf_clm5_inparm" valid_values="" >
 Limit applied to integrated snowfall when determining changes in snow-covered fraction during melt
 (mm H2O)
+Only applies when using the CLM5 snow cover fraction method
 </entry>
 
 <entry id="n_melt_glcmec" type="real" category="clm_physics"
-       group="clm_inparm" valid_values="" >
+       group="scf_clm5_inparm" valid_values="" >
 SCA shape parameter for glc_mec (glacier multiple elevation class) columns
 For most columns, n_melt is based on the standard deviation of 1km topography in the grid cell;
 but glc_mec columns already account for subgrid topographic variability through their use of
 multiple elevation classes; thus, to avoid double-accounting for topographic variability
 in these columns, we use a fixed value of n_melt.
+Only applies when using the CLM5 snow cover fraction method
 </entry>
 
 <entry id="wind_dependent_snow_density" type="logical" category="clm_physics"
@@ -2549,10 +2551,11 @@ below a given elevation, you risk forcing the system to evolve
 differently in areas below and above reset_snow_glc_ela.
 </entry>
 
-<entry id="oldfflag" type="integer" category="clm_physics"
-       group="clm_snowhydrology_inparm" valid_values="0,1" >
-Use old snow cover fraction from Niu et al. 2007
-(deprecated -- will be removed)
+<entry id="snow_cover_fraction_method" type="char*64" category="clm_physics"
+       group="clm_inparm" valid_values="ny07,clm5" >
+Parameterization to use for snow cover fraction
+ny07: Niu and Yang 2007
+clm5: Newer method, default in CLM5
 </entry>
 
 <!-- ========================================================================================  -->

diff --git a/cime_config/SystemTests/lwiso.py b/cime_config/SystemTests/lwiso.py
@@ -0,0 +1,44 @@
+"""
+Implementation of the CIME LWISO (Land Water Isotope) test.
+
+This is a CTSM specific test:
+Verifies turning on water isotopes doesn't change answers
+(1) do a run with water isotopes on
+(2) do a run with water isotopes off
+
+The reason for having water isotopes on in the first case (rather than the other way
+around) is that this results in having water isotope quantities present in the baselines.
+
+"""
+
+from CIME.SystemTests.system_tests_compare_two import SystemTestsCompareTwo
+from CIME.XML.standard_module_setup import *
+from CIME.SystemTests.test_utils.user_nl_utils import append_to_user_nl_files
+
+logger = logging.getLogger(__name__)
+
+class LWISO(SystemTestsCompareTwo):
+
+    def __init__(self, case):
+        SystemTestsCompareTwo.__init__(self, case,
+                                       separate_builds = False,
+                                       run_two_suffix = 'nowiso',
+                                       run_one_description = 'water isotopes on',
+                                       run_two_description = 'water isotopes off',
+                                       ignore_fieldlist_diffs = True)
+
+    def _case_one_setup(self):
+        # BUG(wjs, 2019-07-30, ESCOMP/ctsm#495) We currently can't turn on actual water
+        # isotopes in a multi-timestep test, so we're setting
+        # enable_water_tracer_consistency_checks rather than enable_water_isotopes;
+        # eventually, though, we should change this to the latter. (See
+        # <https://github.com/ESCOMP/ctsm/issues/495#issuecomment-516619853>.)
+        append_to_user_nl_files(caseroot = self._get_caseroot(),
+                                component = "clm",
+                                contents = "enable_water_tracer_consistency_checks=.true.")
+
+    def _case_two_setup(self):
+        append_to_user_nl_files(caseroot = self._get_caseroot(),
+                                component = "clm",
+                                contents = "enable_water_tracer_consistency_checks=.false.")
+
diff --git a/cime_config/config_compsets.xml b/cime_config/config_compsets.xml
@@ -143,6 +143,10 @@
     <alias>I2000Clm50BgcCropQianRsGs</alias>
     <lname>2000_DATM%QIA_CLM50%BGC-CROP_SICE_SOCN_SROF_SGLC_SWAV</lname>
   </compset>
+  <compset>
+    <alias>I2000Clm45BgcCropQianRsGs</alias>
+    <lname>2000_DATM%QIA_CLM45%BGC-CROP_SICE_SOCN_SROF_SGLC_SWAV</lname>
+  </compset>
 
   <compset>
     <alias>I2000Clm50BgcCruGs</alias>

diff --git a/cime_config/config_tests.xml b/cime_config/config_tests.xml
@@ -75,6 +75,16 @@ SSP    smoke CLM spinup test (only valid for CLM compsets with CLM45)
     <HIST_N>$STOP_N</HIST_N>
   </test>
 
+  <test NAME="LWISO">
+    <DESC>CTSM Land model test to ensure Water isotopes on or off does NOT change answers</DESC>
+    <INFO_DBUG>1</INFO_DBUG>
+    <DOUT_S>FALSE</DOUT_S>
+    <CONTINUE_RUN>FALSE</CONTINUE_RUN>
+    <REST_OPTION>never</REST_OPTION>
+    <HIST_OPTION>$STOP_OPTION</HIST_OPTION>
+    <HIST_N>$STOP_N</HIST_N>
+  </test>
+
   <test NAME="SSP">
     <DESC>smoke CLM spinup test</DESC>
     <INFO_DBUG>1</INFO_DBUG>

diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml
@@ -1186,6 +1186,16 @@
       <option name="wallclock">00:20:00</option>
     </options>
   </test>
+  <test name="ERP_D_P36x2_Ld3" grid="f10_f10_musgs" compset="I2000Clm45BgcCrop" testmods="clm/no_subgrid_fluxes">
+    <machines>
+      <machine name="cheyenne" compiler="gnu" category="aux_clm">
+        <options>
+          <option name="wallclock">00:20:00</option>
+          <option name="comment">This covers some code that isn't covered by any existing tests (such as the oldhyd test), though the amount of additional code coverage is small, so we don't necessarily need to keep this test long-term.</option>
+        </options>
+      </machine>
+    </machines>
+  </test>
   <test name="SMS_D_Ld1_P48x1" grid="f10_f10_musgs" compset="I2000Clm50BgcCruGs" testmods="clm/af_bias_v7">
     <machines>
       <machine name="hobart" compiler="nag" category="aux_clm"/>
@@ -1674,6 +1684,16 @@
 
     </machines>
   </test>
+  <test name="LWISO_Ld10" grid="f10_f10_musgs" compset="I2000Clm50BgcCropGs" testmods="clm/coldStart">
+    <machines>
+      <machine name="cheyenne" compiler="gnu" category="aux_clm">
+        <options>
+          <option name="wallclock">0:30:00</option>
+          <option name="comment">Ensure that turning on water tracers doesn't change answers. Cold start for now, until we can use initial conditions from a non-isotope case in an isotope case; once we can do that, this should be changed to not be cold start (e.g., 5-day decStart transient test: see also https://github.com/ESCOMP/ctsm/issues/495#issuecomment-516619853).</option>
+        </options>
+      </machine>
+    </machines>
+  </test>
   <test name="ERP_P36x2_D_Ld5" grid="f10_f10_musgs" compset="I2000Ctsm50NwpSpGswpGs" testmods="clm/default">
     <machines>
 

diff --git a/cime_config/testdefs/testmods_dirs/clm/no_subgrid_fluxes/README b/cime_config/testdefs/testmods_dirs/clm/no_subgrid_fluxes/README
@@ -0,0 +1,8 @@
+The purpose of this testmod is to exercise the newer snow cover fraction
+method along with having use_subgrid_fluxes = .false.
+
+Note that we also need h2osfcflag = 0 in this testmod, because
+h2osfcflag = 1 is incompatible with use_subgrid_fluxes = .false. Also,
+this needs to be used in a clm45 configuration for compatibility with
+other options (or else would need to set some other flags:
+lower_boundary_condition and possibly others)
diff --git a/cime_config/testdefs/testmods_dirs/clm/no_subgrid_fluxes/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/no_subgrid_fluxes/include_user_mods
@@ -0,0 +1 @@
+../default
diff --git a/cime_config/testdefs/testmods_dirs/clm/no_subgrid_fluxes/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/no_subgrid_fluxes/user_nl_clm
@@ -0,0 +1,2 @@
+ h2osfcflag  = 0
+ use_subgrid_fluxes = .false.
diff --git a/cime_config/testdefs/testmods_dirs/clm/oldhyd/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/oldhyd/include_user_mods
@@ -0,0 +1 @@
+../default
diff --git a/cime_config/testdefs/testmods_dirs/clm/oldhyd/shell_commands b/cime_config/testdefs/testmods_dirs/clm/oldhyd/shell_commands
diff --git a/cime_config/testdefs/testmods_dirs/clm/oldhyd/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/oldhyd/user_nl_clm
@@ -1,5 +1,4 @@
- hist_ndens  = 1,1
- oldfflag    = 1
+ snow_cover_fraction_method = 'ny07'
  h2osfcflag  = 0
  origflag    = 1
- subgridflag = 0
+ use_subgrid_fluxes = .false.
diff --git a/src/biogeophys/CMakeLists.txt b/src/biogeophys/CMakeLists.txt
@@ -15,6 +15,7 @@ list(APPEND clm_sources
   QSatMod.F90
   RootBiophysMod.F90
   SaturatedExcessRunoffMod.F90
+  SnowCoverFractionMod.F90
   SnowHydrologyMod.F90
   SnowSnicarMod.F90
   SoilHydrologyMod.F90

diff --git a/src/biogeophys/HydrologyNoDrainageMod.F90 b/src/biogeophys/HydrologyNoDrainageMod.F90
@@ -100,7 +100,7 @@ subroutine HydrologyNoDrainage(bounds, &
        wateratm2lndbulk_inst, &
        waterfluxbulk_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, &
        soilhydrology_inst, saturated_excess_runoff_inst, infiltration_excess_runoff_inst, &
-       aerosol_inst, canopystate_inst, soil_water_retention_curve, topo_inst)
+       aerosol_inst, canopystate_inst, scf_method, soil_water_retention_curve, topo_inst)
     !
     ! !DESCRIPTION:
     ! This is the main subroutine to execute the calculation of soil/snow
@@ -122,7 +122,8 @@ subroutine HydrologyNoDrainage(bounds, &
     use SoilHydrologyMod     , only : UpdateUrbanPonding
     use SoilHydrologyMod     , only : WaterTable, PerchedWaterTable
     use SoilHydrologyMod     , only : ThetaBasedWaterTable, RenewCondensation
-    use SoilWaterMovementMod , only : SoilWater 
+    use SoilWaterMovementMod , only : SoilWater
+    use SnowCoverFractionMod , only : snow_cover_fraction_type
     use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type
     use SoilWaterMovementMod , only : use_aquifer_layer
     use SoilWaterPlantSinkMod , only : Compute_EffecRootFrac_And_VertTranSink
@@ -155,6 +156,7 @@ subroutine HydrologyNoDrainage(bounds, &
     type(saturated_excess_runoff_type), intent(inout) :: saturated_excess_runoff_inst
     type(infiltration_excess_runoff_type), intent(inout) :: infiltration_excess_runoff_inst
     type(canopystate_type)   , intent(inout) :: canopystate_inst
+    class(snow_cover_fraction_type), intent(in) :: scf_method
     class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve
     class(topo_type)   , intent(in)    :: topo_inst
     !
@@ -322,6 +324,7 @@ subroutine HydrologyNoDrainage(bounds, &
 
       ! Natural compaction and metamorphosis.
       call SnowCompaction(bounds, num_snowc, filter_snowc, &
+           scf_method, &
            temperature_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, atm2lnd_inst)
 
       ! Combine thin snow elements

diff --git a/src/biogeophys/LakeHydrologyMod.F90 b/src/biogeophys/LakeHydrologyMod.F90
@@ -50,6 +50,7 @@ module LakeHydrologyMod
   subroutine LakeHydrology(bounds, &
        num_lakec, filter_lakec, num_lakep, filter_lakep, &
        num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc, &
+       scf_method, &
        atm2lnd_inst, temperature_inst, soilstate_inst, waterstatebulk_inst, &
        waterdiagnosticbulk_inst, waterbalancebulk_inst, waterfluxbulk_inst, &
        wateratm2lndbulk_inst, energyflux_inst, aerosol_inst, lakestate_inst, topo_inst)
@@ -80,7 +81,8 @@ subroutine LakeHydrology(bounds, &
     use SnowHydrologyMod, only : ZeroEmptySnowLayers, BuildSnowFilter, SnowCapping
     use SnowHydrologyMod, only : DivideSnowLayers, NewSnowBulkDensity
     use LakeCon         , only : lsadz
-    use TopoMod, only : topo_type
+    use TopoMod         , only : topo_type
+    use SnowCoverFractionMod, only : snow_cover_fraction_type
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds  
@@ -92,6 +94,7 @@ subroutine LakeHydrology(bounds, &
     integer                , intent(out)   :: filter_shlakesnowc(:)   ! column filter for snow points
     integer                , intent(out)   :: num_shlakenosnowc       ! number of column non-snow points
     integer                , intent(out)   :: filter_shlakenosnowc(:) ! column filter for non-snow points
+    class(snow_cover_fraction_type), intent(in) :: scf_method
     type(atm2lnd_type)     , intent(in)    :: atm2lnd_inst
     type(temperature_type) , intent(inout) :: temperature_inst
     type(soilstate_type)   , intent(in)    :: soilstate_inst
@@ -435,6 +438,7 @@ subroutine LakeHydrology(bounds, &
       ! Natural compaction and metamorphosis.
 
       call SnowCompaction(bounds, num_shlakesnowc, filter_shlakesnowc, &
+           scf_method, &
            temperature_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, atm2lnd_inst)
 
       ! Combine thin snow elements