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:
We set
In the hydrostatic case \(\epsilon_{\rm nh}=0\), allowing us to find \({\eta}^{n+1}\), thus:
and then to compute (CORRECTION_STEP):
Notes:
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 filedata
). 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} )\)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} }\)
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
, namelistPARM01
).