Split Time Stepping

Table of Contents

  1. Overview

  2. Requirements

  3. Algorithmic Formulation

  4. Design

  5. Verification and Testing

1. Overview

To enhance computational efficiency by allowing longer timesteps, ocean models require split barotropic-baroclinic time stepping methods. The implementation described here is based on the approach of Higdon (2005) and the MPAS-Ocean split_explicit scheme, with modifications for Omega’s non-Boussinesq pseudo-height \(\tilde{z}\) vertical coordinate.

Omega currently implements two closely related steppers: SplitExplicitRK2, which explicitly subcycles the barotropic mode, and UnsplitRK2, which uses the same outer predictor-corrector framework without a barotropic split. The only implemented barotropic algorithm is Predictor-Corrector; implicit barotropic stepping and an Adams–Bashforth split-explicit method are future extensions.

The split-explicit method involves the following sequence:

  • Decompose velocity into barotropic and baroclinic components.

  • Advance the baroclinic velocities using a large timestep, and compute the vertically averaged forcing, the \(\overline{G}\) term.

  • Subcycle the barotropic velocity using small explicit timesteps.

  • Construct a corrected transport velocity, update pseudo thickness and tracers, and reconstruct the physical velocity.

2. Requirements

2.1 Requirement: A split time-stepping method in the pseudo-height coordinate

The algorithm is based on Section 2.3 of Higdon (2005), with modifications to accommodate the \(\tilde{z}\)-coordinate variables in the non-Boussinesq framework. It accepts the outer time step, requested barotropic time step, barotropic algorithm, number of outer predictor-corrector iterations, number of baroclinic Coriolis iterations, and an option to recompute the velocity split. The UnsplitRK2 variant advances the full velocity through the baroclinic path and skips all barotropic operations.

2.2 Requirement: Stable time integration for long-term high-resolution simulations

Stability constrains the maximum allowable timestep, which in turn affects computational cost. The implemented split-explicit time stepping methods must allow for reasonably long timesteps while preventing numerical instabilities, such as those arising from internal gravity waves and barotropic modes, which are particularly important in global-scale and high-resolution ocean modeling. At a minimum, the time-stepping approach used in Omega V2 should accommodate the same timestep sizes as MPAS-Ocean for both the baroclinic and barotropic subsystems since Omega V2 is non-Boussinesq but hydrostatic.

2.3 Requirement: Modularization of baroclinic and barotropic time-stepping methods

Modularity ensures ease of testing and future-proofing of the Omega V2 codebase. Implementing a modular design enables mix-and-match time-stepping schemes of the baroclinic and barotropic subsystems, straightforward integration of alternative time-stepping schemes, and easier maintenance by separating the baroclinic and barotropic time-stepping codes, thereby enhancing flexibility.

3. Algorithmic Formulation

3.1 Barotropic (external) and baroclinic (internal) mode splitting

The split-explicit method separates ocean velocity into depth-integrated barotropic and depth-dependent baroclinic components. This separation allows computationally expensive baroclinic modes to run at longer timesteps and computationally efficient barotropic modes to run at shorter timesteps, enhancing computational efficiency and accuracy.

The layered discrete governing equations for Omega V2 are described in the Omega V1 governing equations design document. The mass, tracer, and velocity equations used as the starting point for the mode split are summarized below.

Mass:

(1)\[ \frac{\partial \tilde{h}_{i,k}}{\partial t} + \nabla \cdot \left( [\tilde{h}_{k}]_e {\bf u}_{e,k} \right) + \left[ \tilde{W}_{tr} \right]^{\text{top}}_{k} - \left[ \tilde{W}_{tr} \right]^{\text{top}}_{k+1} = 0 . \]

Tracer:

(2)\[ \frac{\partial \tilde{h}_{i,k}\varphi_{i,k}}{\partial t} + \nabla \cdot \left( [\tilde{h}_{i,k}\varphi_{i,k}]_e {\bf u}_{e,k} \right) + \left\{ \left[\varphi \tilde{W}_{tr}\right]^{\text{top}}_{k} - \left[\varphi \tilde{W}_{tr}\right]^{\text{top}}_{k+1} \right\} = [D_h^\varphi]_{i,k} - [D_v^\varphi]_{i,k}, \]

where

(3)\[ [D_h^\varphi]_{i,k} = \kappa_{2,e}\nabla^2\varphi_{i,k} - \kappa_{4,e}\nabla^4\varphi_{i,k}, \]

and

(4)\[ [D_v^\varphi]_{i,k} = [\tilde{\kappa}_{v}]_{i,k} \left[\frac{\partial h\varphi}{\partial \tilde{z}}\right]_{i,k} - [\tilde{\kappa}_{v}]_{i,k+1} \left[\frac{\partial h\varphi}{\partial \tilde{z}}\right]_{i,k+1} . \]

Velocity:

(5)\[ \frac{\partial {\bf u}_{e,k}}{\partial t} + \left[ {\bf k} \cdot \nabla \times {\bf u}_{e,k} + f_v \right]_e {\bf u}^{\perp}_{e,k} + [\nabla K]_e + \frac{1}{[\tilde{h}_{i,k}]_e} \left[\tilde{W}_{tr}\frac{\partial U}{\partial \tilde{z}}\right]_{e,k} = - [(\alpha \nabla p + \nabla \Phi)]_{e,k} + [D_h^{\bf u}]_{e,k} - [D_v^{\bf u}]_{e,k}, \]

where

(6)\[ [D_h^{\bf u}]_{e,k} = \nu_{2,e}\nabla^2 {\bf u}_{e,k} - \nu_{4,e}\nabla^4 {\bf u}_{e,k}, \]

and

(7)\[ [D_v^{\bf u}]_{e,k} = \frac{1}{[\tilde{h}_{i,k}]_e} \left\{ [\tilde{\nu}_{v}]^{\text{top}}_{e,k} \left[\frac{\partial {\bf u}}{\partial \tilde{z}}\right]_{e,k} - [\tilde{\nu}_{v}]^{\text{top}}_{e,k+1} \left[\frac{\partial {\bf u}}{\partial \tilde{z}}\right]_{e,k+1} \right\} . \]

Define the barotropic velocity, baroclinic velocity, barotropic pressure, and barotropic pressure anomaly as follows.

Barotropic velocity:

(8)\[ \overline{{\bf u}} \equiv \frac{1}{\tilde{H}} \sum_{k=0}^{K_{\max}} \tilde{h}_{k}{\bf u}_{k}. \]

Baroclinic velocity:

