Update MALI version to include higher-order advection and time integration

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

@@ -456,6 +456,8 @@ if ($MALI_DYNAMIC eq 'TRUE') {
 } else {
     add_default($nl, 'config_tracer_advection', 'val'=>"none");
 }
+add_default($nl, 'config_horiz_tracer_adv_order');
+add_default($nl, 'config_advection_coef_3rd_order');
 add_default($nl, 'config_restore_thickness_after_advection');
 add_default($nl, 'config_zero_sfcMassBalApplied_over_bare_land');
 
@@ -509,6 +511,7 @@ add_default($nl, 'config_damagecalvingParameter');
 add_default($nl, 'config_ismip6_retreat_k');
 add_default($nl, 'config_calving_error_threshold');
 add_default($nl, 'config_distribute_unablatedVolumeDynCell');
+add_default($nl, 'config_update_velocity_before_calving');
 
 ##################################
 # Namelist group: thermal_solver #
@@ -574,6 +577,8 @@ add_default($nl, 'config_dynamic_thickness');
 
 add_default($nl, 'config_dt');
 add_default($nl, 'config_time_integration');
+add_default($nl, 'config_rk_order');
+add_default($nl, 'config_rk3_stages');
 add_default($nl, 'config_adaptive_timestep');
 add_default($nl, 'config_min_adaptive_timestep');
 add_default($nl, 'config_max_adaptive_timestep');
@@ -635,6 +640,7 @@ add_default($nl, 'config_print_calving_info');
 add_default($nl, 'config_print_thermal_info');
 add_default($nl, 'config_always_compute_fem_grid');
 add_default($nl, 'config_print_velocity_cleanup_details');
+add_default($nl, 'config_check_tracer_monotonicity');
 
 ####################################
 # Namelist group: subglacial_hydro #

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

@@ -22,6 +22,8 @@ add_default($nl, 'config_effective_pressure_max');
 
 add_default($nl, 'config_thickness_advection');
 add_default($nl, 'config_tracer_advection');
+add_default($nl, 'config_horiz_tracer_adv_order');
+add_default($nl, 'config_advection_coef_3rd_order');
 add_default($nl, 'config_restore_thickness_after_advection');
 add_default($nl, 'config_zero_sfcMassBalApplied_over_bare_land');
 
@@ -71,6 +73,7 @@ add_default($nl, 'config_damagecalvingParameter');
 add_default($nl, 'config_ismip6_retreat_k');
 add_default($nl, 'config_calving_error_threshold');
 add_default($nl, 'config_distribute_unablatedVolumeDynCell');
+add_default($nl, 'config_update_velocity_before_calving');
 
 ##################################
 # Namelist group: thermal_solver #
@@ -136,6 +139,8 @@ add_default($nl, 'config_dynamic_thickness');
 
 add_default($nl, 'config_dt');
 add_default($nl, 'config_time_integration');
+add_default($nl, 'config_rk_order');
+add_default($nl, 'config_rk3_stages');
 add_default($nl, 'config_adaptive_timestep');
 add_default($nl, 'config_min_adaptive_timestep');
 add_default($nl, 'config_max_adaptive_timestep');
@@ -197,6 +202,7 @@ add_default($nl, 'config_print_calving_info');
 add_default($nl, 'config_print_thermal_info');
 add_default($nl, 'config_always_compute_fem_grid');
 add_default($nl, 'config_print_velocity_cleanup_details');
+add_default($nl, 'config_check_tracer_monotonicity');
 
 ####################################
 # Namelist group: subglacial_hydro #

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

@@ -19,6 +19,8 @@
 <!-- advection -->
 <config_thickness_advection>'fo'</config_thickness_advection>
 <config_tracer_advection>'none'</config_tracer_advection>
+<config_horiz_tracer_adv_order>3</config_horiz_tracer_adv_order>
+<config_advection_coef_3rd_order>0.25</config_advection_coef_3rd_order>
 <config_restore_thickness_after_advection>.false.</config_restore_thickness_after_advection>
 <config_zero_sfcMassBalApplied_over_bare_land>.true.</config_zero_sfcMassBalApplied_over_bare_land>
 
@@ -62,6 +64,7 @@
 <config_ismip6_retreat_k>-170.0</config_ismip6_retreat_k>
 <config_calving_error_threshold>0.1</config_calving_error_threshold>
 <config_distribute_unablatedVolumeDynCell>.false.</config_distribute_unablatedVolumeDynCell>
+<config_update_velocity_before_calving>.false.</config_update_velocity_before_calving>
 
 <!-- thermal_solver -->
 <config_thermal_solver>'none'</config_thermal_solver>
