Adding excess ground ice to CTSM

bld/CLMBuildNamelist.pm

@@ -1732,6 +1732,11 @@ sub process_namelist_inline_logic {
   # namelist group: clm_initinterp_inparm #
   #########################################
   setup_logic_initinterp($opts, $nl_flags, $definition, $defaults, $nl);
+
+  ###############################
+  # namelist group: exice_streams   #
+  ###############################
+  setup_logic_exice($opts, $nl_flags, $definition, $defaults, $nl);
 }
 
 #-------------------------------------------------------------------------------
@@ -2132,6 +2137,7 @@ sub setup_logic_soilstate {
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'organic_frac_squared' );
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_bedrock',
               'use_fates'=>$nl_flags->{'use_fates'}, 'vichydro'=>$nl_flags->{'vichydro'} );
+  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_excess_ice'); # excess ice flag should be read before stream vars
 
   my $var1 = "soil_layerstruct_predefined";
   my $var2 = "soil_layerstruct_userdefined";
@@ -4137,6 +4143,24 @@ sub setup_logic_fates {
     }
 }
 
+#-------------------------------------------------------------------------------
+sub setup_logic_exice {
+  #
+  # excess ice streams
+  #
+  my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
+  my $use_exice = $nl->get_value( 'use_excess_ice'               );
+  if (defined($use_exice) && $opts->{'driver'} eq "nuopc" && value_is_true($use_exice)) {
+    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_exice');
+    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_exice');
+    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_mapalgo_exice');
+  } elsif (defined($use_exice) && (! $opts->{'driver'} eq "nuopc") && value_is_true($use_exice)) {
+       $log->fatal_error("nuopc driver is required when use_excess_ice is set to true" );
+}
+
+
+} # end exice streams
+
 #-------------------------------------------------------------------------------
 
 sub setup_logic_misc {
@@ -4204,6 +4228,7 @@ sub write_output_files {
   push @groups, "nitrif_inparm";
   push @groups, "lifire_inparm";
   push @groups, "ch4finundated";
+  push @groups, "exice_streams";
   push @groups, "soilbgc_decomp";
   push @groups, "clm_canopy_inparm";
   if (remove_leading_and_trailing_quotes($nl->get_value('snow_cover_fraction_method')) eq 'SwensonLawrence2012') {

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -2586,4 +2586,14 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts
 
 <init_interp_method>general</init_interp_method>
 
+<!-- =========================================  -->
+<!-- Defaults for excess_ice                    -->
+<!-- =========================================  -->
+
+<use_excess_ice>.false.</use_excess_ice>
+
+<stream_fldfilename_exice use_excess_ice=".true.">lnd/clm2/paramdata/exice_init_0.125x0.125_c20220516.nc</stream_fldfilename_exice>
+<stream_meshfile_exice    use_excess_ice=".true.">lnd/clm2/paramdata/exice_init_0.125x0.125_ESMFmesh_c20220516.nc</stream_meshfile_exice>
+<stream_mapalgo_exice     use_excess_ice=".true.">bilinear</stream_mapalgo_exice>
+
 </namelist_defaults>

bld/namelist_files/namelist_definition_ctsm.xml

@@ -2818,4 +2818,30 @@ use case.)
 
 </entry>
 
+<!-- ========================================================================================  -->
+<!-- Namelist options related to excess ice                                                    -->
+<!-- ========================================================================================  -->
+<entry id="use_excess_ice" type="logical" category="clm_physics"
+        group="clm_inparm" valid_values="" >
+If TRUE turn on the excess ice physics, (Lee et al., 2014; Cai et al., 2020)
+</entry>
+
+<entry id="stream_fldfilename_exice" type="char*256(30)" category="datasets"
+       input_pathname="abs" group="exice_streams" valid_values="" >
+Filename of input stream data for excess ice data
+</entry>
+
+<entry id="stream_meshfile_exice" type="char*256" category="datasets"
+       input_pathname="abs" group="exice_streams" valid_values="" >
+mesh filename of input stream data for excess ice
+</entry>
+
+<entry id="stream_mapalgo_exice" type="char*256" category="datasets"
+       group="exice_streams" valid_values="bilinear,nn" >
+Mapping method from excess ice input stream data to the model resolution
+    bilinear = bilinear interpolation
+    nn       = nearest neighbor
+</entry>
+
+
 </namelist_definition>

src/biogeophys/SoilHydrologyMod.F90

@@ -179,6 +179,7 @@ subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, &
     integer :: j, fc, c
     real(r8) :: vol_ice(bounds%begc:bounds%endc,1:nlevsoi) !partial volume of ice lens in layer
     real(r8) :: icefrac_orig ! original formulation for icefrac
+    real(r8) :: dz_ext(bounds%begc:bounds%endc,1:nlevsoi)
 
     character(len=*), parameter :: subname = 'SetSoilWaterFractions'
     !-----------------------------------------------------------------------
@@ -191,6 +192,7 @@ subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, &
 
          h2osoi_liq       =>    waterstatebulk_inst%h2osoi_liq_col      , & ! Input:  [real(r8) (:,:) ]  liquid water (kg/m2)
          h2osoi_ice       =>    waterstatebulk_inst%h2osoi_ice_col      , & ! Input:  [real(r8) (:,:) ]  ice water (kg/m2)
+         excess_ice       =>    waterstatebulk_inst%excess_ice_col      , & ! Input:  [real(r8) (:,:) ]  excess ice (kg/m2)
 
          origflag         =>    soilhydrology_inst%origflag         , & ! Input:  logical
          icefrac          =>    soilhydrology_inst%icefrac_col      , & ! Output: [real(r8) (:,:) ]                                                  
@@ -203,7 +205,8 @@ subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, &
 
           ! Porosity of soil, partial volume of ice and liquid, fraction of ice in each layer,
           ! fractional impermeability
-          vol_ice(c,j) = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))
+          dz_ext(c,j) = dz(c,j) + excess_ice(c,j)/denice ! extended layer thickness, should be good for all the columns
+          vol_ice(c,j) = min(watsat(c,j), (h2osoi_ice(c,j) + excess_ice(c,j))/(dz_ext(c,j)*denice))
           eff_porosity(c,j) = max(0.01_r8,watsat(c,j)-vol_ice(c,j))
           icefrac(c,j) = min(1._r8,vol_ice(c,j)/watsat(c,j))