Add flux-corrected transport scheme to MALI

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."
@@ -638,6 +646,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 +753,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 +863,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"/>
@@ -1270,6 +1284,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"
                 />

components/mpas-albany-landice/src/mode_forward/Makefile

@@ -7,6 +7,8 @@ OBJS =  mpas_li_core.o \
     mpas_li_time_integration_fe.o \
     mpas_li_diagnostic_vars.o \
     mpas_li_advection.o \
+    mpas_li_advection_fct_shared.o \
+    mpas_li_advection_fct.o \
     mpas_li_calving.o \
     mpas_li_statistics.o \
     mpas_li_velocity.o \
@@ -45,8 +47,12 @@ mpas_li_time_integration_fe.o: mpas_li_advection.o \
 			       mpas_li_bedtopo.o
 
 mpas_li_advection.o: mpas_li_thermal.o \
+                     mpas_li_advection_fct.o \
+                     mpas_li_advection_fct_shared.o \
 	             mpas_li_diagnostic_vars.o
 
+mpas_li_advection_fct.o:  mpas_li_advection_fct_shared.o
+
 mpas_li_calving.o: mpas_li_thermal.o \
 		   mpas_li_advection.o