2.10.1. Crank-Nicolson barotropic time stepping

The full implicit time stepping described previously is unconditionally stable but damps the fast gravity waves, resulting in a loss of potential energy. The modification presented now allows one to combine an implicit part (\(\gamma,\beta\)) and an explicit part (\(1-\gamma,1-\beta\)) for the surface pressure gradient (\(\gamma\)) and for the barotropic flow divergence (\(\beta\)). For instance, \(\gamma=\beta=1\) is the previous fully implicit scheme; \(\gamma=\beta=1/2\) is the non damping (energy conserving), unconditionally stable, Crank-Nicolson scheme; \((\gamma,\beta)=(1,0)\) or \(=(0,1)\) corresponds to the forward - backward scheme that conserves energy but is only stable for small time steps. In the code, \(\gamma,\beta\) are defined as parameters, respectively implicSurfPress, implicDiv2DFlow. They are read from the main parameter file data (namelist PARM01) and are set by default to 1,1.

Equations (2.12)(2.17) are modified as follows:

\[\frac{ \vec{\bf v}^{n+1} }{ \Delta t } + {\bf \nabla}_h b_s [ \gamma {\eta}^{n+1} + (1-\gamma) {\eta}^{n} ] + \epsilon_{\rm nh} {\bf \nabla}_h {\phi'_{\rm nh}}^{n+1} = \frac{ \vec{\bf v}^{n} }{ \Delta t } + \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)} + {\bf \nabla}_h {\phi'_{\rm hyd}}^{(n+1/2)}\]
(2.69)\[\epsilon_{\rm fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t} + {\bf \nabla}_h \cdot \int_{R_{\rm fixed}}^{R_o} [ \beta \vec{\bf v}^{n+1} + (1-\beta) \vec{\bf v}^{n}] dr = \epsilon_{\rm fw} ({\mathcal{P-E}})\]

We set

\[\begin{split}\begin{aligned} \vec{\bf v}^* & = & \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)} + (\gamma-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n} + \Delta t {\bf \nabla}_h {\phi'_{\rm hyd}}^{(n+1/2)} \\ {\eta}^* & = & \epsilon_{fs} {\eta}^{n} + \epsilon_{\rm fw} \Delta t ({\mathcal{P-E}}) - \Delta t {\bf \nabla}_h \cdot \int_{R_{\rm fixed}}^{R_o} [ \beta \vec{\bf v}^* + (1-\beta) \vec{\bf v}^{n}] dr\end{aligned}\end{split}\]

In the hydrostatic case \(\epsilon_{\rm nh}=0\), allowing us to find \({\eta}^{n+1}\), thus:

\[\epsilon_{\rm fs} {\eta}^{n+1} - {\bf \nabla}_h \cdot \gamma\beta \Delta t^2 b_s (R_o - R_{\rm fixed}) {\bf \nabla}_h {\eta}^{n+1} = {\eta}^*\]

and then to compute (CORRECTION_STEP):

\[\vec{\bf v}^{n+1} = \vec{\bf v}^{*} - \gamma \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}\]

Notes:

  1. The RHS term of equation (2.69) corresponds the contribution of fresh water flux ({mathcal{P-E}}) to the free-surface variations (\(\epsilon_{fw}=1\), useRealFreshWaterFlux =.TRUE. in parameter file data). In order to remain consistent with the tracer equation, specially in the non-linear free-surface formulation, this term is also affected by the Crank-Nicolson time stepping. The RHS reads: \(\epsilon_{\rm fw} ( \beta ({\mathcal{P-E}})^{n+1/2} + (1-\beta) ({\mathcal{P-E}})^{n-1/2} )\)

  2. The stability criteria with Crank-Nicolson time stepping for the pure linear gravity wave problem in cartesian coordinates is:

    • \(\gamma + \beta < 1\) : unstable

    • \(\gamma \geq 1/2\) and \(\beta \geq 1/2\) : stable

    • \(\gamma + \beta \geq 1\) : stable if \(c_{\rm max}^2 (\gamma - 1/2)(\beta - 1/2) + 1 \geq 0\) with \(c_{\rm max} = 2 \Delta t \sqrt{gH} \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }\)

  3. A similar mixed forward/backward time-stepping is also available for the non-hydrostatic algorithm, with a fraction \(\gamma_{\rm nh}\) (\(0 < \gamma_{\rm nh} \leq 1\)) of the non-hydrostatic pressure gradient being evaluated at time step \(n+1\) (backward in time) and the remaining part (\(1 - \gamma_{\rm nh}\)) being evaluated at time step \(n\) (forward in time). The run-time parameter implicitNHPress corresponding to the implicit fraction \(\gamma_{\rm nh}\) of the non-hydrostatic pressure is set by default to the implicit fraction \(\gamma\) of surface pressure (implicSurfPress), but can also be specified independently (in main parameter file data, namelist PARM01).