crank-nicol.rst

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 (β, γ) and an explicit part (1 − β, 1 − γ) for the surface pressure gradient (β) and for the barotropic flow divergence (γ). For instance, β = γ = 1 is the previous fully implicit scheme; β = γ = 1/2 is the non damping (energy conserving), unconditionally stable, Crank-Nicolson scheme; (β, γ) = (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, β, γ 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 ustar-backward-free-surface – vn+1-backward-free-surface are modified as follows:

$$\frac{ \vec{\bf v}^{n+1} }{ \Delta t } + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ] + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1} = \frac{ \vec{\bf v}^{n} }{ \Delta t } + \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)} + {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}$$

$$\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t} + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr = \epsilon_{fw} (P-E)$$

We set

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

In the hydrostatic case ϵ_nh = 0, allowing us to find η^n + 1, thus:

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

and then to compute (CORRECTION_STEP <model/src/correction_step.F>):

$$\vec{\bf v}^{n+1} = \vec{\bf v}^{*} - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}$$

Notes:

1. The RHS term of equation eta-n+1-CrankNick corresponds the contribution of fresh water flux (P-E) to the free-surface variations (ϵ_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: ϵ_fw(γ(P − E)^n + 1/2 + (1 − γ)(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:
   • β + γ < 1 : unstable
   • β ≥ 1/2 and γ ≥ 1/2 : stable
   • β + γ ≥ 1 : stable if c_max²(β − 1/2)(γ − 1/2) + 1 ≥ 0 with $c_{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 β_nh (0 < β_nh ≤ 1) of the non-hydrostatic pressure gradient being evaluated at time step n + 1 (backward in time) and the remaining part (1 − β_nh) being evaluated at time step n (forward in time). The run-time parameter implicitNHPress corresponding to the implicit fraction β_nh of the non-hydrostatic pressure is set by default to the implicit fraction β of surface pressure (implicSurfPress), but can also be specified independently (in main parameter file data, namelist PARM01).