@@ -115,6 +118,8 @@
 <!-- time_integration -->
 <config_dt>'0000-00-01_00:00:00'</config_dt>
 <config_time_integration>'forward_euler'</config_time_integration>
+<config_rk_order>2</config_rk_order>
+<config_rk3_stages>3</config_rk3_stages>
 <config_adaptive_timestep>.false.</config_adaptive_timestep>
 <config_min_adaptive_timestep>3600.0</config_min_adaptive_timestep>
 <config_max_adaptive_timestep>3.15e9</config_max_adaptive_timestep>
@@ -129,8 +134,8 @@
 <!-- time_management -->
 <config_do_restart>.false.</config_do_restart>
 <config_restart_timestamp_name>'rpointer.glc'</config_restart_timestamp_name>
-<config_start_time>'0000-01-01_00:00:00'</config_start_time>
-<config_stop_time>'0000-01-01_00:00:00'</config_stop_time>
+<config_start_time>'0001-01-01_00:00:00'</config_start_time>
+<config_stop_time>'0001-01-01_00:00:00'</config_stop_time>
 <config_run_duration>'none'</config_run_duration>
 <config_calendar_type CALENDAR="NO_LEAP">'noleap'</config_calendar_type>
 <config_calendar_type>'gregorian'</config_calendar_type>
@@ -159,6 +164,7 @@
 <config_print_thermal_info>.false.</config_print_thermal_info>
 <config_always_compute_fem_grid>.true.</config_always_compute_fem_grid>
 <config_print_velocity_cleanup_details>.false.</config_print_velocity_cleanup_details>
+<config_check_tracer_monotonicity>.false.</config_check_tracer_monotonicity>
 
 <!-- subglacial_hydro -->
 <config_SGH>.false.</config_SGH>

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

@@ -139,15 +139,31 @@ Default: Defined in namelist_defaults.xml
 	category="advection" group="advection">
 Selection of the method for advecting thickness ('fo' = first-order upwinding).
 
-Valid values: 'fo', 'none'
+Valid values: 'fo', 'fct', 'none'
 Default: Defined in namelist_defaults.xml
 </entry>
 
 <entry id="config_tracer_advection" type="char*1024"
 	category="advection" group="advection">
 Selection of the method for advecting tracers.
 
-Valid values: 'fo', 'none'
+Valid values: 'fo', 'fct', 'none'
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_horiz_tracer_adv_order" type="integer"
+	category="advection" group="advection">
+Order of polynomial used for tracer reconstruction at cell edges
+
+Valid values: 2, 3 and 4
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_advection_coef_3rd_order" type="real"
+	category="advection" group="advection">
+Reconstruction of 3rd-order reconstruction to blend with 4th-order reconstuction. Equivalent to beta in Skamarock and Gassmann (2011) eq. 7. 0 is fully 4th order, 1 is fully 3rd order.
+
+Valid values: any real between 0 and 1
 Default: Defined in namelist_defaults.xml
 </entry>
 
@@ -469,6 +485,14 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_update_velocity_before_calving" type="logical"
+	category="calving" group="calving">
+If true, add an additional velocity solve between advection and calving. If false, use velocity field from beginning of time step to calculate calving rate. For certain calving laws, like damage threshold calving, it is not necessary to update the velocity before calving, while for von Mises stress and eigencalving, it is more accurate to have an updated velocity state before solving for calvingThickness.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
 
 <!-- thermal_solver -->
 
@@ -835,9 +859,25 @@ Default: Defined in namelist_defaults.xml
 
 <entry id="config_time_integration" type="char*1024"
 	category="time_integration" group="time_integration">
-Time integration method (currently, only forward Euler is supported).
+Time integration method.
 
-Valid values: 'forward_euler'
+Valid values: 'forward_euler' or 'runge_kutta'
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_rk_order" type="integer"
+	category="time_integration" group="time_integration">
+Order of Runge-Kutta time integration to use. A value of 1 would be equivalent to forward euler, but will cause an error to avoid unnecessary redundancy. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods.
+
+Valid values: 2 or 3
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_rk3_stages" type="integer"
+	category="time_integration" group="time_integration">
+Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2.
+
+Valid values: 3 or 4
 Default: Defined in namelist_defaults.xml
 </entry>
 
@@ -1133,6 +1173,14 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_check_tracer_monotonicity" type="logical"
+	category="debug" group="debug">
+Check tracer monotonicity at the end of the monotonic advection routine and write warnings to log file if not monotonic.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
 
 <!-- subglacial_hydro -->
 

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

