Merge branch 'matthewhoffman/mali/malidev-update-20240814' (PR #6547)

MALI thermal solver bug fix and improved iceberg handling

This merge updates the E3SM-Project/E3SM repo with MALI improvements
from MALI-Dev/develop since February 2024. The primary update needed is
a bug fix to the thermal solver that prevents model crashes in some
situations when ice thickness is very small. Bringing this to master is
in preparation for enabling the MALI thermal solver in E3SM, which is
required for completing ice-sheet/ocean coupling.

Additionally, this brings in improvements to the detection of icebergs
and small islands that makes the velocity solver more robust. There are
a number of other changes related to the subglacial hydrology model and
sea-level model in MALI, both of which are still in development for
coupled E3SM applications.

[NML] for configurations with MALI
[non-BFB] only for configurations with MALI

components/mpas-albany-landice/bld/build-namelist

@@ -663,6 +663,7 @@ add_default($nl, 'config_SGH_bed_roughness');
 add_default($nl, 'config_SGH_bed_roughness_max');
 add_default($nl, 'config_SGH_creep_coefficient');
 add_default($nl, 'config_SGH_englacial_porosity');
+add_default($nl, 'config_SGH_use_iceThicknessHydro');
 add_default($nl, 'config_SGH_chnl_active');
 add_default($nl, 'config_SGH_chnl_include_DCFL');
 add_default($nl, 'config_SGH_chnl_alpha');
@@ -673,6 +674,7 @@ add_default($nl, 'config_SGH_incipient_channel_width');
 add_default($nl, 'config_SGH_include_pressure_melt');
 add_default($nl, 'config_SGH_shmip_forcing');
 add_default($nl, 'config_SGH_basal_melt');
+add_default($nl, 'config_SGH_iter_smooth_waterPressureSlopeNormal');
 
 ##################################
 # Namelist group: AM_globalStats #

components/mpas-albany-landice/bld/build-namelist-section

@@ -225,6 +225,7 @@ add_default($nl, 'config_SGH_bed_roughness');
 add_default($nl, 'config_SGH_bed_roughness_max');
 add_default($nl, 'config_SGH_creep_coefficient');
 add_default($nl, 'config_SGH_englacial_porosity');
+add_default($nl, 'config_SGH_use_iceThicknessHydro');
 add_default($nl, 'config_SGH_chnl_active');
 add_default($nl, 'config_SGH_chnl_include_DCFL');
 add_default($nl, 'config_SGH_chnl_alpha');
@@ -235,6 +236,7 @@ add_default($nl, 'config_SGH_incipient_channel_width');
 add_default($nl, 'config_SGH_include_pressure_melt');
 add_default($nl, 'config_SGH_shmip_forcing');
 add_default($nl, 'config_SGH_basal_melt');
+add_default($nl, 'config_SGH_iter_smooth_waterPressureSlopeNormal');
 
 ##################################
 # Namelist group: AM_globalStats #

components/mpas-albany-landice/bld/namelist_files/albany_input.yaml

@@ -17,6 +17,7 @@ ANONYMOUS:
           Zero Effective Pressure On Floating Ice At Nodes: true
           Zero Beta On Floating Ice: false
           Use Pressurized Bed Above Sea Level: false
+
 # Discretization Description
   Discretization:
     #Exodus Output File Name: albany_output.exo

components/mpas-albany-landice/bld/namelist_files/namelist_defaults_mali.xml

@@ -26,7 +26,7 @@
 
 <!-- solidearth -->
 <config_uplift_method>'none'</config_uplift_method>
-<config_slm_coupling_interval>'0002-00-00_00:00:00'</config_slm_coupling_interval>
+<config_slm_coupling_interval>2</config_slm_coupling_interval>
 <config_MALI_to_SLM_weights_file>'mpas_to_grid.nc'</config_MALI_to_SLM_weights_file>
 <config_SLM_to_MALI_weights_file>'grid_to_mpas.nc'</config_SLM_to_MALI_weights_file>
 
@@ -184,6 +184,7 @@
 <config_SGH_bed_roughness_max>0.1</config_SGH_bed_roughness_max>
 <config_SGH_creep_coefficient>0.04</config_SGH_creep_coefficient>
 <config_SGH_englacial_porosity>0.01</config_SGH_englacial_porosity>
+<config_SGH_use_iceThicknessHydro>.true.</config_SGH_use_iceThicknessHydro>
 <config_SGH_chnl_active>.false.</config_SGH_chnl_active>
 <config_SGH_chnl_include_DCFL>.false.</config_SGH_chnl_include_DCFL>
 <config_SGH_chnl_alpha>1.25</config_SGH_chnl_alpha>
@@ -194,6 +195,7 @@
 <config_SGH_include_pressure_melt>.true.</config_SGH_include_pressure_melt>
 <config_SGH_shmip_forcing>'none'</config_SGH_shmip_forcing>
 <config_SGH_basal_melt>'file'</config_SGH_basal_melt>
+<config_SGH_iter_smooth_waterPressureSlopeNormal>1</config_SGH_iter_smooth_waterPressureSlopeNormal>
 
 <!-- AM_globalStats -->
 <config_AM_globalStats_enable>.true.</config_AM_globalStats_enable>

components/mpas-albany-landice/bld/namelist_files/namelist_definition_mali.xml

@@ -194,11 +194,11 @@ Valid values: 'none', 'data', 'sealevelmodel'
 Default: Defined in namelist_defaults.xml
 </entry>
 
-<entry id="config_slm_coupling_interval" type="char*1024"
+<entry id="config_slm_coupling_interval" type="integer"
 	category="solidearth" group="solidearth">
-Time interval at which the sea-level model is called by MALI. The interval has to be an even multiple of the option 'config_adaptive_timestep_force_interval
+Time interval at which the sea-level model is called by MALI. Only integer years are currently supported. The user must set 'dt1' in namelist.sealevel to match this value.  Note that the user needs to set config_adaptive_timestep_force_interval to divide evenly into config_slm_coupling_interval.  Also, restart file output_interval needs to be a multiple of config_slm_coupling_interval.
 
-Valid values: Any time interval of the format 'YYYY-MM-DD_HH:MM:SS'
+Valid values: Any positive integer
 Default: Defined in namelist_defaults.xml
 </entry>
 
@@ -1218,7 +1218,7 @@ Default: Defined in namelist_defaults.xml
 
 <entry id="config_SGH_tangent_slope_calculation" type="char*1024"
 	category="subglacial_hydro" group="subglacial_hydro">
-Selection of the method for calculating the tangent component of slope at edges. 'from_vertex_barycentric' interpolates scalar values from cell centers to vertices using the barycentric interpolation routine in operators (mpas_cells_to_points_using_baryweights) and then calculates the slope between vertices.  It works for obtuse triangles, but will not work correctly across the edges of periodic meshes. 'from_vertex_barycentric_kiteareas' interpolates scalar values from cell centers to vertices using barycentric interpolation based on kiterea values and then calculates the slope between vertices.  It will work across the edges of periodic meshes, but will not work correctly for obtuse triangles. 'from_normal_slope' uses the vector operator mpas_tangential_vector_1d to calculate the tangent slopes from the normal slopes on the edges of the adjacent cells.  It will work for any mesh configuration, but is the least accurate method.
+Selection of the method for calculating the tangent component of slope at edges. 'from_vertex_barycentric' interpolates scalar values from cell centers to vertices using the barycentric interpolation routine in operators (mpas_cells_to_points_using_baryweights) and then calculates the slope between vertices.  It works for obtuse triangles, but will not work correctly across the edges of periodic meshes. 'from_vertex_barycentric_kiteareas' interpolates scalar values from cell centers to vertices using barycentric interpolation based on kiterea values and then calculates the slope between vertices.  It will work across the edges of periodic meshes, but will not work correctly for obtuse triangles. 'from_normal_slope' uses the vector operator mpas_tangential_vector_1d to calculate the tangent slopes from the normal slopes on the edges of the adjacent cells.  It will work for any mesh configuration. 'from_normal_slope' uses a larger stencil, so may therefore produce a smoother 'gradMagPhiEdge' field. Detailed testing yielded nearly identical results between 'from_normal_slope' and 'from_vertex_barycentric' methods, but 'from_normal_slope' seemed to produce slightly more stable results at the grounding line.
 
 Valid values: 'from_vertex_barycentric', 'from_vertex_barycentric_kiteareas', 'from_normal_slope'
 Default: Defined in namelist_defaults.xml
@@ -1320,6 +1320,14 @@ Valid values: positive real number
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_SGH_use_iceThicknessHydro" type="logical"
+	category="subglacial_hydro" group="subglacial_hydro">
+Option to use an altered ice thickness field called iceThicknessHydro that replaces local maxima/minima in upperSurface with a mean of the cells neighbors. This option has no significant effect on the behavior of the model but makes it more stable.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
 <entry id="config_SGH_chnl_active" type="logical"
 	category="subglacial_hydro" group="subglacial_hydro">
 activate channels in subglacial hydrology model
@@ -1400,6 +1408,14 @@ Valid values: 'file', 'thermal', 'basal_heat'
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_SGH_iter_smooth_waterPressureSlopeNormal" type="integer"
+	category="subglacial_hydro" group="subglacial_hydro">
+number of iterations to smooth waterPressure over when calculating waterPressureSlopeNormal. Used only to keep channelPressureFreeze stable and will not affect other aspects of the model that rely on waterPressure.
+
+Valid values: positive integer or zero
+Default: Defined in namelist_defaults.xml
+</entry>
+
 
 <!-- AM_globalStats -->
 

components/mpas-albany-landice/src/Registry.xml

@@ -145,9 +145,9 @@
 		            description="Selection of the method for bedrock uplift calculation."
 		            possible_values="'none', 'data', 'sealevelmodel'"
         />
-        <nml_option name="config_slm_coupling_interval" type="character" default_value="0002-00-00_00:00:00" units="time"
-            description="Time interval at which the sea-level model is called by MALI. The interval has to be an even multiple of the option 'config_adaptive_timestep_force_interval"
-            possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS'"
+        <nml_option name="config_slm_coupling_interval" type="integer" default_value="2" units="years"
+            description="Time interval at which the sea-level model is called by MALI. Only integer years are currently supported. The user must set 'dt1' in namelist.sealevel to match this value.  Note that the user needs to set config_adaptive_timestep_force_interval to divide evenly into config_slm_coupling_interval.  Also, restart file output_interval needs to be a multiple of config_slm_coupling_interval."
+            possible_values="Any positive integer"
 		/>
         <nml_option name="config_MALI_to_SLM_weights_file" type="character" default_value="mpas_to_grid.nc" units="unitless"
             description="File containing the interpolation weights for regridding from MPAS mesh to the Gaussian grid used by the Sea Level Model."
@@ -842,6 +842,7 @@
 			<var name="sfcMassBal"/>
                         <var name="floatingBasalMassBal"/>
                         <var name="regionCellMasks"/>
+                        <var name="regionEdgeMasks"/>
                         <!-- If restart file sizes become too large, eigencalvingParameter,
                              groundedVonMisesThresholdStress, and floatingVonMisesThresholdStress
                              could be moved to packages. -->
@@ -869,7 +870,8 @@
                         <var name="groundedBasalMassBal" packages="hydro"/>  <!-- This is only needed if config_SGH_basal_melt='thermal', but adding it always-->
 			<var name="basalMeltInput" packages="hydro"/>  <!-- This is NOT needed if config_SGH_basal_melt='thermal', but adding it always-->
                         <var name="externalWaterInput" packages="hydro"/>
-                        <!-- SGH state variables: -->
+                        <!-- SGH variables: -->
+                        <var name="deltatSGH" packages="hydro"/>
 			<var name="waterThickness" packages="hydro"/>
 			<var name="tillWaterThickness" packages="hydro"/>
 			<var name="waterPressure" packages="hydro"/>
@@ -935,6 +937,7 @@
                                 input_interval="initial_only"
                                 runtime_format="single_file">
                         <var name="regionCellMasks"/>
+                        <var name="regionEdgeMasks"/>
                         <!--var name="nCellsInRegion"/>
                         <var name="regionVertexMasks"/>
                         <var name="nVerticesInRegion"/>
@@ -1125,10 +1128,10 @@ is the value of that variable from the *previous* time level!
                 />
 		<!-- Variables only needed by SIA solver -->
 		<var name="baryCellsOnVertex" type="integer" dimensions="R3 nVertices" units="unitless"
-			 description="Cell center indices to use for interpolating from cell centers to vertex locations.  Note these are local indices!" packages="SIAvelocity"
+			 description="Cell center indices to use for interpolating from cell centers to vertex locations.  Note these are local indices!" packages="SIAvelocity;hydro"
 		/>
 		<var name="baryWeightsOnVertex" type="real" dimensions="R3 nVertices" units="unitless"
-			 description="Weights to interpolate from cell centers to vertex locations.  Each weight is used with the corresponding cell center index indentified by baryCellsOnVertex." packages="SIAvelocity"
+			 description="Weights to interpolate from cell centers to vertex locations.  Each weight is used with the corresponding cell center index indentified by baryCellsOnVertex." packages="SIAvelocity;hydro"
 		/>
 		<!-- Variables only needed by HO solver -->
 		<var name="wachspressWeightVertex" type="real" dimensions="maxEdges nCells" units="unitless"
@@ -1734,6 +1737,9 @@ is the value of that variable from the *previous* time level!
 	  <var name="regionCellMasks" type="integer" dimensions="nRegions nCells" units="unitless" default_value="1"
 	       description="masks set to 1 in cells that fall within a given region and 0 otherwise"
 	       />
+	  <var name="regionEdgeMasks" type="integer" dimensions="nRegions nEdges" units="unitless" default_value="1"
+	       description="masks set to 1 in edges that fall within a given region and 0 otherwise"
+	       />
 	  <!--var name="nCellsInRegion" type="integer" dimensions="nRegions" units="unitless"
 	       description="number of cells that fall within a given region"
 	       />