# 4.2. Time Stepping¶

CROCO is discretized in time using a third-order predictor-corrector scheme (referred to as LFAM3) for tracers and baroclinic momentum. It is a split-explicit, free-surface ocean model, where short time steps are used to advance the surface elevation and barotropic momentum, with a much larger time step used for tracers, and baroclinic momentum. The model has a 2-way time-averaging procedure for the barotropic mode, which satisfies the 3D continuity equation. The specially designed 3rd order predictor-corrector time step algorithm is described in Shchepetkin and McWilliams (2005) and is summarized in this subsection.

General structure of the time-stepping:

call prestep3D_thread()    ! Predictor step for 3D momentum and tracers
call step3D_uv_thread()   ! Corrector step for momentum
call step3D_t_thread()    ! Corrector step for tracers


## 4.2.1. 3D momentum and tracers¶

Predictor-corrector approach : Leapfrog (LF) predictor with 3rd-order Adams-Moulton (AM) interpolation (LFAM3 timestepping). This scheme is used to integrate 3D advection, the pressure gradient term, the continuity equation and the Coriolis term which are all contained in the $${\rm RHS}$$ operator.

For a given quantity $$q$$

$\begin{split}\left\{ \begin{array}{rcll} \displaystyle q^{n+1,\star} & = & \displaystyle q^{n-1} + 2 \varDelta t \; {\rm RHS} \left\{ q^n \right\} & \qquad \qquad {\rm (LF)}\\ \displaystyle q^{n+\frac{1}{2}} & = & \displaystyle \frac{5}{12} \; q^{n+1,\star} + \frac{2}{3} \; q^n - \frac{1}{12} \; q^{n-1} & \qquad \qquad {\rm (AM3)} \\ \displaystyle q^{n+1} & = & \displaystyle q^n + \varDelta t \; {\rm RHS} \left\{ q^{n+\frac{1}{2}} \right\}& \qquad \qquad {\rm (corrector)}\\ \end{array} \right.\end{split}$

which can be rewritten in a compact way as used in the Croco code :

\begin{eqnarray*} q^{n+\frac{1}{2}} & = & \left( \frac{1}{2} - \gamma \right) q^{n-1} + \left( \frac{1}{2} + \gamma \right) q^{n} + (1 - \gamma)\varDelta t \; {\rm RHS} \left\{ q^n \right\} \\ q^{n+1} & = & q^n + \varDelta t \; {\rm RHS} \left\{ q^{n+\frac{1}{2}} \right\} \end{eqnarray*}

with $$\gamma = \frac{1}{6}$$.

Physical parameterizations for vertical mixing, rotated diffusion and viscous/diffusion terms are computed once per time-step using an Euler step.

## 4.2.2. Tracers-momentum coupling¶

The numerical integration of internal waves can be studied using the following subsystem of equations

$\begin{split}\left\{ \begin{array}{rcl} \displaystyle \partial_z w + \partial_x u &=& 0 \\ \displaystyle \partial_z p + \rho g &=& 0 \\ \displaystyle \partial_t u + \frac{1}{\rho_0} \partial_x p &=& 0 \\ \displaystyle \partial_t \rho + \partial_z (w \rho) &=& 0 \end{array} \right.\end{split}$

Predictor step:

\begin{eqnarray*} \partial_x p^n = g \; \partial_x \left( \int_{z}^0 \rho^n dz \right) & \qquad \rightarrow \qquad & u^{n+\frac{1}{2}} = \left( \frac{1}{2} - \gamma \right) u^{n-1} + \left( \frac{1}{2} + \gamma \right) u^{n} + (1 - \gamma) \frac{\Delta t}{\rho_0} \; (\partial_x p^{n}) \\ w^n = - \int_{-H}^{z} \partial_x u^n dz' &\qquad \rightarrow \qquad& \rho^{n+\frac{1}{2}} = \left( \frac{1}{2} - \gamma \right) \rho^{n-1} + \left( \frac{1}{2} + \gamma \right) \rho^{n} + (1 - \gamma) \Delta t\; \partial_z (w^{n} \rho^{n}) \end{eqnarray*}

Corrector step:

\begin{eqnarray*} \partial_x p^{n+\frac{1}{2}} = g \; \partial_x \left( \int_{z}^0 \rho^{n+\frac{1}{2}} dz \right) & \qquad \rightarrow \qquad & u^{n+1} = u^{n} + \frac{\Delta t}{\rho_0} \; (\partial_x p^{n+\frac{1}{2}}) \\ w^{n+\frac{1}{2}} = - \int_{-H}^{z} \partial_x \left\{ \frac{3 u^{n+\frac{1}{2}}}{4} + \frac{u^n + u^{n+1}}{8} \right\} dz' & \qquad \rightarrow \qquad & \rho^{n+1} = \rho^{n} + \Delta t\; \partial_z (w^{n+\frac{1}{2}} \rho^{n+\frac{1}{2}}) \end{eqnarray*}

Consequences:

• 3D-momentum integrated before the tracers in the corrector

• 2 evaluations of the pressure gradient per time-step

• 3 evaluations of the continuity equation per time-step

## 4.2.3. Barotropic mode¶

Generalized forward-backward (predictor-corrector)

1. AB3-type extrapolation

\begin{eqnarray*} D^{m+\frac{1}{2}} &=& H + \left( \frac{3}{2} + \beta \right) \zeta^{m} - \left( \frac{1}{2} + 2\beta \right) \zeta^{m-1} + \beta \zeta^{m-2} \\ \overline{u}^{m+\frac{1}{2}} &=& \left( \frac{3}{2} + \beta \right) \overline{u}^{m} - \left( \frac{1}{2} + 2\beta \right) \overline{u}^{m-1} + \beta \overline{u}^{m-2} \end{eqnarray*}
1. Integration of $$\zeta$$

$\zeta^{m+1} = \zeta^{m} - \Delta \tau \; \partial_x( D^{m+\frac{1}{2}} \overline{u}^{m+\frac{1}{2}} )$
1. AM4 interpolation

$\zeta^{\star} = \left(\frac{1}{2}+\gamma+2\varepsilon\right) \zeta^{m+1} + \left(\frac{1}{2}-2\gamma-3\varepsilon\right) \zeta^{m} + \gamma \zeta^{m-1} + \varepsilon \zeta^{m-2}$
1. Integration of $$\bar{u}$$

$\overline{u}^{m+1} = \frac{1}{D^{m+1}} \left[ D^{m} \overline{u}^{m} + \Delta \tau \; {\rm RHS2D}(D^{m+\frac{1}{2}}, \overline{u}^{m+\frac{1}{2}}, \zeta^{\star}) \right]$

where the parameter values are $$(\beta,\gamma,\varepsilon)=(0.281105,0.088,0.013)$$ except when the filter_none option is activated (see below).

## 4.2.4. Baroclinic-barotropic coupling¶

Slow forcing term of the barotrope by the barocline is extrapolated

${\cal F}_{3D}^{n+\frac{1}{2}} = \left\{ \int {\rm rhs}(u,v) dz - {\rm rhs2D}(\bar{u},\bar{v}) \right\}^{n+\frac{1}{2}} = {\rm Extrap}( {\cal F}_{3D}^{n}, {\cal F}_{3D}^{n-1}, {\cal F}_{3D}^{n-2})$

### 4.2.4.1. M2_FILTER_POWER option¶

Barotropic integration from $$n$$ to $$n+M^{*} \Delta \tau$$ ($$M^{*} \le 1.5 M$$)

Because of predictor-corrector integration two barotropic filters are needed

• $$\left< \zeta \right>^{n+1} \rightarrow \;$$ update of the vertical grid

• $$\left< U \right>^{n+1} \rightarrow \;$$ correction of baroclinic velocities at time $$n+1$$

• $$\left<\left< U \right>\right>^{n+\frac{1}{2}} \rightarrow \;$$ correction of baroclinic velocities at time $$n+\frac{1}{2}$$

### 4.2.4.2. M2_FILTER_NONE option¶

Motivation: averaging filters can lead to excessive dissipation in the barotropic mode

Objective: put the minimum amount of dissipation to stabilize the splitting

Diffusion is introduced within the barotropic time-stepping rather than averaging filters by adapting the parameters in the generalized forward-backward scheme

$(\beta,\gamma,\varepsilon)=(0.281105,0.08344500 - 0.51358400 \alpha_{\rm d},0.00976186 - 0.13451357 \alpha_{\rm d})$

with $$\alpha_{\rm d} \approx 0.5$$.

Remarks:

• This option may require to increase $${\rm NDTFAST} = \varDelta t_{\rm 3D} / \varDelta t_{\rm 2D}$$ because the stability constraint of the modified generalized forward-backward scheme is less than the one of the original generalized forward-backward scheme.

• The filter_none approach is systematically more efficient than averaging filters

## 4.2.5. Stability constraints¶

• Barotropic mode (note that considering an Arakawa C-grid divides the theoretical stability limit by a factor of 2)

$\Delta t \sqrt{ g H \left( \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} \right) } \le 0.89$

$\alpha_{\rm adv}^x + \alpha_{\rm adv}^y + \beta \alpha_{\rm adv}^z \le \alpha^\star_{\rm horiz}$

where $$\alpha_{\rm adv}^x$$, $$\alpha_{\rm adv}^y$$, and $$\alpha_{\rm adv}^z$$ are the Courant numbers in each direction and $$\beta = \alpha^\star_{\rm horiz} / \alpha^\star_{\rm vert}$$ a coefficient arising from the fact that different advection schemes with different stability criteria may be used in the horizontal and vertical directions. Typical CFL values for $$\alpha^\star_{\rm horiz}$$ and $$\alpha^\star_{\rm vert}$$ with Croco time-stepping algorithm are

Max Courant number ($$\alpha^\star$$)

C2

1.587

UP3

0.871

SPLINES

0.916

C4

1.15

UP5

0.89

C6

1.00

• Internal waves

$\Delta t c_1 \sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} } \le 0.843686$

where $$c_1$$ the phase speed associated with the first (fastest) baroclinic mode.

• Coriolis

$f \Delta t \le 1.58$