@@ -42,7 +42,7 @@
 		<dim name="nISMIP6OceanLayers" units="unitless"
 		     description="The number of layers in the ISMIP6 ocean temperature dataset."
 		/>
-		<dim name="maxTracersAdvect" definition="2" units="unitless"
+		<dim name="maxTracersAdvect" definition="4" units="unitless"
 		     description="The maximum number of tracers to be advected."
 		/>
 		<dim name="nRegions" definition="1" units="unitless"
@@ -110,12 +110,20 @@
 	<nml_record name="advection" in_defaults="true">
 		<nml_option name="config_thickness_advection" type="character" default_value="fo" units="unitless"
 		            description="Selection of the method for advecting thickness ('fo' = first-order upwinding)."
-		            possible_values="'fo', 'none'"
+		            possible_values="'fo', 'fct', 'none'"
 		/>
 		<nml_option name="config_tracer_advection" type="character" default_value="none" units="unitless"
 		            description="Selection of the method for advecting tracers."
-		            possible_values="'fo', 'none'"
-                />
+		            possible_values="'fo', 'fct', 'none'"
+		/>
+		<nml_option name="config_horiz_tracer_adv_order" type="integer" default_value="3" units="unitless"
+			    description="Order of polynomial used for tracer reconstruction at cell edges"
+			    possible_values="2, 3 and 4"
+		/>
+		<nml_option name="config_advection_coef_3rd_order" type="real" default_value="0.25" units="unitless"
+					description="Reconstruction of 3rd-order reconstruction to blend with 4th-order reconstuction. Equivalent to beta in Skamarock and Gassmann (2011) eq. 7. 0 is fully 4th order, 1 is fully 3rd order."
+					possible_values="any real between 0 and 1"
+		/>
 		<nml_option name="config_restore_thickness_after_advection" type="logical" default_value=".false." units="unitless"
 			    description="If true, reset thickness to values at previous timestep after advection occurs. This is used for spinning up tracer fields such as damage.  When this is true, geometry changes from surface and basal mass balance (grounded or floating) and facemelting are not retained, but changes from calving are."
 			    possible_values=".true. or .false."
@@ -284,6 +292,10 @@
                         description="If true, then distribute unablatedVolumeDynCell among dynamic neighbors when converting ablation velocity to ablation thickness. This should only be used as a clean-up measure, while limiting the timestep based on ablation velocity should be used as the primary method of getting accurate ablation thickness from ablation velocity.  If you choose to set config_adaptive_timestep_calvingCFL_fraction much larger than 1.0 (which is not recommended), setting this option to true usually results in more accurate calving behavior. "
                         possible_values=".true. or .false."
                 />
+                <nml_option name="config_update_velocity_before_calving" type="logical" default_value=".false." units="unitless"
+                            description="If true, add an additional velocity solve between advection and calving. If false, use velocity field from beginning of time step to calculate calving rate. For certain calving laws, like damage threshold calving, it is not necessary to update the velocity before calving, while for von Mises stress and eigencalving, it is more accurate to have an updated velocity state before solving for calvingThickness."
+                            possible_values=".true. or .false."
+                />
 	</nml_record>
 
 	<nml_record name="thermal_solver" in_defaults="true">
@@ -477,8 +489,16 @@
 		            possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS', but limited by CFL condition."
 		/>
 		<nml_option name="config_time_integration" type="character" default_value="forward_euler" units="unitless"
-		            description="Time integration method (currently, only forward Euler is supported)."
-		            possible_values="'forward_euler'"
+		            description="Time integration method."
+		            possible_values="'forward_euler' or 'runge_kutta'"
+		/>
+                <nml_option name="config_rk_order" type="integer" default_value="2" units="unitless"
+                            description="Order of Runge-Kutta time integration to use. A value of 1 would be equivalent to forward euler, but will cause an error to avoid unnecessary redundancy. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods."
+                            possible_values="2 or 3"
+                />
+		<nml_option name="config_rk3_stages" type="integer" default_value="3" units="stages"
+                        description="Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2."
+                        possible_values="3 or 4"
 		/>
 		<nml_option name="config_adaptive_timestep" type="logical" default_value=".false." units="unitless"
 			description="Determines if the time step should be adjusted based on the CFL condition or should be steady in time. If true, the config_dt_* options are ignored."
@@ -533,11 +553,11 @@
 		description="Path to the filename for restart timestamps to be read and written from."
 		possible_values="Path to a file."
 		/>
