Merge branch 'trhille/add_higher_order_advection' into MALI-Dev/develop

This merge adds a flux-correct transport (FCT) scheme to MALI for thickness and
tracer advection, ported over with MALI-relevant modification from MPAS-Ocean's
routines, which are based on Skamarock and Gassmann (2011)
(https://doi.org/10.1175/MWR-D-10-05056.1). This uses a blend of 3rd and 4th
order fluxes to achieve monotonicity. The FCT routine is only used for tracers
in MPAS-Ocean, whereas here it is modified for use with both thickness and
tracers. The user can specify 2nd, 3rd, or 4th order advection with the
config_horiz_tracer_adv_order option, but only config_horiz_tracer_adv_order = 3
with 0 < config_advection_coef_3rd_order < 1 is truly FCT. The
config_advection_coef_3rd_order option specifies the blend between 3rd and 4th
order fluxes used in the flux correction. config_advection_coef_3rd_order = 1.0
is purely 3rd order, while config_advection_coef_3rd_order = 0.0 is purely 4th
order. The default value of 0.25 is taken from the MPAS-Ocean default and may
not be appropriate for all situations. Note that all higher-order advection must
reduce to 1st order at the boundaries.

This also adds a new variable passiveTracer2d, that can be used to verify
advection schemes.

Currently supported combinations of thickness and tracer advection with fct
include:

1. config_thickness_advection = 'fo'; config_tracer_advection = 'fct'
2. config_thickness_advection = 'fct'; config_tracer_advection = 'fct'
3. config_thickness_advection = 'fct'; config_tracer_advection = 'none'

FCT tracer advection with no thickness advection and FCT thickness advection
with FO tracer advection could be added, but we do not currently anticipate
using them.  Therefore, we have left them out of this PR as they would add
unnecessary complexity to the code.  When used with forward Euler time
integration, FCT requires a severely reduced time step relative to first order
advection. Testing shows that CFL fraction on the order of 0.1 is likely
sufficiently small, but this is probably case-dependent. Once Runge-Kutta time
integration is operational, it is recommended to use that instead of forward
Euler.

At the time of merging, there is a very slight conservation error (<0.1% in the
grounded slab test shown below) for tracers when using FCT for both thickness
and tracer advection, while mass is fully conserved. Tracers are conserved when
using FO thickness advection with FCT tracer advection. The convergence order
has not yet been verified, but testing is under way using a new mesh convergence
COMPASS test: MPAS-Dev/compass#690.

* trhille/add_higher_order_advection: (46 commits)
  Enable fct thickness advection without tracer advection
  Change passiveTracer to passiveTracer2d
  Only pass layer thickness tracer for first call to fct
  Fix the case of fo thickness, none tracer advection
  Fix bug with activeTracerHorizontalAdvectionEdgeFlux
  Make activeTracerHorizontalAdvectionEdgeFlux optional
  Throw error forr fct thickness with fo tracer
  Simple cleanup after code review
  Make conserve tracer volume
  Fix first order flux at ice edge
  Make fct conserve mass
  Call li_tracer_advection_fct_tend for thickness and tracers separately
  Clean up 2nd order mask
  Make 2nd order and 3rd-4th order masks mutually exclusive
  Try new mask for 2nd order terms
  Change normalThicknessFlux, layerThickness, and tracer definitions
  Pass layerThickness as tracer instead of array of 1s
  Increase max number of tracers to accomodate passiveTracer
  Add passiveTracer to help verify advection schemes
  Fix bug in marking boundaryCell
  ...

Author: matthewhoffman

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"/>
@@ -1273,6 +1287,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