Ctrl+K

Time Stepping

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.

../_images/timestepping.png

Fig: schematic view of the Croco predictor-corrector hydrostatic kernel#

General structure of the time-stepping:

call prestep3D_thread()    ! Predictor step for 3D momentum and tracers
call step2d_thread()       ! Barotropic mode
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 (the time tendencies).

../_images/timestepping_3Dequations.png

For a given quantity \(q\) with \(q_t = {\rm RHS}(q)\) we can write

\[\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#

../_images/timestepping_veltra_cpl.png

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#

../_images/timestepping_barotropic.png

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#

../_images/timestepping_baro_cpl.png

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\))

../_images/Shortstep.png

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}\)

../_images/filters.png

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\]

  • 3D advection

\[\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

Advection scheme

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\]
Contents