-		<nml_option name="config_start_time" type="character" default_value="0000-01-01_00:00:00" units="unitless"
+		<nml_option name="config_start_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
 		            description="Timestamp describing the initial time of the simulation.  If it is set to 'file', the initial time is read from the filename specified by config_restart_timestamp_name (defaults to 'restart_timestamp')."
 		            possible_values="'YYYY-MM-DD_HH:MM:SS' (items in the format string may be dropped from the left if not needed, and the components on either side of the underscore may be replaced with a single integer representing the rightmost unit)"
 		/>
-		<nml_option name="config_stop_time" type="character" default_value="0000-01-01_00:00:00" units="unitless"
+		<nml_option name="config_stop_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
 		            description="Timestamp describing the final time of the simulation. If it is set to 'none' the final time is determined from config_start_time and config_run_duration.  If config_run_duration is also specified, it takes precedence over config_stop_time.  Set config_stop_time to be equal to config_start_time (and config_run_duration to 'none') to perform a diagnostic solve only."
 		            possible_values="'YYYY-MM-DD_HH:MM:SS' or 'none' (items in the format string may be dropped from the left if not needed, and the components on either side of the underscore may be replaced with a single integer representing the rightmost unit)"
 		/>
@@ -638,6 +658,10 @@
 		            description="After velocity is calculated there are a few checks for appropriate values in certain geometric configurations.  Setting this option to .true. will cause detailed information about those adjustments to be printed."
 		            possible_values=".true. or .false."
 		/>
+		<nml_option name="config_check_tracer_monotonicity" type="logical" default_value=".false."
+                            description="Check tracer monotonicity at the end of the monotonic advection routine and write warnings to log file if not monotonic."
+                            possible_values=".true. or .false."
+		/>
 	</nml_record>
 
 
@@ -741,6 +765,7 @@
                         <var name="uReconstructX" packages="higherOrderVelocity"/>
 			<var name="uReconstructY" packages="higherOrderVelocity"/>
                         <var name="basalFrictionFlux"/>
+                        <var name="passiveTracer2d"/>
 			<var name="damage"/>
 			<var name="calvingVelocityData"/>
 			<var name="basalMeltInput" packages="hydro"/>
@@ -850,6 +875,7 @@
 			<var name="waterPressure" packages="hydro"/>
                         <var name="channelArea" packages="hydro"/>  <!-- this is only needed if running with channels - could be package-controlled -->
 			<var name="waterFluxMask" packages="hydro"/>
+                        <var name="passiveTracer2d"/>
 			<var name="damage"/>
 			<var name="calvingVelocityData"/>
 			<var name="damageMax"/>
@@ -1207,9 +1233,12 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="Applied basal mass balance on floating regions"
                 />
-                <var name="calvingMask" type="integer" dimensions="nCells Time" units="m" time_levs="1"
+                <var name="calvingMask" type="integer" dimensions="nCells Time" units="" time_levs="1"
                         description="mask of grid cells that should be calved.  0=no calving, 1=should be calved"
                 />
+                <var name="openOceanMask" type="integer" dimensions="nCells Time" units="" time_levs="1"
+                        description="mask of open ocean cells, including floating non-dynamic cells.  Used to eliminate holes in ice shelf from calving"
+                />
                 <var name="calvingFrontMask" type="integer" dimensions="nCells Time" units="none" time_levs="1"
                         description="mask of the calving front position on cells"
                 />
@@ -1270,6 +1299,9 @@ is the value of that variable from the *previous* time level!
                 <var name="calvingVelocityData" type="real" dimensions="nCells Time" units="m s^{-1}" time_levs="1"
                      description="rate of calving front retreat due to calving, represented as a velocity normal to the calving front (in the x-y plane), given by the input netCDF file."
                 />
+                <var name="passiveTracer2d" type="real" dimensions="nCells Time" units="none" time_levs="1"
+                     description="passive tracer used to verify advection routines"
+                />
                 <var name="damage" type="real" dimensions="nCells Time" units="none" time_levs="1"
                      description="Damage is parameterized as the local ice thickness divided by the fracture depth, such that unfractured ice has damage=0 and ice fractured through its full thickness has damage=1"
                 />
@@ -1413,6 +1445,9 @@ is the value of that variable from the *previous* time level!
                 <var name="layerNormalVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
                      description="horizonal velocity, normal component to an edge, layer midpoint"
                 />
+		<var name="layerThicknessEdgeFlux" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-1}"
+                     description="layer-normal thickness flux on edges"
+		/>
                 <var name="normalVelocityInitial" type="real" dimensions="nVertInterfaces nEdges Time" units="m s^{-1}"
                      description="horizonal velocity, normal component to an edge, computed at initialization"
                 />