(9)\[ {\bf u}'_k \equiv {\bf u}_k - \overline{{\bf u}}. \]

Barotropic pressure:

(10)\[ B \equiv p^{\text{floor}} - p^{\text{surf}}. \]

Barotropic pressure anomaly:

(11)\[ B' \equiv B - \rho_0 g b, \]

so that

(12)\[ B = B' + \rho_0 g b, \]

and

(13)\[ p^{\text{floor}} = p^{\text{surf}} + B' + \rho_0 g b, \]

where \(b \equiv -z^{\text{floor}}\) is the geometric bottom depth.

Define the column-integrated pseudo thickness and geometric thickness as follows.

Column-integrated pseudo thickness, or mass-equivalent column depth:

(14)\[ \tilde{H} \equiv \sum_{k=0}^{K_{\max}} \tilde{h}_{k} = \frac{1}{\rho_0 g}\sum_{k=0}^{K_{\max}} \left(p_{k+1}^{\text{top}} - p_k^{\text{top}}\right) = \frac{1}{\rho_0 g}\left(p^{\text{floor}}-p^{\text{surf}}\right) = \frac{B}{\rho_0 g} = \frac{1}{\rho_0 g}\left(B' + \rho_0 g b\right). \]

Column-integrated geometric thickness, or geometric column depth:

(15)\[ H = \sum_{k=0}^{K_{\max}} h_k = \sum_{k=0}^{K_{\max}} \frac{\rho_0}{\rho_k}\tilde{h}_k = \rho_0 \sum_{k=0}^{K_{\max}} \alpha_k \tilde{h}_k = \rho_0 S, \]

where

(16)\[ S = \sum_{k=0}^{K_{\max}} \alpha_k \tilde{h}_k \]

is the column-integrated specific volume.

Sea-surface height can be diagnosed by

(17)\[ \eta = \rho_0 S - b. \]

The sea-surface height is decomposed into two components:

(18)\[ \eta = \eta^{\text{mass}} + \eta^{\text{steric}}, \]

where

(19)\[ \eta^{\text{mass}} = \frac{B'}{\rho_0 g}, \qquad \eta^{\text{steric}} = \rho_0 S - \frac{B}{\rho_0 g}. \]

This expression is equivalent to Eq. (10.60) in Griffies (2012) and Madec et al. (2015), but formulated using \(S\) and \(B\). The term \(\eta^{\text{steric}}\) represents the change in sea-surface height resulting from the non-Boussinesq steric effect, while \(\eta^{\text{mass}}\) is the mass-related component of sea-surface height that is advanced in time in the barotropic system as the fast process.

Under the Boussinesq approximation, the steric height \(\eta^{\text{steric}}\) reduces to zero:

(20)\[ \eta^{\text{steric}} = H - \tilde{H}; \qquad \rho \rightarrow \rho_0; \qquad H - H = 0. \]

A common expression for steric height, following Griffies (2012) and Madec et al. (2015), is

(21)\[ \eta^{\text{steric}} = -\int_{-H}^{\eta} \left(\frac{\rho - \rho_0}{\rho_0}\right) dz . \]

The derivation below demonstrates that this definition of \(\eta^{\text{steric}}\) is identical to the definition employed in the split system of Omega:

(22)\[\begin{split} \begin{aligned} \eta^{\text{steric}} &= -\int_{-H}^{\eta} \left(\frac{\rho - \rho_0}{\rho_0}\right) dz \\ &= \int_{-H}^{\eta} \left(1 - \frac{\rho}{\rho_0}\right) dz \\ &= \int_{-H}^{\eta} dz - \int_{-H}^{\eta}\frac{\rho}{\rho_0} dz \\ &= \int_{-H}^{\eta} dz - \int_{-H}^{\eta} d\tilde{z} \\ &= H - \tilde{H} \\ &= \rho_0 S - \frac{B}{\rho_0 g} . \end{aligned} \end{split}\]

The barotropic continuity and momentum equations are written as follows.

Barotropic continuity equation:

(23)\[ \frac{\partial B_i'}{\partial t} + \left[ \nabla \cdot \left( [(B_i' + \rho_0 g b_i)]_e \overline{{\bf u}}_e \right) \right]_i = -\rho_0 g Q_i . \]

Barotropic momentum equation:

(24)\[ \frac{\partial \overline{{\bf u}}_e}{\partial t} + f_e \overline{{\bf u}}^{\perp}_e = -[\overline{\alpha}_i]_e [\nabla B_i']_e + \overline{G}_e . \]

Baroclinic momentum equation:

(25)\[ \frac{\partial {\bf u}'_{e,k}}{\partial t} = -[f_v]_e {\bf u}_{e,k}^{\prime\perp} + \Gamma_{e,k} + [\overline{\alpha}_i]_e [\nabla B_i']_e - \overline{G}_e, \]

where

(26)\[ \Gamma_{e,k} \equiv -[\nabla K]_e -(\alpha \nabla p + \nabla \Phi)_{e,k} - [ {\bf k} \cdot \nabla \times {\bf u}_{e,k}]_e {\bf u}^{\perp}_{e,k} - \frac{1}{[\tilde{h}_{i,k}]_e} \left[\tilde{W}_{tr}\frac{\partial U}{\partial \tilde{z}}\right]_{e,k} + [D_h^{\bf u}]_{e,k} + [F_s^{\bf u}]_{e,k} + [F_b^{\bf u}]_{e,k}. \]

Here, \(\overline{G}_e\) includes all remaining terms in the barotropic equation.

Each term of \(\Gamma\) is applied only when its own tendency switch is enabled, and each is evaluated from the full normal velocity of the working state rather than from the baroclinic velocity alone. \([F_s^{\bf u}]\) and \([F_b^{\bf u}]\) are the surface stress forcing and the explicit bottom drag, which the implementation includes in the baroclinic tendency alongside the terms of Eq. (5). A user-registered custom velocity tendency, when present, is also added here. \(Q\) is the surface freshwater flux.

The vertical mixing term \([D_v^{\bf u}]\) does not appear in \(\Gamma\). Vertical mixing of momentum and tracers is applied implicitly once per slow step, after the outer iterations and the time-level rotation have completed; see Section 3.2.6. Consequently \(\Gamma\) contains only the terms integrated explicitly by the outer scheme.

3.2 Split-explicit time stepping algorithm

The mode-splitting algorithm follows the MPAS-Ocean split-explicit framework. It advances baroclinic velocity over the slow time step, explicitly subcycles the barotropic system, and uses the time-averaged barotropic transport to correct the velocity used for pseudo-thickness and tracer transport. The implementation calls this option SplitExplicitRK2.

SplitExplicitRK2 uses two outer iterations by default. The first evaluates momentum at the time-\(n\) state and constructs a midpoint state. The second evaluates momentum at that midpoint and retains a full-step conservative pseudo-thickness and tracer update. Each outer iteration passes the full slow time step \(\Delta t\) to all three stages; it does not use \(\Delta t/\text{NTimeStepIteration}\). The midpoint is produced by the one-half coefficients in the baroclinic, pseudo-thickness, and tracer updates. Setting NTimeStepIteration to one is permitted but leaves only the predictor and reduces the outer method to forward Euler.

An SE-AB2 method is not implemented. It can be added later using the same barotropic-stepper interface and additional previous-tendency storage.

3.2.1 Initialization

The default outer-iteration count is

(27)\[ \text{NTimeStepIteration} = 2. \]

Compute NBtrSubcycle as

(28)\[ \text{NBtrSubcycles} = \max\left(1, \left\lceil \frac{\text{TimeStep}}{\text{BtrTimeStep}} \right\rceil\right). \]

BtrTimeStep is therefore an upper bound requested by the user. The effective barotropic step is

(29)\[ \Delta t_{\mathrm{bt}} = \frac{\Delta t}{\text{NBtrSubcycles}}. \]

Compute the barotropic velocity:

(30)\[ \overline{{\bf u}}_e \equiv \frac{1}{[\tilde{H}_i]_e} \sum_{k=0}^{K_{\max}} [\tilde{h}_{i,k}]_e {\bf u}_{e,k}. \]

For a startup from initial conditions, the code computes this split using the arithmetic mean of the two neighboring cell pseudo thicknesses at every edge, and masks the result with the edge mask of the top active layer so that dry edges carry no barotropic velocity. For a restart, NormalBarotropicVelocity and BarotropicPressureAnomaly are preserved, while NormalBaroclinicVelocity is made consistent with the read full and barotropic velocities through \({\bf u}'={\bf u}-\overline{\bf u}\).

Compute the baroclinic velocity:

(31)\[ {\bf u}'_{e,k} = {\bf u}_{e,k} - \overline{{\bf u}}_e. \]

Before initializing the pressure split, the equation of state and the momentum vertical auxiliary fields are computed from the time-level-0 state and tracers through AuxiliaryState::computeMomVertAux. That call computes, in order, the interface and midpoint pressures together with the column-integrated pseudo thickness \(\tilde H\), the specific volume, the depth-integrated and depth-mean specific volume, the column-integrated geometric thickness \(H=\rho_0 S\), the geometric height, and the target thickness.

Compute the barotropic pressure:

(32)\[ B_i = p_{i,K_{\max}+1} - p_{i}^{\text {surf}}. \]

Here \(p_{i,0}\) is the surface pressure and \(p_{i,K_{\max}+1}\) is the pressure at the interface below the bottom-most active layer. \(B\) is stored once as a diagnostic in the split-explicit scratch space and is not used again during time stepping; the prognostic quantity is the anomaly \(B'\).

Compute the barotropic pressure anomaly:

(33)\[ B_i' = B_i - \rho_0 g b_i. \]

Time level 1 is then initialized from time level 0 for normal velocity, baroclinic velocity, barotropic velocity, pseudo thickness, and barotropic pressure anomaly; tracers are copied at the start of doStep. If ReinitSplitVelocity is true, the velocity split is recomputed from NormalVelocity at the beginning of each slow step, including a new depth mean for NormalBarotropicVelocity. When it is false, the existing barotropic velocity is retained and only baroclinic velocity is reconstructed from the other two velocity fields. In conceptual notation, the working state initially satisfies

(34)\[ {\bf u}^{*}_{e,k} = {\bf u}^n_{e,k}, \qquad \tilde{W}^{*}_{i,k} = \tilde{W}^{n}_{i,k}, \qquad \tilde{h}^{*}_{i,k} = \tilde{h}^{n}_{i,k}, \qquad \varphi^{*}_{i,k} = \varphi^{n}_{i,k}, \qquad p^{*}_{i,k} = p^{n}_{i,k}, \qquad \tilde{z}^{*}_{i,k} = \tilde{z}^{n}_{i,k}, \quad \text{etc.} \]

3.2.2 Stage 1: Baroclinic velocity advance with long time step

This stage advances the baroclinic velocity \(u'\) with the long time step and computes the barotropic forcing term \(\overline{G}\).

Compute the baroclinic forcing plus the barotropic pressure-gradient contribution:

(35)\[ \Gamma^*_{e,k} + [\overline{\alpha}^{*}_{i}]_e [\nabla B_i^{\prime *}]_e . \]

Both terms are evaluated from the working state: the pseudo thickness, full normal velocity, barotropic velocity, and barotropic pressure anomaly are all taken at the working time level, and the momentum auxiliary variables are recomputed from that state first. The barotropic pressure-anomaly gradient term is gated on SplitFactor, because it exists solely to cancel the barotropic part of the full pressure gradient and must be present whenever the mode split is active. The horizontal-advection term uses the relative vorticity alone, with the planetary contribution removed, so that the Coriolis acceleration can be iterated separately below.

Compute the column-integrated pseudo thickness:

(36)\[ \tilde{H}^{*}_i = \sum_{k=0}^{K} \tilde{h}^{*}_{i,k}. \]

The non-Coriolis tendency is saved in BaseVelocityTend. Compute the Coriolis term using a fixed-point treatment. For \(j = 0, \ldots, \text{NBclCoriolisIteration}-1\), with default value 2, restore the saved tendency after the first iteration, add the Coriolis acceleration computed from the latest working baroclinic velocity, recompute \(\overline G\), and update the baroclinic velocity from the fixed time-\(n\) state:

(37)\[ {\bf u}_{e,k}^{\prime\perp *} = \sum_{e'\in ECP(e)} \tilde{E}_{e,e'} f_{e'} {\bf u}_{e',k}^{\prime *}. \]

The Coriolis kernel shares its enable switch with the potential-vorticity horizontal-advection term: when that tendency is disabled, the Coriolis acceleration is skipped as well, in both the baroclinic and barotropic systems.

Before each update, form the provisional full-step baroclinic velocity

(38)\[ {\bf u}_{e,k}^{\prime\mathrm{prov}} = {\bf u}_{e,k}^{\prime n} + \Delta t\,{\cal R}_{e,k}^{*}, \]

where \({\cal R}^{*}\) is the current velocity tendency including the Coriolis term. The implemented barotropic forcing is the pseudo-thickness-weighted column mean of this provisional velocity divided by \(\Delta t\):

(39)\[ \overline G_e^{*} = \frac{1}{\Delta t} \frac{\displaystyle\sum_k [\tilde h_k^*]_e {\bf u}_{e,k}^{\prime\mathrm{prov}}} {\displaystyle\sum_k [\tilde h_k^*]_e}. \]

Here \([\tilde h_k^*]_e\) is the arithmetic mean of the two neighboring cell pseudo thicknesses at the working time level, and \({\bf u}^{\prime n}\) is the baroclinic velocity at time level 0, held fixed through the Coriolis iteration. The result is stored in BarotropicForcing and is the \(\overline G\) used unchanged by every barotropic subcycle pass in Stage 2.

The forcing is subtracted from every active layer of \({\cal R}^{*}\), ensuring that the tendency passed to the baroclinic update has its depth-mean component removed. The working baroclinic velocity is then advanced to the midpoint:

(40)\[ {\bf u}_{e,k}^{\prime n+1/2} = {\bf u}_{e,k}^{\prime n} + \frac{\Delta t}{2} \left( -[f_v]_e {\bf u}_{e,k}^{\prime\perp *} + \Gamma^*_{e,k} + [\overline{\alpha}^{*}_{i}]_e [\nabla B_i^{\prime *}]_e - \overline G_e^{*} \right). \]

For UnsplitRK2, set

(41)\[ \overline{G}^{*}_e = 0 \]

and use one Coriolis iteration. After Stage 1, the implementation exchanges only the working baroclinic-velocity halo.

3.2.3 Stage 2: Barotropic velocity advance, explicitly subcycled

This stage advances \(B'\) and \(\overline{\bf {u}}\) as a coupled system through \(2M\) predictor-corrector passes, where \(M=\text{NBtrSubcycles}\) and \(\Delta t_{\mathrm{bt}}=\Delta t/M\). As in the MPAS-Ocean algorithm, the passes span an extended \(2\Delta t\) barotropic averaging window. For UnsplitRK2, this stage is skipped, \(\overline u=0\), and \(u=u'\).

The discrete barotropic continuity equation is

(42)\[ B_i^{\prime n+1} = B_i^{\prime n} - \Delta t \left[ \nabla \cdot \left( [(B_i^{\prime n}+\rho_0 g b_i)]_e \overline{{\bf u}}_e^n \right) \right]_i - \Delta t\,\rho_0 g Q_i^n. \]

The discrete barotropic momentum update is

(43)\[ \overline{{\bf u}}^{n+1}_e = \overline{{\bf u}}^{n}_e + \Delta t \left( - f_e \overline{{\bf u}}^{\perp n}_e - [\overline{\alpha}_i^{*}\nabla B_i^{\prime n}]_e + \overline{G}^{*}_e \right). \]

Throughout this section a hat marks a quantity that lives inside the barotropic subcycle and changes from pass to pass, as distinct from the unhatted outer-iteration quantities carried by Stages 1 and 3.

At the start of an outer iteration, initialize the subcycle variables and accumulators:

(44)\[ \hat{\overline{{\bf u}}}^{n}_e = \overline{{\bf u}}^{n}_e, \qquad \hat{B}^{\prime n}_i = B_i^{\prime n}, \qquad \overline{{\bf u}}_e^{\mathrm{sum}}=\overline{{\bf u}}_e^n, \qquad F_e^{\mathrm{sum}}=0. \]

The barotropic mass flux is not formed directly from \((B'+\rho_0gb){\overline {\bf{u}}}\). Instead, the implementation uses the most recent baroclinic edge pseudo thickness and a pressure-anomaly correction relative to the provisional outer-iteration state \(B'^{*}\):

(45)\[ {\cal P}_e({\cal B}') = \rho_0g\sum_k[\tilde h_k^*]_e + [{\cal B}'-B'^{*}]_e. \]

Here \({\cal B}'\) is a formal argument rather than a field: it stands for whichever subcycle-level pressure anomaly the caller supplies. It is set calligraphic rather than hatted so that it cannot be read as any one particular buffer, following the same convention as the other calligraphic symbols in this document, \({\cal P}_e\), \({\cal R}\), and \({\cal T}\), none of which is a prognostic field. Each pass forms a barotropic mass flux in three of its kernels, the \(B'\) predictor, the \(B'\) corrector, and the transport accumulation, and each supplies a different blend of the subcycle buffers, given with the equations below.

The two velocity kernels are not among them: they need the gradient of the pressure anomaly across an edge, not the edge pressure \({\cal P}_e\) that multiplies a velocity to make a flux. Only the bracketed difference in Eq. (45) depends on \({\cal B}'\); the column sum is common to all three.

The column sum \(\sum_k[\tilde h_k^*]_e\) uses the flux pseudo thickness on edges produced by the pseudo-thickness auxiliary state, so it already carries the configured centered or upwind edge reconstruction. It is evaluated once, before the subcycle loop begins, and is held fixed for all \(2M\) passes.

The reference \(B'^{*}\) is the provisional barotropic pressure anomaly at the working time level on entry to Stage 2, in the same sense the asterisk carries elsewhere in Section 3.2: \(B'^n\) on the first outer iteration and the \(\rho_0 g(\tilde H^{*}-b)\) of Eq. (70) left by the previous outer iteration afterward. It is likewise held fixed and is only overwritten by the final corrected anomaly after the last pass, so \({\cal P}_e\) measures the departure of the subcycled anomaly from the baroclinic column mass. This is what enforces the Hallberg and Adcroft (2009) consistency described in Section 3.2.4. Note that \(B'^{*}\) is unhatted and so is an outer-iteration quantity; the hatted \(\hat B^{\prime *}\) appearing in the predictor-corrector equations below is a different thing, the subcycle predictor output.

The edge value of the anomaly correction uses the same centered or upwind choice selected for pseudo-thickness fluxes. With the upwind choice, ties at zero normal velocity take the larger of the two neighboring cell values.

The depth-mean specific volume \([\overline{\alpha}_i]_e\) is Eos::DepthMeanSpecificVolume averaged to the edge. It is computed from the working state during Stage 1 and is not updated during subcycling.

Predictor-Corrector is the only barotropic time stepper currently implemented, so the passes below are specific to it. Section 4.2.3 gives the interface an alternative barotropic stepper would implement.

For each predictor-corrector pass \(m=0,\ldots,2M-1\), use the following steps. Both velocity updates are multiplied by the edge mask of the top active layer, and edges with no active layers are set to zero.

\(\overline{\bf {u}}\) predictor:

(46)\[ [\hat{\overline{{\bf u}}}^{*}_e]^{n+(m+1)/M} = [\hat{\overline{{\bf u}}}_e]^{n+m/M} + \frac{\Delta t}{M} \left( -f_e[\hat{\overline{{\bf u}}}_e^{\perp}]^{n+m/M} -[\overline{\alpha}_i]^{*}_e[\nabla \hat{B}_i']^{n+m/M}_e + \overline{G}_e \right). \]

\(B'\) predictor:

(47)\[ [F^{*}_e]^{m+1} = {\cal P}_e\left([\hat{B}_i']^{n+m/M}\right) \left( (1-\gamma_1)[\hat{\overline{{\bf u}}}_e]^{n+m/M} + \gamma_1[\hat{\overline{{\bf u}}}^{*}_e]^{n+(m+1)/M} \right), \]

and

(48)\[ [\hat{B}_i^{\prime *}]^{n+(m+1)/M} = [\hat{B}_i']^{n+m/M} + \frac{\Delta t}{M} [-\nabla \cdot [F^{*}_e]^{m+1}]_i. \]

\(\overline{\bf {u}}\) corrector:

(49)\[ [\hat{\overline{{\bf u}}}_e]^{n+(m+1)/M} = [\hat{\overline{{\bf u}}}_e]^{n+m/M} + \frac{\Delta t}{M} \left( -f_e[\hat{\overline{{\bf u}}}_e^{\perp *}]^{n+(m+1)/M} - [\overline{\alpha}_i]^{*}_e \nabla \left( (1-\gamma_2)[\hat{B}_i']^{n+m/M} + \gamma_2[\hat{B}_i^{\prime *}]^{n+(m+1)/M} \right)_e + \overline{G}_e \right). \]

\(B'\) corrector:

(50)\[ [F_e]^{m+1} = {\cal P}_e\left( (1-\gamma_2)[\hat{B}_i']^{n+m/M} + \gamma_2[\hat{B}_i^{\prime *}]^{n+(m+1)/M} \right) \left( (1-\gamma_3)[\hat{\overline{{\bf u}}}_e]^{n+m/M} + \gamma_3[\hat{\overline{{\bf u}}}^{*}_e]^{n+(m+1)/M} \right), \]

and

(51)\[ [\hat{B}_i^{\prime \mathrm{cor}}]^{n+(m+1)/M} = [\hat{B}_i']^{n+m/M} + \frac{\Delta t}{M} [-\nabla \cdot [F_e]^{m+1}]_i. \]

On hatted subcycle quantities the asterisk denotes a different buffer depending on where it appears, so the mapping to the code is worth stating explicitly. Starred hatted pressures are always the predictor output Pre. Starred hatted velocities are the predictor output Pre in the \(\overline {\bf{u}}\) predictor, in the \(B'\) predictor flux, and in the Coriolis term of the \(\overline {\bf {u}}\) corrector; they are the corrector output Cor in the \(B'\) corrector flux and in the transport accumulation. On unhatted quantities such as \(B'^{*}\), \(\tilde h^{*}\), and \(\overline{\alpha}^{*}\), the asterisk keeps its Section 3.2 meaning of a provisional outer-iteration value, which is frozen for the whole of Stage 2. The coefficients are

(52)\[ \gamma_1=0.5333,\qquad \gamma_2=0.5333,\qquad \gamma_3=1. \]

Accumulate the barotropic velocity and flux during subcycling:

(53)\[ \overline{{\bf u}}_e^{\mathrm{sum}} = \sum_{m=0}^{2M} [\hat{\overline{{\bf u}}}_e]^{n+m/M}, \]

and

(54)\[ F_e^{\mathrm{sum}} = \sum_{m=0}^{2M-1} [F_e]^{(m+1)/M}. \]

Compute the time average after subcycling:

(55)\[ \overline{ {\bar {\bf u}}}_e^{\text{bt}} = \frac{1}{2M+1} \overline{{\bf u}}_e^{\mathrm{sum}} \]

and

(56)\[ \overline{F}_e^{\text{bt}} = \frac{1}{2M}F_e^{\mathrm{sum}} \]

The flux accumulated for transport uses the corrected velocity and evaluates \({\cal P}_e\) with \((1-\gamma_2)\hat B^{\prime n+m/M} +\gamma_2\hat B^{\prime\mathrm{cor}}\), which differs from the pressure used in the \(B'\) corrector flux above, where the second term is the predictor anomaly. Both accumulations run over owned edges only.

3.2.4 Barotropic-baroclinic coupling and barotropic pressure consistency

For the mode-split consistency of the barotropic pressure anomaly between \(B'\) from the barotropic mode and \(\rho_0 g(\tilde{H}-b)\) from the baroclinic mode, Omega follows the scheme from Hallberg and Adcroft (2009), as implemented in MPAS-Ocean.

The barotropic update of \(B'\) is given by

(57)\[ \frac{B_i^{\prime n+1}-B_i^{\prime n}}{\Delta t} + \left[ \nabla \cdot \overline{\left([(B_i' + \rho_0 g b_i)]_e \overline{{\bf u}}_e\right)}^{\text{bt}} \right]_i = -\rho_0 g Q_i . \]

Here, \(\overline{\varphi}^{\text{bt}}\) denotes a time-averaged quantity from the barotropic subcycles, and \(n\) indicates the baroclinic time step.

The velocity correction \({\bf u}^{\text{co}}\) is written as

(58)\[ {\bf u}_{e,k}^{\text{co}} = \left\{ \frac{\overline{F}_e^{\text{bt}}}{\rho_0g} - \sum_{k=0}^{K} [\tilde{h}_i^{*}]_{e,k} \left( \overline{{\bar {\bf u}}}_e^{\text{bt}} +{\bf u}_{e,k}^{\prime n+0.5} \right) \right\} \bigg/ [\tilde{H}_i^{*}]_e , \]

where \(\overline{F}_e^{\text{bt}}\) is the time-averaged pressure transport from Stage 2. The same correction is added uniformly to every active layer on an edge. Here \([\tilde h_i^*]_{e,k}\) is the arithmetic mean of the two neighboring cell pseudo thicknesses at the working time level, not the flux pseudo thickness used by \({\cal P}_e\) in Stage 2.

The asterisk indicates a provisional variable updated during the outer time-step iteration; the most recent available value is used for forcing terms.

The same kernel that forms \({\bf u}^{\text{co}}\) also writes the uncorrected sum \(\overline{{\bar{\bf u}}}^{\text{bt}}_e + {\bf u}^{\prime n+0.5}_{e,k}\) into the working NormalVelocity, which is what the Stage-3 auxiliary variables and tendencies see. On the final outer iteration this value is replaced later by the \(n+1\) reconstruction of Eq. (73).

The transport velocity \(u^{\text{tr}}\) is defined as

(59)\[ {\bf u}^{\text{tr}}_{e,k} = \overline{{\bar {\bf u}}}_e^{\text{bt}} +{\bf u}_{e,k}^{\prime n+0.5} +{\bf u}_{e,k}^{\text{co}}. \]

The transport velocity is used to compute vertical transport velocity and horizontal transport for both pseudo thickness and tracers.

For the unsplit algorithm, the above processes are skipped except that

(60)\[ {\bf u}_{e,k}^{\text{tr}} = {\bf u}_{e,k}^{\prime n+0.5}, \]

where

(61)\[ {\bf u}_{e,k}^{\prime n+0.5} = {\bf u}_{e,k}^{n+0.5}. \]

3.2.5 Stage 3: Update tracers and diagnostics

Compute the pseudo-thickness and tracer auxiliary variables from \({\bf u}_{e,k}^{\text{tr}}\) through AuxiliaryState::computePseudoThicknessTracerAux, which also refreshes the vertical momentum auxiliary variables, and therefore the pressure, specific volume, column thicknesses, geometric height, and target thickness, from the working pseudo thickness and tracers.

Compute \(\tilde{W}_{i,k}^{*}\) using \({\bf u}_{e,k}^{\text{tr}}\).

Compute pseudo thickness tendencies using \({\bf u}_{e,k}^{\text{tr}}\):

(62)\[ {\cal T}_{\tilde h} = - \nabla \cdot \left([\tilde{h}_{k}^{*}]_e {\bf u}_{e,k}^{\text{tr}}\right) - \left( [\tilde{W}_{tr}]^{\text{top}}_{k} - [\tilde{W}_{tr}]^{\text{top}}_{k+1} \right). \]

Compute tracer tendencies using \({\bf u}_{e,k}^{\text{tr}}\):

(63)\[ {\cal T}_{\varphi} = - \nabla \cdot \left([ \tilde{h}_{i,k}^{*}]_e [\varphi_{i,k}^{*}]_e {\bf u}_{e,k}^{\text{tr}} \right) - \left( [\varphi^{*}\tilde{W}_{tr}]^{\text{top}}_{i,k} - [\varphi^{*}\tilde{W}_{tr}]^{\text{top}}_{i,k+1} \right). \]

On a non-final outer iteration, pseudo thickness is updated with \(\Delta t/2\):

(64)\[ \tilde h^{*} = \tilde h^n+\frac{\Delta t}{2}{\cal T}_{\tilde h}. \]

For each tracer, the implementation first forms the conservative provisional end concentration:

(65)\[ \varphi^{\mathrm{end}} =\frac{\tilde h^n\varphi^n+\Delta t\,{\cal T}_\varphi} {\tilde h^n+\Delta t\,{\cal T}_{\tilde h}}, \]

and stores \(\varphi^*=(\varphi^n+\varphi^{\mathrm{end}})/2\). On the final outer iteration, it retains the full-step conservative updates shown above.

3.2.6 Reset variables

The reconstruction of the full velocity and the reset of \(\tilde H\) and \(B'\) are performed by a single routine at the end of Stage 3 on every outer iteration. Only the velocity reconstruction differs between an intermediate iteration and the final one; the column pseudo thickness and the barotropic pressure anomaly are recomputed the same way in both cases. For unsplit stepping the routine reconstructs the velocity and returns without touching \(B'\).

If iterating, reset the provisional variables as follows:

(66)\[ {\bf u}^{\prime *} = {\bf u}^{\prime n+0.5} \quad \text{from Stage 1}, \]
(67)\[ \overline{{\bf u}}^{*} = \overline{{\bar {\bf u}}}^{\text{bt}} \quad \text{from Stage 2}, \]
(68)\[ {\bf u}^{*} = \overline{{\bf u}}^{*} + u^{\prime *}, \]

with pseudo thickness and tracer concentration reset by Eqs. (64) and (65) of Section 3.2.5. The column totals then follow:

(69)\[ \tilde{H}^{*} = \sum_{k=0}^{K} \tilde{h}_k^{*}, \]

and

(70)\[ B^{\prime *} = \rho_0 g\left(\tilde{H}^{*}-b\right). \]

The full working velocity is the sum of the working split velocities. The state and tracer halos are exchanged before the next outer iteration.

After the final iteration, Stage 1 still stores the working baroclinic velocity at the midpoint:

(71)\[ {\bf u}^{\prime *}={\bf u}^{\prime n+1/2}. \]
(72)\[ \overline{{\bf u}}^{n+1} = \overline{{\bar {\bf u}}}^{\text{bt}} \quad \text{from Stage 2}, \]
(73)\[ {\bf u}^{n+1} = \overline{{\bf u}}^{n+1} +2{\bf u}^{\prime n+1/2}-{\bf u}^{\prime n}, \]
(74)\[ \psi^{n+1} \quad \text{is retained for pseudo thickness and tracers}, \]
(75)\[ \tilde{H}^{n+1} = \sum_{k=0}^{K} \tilde{h}_k^{n+1}, \]

and

(76)\[ B^{\prime n+1} = \rho_0 g\left(\tilde{H}^{n+1}-b\right). \]

The state and tracer time levels are rotated once, making the completed working state the new time level 0. Kinetic energy, velocity divergence, and the other kinetic auxiliary variables are then recomputed from the rotated \(n+1\) velocity so that history output and state validation see diagnostics consistent with the reconstructed velocity. Implicit vertical mixing is applied next if either velocity or tracer vertical mixing is enabled, followed by a state and tracer halo exchange. Finally the state is validated and the clock and persistent step count advance.

3.3 Unsplit time stepping algorithm

UnsplitRK2 uses the same SplitExplicitRK2Stepper implementation with SplitFactor=0. The stored baroclinic velocity is set equal to full velocity, barotropic velocity is zero, barotropic forcing is zero, and Stage 2 is skipped. Barotropic configuration keys are ignored.

3.3.1 Initialization

Compute the pressure \(p\) and the other vertical momentum auxiliary variables.

Set the velocity split trivially, at both time levels and for both a startup and a restart:

(77)\[ \overline{{\bf u}}_e = 0, \qquad {\bf u}'_{e,k} = {\bf u}_{e,k}. \]

Barotropic pressure and pressure anomaly are not initialized, since \(B'\) is never referenced when SplitFactor is zero.

Prepare variables before the first iteration:

(78)\[ {\bf u}^{*}_{e,k} = {\bf u}^n_{e,k}, \qquad \tilde{W}^{*}_{i,k} = \tilde{W}^{n}_{i,k}, \qquad \tilde{h}^{*}_{i,k} = \tilde{h}^{n}_{i,k}, \qquad \varphi^{*}_{i,k} = \varphi^{n}_{i,k}, \qquad p^{*}_{i,k} = p^{n}_{i,k}, \qquad \tilde{z}^{*}_{i,k} = \tilde{z}^{n}_{i,k}, \quad \text{etc.} \]

3.3.2 Stage 1: Velocity advance

This stage advances the full velocity \(u\).

Compute \(\Gamma_{e,k}^{*}\):

(79)\[ \Gamma_{e,k}^{*} = -[{\bf k}\cdot\nabla\times {\bf u}_{e,k}]_e {\bf u}_{e,k}^{\perp} -[\nabla K]_e - \frac{1}{[\tilde{h}_{i,k}]_e} \left[\tilde{W}_{tr}\frac{\partial U}{\partial \tilde{z}}\right]_{e,k} -(\alpha\nabla p+\nabla\Phi)_{e,k} +[D_h^{\bf u}]_{e,k} +[F_s^{\bf u}]_{e,k} +[F_b^{\bf u}]_{e,k}. \]

As in Section 3.1, the vertical mixing term is absent because it is applied implicitly after the outer iterations, and the barotropic pressure-anomaly gradient term of Eq. (35) is absent because SplitFactor is zero.

Compute the column-integrated pseudo thickness:

(80)\[ \tilde{H}_i^{*} = \sum_{k=0}^{K} \tilde{h}^{*}_{i,k}. \]

The unsplit configuration forces NBclCoriolisIteration=1. Compute the Coriolis term from the current working full velocity:

(81)\[ {\bf u}_{e,k}^{\perp *} = \sum_{e'\in ECP(e)} \tilde{E}_{e,e'} f_{e'} {\bf u}_{e',k}^{*}. \]

Set

(82)\[ \overline{G}_e^{*} = 0, \]

so the barotropic forcing is neither computed nor subtracted from the tendency. The shared update routine then advances the working velocity directly to the midpoint with the half coefficient:

(83)\[ {\bf u}_{e,k}^{n+0.5} = {\bf u}_{e,k}^{n} + \frac{\Delta t}{2} \left( -[f_v]_e {\bf u}_{e,k}^{\perp *} + \Gamma_{e,k}^{*} \right). \]

The full-step velocity \({\bf u}^{n+1}\) is not formed here; it is recovered from the midpoint at the end of the final outer iteration by Eq. (86).

3.3.3 Stage 2: Barotropic velocity advance, explicitly subcycled

For the unsplit time stepper, \(\overline{{\bf u}}=0\). This stage is skipped.

3.3.4 Stage 3: Update tracers and diagnostics

Compute \({\bf u}_{e,k}^{\text{tr}}\):

(84)\[ {\bf u}_{e,k}^{\text{tr}} = {\bf u}_{e,k}^{n+0.5}. \]

There is no barotropic mode to reconcile, so the velocity correction \({\bf u}^{\text{co}}\) of Eq. (58) is identically zero and the transport velocity is simply the midpoint velocity from Stage 1.

From this point Stage 3 is identical to the split case, because it depends on the mode split only through \({\bf u}^{\text{tr}}\). The same routine computes the pseudo-thickness and tracer auxiliary variables, diagnoses \(\tilde{W}_{tr}\), and advances pseudo thickness and tracers, applying the half-step midpoint construction on a non-final outer iteration and retaining the full-step conservative update on the final one. Section 3.2.5 gives the equations.

3.3.5 Reset variables

Only the velocity handling differs from the split case, since the mode split enters Section 3.2.6 solely through the velocity. If iterating, the working full velocity is simply the Stage 1 midpoint,

(85)\[ {\bf u}^{*} = {\bf u}^{n+0.5} \quad \text{from Stage 1}, \]

and after the final iteration the physical full-step velocity is reconstructed from it with no barotropic contribution:

(86)\[ {\bf u}^{n+1}=2{\bf u}^{n+1/2}-{\bf u}^{n}. \]

Everything else follows Section 3.2.6 unchanged: the pseudo-thickness and tracer resets, the column totals \(\tilde H^{*}\) and \(\tilde H^{n+1}\), the working-state and tracer halo exchange between outer iterations, and the closing sequence of time-level rotation, kinetic diagnostics, implicit vertical mixing, state validation, and clock advancement. Because SplitFactor is zero, the routine skips the \(B'\) reset of Eq. (70).

4. Design

The split-explicit time stepping implementation consists of a baroclinic stepper with a barotropic subcycle stage. The same stepper also supports an unsplit mode by setting the split forcing factor to zero and skipping the barotropic subcycling stage.

The implementation is spread across the following files.

File

Contents

timeStepping/SplitExplicitTypes.h

SplitExplicitConfig and SplitExplicitScratch

timeStepping/SplitExplicitInit.{h,cpp}

Config reading, scratch allocation, velocity split and recombination, barotropic pressure initialization

timeStepping/SplitExplicitRK2Stepper.{h,cpp}

Outer iteration, Stage 1, Stage 3, transport velocity, time-level management

timeStepping/SplitExplicitBarotropicPCStepper.{h,cpp}

Stage 2 predictor-corrector barotropic subcycle

timeStepping/TimeStepper.{h,cpp}

New stepper types and the initializeStateFromInput hook

ocn/OceanState.{h,cpp}

The three new split prognostic fields

ocn/Tendencies.{h,cpp}

Baroclinic velocity tendency and standalone Coriolis acceleration

ocn/TendencyTerms.h

CoriolisAccelerationOnEdge and the split-specific operator overloads

ocn/AuxiliaryState.{h,cpp}

computePseudoThicknessTracerAux and transport-velocity overloads

ocn/Eos.{h,cpp}

Depth-integrated and depth-mean specific volume

ocn/VertCoord.{h,cpp}

Column-integrated pseudo and geometric thickness

ocn/OceanInit.cpp

Calls initializeStateFromInput after input is read

4.1 Data types and parameters

The slow TimeStepper and TimeStep options are read from TimeIntegration. Split-specific options are read from its optional ModeSplitShare subgroup. Missing keys retain the defaults in SplitExplicitConfig.

4.1.1 Parameters

  • TimeStepper: SplitExplicitRK2 selects mode-split stepping and UnsplitRK2 selects its unsplit counterpart. These are the only recognized split-stepper names.

  • NTimeStepIteration: Number of outer predictor-corrector iterations. The default is 2 and values must be positive. Every iteration uses the full TimeStep; one iteration reduces the outer scheme to forward Euler.

  • NBclCoriolisIteration: Number of centered Coriolis iterations used in the baroclinic velocity update. The default is 2 for split-explicit stepping and 1 for unsplit stepping.

  • BtrTimeStep: Requested maximum barotropic subcycle time step. It must be positive and is ignored by UnsplitRK2.

  • BtrTimeStepper: Barotropic subcycle algorithm. The initial implementation supports Predictor-Corrector.

  • SplitFactor: Internal factor multiplying the barotropic pressure-anomaly contribution to the baroclinic velocity tendency. It is 1 for SplitExplicitRK2 and 0 for unsplit stepping.

  • ReinitSplitVelocity: If true, recompute barotropic and baroclinic velocity from NormalVelocity at the beginning of every slow step. The default is false.

  • NBtrSubcycles: Internal count computed as max(1, ceil(TimeStep / BtrTimeStep)) when the stepper is constructed, and only for SplitExplicitRK2. TimeStepper::changeTimeStep does not currently recompute this count, so a stepper whose time step is changed after construction keeps its original subcycle count and shrinks \(\Delta t_{\mathrm{bt}}\) proportionally.

If ModeSplitShare or BtrTimeStep is omitted from a user configuration, the requested barotropic step defaults to the slow TimeStep, giving one nominal subcycle and two predictor-corrector passes. configs/Default.yml supplies the whole subgroup, so the fallback applies only when a run configuration replaces rather than extends the defaults.

The ModeSplitShare settings shipped in configs/Default.yml, together with a SplitExplicitRK2 selection, are:

TimeIntegration:
  TimeStepper: SplitExplicitRK2
  TimeStep: 0000_00:10:00
  ModeSplitShare:
    BtrTimeStepper: Predictor-Corrector
    BtrTimeStep: 0000_00:00:20
    NTimeStepIteration: 2
    NBclCoriolisIteration: 2
    ReinitSplitVelocity: false

4.1.2 Class/structs/data types

  • SplitExplicitConfig: Holds split-explicit configuration, including barotropic stepper selection, subcycle count, time-step iteration count, Coriolis iteration count, and SplitFactor.

  • SplitExplicitScratch: Holds temporary arrays for barotropic subcycling and baroclinic updates:

    • NormalBarotropicVelocitySubcycle{Cur,Pre,Cor} and BarotropicPressureAnomalySubcycle{Cur,Pre,Cor}, the three subcycle buffers per field described in Section 3.2.3.

    • BarotropicPressure, the diagnostic \(B\) formed at initialization.

    • BarotropicForcing, the \(\overline G\) produced by Stage 1 and consumed by Stage 2.

    • BarotropicFlux, the time-averaged pressure transport \(\overline F^{\mathrm{bt}}\) produced by Stage 2.

    • BaroclinicPseudoThicknessEdge, the column sum of flux pseudo thickness on edges used by \({\cal P}_e\).

    • BaseVelocityTend, the saved non-Coriolis baroclinic velocity tendency.

    • NormalTransportVelocity, the corrected transport velocity \({\bf u}^{\mathrm{tr}}\) used by Stage 3.

    All arrays are sized over the full halo extent and zeroed at allocation. The scratch struct is a mutable member of the stepper, since doStep is const.

  • SplitExplicitInit: Provides initialization utilities for reading split-explicit options, allocating scratch arrays, splitting full normal velocity into barotropic and baroclinic parts, combining the split velocity, and initializing barotropic pressure and pressure anomaly.

  • SplitExplicitRK2Stepper: Implements the RK2 time step, including the baroclinic velocity stage, optional barotropic subcycle stage, thickness and tracer stage, halo exchanges, and time-level rotation.

  • SplitExplicitBarotropicPCStepper: Implements the predictor-corrector barotropic stage interface. This class owns the details of the barotropic algorithm so that SplitExplicitRK2Stepper can call a generic stage-2 function.

  • TimeStepper: Gains the SplitExplicitRK2 and UnsplitRK2 enumerators, the matching strings in getTimeStepperFromStr, the construction cases in TimeStepper::create, and a virtual initializeStateFromInput whose base implementation is an empty no-op so that the other steppers are unaffected.

  • OceanState: Stores the additional split-explicit prognostic fields: NormalBaroclinicVelocity (edges, layered), NormalBarotropicVelocity (edges, single level), and BarotropicPressureAnomaly (cells, single level). All three are added to the Restart field group but not to the State group, so they are written to and read from restart files without appearing in default history output.

  • Tendencies: Exposes computeCoriolisAccelerationOnEdge in layered and single-layer forms, so that the Coriolis term can be applied on its own inside the Stage-1 iteration and by the barotropic subcycle. Both forms return without acting when the potential-vorticity tendency is disabled.

  • TendencyTerms: Adds a CoriolisAccelerationOnEdge operator with layered and single-layer call operators; an overload of PotentialVortHAdvOnEdge that advects with relative vorticity only, leaving the planetary part to the Coriolis operator.

  • AuxiliaryState: Adds computePseudoThicknessTracerAux, which computes just the pseudo-thickness and tracer auxiliary variables needed by Stage 3, in time-level and explicit-velocity-array forms.

  • Eos: Stores DepthIntegSpecificVolume, the column integral of specific volume weighted by pseudo thickness, and DepthMeanSpecificVolume, that integral divided by the column-integrated pseudo thickness. Both are computed by computeDepthIntegratedSpecificVolume and are registered as Eos fields. DepthMeanSpecificVolume is the \(\overline{\alpha}\) used by both the barotropic momentum equation and the barotropic pressure-anomaly contribution to the baroclinic velocity tendency.

  • VertCoord: Stores TotalPseudoThickness, the column sum \(\tilde H\), and TotalGeometricThickness, the column geometric depth \(H=\rho_0 S\) of Eq. (15).

  • OceanInit: initStateForTimeStepper calls TimeStepper::initializeStateFromInput after the initial-state or restart read and after initUpdateHaloAndHostArrays, then exchanges the state halo and copies to the host.

4.2 Methods

The implemented doStep control flow is shown below. Pacer instrumentation and the communicator used only by its timing barriers are elided.

void SplitExplicitRK2Stepper::doStep(OceanState *State,
                                     TimeInstant &SimTime) const {

   const int CurLevel  = 0;
   const int NextLevel = 1;

   Array3DReal CurTracerArray  = Tracers::getAll(CurLevel);
   Array3DReal NextTracerArray = Tracers::getAll(NextLevel);

   // Initialize NextLevel from CurLevel, optionally recomputing the split
   initializeNextState(State, CurLevel, NextLevel, SEConfig.SplitFactor,
                       SEConfig.ReinitSplitVelocity);
   deepCopy(NextTracerArray, CurTracerArray);

   const TimeInstant StageTime = SimTime;
   for (I4 TimeStepIteration = 0;
        TimeStepIteration < SEConfig.NTimeStepIteration; ++TimeStepIteration) {

      const bool FinalIteration =
          TimeStepIteration + 1 == SEConfig.NTimeStepIteration;

      // The first iteration evaluates the momentum right-hand side at the n
      // state copied into NextLevel; every later iteration sees the midpoint
      // state left by its predecessor, so time-dependent terms are sampled at
      // n+1/2 to keep the predictor-corrector second order.
      const TimeInstant VelStageTime =
          TimeStepIteration == 0 ? StageTime : StageTime + 0.5 * TimeStep;

      // Stage 1: Baroclinic velocity advance, with long time step
      doBaroclinicVelocityUpdate(State, NextTracerArray, CurLevel, NextLevel,
                                 VelStageTime, TimeStep);

      MeshHalo->exchangeFullArrayHalo(
          State->getNormalBaroclinicVelocity(NextLevel), OnEdge);

      if (SEConfig.SplitFactor != 0._Real) {
         // Stage 2: Barotropic velocity advance, explicitly subcycled
         doBarotropicVelocityUpdate(State, CurLevel, NextLevel,
                                    StageTime + 0.5 * TimeStep, TimeStep);
      }

      // Physical total velocity and the corrected transport velocity
      computeTransportVelocity(State, NextLevel);

      // Stage 3: Update thickness, tracers, other diagnostics
      doThicknessTracerUpdate(State, CurTracerArray, NextTracerArray, CurLevel,
                              NextLevel, StageTime, TimeStep, FinalIteration);

      if (!FinalIteration) {
         State->exchangeHalo(NextLevel);
         MeshHalo->exchangeFullArrayHalo(NextTracerArray, OnCell);
      }
   }

   State->updateTimeLevels();
   Tracers::updateTimeLevels();

   // Refresh the kinetic auxiliary variables from the completed n+1 velocity,
   // an inline parallelForOuter over KineticAux::computeVarsOnCell

   // Apply implicit vertical mixing, then re-exchange the affected halos
   VertMix *VMix = VertMix::getInstance();
   if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) {
      VMix->VertMixImplicit(State, AuxState, Tracers::getAll(CurLevel),
                            Tracers::getNumTracers(), CurLevel);
      State->exchangeHalo(CurLevel);
      Tracers::exchangeHalo(CurLevel);
   }

   validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel);

   StepClock->advance();
   SimTime = StepClock->getCurrentTime();
   ++StepCount;
}

The stepper also implements finalizeInit, which validates the tendency, mesh, vertical coordinate, and halo pointers, allocates the scratch arrays, and binds the barotropic stage callback for the configured BtrTimeStepper.

4.2.1 Initialization

After the initial condition or restart state has been read, OceanInit initializes the split-explicit state through TimeStepper::initializeStateFromInput. For SplitExplicitRK2Stepper, this initialization performs the following operations:

  1. Compute momentum vertical auxiliary variables needed for pressure.

  2. For a non-restart SplitExplicitRK2 run, compute barotropic pressure and pressure anomaly from pressure-interface, surface-pressure, and bottom-depth fields, then initialize all three barotropic subcycle pressure buffers.

  3. Initialize velocity split fields. For a non-restart run with SplitExplicitRK2, full normal velocity is split into barotropic and baroclinic components. For a restart run, barotropic velocity is preserved and baroclinic velocity is reconstructed as full minus barotropic velocity. For unsplit stepping, barotropic velocity is set to zero and baroclinic velocity is set equal to full normal velocity.

  4. Initialize time level 1 from time level 0 so that the first time step starts from a consistent pair of time levels.

OceanInit then exchanges the time-level-0 state halo and copies it to the host, so that the newly split velocity fields are consistent across ranks and on the host before any output or time stepping.

At the beginning of every full time step, time level 1 is refreshed from time level 0 for the state fields advanced by the split-explicit scheme, and tracer time level 1 is copied from tracer time level 0. The state at time level 0 is kept fixed during the internal time-step iterations. Time level 1 is updated repeatedly and provides the most recent state for tendency computations.

This refresh is also where ReinitSplitVelocity acts. When it is true and the mode split is active, the velocity split is recomputed from the time-level-0 NormalVelocity, including a new depth mean for NormalBarotropicVelocity. When it is false, the stored barotropic velocity is kept and only the baroclinic velocity is rebuilt as \({\bf u}-\overline{\bf u}\), which preserves the barotropic state carried forward from the previous step.

4.2.2 Stage 1: baroclinic velocity

Stage 1 computes the baroclinic velocity tendency at the current iteration state and advances NormalBaroclinicVelocity by half of the full slow time step. The first outer iteration samples at \(t^n\) and every later iteration samples time-dependent forcing at \(t^{n+1/2}\). The non-Coriolis tendency is saved in SplitExplicitScratch before the centered Coriolis iteration begins.

The centered Coriolis treatment repeats NBclCoriolisIteration times. Each iteration restores the saved non-Coriolis tendency, adds the Coriolis acceleration computed from the updated baroclinic velocity, and updates NormalBaroclinicVelocity again from the fixed time-level-0 base state. After stage 1, only the updated baroclinic velocity halo is exchanged.

4.2.3 Stage 2: barotropic subcycling

Stage 2 advances the barotropic velocity and barotropic pressure anomaly using the configured barotropic stepper. The RK2 stepper delegates this work through a generic barotropic-stage callback so that additional barotropic algorithms can be added without complicating the RK2 stage logic.

For SplitExplicitRK2, the barotropic stage is active and currently uses the predictor-corrector barotropic stepper. For unsplit stepping, SplitFactor is zero and stage 2 is skipped. The predictor-corrector implementation owns three buffers (Cur, Pre, and Cor) for both barotropic velocity and pressure anomaly. It exchanges the current buffers once per pass, runs communication- avoiding kernels over shrinking halo ranges, accumulates corrected velocity and pressure transport, and exchanges the resulting averaged fields after the last pass.

The callback signature takes only the state, the two time levels, the stage time, and the stage time step; the barotropic stepper reaches everything else it needs, including the scratch arrays, the configuration, the mesh, the halo, the vertical coordinate, and the equation of state, through the arguments the RK2 stepper binds into the lambda. Adding an implicit or Adams–Bashforth barotropic algorithm therefore means adding an enumerator, a class implementing the same doBarotropicVelocityUpdate entry point, and a branch in initBarotropicStepper, with no change to the RK2 stage logic.

4.2.4 Stage 3: thickness and tracers

Before Stage 3, computeTransportVelocity reconstructs physical normal velocity from the split fields and adds a depth-uniform correction to a separate NormalTransportVelocity scratch field. Stage 3 computes pseudo-thickness and tracer auxiliary variables with this corrected transport, diagnoses vertical pseudo velocity, and advances pseudo thickness and tracers. These Stage-3 tendencies are sampled at \(t^{n+1/2}\) on every outer iteration.

During intermediate time-step iterations, full normal velocity is reconstructed as the sum of barotropic and baroclinic velocity at time level 1. On the final iteration, full normal velocity is reconstructed at time level \(n+1\) using

\[ {\bf u}^{n+1} = {\overline {\bar {\bf u}}}^{\mathrm{bt}} + 2 {\bf u}'^{n+1/2} - {\bf u}'^{n}. \]

This final reconstruction follows the MPAS-Ocean split-explicit update: the baroclinic velocity stored at time level 1 is the midpoint value, but the full normal velocity must be available at the final time level for output, diagnostics, restart, and the next time step.

4.2.5 Time-level management

The internal time-step iterations do not rotate Ocean state or tracer time levels. Instead, time level 0 remains the base state for RK2 updates and time level 1 is updated in place. Halo exchanges are performed on time level 1 between internal iterations. After the final iteration, updateTimeLevels is called once for the Ocean state and tracers, making the completed time-level-1 state the new time-level-0 state. NormalBaroclinicVelocity, NormalBarotropicVelocity, and BarotropicPressureAnomaly participate in OceanState::exchangeHalo and in restart I/O, but are registered only in the Restart field group. The physical NormalVelocity and PseudoThickness remain in the normal state/history group.

Because the three split fields are in the restart group, a restart carries the barotropic velocity and pressure anomaly forward exactly, and Section 4.2.1 rebuilds only the baroclinic velocity from them. The alternative, recomputing the split from NormalVelocity alone, would discard the barotropic pressure anomaly and break restart reproducibility.

5. Verification and testing

5.1 Unit testing

TimeStepperTest exercises both new time-stepper types with a velocity-decay tendency. It initializes each stepper through initializeStateFromInput, performs time-step refinement, and checks the observed convergence rate. The current expected rates in the implementation are first order for SplitExplicitRK2 (tolerance 0.15) and second order for UnsplitRK2 (tolerance 0.1). This records the behavior of the current code and does not claim second-order convergence for the complete split algorithm.

The same test creates SplitExplicitRK2 and UnsplitRK2 without a stop time, verifies that no end alarm is present, calls doStep twice, and checks that the persistent step count increments across calls. Test setup overrides BtrTimeStep in the in-memory configuration so that the coarsest refinement level runs ten barotropic subcycles per baroclinic step, exercising the subcycle loop rather than the degenerate single-subcycle case. Because changeTimeStep does not recompute NBtrSubcycles, the finer refinement level runs the same ten subcycles at half the barotropic step.

The supporting operators added for the split scheme have their own unit tests:

  • TendencyTermsTest adds testCoriolisAccelerationOnEdge, which checks both the layered and single-layer forms of CoriolisAccelerationOnEdge against an independently computed tangential reconstruction.

  • EosTest adds testDepthIntegratedSpecificVolume, which checks the column integral against a constant specific volume and thickness.

  • VertCoordTest checks TotalPseudoThickness after computePressure and TotalGeometricThickness after computeTotalGeometricThickness against analytic column values.

Additional focused unit tests should cover velocity split and recombination, barotropic-pressure initialization, ceiling-based subcycle selection, centered and upwind anomaly interpolation, the transport-velocity correction, restart preservation of the three split prognostic fields, and rejection of halo widths smaller than three.

5.2 Polaris tests

The following end-to-end verification cases are available in Polaris:

  • Overflow

  • Baroclinic channel

  • Realistic global ocean