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 [ \gamma {\eta}^{n+1} + (1-\gamma) {\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} [ \beta \vec{\bf v}^{n+1} + (1-\beta) \vec{\bf v}^{n}] dr = \epsilon_{fw} ({\mathcal{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)} + (\gamma-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 ({\mathcal{P-E}}) - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} [ \beta \vec{\bf v}^* + (1-\beta) \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 \gamma\beta \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}^{*} - \gamma \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 ({mathcal{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(β(𝒫 − ℰ)^n + 1/2 + (1 − β)(𝒫 − ℰ)^{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).