diff --git a/docs/images/ALE_general_schematic.png b/docs/images/ALE_general_schematic.png
deleted file mode 100644
index 3f492ed56d..0000000000
Binary files a/docs/images/ALE_general_schematic.png and /dev/null differ
diff --git a/docs/ocean.bib b/docs/ocean.bib
index a5d7483f9c..2297f25354 100644
--- a/docs/ocean.bib
+++ b/docs/ocean.bib
@@ -70,63 +70,3 @@ @article{Kasahara1974
title = {Various Vertical Coordinate Systems Used for Numerical Weather Prediction},
journal = {Monthly Weather Rev.}
}
-
-@Article{Griffies_Adcroft_Hallberg2020,
-author = "S.M. Griffies and A. Adcroft and R.W. Hallberg",
-title = "A primer on the vertical Lagrangian-remap method in
- ocean models based on finite volume generalized vertical coordinates",
-journal = "Journal of Advances in Modeling Earth Systems",
-year = "2020",
-volume = "12",
-doi = "10.1029/2019MS001954",
-}
-
-@Article{Shao_etal_2020,
-author = "A. Shao and A.J. Adcroft and R.W. Hallberg and S.M. Griffies",
-title = "A general-coordinate, nonlocal neutral diffusion operator",
-journal = "Journal of Advances in Modeling Earth Systems",
-year = "2020",
-volume = "12",
-doi = "10.1029/2019MS001992",
-}
-
-@Article{GM95,
-author = "P. R. Gent and J. Willebrand and T. J. McDougall and J. C. McWilliams",
-title = "Parameterizing eddy-induced tracer transports in ocean circulation models",
-journal = "Journal of Physical Oceanography",
-year = "1995",
-volume = "25",
-pages = "463--474",
-doi = "10.1175/1520-0485(1995)025<0463:PEITTI>2.0.CO;2",
-}
-
-@Article{foxkemper_etal2008,
-author = "Baylor Fox-Kemper and Raffaele Ferrari and Robert Hallberg",
-title = "Parameterization of mixed layer eddies. {I}: {T}heory and diagnosis",
-journal = "Journal of Physical Oceanography",
-year = "2008",
-volume = "38",
-pages = "1145--1165",
-doi = "10.1175/2007JPO3792.1",
-}
-
-@Article{McDougall_etal_2021,
-author = "T. J. McDougall and P.M.\ Barker and R.M.\ Holmes and R.\ Pawlowicz and S.M.\ Grif\/f\/ies and P.J.\ Durack",
-title = "The interpretation of temperature and salinity variables in numerical ocean model output,
- and the calculation of heat fluxes and heat content",
-journal = "Geoscientific Model Development",
-year = "2021",
-volume = "14",
-pages = "6445--6466",
-doi = "10.5194/gmd-14-6445-2021",
-}
-
-@article{Young2010,
-author = "W. R. Young",
-year = "2010",
-title = "Dynamic Enthalpy, {Conservative Temperature}, and the Seawater {Boussinesq} Approximation",
-journal = "Journal of Physical Oceanography",
-volume = "40",
-pages = "394--400",
-doi = "10.1175/2009JPO4294.1",
-}
diff --git a/src/ALE/_ALE.dox b/src/ALE/_ALE.dox
index 8239a79a0e..9313ed2aa1 100644
--- a/src/ALE/_ALE.dox
+++ b/src/ALE/_ALE.dox
@@ -1,184 +1,90 @@
-/*! \page ALE Vertical Lagrangian method: conceptual
+/*! \page ALE ALE
-\section section_ALE Lagrangian and ALE
+\section section_ALE Basics of the Vertical Lagrangian-Remap Method in MOM6
-As discussed by Adcroft and Hallberg (2008) \cite adcroft2006 and
-Griffies, Adcroft and Hallberg (2020) \cite Griffies_Adcroft_Hallberg2020,
-we can conceive of two general classes
+As discussed by \cite adcroft2006, there are two general classes
of algorithms that frame how hydrostatic ocean models are
formulated. The two classes differ in how they treat the vertical
direction. Quasi-Eulerian methods follow the approach traditionally
-used in geopotential coordinate models, whereby vertical motion is
-diagnosed via the continuity equation. Quasi-Lagrangian methods are
-traditionally used by layered isopycnal models, with the vertical
-Lagrangian approach specifying motion that crosses coordinate
+used in geopotential coordinate models, whereby vertical motion
+is diagnosed via the continuity equation. Quasi-Lagrangian
+methods are traditionally used by layered isopycnal models, with
+the Lagrangian approach specifying motion that crosses coordinate
surfaces. Indeed, such dia-surface flow can be set to zero using
-Lagrangian methods for studies of adiabatic dynamics. MOM6 makes use
-of the vertical Lagrangian remap method, as pioneered for ocean
-modeling by Bleck (2002) \cite bleck2002 and further documented by
-\cite Griffies_Adcroft_Hallberg2020, with this method a limit case of
-the Arbitrary-Lagrangian-Eulerian method (\cite hirt1997). Dia-surface
+Lagrangian methods for studies of adiabatic dynamics. MOM6 makes
+use of the vertical Lagrangian remap method, as pioneered for
+ocean modeling by \cite bleck2002, which is a limit case of the
+Arbitrary-Lagrangian-Eulerian method (\cite hirt1997). Dia-surface
transport is implemented via a remapping so that the method can be
-summarized as the Lagrangian plus remap approach and so it is a
-one-dimensional version of the incremental remapping of Dukowicz (2000)
+summarized as the Lagrangian plus remap approach and is essentially
+a one-dimensional version of the incremental remapping of
\cite dukowicz2000.
-\image html ALE_general_schematic.png "Schematic of the 3d Lagrangian regrid/remap method" width=70%
-\image latex ALE_general_schematic.png "Schematic of the 3d Lagrangian regrid/remap method" width=0.7\textwidth
+The MOM6 implementation of the vertical Lagrangian-remap method makes use
+of two general steps. The first evolves the ocean state forward in
+time according to a vertical Lagrangian limit with \f$\dot{r}=0\f$. Hence,
+the horizontal momentum, thickness, and tracers are time stepped
+with the red terms removed in equations \eqref{eq:h-horz-momentum,h-equations,momentum},
+\eqref{eq:h-thickness-equation,h-equations,thickness}, \eqref{eq:h-temperature-equation,h-equations,potential temperature},
+and \eqref{eq:h-salinity-equation,h-equations,salinity}. All advective transport thus
+occurs within a layer as defined by constant \f$r\f$-surfaces so that
+the volume within each layer is fixed. All other terms are retained in
+their full form, including subgrid scale terms that contribute to
+the transfer of tracer and momentum into distinct \f$r\f$ layers (e.g.,
+dia-surface diffusion of tracer and velocity). Maintaining constant
+volume within a layer yet allowing for tracers to move between layers
+engenders no inconsistency between tracer and thickness evolution. The
+reason is that tracer diffusion, even dia-surface diffusion, does
+not transfer volume.
-Refer to the above figure taken from Griffies, Adcroft, and Hallberg
-(2020) \cite Griffies_Adcroft_Hallberg2020. It shows a schematic of
-the Lagrangian-remap method as well as the Arbitrary
-Lagrangian-Eulerian (ALE) method. The first panel shows a square fluid
-region and square grid used to represent the fluid, along with
-rectangular subregions partitioned by grid lines. The second panel
-shows the result of evolving the fluid region and evolving the
-grid. The grid can evolve according to the fluid flow, as per a
-Lagrangian method, or it can evolve according to some specified grid
-evolution, as per an ALE method. The right panel depicts the grid
-reinitialization onto a target grid (the regrid step). A regrid step
-necessitates a corresponding remap step to estimate the ocean state on
-the target grid, with conservative remapping required to preserve
-integrated scalar contents (e.g., potential enthalpy, salt mass, and
-seawater mass). The regrid/remap steps are needed for Lagrangian
-methods in order for the grid to retain an accurate representation of
-the ocean state. Ideally, the remap step does not affect any changes
-to the fluid state; rather, it only modifies where in space the fluid
-state is represented. However, any numerical realization incurs
-interpolation inaccuracies that lead to unphysical (spurious) state
-changes.
-
-\section section_ALE_MOM Vertical Lagrangian regrid/remap method
-
-We now get a bit more specific to the vertical Lagrangian method.
-For this purpose, recall recall the basic dynamical equations (those
-equations with a time derivative) of MOM6 discussed in
-\ref General_Coordinate
-\f{align}
-\rho_0
-\left[ \frac{\partial \mathbf{u}}{\partial t} + \frac{( f + \zeta )}{h} \,
-\hat{\mathbf{z}} \times h \, \mathbf{u} + \underbrace{ \dot{r} \,
-\frac{\partial \mathbf{u}}{\partial r} }
-\right]
-&= -\nabla_r \, (p + \rho_{0} \, K) -
-\rho \nabla_r \, \Phi + \mathbf{\mathcal{F}}
-&\mbox{horizontal momentum}
-\label{eq:h-horz-momentum-vlm}
-\\
-\frac{\partial h}{\partial t} + \nabla_r \cdot \left( h \, \mathbf{u} \right) +
-\underbrace{ \delta_r ( z_r \dot{r} ) }
- &= 0
-&\mbox{thickness}
-\label{eq:h-thickness-equation-vlm}
-\\
-\frac{\partial ( \theta \, h )}{\partial t} + \nabla_r \cdot \left( \theta h \,
-\mathbf{u} \right) + \underbrace{ \delta_r ( \theta \, z_r \dot{r} ) }
-&=
-h \mathbf{\mathcal{N}}_\theta^\gamma - \delta_r J_\theta^{(z)}
-&\mbox{potential/Conservative temp}
-\label{eq:h-temperature-equation-vlm}
- \\
-\frac{\partial ( S \, h )}{\partial t} + \nabla_r \cdot \left( S \, h \,
-\mathbf{u} \right) + \underbrace{ \delta_r ( S \, z_r \dot{r} ) }
- &=
-h \mathbf{\mathcal{N}}_S^\gamma - \delta_r J_S^{(z)}
-&\mbox{salinity}
-\label{eq:h-salinity-equation-vlm}
-\f}
-The MOM6 implementation of the vertical Lagrangian method makes
-use of two general steps. The first evolves the ocean state forward in
-time according to a vertical Lagrangian approach with with
-\f$\dot{r}=0\f$. Hence, the horizontal momentum, thickness, and
-tracers are time stepped with the underbraced terms removed in the
-above equations. All advective transport occurs within a layer as
-defined by constant \f$r\f$-surfaces so that the volume within each
-layer is fixed. All other terms are retained in their full form,
-including subgrid scale terms that contribute to the transfer of
-tracer and momentum into distinct \f$r\f$ layers (e.g., dia-surface
-diffusion of tracer and velocity). Maintaining constant volume within
-a layer yet allowing for tracers to move between layers engenders no
-inconsistency between tracer and thickness evolution. The reason is
-that tracer diffusion, even dia-surface diffusion, does not transfer
-volume.
-
-The second step in the method comprises the generation of a new
+The second step in the algorithm comprises the generation of a new
vertical grid following a prescription, such as whether the grid
-should align with isopcynals or constant \f$z^{*}\f$ or a combination.
-This second step is known as the regrid step. The ocean state is then
-vertically remapped to the newly generated vertical grid. This
-remapping step incorporates dia-surface transfer of properties, with
-such transfer depending on the prescription given for the vertical
+should align with isopcynals or constant \f$z^{*}\f$ or a combination. The
+ocean state is then vertically remapped to the newly generated vertical
+grid. The remapping step incorporates dia-surface transfer of properties,
+with such transfer depending on the prescription given for the vertical
grid generation. To minimize discretization errors and the associated
-spurious mixing, the remapping step makes use of the high order
-accurate methods developed by \cite white2008 and \cite white2009.
+spurious mixing, the remapping step makes use of the high order accurate
+methods developed by \cite white2008 and \cite white2009.
+The underlying algorithm for treatment of the vertical can
+be related to operator-splitting of the red terms in equations
+\eqref{eq:h-thickness-equation,h-equations,thickness}--\eqref{eq:h-temperature-equation,h-equations,potential temperature}. If we
+consider, for simplicity, an Euler-forward update for a time-step \f$\Delta
+t\f$, the time-stepping for the continuity and temperature equation can
+be summarized as
-\section section_ALE_MOM_numerics Outlining the numerical algorithm
-
-The underlying algorithm for treatment of the vertical can be related
-to operator-splitting of the underbraced terms in the above equations.
-If we consider, for simplicity, an Euler-forward update for a
-time-step \f$\Delta t\f$, the time-stepping for the thickness and
-tracer equation (\f$C\f$ is an arbitrary tracer) can be summarized as
-(from Table 1 in Griffies, Adcroft and Hallberg (2020)
-\cite Griffies_Adcroft_Hallberg2020)
-\f{align}
-\label{html:ale-equations}\notag
-\\
- \delta_{r} w^{\scriptstyle{\mathrm{grid}}}
- &= -\nabla_{r} \cdot [h \, \mathbf{u}]^{(n)}
- &\mbox{layer motion via convergence of horz advection}
-\\
- h^{\dagger} &= h^{(n)} + \Delta t \, \delta_{r} w^{\scriptstyle{\mathrm{grid}}}
-= h^{(n)} - \Delta t \, \nabla_{r} \cdot [h \, \mathbf{u}]^{(n)}
- &\mbox{horz advection thickness update}
-\\
- [h \, C]^{\dagger} &= [h \, C]^{(n)} -\Delta t \, \nabla_{r} \cdot [ h \, C \, \mathbf{u} ]^{(n)}
- &\mbox{horz advection tracer update}
-\\
- h^{(n+1)} &= h^{\scriptstyle{\mathrm{target}}}
- &\mbox{regrid to the target grid}
-\\
- \delta_{r} w^{(\dot{r})} &= -(h^{\scriptstyle{\mathrm{target}}} - h^{\dagger})/\Delta t
- &\mbox{diagnose dia-surface transport}
-\\
- [h \, C]^{(n+1)} &= [h \, C]^{\dagger} - \Delta t \, \delta_{r} ( w^{(\dot{r})} \, C^{\dagger})
- &\mbox{remap tracer using dia-surface transport}
-\f}
-The first three equations constitute the Lagrangian portion of the
-algorithm. In particular, the second equation provides an
-intermediate or ``predictor'' value for the updated thickness,
-\f$h^{\dagger}\f$, resulting from the vertical Lagrangian update.
-Similarly, the third equation performs a Lagrangian update of the
-thickness-weighted tracer to intermediate values, again operationally
-realized by dropping the \f$w^{(\dot{r})}\f$ contribution.
-The fourth equation is the regrid step, which is the key step in the
-algorithm with the new grid defined by the new thickness
-\f$h^{(n+1)}\f$. The new thickness is prescribed by the target values
-for the vertical grid,
-\f{align}
- h^{(n+1)} = h^{\scriptstyle{\mathrm{target}}}.
-\f}
-The prescribed target grid thicknesses are then used to diagnose the
-dia-surface velocity according to
-\f{align}
- \delta_{r} w^{(\dot{r})} = -(h^{\scriptstyle{\mathrm{target}}} - h^{\dagger})/\Delta t.
+\f{eqnarray}
+\label{html:ale-equations}\notag \\
+h^\dagger &= h^{(n)} - \Delta t \left[ \nabla_r \cdot \left( h \, \mathbf{u} \right) \right]
+&\mbox{thickness} \label{eq:ale-thickness-equation} \\
+\theta^\dagger \, h^\dagger &= \theta^{(n)} \, h^{(n)} - \Delta t \left[ \nabla_r \cdot \left( \theta h \, \mathbf{u} \right) - h \boldsymbol{\mathcal{N}}_\theta^\gamma + \delta_r J_\theta^{(z)} \right]
+&\;\;\;\;\mbox{potential temp} \label{eq:ale-temperature-equation} \\
+h^{(n+1)} &= h^\dagger - \Delta t \, \delta_r \left( z_r \dot{r} \right)
+&\mbox{move grid} \label{eq:ale-new-grid} \\
+\theta^{(n+1)} h^{(n+1)} &= \theta^\dagger h^\dagger - \Delta t \, \delta_r \left( z_r \dot{r} \, \theta^\dagger \right)
+&\mbox{remap temperature.} \label{eq:ale-remap-temperature}
\f}
-This step, and the remaining step for tracers, constitute the
-remapping portion of the algorithm. For example, if the prescribed
-coordinate surfaces are geopotentials, then \f$w^{(\dot{r})} = w\f$
-and \f$h^{\scriptstyle{\mathrm{target}}} = h^{(n)}\f$, in which case the
-remap step reduces to Cartesian vertical advection.
+
+Substituting \eqref{eq:ale-thickness-equation,ale-equations,thickness} into \eqref{eq:ale-new-grid,ale-equations,move grid}
+recovers a time-discrete form of \eqref{eq:h-thickness-equation,h-equations,thickness}. The
+intermediate quantities indicated by \f$^\dagger\f$-symbols are the result of
+the vertical Lagrangian step of the algorithm. What were the red terms in
+the continuous-in-time equations are used to evolve the the intermediate
+quantities to the final updated quantities each step. In MOM6, equation
+\eqref{eq:ale-new-grid,ale-equations,move grid} is essentially used to define the dia-surface
+transport \f$z_r \dot{r}\f$ by prescribing \f$h^{(n+1)}\f$. For example, to
+recover a z-coordinate model, \f$h^{(n+1)}=\Delta z\f$, and \f$z_r \dot{r}\f$
+becomes the Eulerian vertical velocity, \f$w\f$.
Within the above framework for evolving the ocean state, we make use
of a standard split-explicit time stepping method by decomposing the
horizontal momentum equation into its fast (depth integrated) and slow
-(deviation from depth integrated) components. Furthermore, we follow
-the methods of Hallberg and Adcroft (2009) \cite hallberg2009 to
-ensure that the free surface resulting from time stepping the depth
-integrated thickness equation (i.e., the free surface equation) is
-consistent with the sum of the thicknesses that result from time
-stepping the layer thickness equations for each of the discretized
-layers; i.e., \f$\sum_{k} h = H + \eta\f$.
+(deviation from depth integrated) components. Furthermore, we follow the
+methods of \cite hallberg2009 to ensure that the free surface resulting
+from time stepping the depth integrated thickness equation (i.e., the
+free surface equation) is consistent with the sum of the thicknesses
+that result from time stepping the layer thickness equations for each
+of the discretized layers; i.e., \f$\sum_{k} h = H + \eta\f$.
*/
diff --git a/src/ALE/_ALE_timestep.dox b/src/ALE/_ALE_timestep.dox
index 6aa1ff8f00..e6da55fda9 100644
--- a/src/ALE/_ALE_timestep.dox
+++ b/src/ALE/_ALE_timestep.dox
@@ -1,63 +1,50 @@
-/*! \page ALE_Timestep Vertical Lagrangian method in pictures
-
-\section section_ALE_remap Graphical explanation of vertical Lagrangian method
-
-Vertical Lagrangian regridding/remapping is not a timestep method in
-the traditional sense. Rather, it is a sequence of operations
-performed to bring the vertical grid back to a target specification
-(the regrid step), and then to remap the ocean state onto this new
-grid (the remap step). This regrid/remap process can be chosen to be
-less frequent than the momentum or thermodynamic timesteps. We are
-motivated to choose less frequent regrid/remap steps to save
-computational time and to reduce spurious mixing that occurs due to
-truncation errors in the remap step. However, there is a downside to
-delaying the regrid/remap. Namely, if delayed too long then the layer
-interfaces can become entangled (i.e., no longer monotonic in the
-vertical), which is a common problem with purely Lagrangian methods.
-On this page we illustrate the regrid/remap steps by making use of
-Figure 3 from Griffies, Adcroft, and Hallberg (2020)
-\cite Griffies_Adcroft_Hallberg2020.
-
-For purposes of this example, assume that the target vertical grid is
-comprised of geopotential \f$z\f$-surfaces, with the initial ocean
-state (e.g., the temperature field) shown on the left in the following
-figure.
-
-\image html remapping1.png "Initial state with level surface (left) and perturbed state after a wave has come through (right)" width=60%
-\image latex remapping1.png "Initial state with level surface (left) and perturbed state after a wave has come through (right)" width=0.6\textwidth
-
-Some time later, assume a wave has perturbed the ocean state. During
-the Lagrangian portion of the algorithm, the coordinate surfaces move
-vertically with the ocean fluid according to \f$\dot{r}=0\f$. Assume
-now that the algorithm has determined that a regrid step is needed,
-with the target vertical grid still geopotential \f$z\f$-surfaces, so
-this new target grid is shown overlaid on the left as a regrid.
-
-\image html remapping2.png "The regrid operation (left) and the remap operation (right)" width=60%
-\image latex remapping2.png "The regrid operation (left) and the remap operation (right)" width=0.6\textwidth
-
-The most complex part of the method involves remapping the wavy ocean
-field onto the new grid. This step also incurs truncation errors that
-are a function of the vertical grid spacing and the numerical method
-used to perform the remapping. We illustrate this remap step in the
-figure above, as well as in the frame below shown after the old
-deformed coordinate grid has been deleted:
-
-\image html remapping3.png "The final state after regriddinig and remapping" width=30%
-\image latex remapping3.png "The final state after regridding and remapping" width=0.3\textwidth
-
-The new layer thicknesses, \f$h_k\f$, are computed and then the layers
-are populated with the new velocities and tracers
-\f{align}
- % h_k^{\mbox{new}} %= \nabla_k z_{\mbox{coord}} \\
- \sum h_k^{\mbox{new}} &= \sum h_k^{\mbox{old}}
-\\
- \mathbf{u}_k^{\mbox{new}}
- &= \frac{1}{h_k}
- \int_{z_{k + \frac{1}{2}}}^{z_{k + \frac{1}{2}} + h_k} \mathbf{u}^{\mbox{old}}(z') \, \mathrm{d}z'
-\\
- \theta^{\mbox{new}} &= \frac{1}{h_k}
- \int_{z_{k + \frac{1}{2}}}^{z_{k + \frac{1}{2}} + h_k} \theta^{\mbox{old}}(z') \, \mathrm{d}z'
-\f}
+/*! \page ALE_Timestep ALE Timestep
+
+\section section_ALE_remap Explanation of ALE remapping
+
+The Arbitrary Lagrangian-Eulerian (ALE) remapping is not a timestep in the traditional
+sense, but rather an operation performed to bring the vertical coordinate back to the target
+specification. This remapping can be less frequent than the momentum or
+thermodynamic timesteps, but must be done before the layer interfaces become entangled
+with each other.
+
+Assuming the target vertical grid is level \f$z\f$-surfaces, the initial state is
+shown on the left in the following figure:
+
+\image html remapping1.png "The initial state with level surface (left) and the perturbed state after a wave has come through (right)."
+\image latex remapping1.png "The initial state with level surface (left) and the perturbed state after a wave has come through (right)."
+
+Some time later, a wave has perturbed the surfaces which move with the
+fluid and it has been determined that a remapping operation is needed. The
+target vertical grid is still level \f$z\f$-surfaces, so this new target
+grid is shown overlaid on the left as regrid:
+
+\image html remapping2.png "The regrid operation (left) and the remap operation (right)."
+\image latex remapping2.png "The regrid operation (left) and the remap operation (right)."
+
+The complex part of the operation is remapping the wavy field onto the new grid as
+shown on the right and again in the final frame after the old deformed coordinate
+system has been deleted:
+
+\image html remapping3.png "The final state after remapping."
+\image latex remapping3.png "The final state after remapping."
+
+Mathematically, the new layer thicknesses, \f$h_k\f$, are computed and then populated
+with the new velocities and tracers:
+
+\f[
+ h_k^{\mbox{new}} = \nabla_k z_{\mbox{coord}}
+\f]
+\f[
+ \sum h_k^{\mbox{new}} = \sum h_k^{\mbox{old}}
+\f]
+\f[
+ \vec{u}_k^{\mbox{new}} = \frac{1}{h_k} \int_{z_{k + \frac{1}{2}}}^{z_{k +
+ \frac{1}{2}} + h_k} \vec{u}^{\mbox{old}}(z')dz'
+\f]
+\f[
+ \theta^{\mbox{new}} = \frac{1}{h_k} \int_{z_{k + \frac{1}{2}}}^{z_{k +
+ \frac{1}{2}} + h_k} \theta^{\mbox{old}}(z')dz'
+\f]
*/
diff --git a/src/core/_General_coordinate.dox b/src/core/_General_coordinate.dox
index d8d7d393cd..cdaf8a34ea 100644
--- a/src/core/_General_coordinate.dox
+++ b/src/core/_General_coordinate.dox
@@ -1,158 +1,76 @@
-/*! \page General_Coordinate Generalized vertical coordinate equations
+/*! \page General_Coordinate General coordinate equations
-The ocean equations discretized by MOM6 are formulated using
-generalized vertical coordinates. Motivation for using generalized
-vertical coordinates, and a full accounting of the ocean equations
-written using these coordinates, can be found in Griffies, Adcroft and
-Hallberg (2020) \cite Griffies_Adcroft_Hallberg2020. Here we provide
-a brief summary.
+Transforming to a vertical coordinate \f$r(z,x,y,t)\f$, with \f$\dot{r} = \frac{\partial r}{\partial t}\f$ ...
-Consider a smooth function of space and time, \f$r(x,y,z,t)\f$, that
-has a single-signed and non-zero vertical derivative known as the
-specific thickness
-\f{align}
- \partial z/\partial r = (\partial r/\partial z)^{-1} = \mbox{specific thickness.}
-\f}
-The specific thickness measures the inverse vertical stratification of
-the vertical coordinate surfaces. As so constrained, \f$r\f$ can
-uniquely prescribe a positiion in the vertical. Consequently, the
-ocean equations can be mapped one-to-one from geopotential vertical
-coordinates to generalized vertical coordinate. Upon transforming to
-\f$r\f$-coordinates, the material time derivative of \f$r\f$ appears
-throughout the equations, playing the role of a pseudo-vertical
-velocity, and we make use of the following shorthand for this
-derivative
-\f{align}
-\dot{r} = D_{t} r.
-\f}
+The Boussinesq hydrostatic equations of motion in general-coordinate
+\f$r\f$ are:
-The Boussinesq hydrostatic ocean equations take the following form using
-generalized vertical coordinates (\f$r\f$-coordinates)
-\f{align}
+\f{eqnarray}
\label{html:r-equations}\notag \\
-\rho_o \left[
- \partial_{t} \mathbf{u} + (f + \zeta) \, \hat{\mathbf{z}} \times \mathbf{u}
- + \dot{r} \, \partial_{r} \mathbf{u} \right]
- &= -\nabla_r \, (p + \rho_{o} \, K) -\rho \nabla_r \Phi + \rho_{o} \, \mathbf{\mathcal{F}}
- &\mbox{horizontal momentum}
-\label{eq:r-horz-momentum}
-\\
-\rho \, \partial_{r} \Phi + \partial_{r}p
- &= 0
-&\mbox{hydrostatic}
-\label{eq:r-hydrostatic-equation}
-\\
- \partial_{t}( z_r)
-+ \nabla_r \cdot ( z_r \, \mathbf{u} )
-+ \partial_{r} ( z_r \, \dot{r} )
-&= 0
-&\mbox{specific thickness}
-\label{eq:r-non-divergence}
-\\
- \partial_{t} ( \theta \, z_r )
-+ \nabla_r \cdot ( \theta z_r \, \mathbf{u} )
-+ \partial_{r} ( \theta \, z_r \, \dot{r} )
-&= z_r \mathbf{\mathcal{N}}_\theta^\gamma
-- \partial_{r} J_\theta^{(z)}
-&\mbox{potential/Conservative temp}
-\label{eq:r-temperature-equation}
-\\
-\partial_{t} ( S \, z_r)
-+ \nabla_r \cdot ( S \, z_r \, \mathbf{u} )
-+ \partial_{r} ( S \, z_r \, \dot{r} )
-&= z_r \mathbf{\mathcal{N}}_S^\gamma
-- \partial_{r} J_S^{(z)}
-&\mbox{salinity}
-\label{eq:r-salinity-equation}
-\\
-\rho &= \rho( S, \theta, -g \rho_0 z )
+\rho_0 \left( \frac{\partial \mathbf{u}}{\partial t} + ( f + \zeta ) \, \hat{\mathbf{z}} \times \mathbf{u} + \dot{r} \, \frac{\partial \mathbf{u}}{\partial r} + \nabla_r \, K \right) &= -\nabla_r \, p - \rho \nabla_r \, \Phi + \boldsymbol{\mathcal{F}}
+&\mbox{momentum} \label{eq:r-horz-momentum} \\
+\rho \, \frac{\partial \Phi}{\partial r} + \frac{\partial p}{\partial r} &= 0
+&\mbox{hydrostatic} \label{eq:r-hydrostatic-equation} \\
+\frac{\partial z_r }{\partial t} + \nabla_r \cdot \, \left( z_r \, \mathbf{u} \right) + \frac{\partial ( z_r \, \dot{r} ) }{\partial r} &= 0
+&\mbox{thickness} \label{eq:r-non-divergence} \\
+\frac{\partial ( \theta \, z_r ) }{\partial t} + \nabla_r \cdot \left( \theta z_r \, \mathbf{u} \right) + \frac{\partial ( \theta \, z_r \, \dot{r} )}{\partial r} &= z_r \boldsymbol{\mathcal{N}}_\theta^\gamma - \frac{\partial J_\theta^{(z)}}{\partial r}
+&\mbox{potential temp} \label{eq:r-temperature-equation} \\
+\frac{\partial ( S \, z_r) }{\partial t} + \nabla_r \cdot \left( S \, z_r \, \mathbf{u} \right) + \frac{\partial ( S \, z_r \, \dot{r} )}{\partial r} &= z_r \boldsymbol{\mathcal{N}}_S^\gamma - \frac{\partial J_S^{(z)}}{\partial r}
+&\mbox{salinity} \label{eq:r-salinity-equation} \\
+\rho &= \rho\left( S, \theta, -g \rho_0 z(r) \right)
&\mbox{equation of state.}
\f}
-The time derivatives appearing in these equations are computed with
-the generalized vertical coordinate fixed rather than the
-geopotential. It is a common misconception that the horizontal
-velocity, \f$\mathbf{u}\f$, is rotated to align with constant \f$r\f$
-surfaces. Such is not the case. Rather, the horizontal velocity,
-\f$\mathbf{u}\f$, is precisely the same horizontal velocity used with
-geopotential coordinates. However, its evolution has here been
-formulated using generalized vertical coordinates.
-As a finite volume model, MOM6 is discretized in the vertical by
-integrating between surfaces of constant \f$r\f$. The layer thickness
-is a basic term appearing in these equations, which results from
-integrating the specific thickness over a layer
-\f{align}
-h = \int z_r \, \mathrm{d}r.
-\f}
-Correspondingly, the model variables are treated as finite volume
-averages over each layer, with full accounting of this finite volume
-approach presented in Griffies, Adcroft and Hallberg (2020)
-\cite Griffies_Adcroft_Hallberg2020, and with the semi-discrete model
-ocean model equations written as follows.
-\f{align}
-\rho_0
-\left[ \frac{\partial \mathbf{u}}{\partial t} + \frac{( f + \zeta )}{h} \,
+The time derivatives are now computed with the generalized vertical
+coordinate fixed rather than the geopotential. We introduced the
+specific thickness, \f$z_r = \partial z/\partial r\f$, which measures the
+inverse vertical stratification of the vertical coordinate surfaces.
+
+ Similar to \cite bleck2002, MOM6 is discretized in the vertical by
+ integrating between surfaces of \f$r\f$ to yield layer equations where the
+ layer thickness is \f$h = \int z_r dr\f$ and variables are treated as finite
+ volume averages over each layer:
+
+\f{eqnarray}
+\label{html:h-equations}\notag \\
+\rho_0 \left( \frac{\partial \mathbf{u}}{\partial t} + \frac{( f + \zeta )}{h} \,
\hat{\mathbf{z}} \times h \, \mathbf{u} + \underbrace{ \dot{r} \,
-\frac{\partial \mathbf{u}}{\partial r} }
-\right]
-&= -\nabla_r \, (p + \rho_{0} \, K) -
-\rho \nabla_r \, \Phi + \mathbf{\mathcal{F}}
-&\mbox{horizontal momentum}
-\label{eq:h-horz-momentum}
-\\
-\rho \, \delta_r \Phi + \delta_r p
-&= 0
-&\mbox{hydrostatic}
-\label{eq:h-hydrostatic-equation}
-\\
+\frac{\partial \mathbf{u}}{\partial r} } + \nabla_r K \right) &= -\nabla_r \, p -
+\rho \nabla_r \, \Phi + \boldsymbol{\mathcal{F}}
+&\mbox{momentum} \label{eq:h-horz-momentum} \\
+\rho \, \delta_r \Phi + \delta_r p &= 0
+&\mbox{hydrostatic} \label{eq:h-hydrostatic-equation} \\
\frac{\partial h}{\partial t} + \nabla_r \cdot \left( h \, \mathbf{u} \right) +
-\underbrace{ \delta_r ( z_r \dot{r} ) }
- &= 0
-&\mbox{thickness}
-\label{eq:h-thickness-equation}
-\\
+\underbrace{ \delta_r ( z_r \dot{r} ) } &= 0
+&\mbox{thickness} \label{eq:h-thickness-equation} \\
\frac{\partial ( \theta \, h )}{\partial t} + \nabla_r \cdot \left( \theta h \,
-\mathbf{u} \right) + \underbrace{ \delta_r ( \theta \, z_r \dot{r} ) }
-&=
-h \mathbf{\mathcal{N}}_\theta^\gamma - \delta_r J_\theta^{(z)}
-&\mbox{potential/Conservative temp}
-\label{eq:h-temperature-equation}
-\\
+\mathbf{u} \right) + \underbrace{ \delta_r ( \theta \, z_r \dot{r} ) } &=
+h \boldsymbol{\mathcal{N}}_\theta^\gamma - \delta_r J_\theta^{(z)}
+&\mbox{potential temp} \label{eq:h-temperature-equation} \\
\frac{\partial ( S \, h )}{\partial t} + \nabla_r \cdot \left( S \, h \,
-\mathbf{u} \right) + \underbrace{ \delta_r ( S \, z_r \dot{r} ) }
-&=
-h \mathbf{\mathcal{N}}_S^\gamma - \delta_r J_S^{(z)}
-&\mbox{salinity}
-\label{eq:h-salinity-equation}
-\\
+\mathbf{u} \right) + \underbrace{ \delta_r ( S \, z_r \dot{r} ) } &=
+h \boldsymbol{\mathcal{N}}_S^\gamma - \delta_r J_S^{(z)}
+&\mbox{salinity} \label{eq:h-salinity-equation} \\
\rho &= \rho\left( S, \theta, -g \rho_0 z(r) \right)
&\mbox{equation of state,} \label{eq:h-equation-of-state}
\f}
-where
-\f{align}
-\delta_{r} = \mathrm{d}r \, (\partial/\partial r)
-\f}
-is the discrete vertical difference operator. The pressure gradient
-accelerations in the momentum equation are written in
-continuous-in-the-vertical form for brevity; the exact discretization
-is detailed in \cite adcroft2008 and
-\cite Griffies_Adcroft_Hallberg2020. The \f$1/h\f$ and \f$h\f$ appearing in
-the horizontal momentum equation are carefully handled in the code to
-ensure proper cancellation even when the layer thickness goes to zero;
-i.e., l'Hospital's rule is respected.
-The MOM6 time-stepping algorithm integrates the above layer-averaged
-equations forward in time allowing the vertical grid to follow the
-motion, i.e. \f$\dot{r}=0\f$, so that the underbraced terms are
-dropped. This approach is generally known as a Lagrangian method, with
-the Lagrangian approach in MOM6 limited to the vertical
-direction. After each Lagrangian step, a regrid step is applied that
-generates a new vertical grid of the user's choosing. The ocean state
-is then remapped from the old to the new grid. The physical state is
-not meant to change during the remap step, yet truncation errors make
-remapping imperfect. We employ high-order accurate reconstructions to
-minimize errors introduced during the remap step (\cite white2008,
-\cite white2009). The connection between time-stepping and remapping
-is described in section \ref ALE_Timestep.
+where \f$\delta_{r} = \mathrm{d}r \, (\partial/\partial r)\f$ is the discrete
+vertical difference operator. The pressure gradient accelerations
+in the momentum equation \eqref{eq:h-horz-momentum,h-equations,momentum} are written in
+continuous-in-the-vertical form for brevity; the exact discretization
+is detailed in \cite adcroft2008. The MOM6 time-stepping algorithm
+integrates the above layer-averaged equations forward allowing the
+vertical grid to follow the motion, i.e. \f$\dot{r}=0\f$, so that the underbraced
+terms are dropped. This approach is generally known as the Lagrangian
+method but here the Lagrangian method is used only in the vertical
+direction. After each Lagrangian step, a remap step is applied that
+generates a new vertical grid of the user's choosing. The ocean state is
+then mapped from the old to the new grid. The physical state is not meant
+to change during the remap step, yet truncation errors make remapping
+imperfect. We employ high-order accurate reconstructions to minimize
+errors introduced during the remap step (\cite white2008, \cite white2009). The
+connection between time-stepping and remapping is described in
+section \ref ALE_Timestep.
*/
diff --git a/src/core/_Governing.dox b/src/core/_Governing.dox
index d1bdeb6f0e..646ba52c09 100644
--- a/src/core/_Governing.dox
+++ b/src/core/_Governing.dox
@@ -1,176 +1,71 @@
/*! \page Governing_Equations Governing Equations
-MOM6 is a hydrostatic ocean circulation model that time steps either
-the non-Boussinesq ocean equations (where the flow velocity is
-divergent: \f$\nabla \cdot \mathbf{v} \ne 0\f$), or the Boussinesq
-ocean equations (where velocity is non-divergent: \f$\nabla \cdot
-\mathbf{v} = 0\f$). We here display the Boussinesq version since
-it is most commonly used (as of 2022). We start by casting the
-equations in geopotentiial coordinates prior to transforming to the
-generalized vertical coordinates used by MOM6. A more thorough
-discussion of these equations, and their finite volume realization
-appropriate for MOM6, can be found in Griffies, Adcroft and Hallberg (2020)
-\cite Griffies_Adcroft_Hallberg2020.
-
-The hydrostatic Boussinesq ocean equations, written using geopotential
-vertical coordinates, are given by
-\f{align}
- \rho_o \left[
- D_t \mathbf{u} + f \hat{\mathbf{z}} \times \mathbf{u}
- \right]
- &= -\rho \, \nabla_z \Phi - \nabla_z p
- + \rho_o \, \mathbf{\mathcal{F}}
- &\mbox{horizontal momentum}
-\\
- \rho \, \partial_{z} \Phi + \partial_{z} p &= 0 &\mbox{hydrostatic}
-\\
- \nabla_z \cdotp \mathbf{u} + \partial_{z} w
- &= 0
- &\mbox{continuity}
-\\
- D_t \theta &= \mathbf{\mathcal{N}}_\theta^\gamma
- - \partial_{z} J_\theta^{(z)}
- &\mbox{potential or Conservative temp}
- \\
- D_t S &= \mathbf{\mathcal{N}}_S^\gamma
-- \partial_{z} J_S^{(z)}
- &\mbox{salinity}
-\\
- \rho &= \rho(S, \theta, z) &\mbox{ equation of state}
-\\
- \mathbf{v} &= \mathbf{u} + \hat{\mathbf{z}} \, w &\mbox{velocity field.}
-\f}
+The Boussinesq hydrostatic equations of motion in height coordinates are
-The acceleration term, \f$\mathbf{\mathcal{F}}\f$, in the
-horizontal momentum equation includes the acceleration due to the
-divergence of internal frictional stresses as well as from bottom and
-surface boundary stresses. Other notation is described in \ref
-Notation.
-
-The prognostic temperature, \f$\theta\f$, is either potential
-temperature or Conservative Temperature, depending on the chosen
-equation of state, and \f$S\f$ is the salinity. We generally follow
-the discussion of \cite McDougall_etal_2021 for how to interpret the
-prognostic temperature and salinity in ocean models. MOM6 has
-historically used the Wright (1997) \cite wright1997 equation of state
-to compute the in situ density, \f$\rho\f$. However, there
-are other options as documented in \ref Equation_of_State. In the
-potential temperature and salinity equations, fluxes due to diabatic
-processes are indicated by \f$J^{(z)}\f$. Tendencies due to the
-convergence of fluxes oriented along neutral directions are indicated
-by \f$\mathbf{\mathcal{N}}^\gamma\f$, with our implementation of
-neutral diffusion detailed in Shao et al (2020)
-\cite Shao_etal_2020.
-
-The total or material time derivative operator is given by
-\f{align}
- D_t &\equiv \partial_{t} + \mathbf{v} \cdotp \nabla
- \\
- &= \partial_{t} + \mathbf{u} \cdotp \nabla_z + w \, \partial_{z},
-\f}
-where the second equality explosed the horizontal and vertical terms. Using the non-divergence condition
-on the three-dimensional velocity allows us to write the material time derivative of an arbitrary scalar field,
-\f$\psi\f$, into a flux-form equation
-\f{align} D_t \psi &= ( \partial_{t} + \mathbf{u} \cdotp \nabla) \, \psi
- \\
- &= \partial_{t} \psi + \nabla \cdotp (\mathbf{v} \, \psi)
-\\
- &= \partial_{t} \psi + \nabla_z \cdotp ( \mathbf{u} \, \psi) + \partial_{z} ( w \, \psi).
+\f{eqnarray} D_t \boldsymbol{u} + f \widehat{\boldsymbol{k}} \times \boldsymbol{u} + \frac{\rho}{\rho_o} \boldsymbol{\nabla}_z \Phi + \frac{1}{\rho_o} \boldsymbol{\nabla}_z p &= \boldsymbol{\mathcal{F}} &\mbox{ momentum} \\
+ \rho \, \frac{\partial \Phi}{\partial z} + \frac{\partial p}{\partial z} &= 0 &\mbox{ hydrostatic} \\
+ \boldsymbol{\nabla}_z \cdotp \boldsymbol{u} + \frac{\partial w}{\partial z} &= 0 &\mbox{ thickness} \\
+ D_t \theta &= \boldsymbol{\mathcal{N}}_\theta^\gamma - \frac{\partial J_\theta^{(z)}}{\partial z} &\mbox{ potential temp} \\
+ D_t S &= \boldsymbol{\mathcal{N}}_S^\gamma - \frac{\partial J_S^{(z)}}{\partial z} &\mbox{ salinity} \\
+ \rho &= \rho(S, \theta, z) &\mbox{ equation of state.}
\f}
-Discretizing the flux-form scalar equations means that fluxes
-transferring scalars between grid cells act in a conservative manner.
-Consequently, the domain integrated scalar (e.g., total seawater volume, total
-salt content, total potential enthalpy) is affected only via surface and bottom
-boundary transport. Such global conservation properties are
-maintained by MOM6 to within computational roundoff, with this level
-of precision found to be essential for using MOM6 to study
-climate. Making use of the flux-form scalar conservation equations
-brings the model equations to the form
-\f{align}
- \rho_o \left[
- D_t \mathbf{u} + f \hat{\mathbf{z}} \times \mathbf{u}
- \right]
- &= -\rho \, \nabla_z \Phi - \nabla_z p
- + \rho_o \, \mathbf{\mathcal{F}}
- &\mbox{horizontal momentum}
-\\
- \rho \, \partial_{z} \Phi + \partial_{z} p &= 0 &\mbox{hydrostatic}
-\\
- \nabla_z \cdotp \mathbf{u} + \partial_{z} w
- &= 0
- &\mbox{continuity}
-\\
-\partial_{t} \theta + \nabla_z \cdotp (\mathbf{u} \, \theta) + \partial_{z} (w \, \theta)
-&= \mathbf{\mathcal{N}}_\theta^\gamma - \partial_{z} J_\theta^{(z)}
-&\mbox{potential or Conservative temp}
-\\
-\partial_{t} S + \nabla_z \cdotp (\mathbf{u} \, S) + \partial_{z}(w \, S)
-&= \mathbf{\mathcal{N}}_S^\gamma -\partial_{z} J_S^{(z)}
- &\mbox{salinity}
-\\
-\rho &= \rho(S, \theta, z) &\mbox{equation of state.}
+
+where notation is described in \ref Notation, \f$\boldsymbol{\mathcal{F}}\f$ represents the accelerations due to
+the divergence of stresses including those provided through boundary interactions.
+
+The prognostic thermodynamic variables are potential temperature,
+\f$\theta\f$, and salinity \f$S\f$, which are related to in situ density
+\f$\rho\f$ through the \cite wright1997 equation of state. In the potential
+temperature and salinity equations, fluxes due to diabatic, vertically
+oriented processes are indicated by \f$J^{(z)}\f$. The tendency due to the
+convergence of fluxes oriented along neutral directions is indicated by
+\f$\boldsymbol{\mathcal{N}}^\gamma\f$. Our implementation of this neutral
+diffusion parameterization is detailed in Shao et al. (personal comm.)
+
+The total derivative is
+
+\f{eqnarray} D_t & \equiv \frac{\partial}{\partial t} + \boldsymbol{v} \cdotp \boldsymbol{\nabla} \\
+ &= \frac{\partial}{\partial t} + \boldsymbol{u} \cdotp \boldsymbol{\nabla}_z + w \frac{\partial}{\partial z}.
\f}
-\section vector_invariant_eqns Vector invariant velocity equation
-
-MOM6 time steps the horizontal velocity equation in its
-vector-invariant form. To derive this equation we make use of the
-following vector identity
-\f{align}
- D_t \mathbf{u}
- &=
- \partial_t \mathbf{u} + \mathbf{v} \cdotp \nabla \mathbf{u}
- \\
- &=
- \partial_t \mathbf{u} + \mathbf{u} \cdotp \nabla_z \mathbf{u} + w \partial_z \mathbf{u}
- \\
- &=
- \partial_t \mathbf{u} + \left( \nabla \times \mathbf{u} \right) \times \mathbf{v}
- + \nabla \left|\mathbf{u}\right|^2/2
- \\
- &=
- \partial_t \mathbf{u} + w \, \partial_{z} \mathbf{u}
- + \zeta \, \hat{\mathbf{z}} \times \mathbf{u} + \nabla_{z} K,
+The non-divergence of flow allows a total derivative to be re-written in flux form:
+
+\f{eqnarray} D_t \theta &= \frac{\partial}{\partial t} + \boldsymbol{\nabla} \cdotp ( \boldsymbol{v} \theta ) \\
+ &= \frac{\partial}{\partial t} + \boldsymbol{\nabla}_z \cdotp ( \boldsymbol{u} \theta ) + \frac{\partial ( w \theta )}{\partial z}.
\f}
-where we introduced the vertical component to the relative vorticity
-\f{align}
- \zeta = \hat{\mathbf{z}} \cdot (\nabla \times \mathbf{u})
- = \partial_{x}v - \partial_{y} u,
-\label{eq:relative-vorticity-z}
+
+The above equations of motion can thus be written as:
+
+\f{eqnarray} D_t \boldsymbol{u} + f \widehat{\boldsymbol{k}} \times \boldsymbol{u} + \frac{\rho}{\rho_o}\boldsymbol{\nabla}_z \Phi + \frac{1}{\rho_o} \boldsymbol{\nabla}_z p &= \boldsymbol{\mathcal{F}} &\mbox{ momentum}\\
+ \rho \, \frac{\partial \Phi}{\partial z} + \frac{\partial p}{\partial z} &= 0 &\mbox{ hydrostatic} \\
+ \boldsymbol{\nabla}_z \cdotp \boldsymbol{u} + \frac{\partial w}{\partial z} &= 0 &\mbox{ thickness} \\
+ \frac{\partial \theta}{\partial t} + \boldsymbol{\nabla}_z \cdotp ( \boldsymbol{u} \theta ) + \frac{\partial ( w \theta )}{\partial z} &= \boldsymbol{\mathcal{N}}_\theta^\gamma - \frac{\partial J_\theta^{(z)}}{\partial z} &\mbox{ potential temp} \\
+ \frac{\partial S}{\partial t} + \boldsymbol{\nabla}_z \cdotp ( \boldsymbol{u} S ) + \frac{\partial ( w S )}{\partial z} &= \boldsymbol{\mathcal{N}}_S^\gamma - \frac{\partial J_S^{(z)}}{\partial z} &\mbox{ salinity} \\
+ \rho &= \rho(S, \theta, z) &\mbox{ equation of state.}
\f}
-as well as the kinetic energy per mass contained in the horizontal flow
-\f{align}
- K = (u^{2} + v^{2})/2.
-\label{eq:kinetic-energy-per-mass}
+
+\section vector_invariant_eqns Vector Invariant Equations
+
+MOM6 solves the momentum equations written in vector-invariant form.
+
+A vector identity allows the total derivative of velocity to be written in the vector-invariant form:
+
+\f{eqnarray} D_t \boldsymbol{u} &= \partial_t \boldsymbol{u} + \boldsymbol{v} \cdotp \boldsymbol{\nabla} \boldsymbol{u} \\
+ &= \partial_t \boldsymbol{u} + \boldsymbol{u} \cdotp \boldsymbol{\nabla}_z \boldsymbol{u} + w \partial_z \boldsymbol{u} \\
+ &= \partial_t \boldsymbol{u} + \left( \boldsymbol{\nabla}
+ \times \boldsymbol{u} \right) \times \boldsymbol{v} + \boldsymbol{\nabla} \underbrace{\frac{1}{2} \left|\boldsymbol{u}\right|^2}_{\equiv K} .
\f}
-It is just the horizontal kinetic energy per mass that appears when
-making the hydrostatic approximation, whereas a non-hydrostatic fluid
-(such as the MITgcm) includes the contribution from vertical motion. With
-these identities we are led to the MOM6 flux-form equations of motion in
-geopotential coordinates
-\f{align}
- \rho_{o} \left[
- \partial_t \mathbf{u} + w \, \partial_{z} \mathbf{u}
- + (f + \zeta) \hat{\mathbf{z}} \times \mathbf{u}
- \right]
- &= -\nabla_{z} (p + K) - \rho \, \nabla_{z} \Phi + \rho_{o} \, \mathbf{\mathcal{F}}
- &\mbox{vector-invariant horz velocity}
-\\
- \rho \, \partial_{z} \Phi + \partial_{z} p &= 0 &\mbox{hydrostatic}
-\\
- \nabla_z \cdotp \mathbf{u} + \partial_{z} w
- &= 0
- &\mbox{continuity}
- \\
- \partial_t \theta + \nabla_z \cdotp ( \mathbf{u} \, \theta ) + \partial_z ( w \, \theta )
- &= \mathbf{\mathcal{N}}_\theta^\gamma - \partial_{z} J_\theta^{(z)}
- &\mbox{potential/Conservative temp}
- \\
- \partial_t S + \nabla_z \cdotp ( \mathbf{u} \, S ) + \partial_z (w \, S)
- &= \mathbf{\mathcal{N}}_S^\gamma - \partial_{z} J_S^{(z)}
- &\mbox{salinity}
- \\
- \rho &= \rho(S, \theta, z) &\mbox{equation of state.}
+
+The flux-form equations of motion in height coordinates can thus be written succinctly as:
+
+\f{eqnarray} \partial_t \boldsymbol{u} + \left( f \widehat{\boldsymbol{k}} +
+ \boldsymbol{\nabla} \times \boldsymbol{u} \right) \times \boldsymbol{v} + \boldsymbol{\nabla} K
+ + \frac{\rho}{\rho_o} \boldsymbol{\nabla} \Phi + \frac{1}{\rho_o} \boldsymbol{\nabla} p &= \boldsymbol{\mathcal{F}} &\mbox{ momentum} \\
+ \boldsymbol{\nabla}_z \cdotp \boldsymbol{u} + \partial_z w &= 0 &\mbox{ thickness} \\
+ \partial_t \theta + \boldsymbol{\nabla}_z \cdotp ( \boldsymbol{u} \theta ) + \partial_z ( w \theta ) &= \boldsymbol{\mathcal{N}}_\theta^\gamma - \frac{\partial J_\theta^{(z)}}{\partial z} &\mbox{ potential temp} \\
+ \partial_t S + \boldsymbol{\nabla}_z \cdotp ( \boldsymbol{u} S ) + \partial_z ( w S ) &= \boldsymbol{\mathcal{N}}_S^\gamma - \frac{\partial J_S^{(z)}}{\partial z} &\mbox{ salinity} \\
+ \rho &= \rho(S, \theta, z) &\mbox{ equation of state}
\f}
+where the horizontal momentum equations and vertical hydrostatic balance equation have been written as a single three-dimensional equation.
*/
diff --git a/src/core/_Notation.dox b/src/core/_Notation.dox
index c192929abb..faecb3b258 100644
--- a/src/core/_Notation.dox
+++ b/src/core/_Notation.dox
@@ -2,60 +2,34 @@
\section Symbols Symbols for variables
-\f$z\f$ refers to geopotential elevation (or height), increasing
-upward and with \f$z=0\f$ defining the resting ocean surface. Much of
-the ocean has \f$z < 0\f$.
+\f$z\f$ refers to elevation (or height), increasing upward so that for much of the ocean \f$z\f$ is negative.
-\f$x\f$ and \f$y\f$ are the Cartesian horizontal coordinates. MOM6
- uses generalized orthogonal curvilinear horizontal
- coordinates. However, the equations are simpler to write using
- Cartesian coordinates, and it is very straightforward to generalize
- the horizontal coordinates using the methods in Chapters 20 and 21 of
- \cite SMGbook.
+\f$x\f$ and \f$y\f$ are the Cartesian horizontal coordinates.
-\f$\lambda\f$ and \f$\phi\f$ are the geographic coordinates on a
-sphere (longitude and latitude, respectively).
+\f$\lambda\f$ and \f$\phi\f$ are the geographic coordinates on a sphere (longitude and latitude respectively).
-Horizontal components of velocity are indicated by \f$u\f$ and \f$v\f$
-and vertical component by \f$w\f$.
+Horizontal components of velocity are indicated by \f$u\f$ and \f$v\f$ and vertical component by \f$w\f$.
-\f$p\f$ is the hydrostatic pressure.
+\f$p\f$ is pressure and \f$\Phi\f$ is geo-potential:
-\f$\Phi\f$ is the geopotential. In the absence of tides, the
-geopotential is given by \f$\Phi = g z,\f$ whereas more general
-expressions hold when including astronomical tide forcing.
+ \f[ \Phi = g z .\f]
-The thermodynamic state variables can be salinity, \f$S\f$, and
-potential temperature, \f$\theta\f$. Alternatively, one can choose
-the Conservative Temperature if using the TEOS10 equation of state
-from \cite TEOS2010.
+The thermodynamic state variables are usually salinity, \f$S\f$, and potential temperature, \f$\theta\f$ or the absolute salinity and conservative temperature, depending on the equation of state. \f$\rho\f$ is in-situ density.
-\f$\rho\f$ is the in-situ density computed as a function
-\f$\rho(S,\theta,p)\f$ for non-Boussinesq ocean or
-\f$\rho(S,\theta,p=-g \, \rho_o \, z)\f$ for Boussinesq ocean. See
-Young (2010) \cite Young2010 or Section 2.4 of Vallis (2017)
-\cite GVbook for reasoning behind the simplified pressure
-used in the Boussinesq equation of state.
+\section vector_notation Vector notation
+The three-dimensional velocity vector is denoted \f$\boldsymbol{v}\f$
+ \f[\boldsymbol{v} = \boldsymbol{u} + \widehat{\boldsymbol{k}} w ,\f]
-\section vector_notation Vector notation
+where \f$\widehat{\boldsymbol{k}}\f$ is the unit vector pointed in the upward vertical direction and \f$\boldsymbol{u} = (u, v, 0)\f$ is the horizontal
+component of velocity normal to the vertical.
-The three-dimensional velocity vector is denoted \f$\boldsymbol{v}\f$
-and it is decomposed into its horizontal and vertical components according to
- \f[\boldsymbol{v}
- = \boldsymbol{u} + \hat{\boldsymbol{z}} w
- = \hat{\boldsymbol{x}} u + \hat{\boldsymbol{y}} v + \hat{\boldsymbol{z}} w,
- \f]
-where \f$\hat{\boldsymbol{z}}\f$ is the unit vector pointed in the
-upward vertical direction and \f$\boldsymbol{u} = (u, v, 0)\f$ is the
-horizontal component of velocity normal to the vertical.
-
-The three-dimensional gradient operator is denoted \f$\nabla\f$, and it is decomposed into
-its horizontal and vertical components according to
- \f[\nabla
- = \nabla_z + \hat{\boldsymbol{z}} \partial_z
- = \hat{\boldsymbol{x}} \partial_x + \hat{\boldsymbol{y}} \partial_y + \hat{\boldsymbol{z}} \partial_z.
- \f]
+The gradient operator without a suffix is three dimensional:
+
+ \f[\boldsymbol{\nabla} = ( \boldsymbol{\nabla}_z, \partial_z ) .\f]
+
+but a suffix indicates a lateral gradient along a surface of constant property indicated by the suffix:
+ \f[\boldsymbol{\nabla}_z = \left( \left. \partial_x \right|_z, \left. \partial_y \right|_z, 0 \right) .\f]
*/
diff --git a/src/parameterizations/vertical/_CVMix_KPP.dox b/src/parameterizations/vertical/_CVMix_KPP.dox
index 72c166c284..7a65b6a6a3 100644
--- a/src/parameterizations/vertical/_CVMix_KPP.dox
+++ b/src/parameterizations/vertical/_CVMix_KPP.dox
@@ -7,7 +7,7 @@
The formulation and implementation of KPP is described in great detail in the
[CVMix manual](https://github.com/CVMix/CVMix-description/raw/master/cvmix.pdf) (written by our own Steve Griffies).
- \section section_KPP_nutshell KPP in a nutshell
+ \section section_KPP_nutshell KPP in a nutshell
Large et al., \cite large1994, decompose the parameterized boundary layer turbulent flux of a scalar, \f$ s \f$, as
\f[ \overline{w^\prime s^\prime} = -K \partial_z s + K \gamma_s(\sigma), \f]