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 :varlink:`implicSurfPress`, :varlink:`implicDiv2DFlow`. They are read
from the main parameter file data (namelist PARM01) and are set
by default to 1,1.
Equations :eq:`ustar-backward-free-surface` – :eq:`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}
\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}
In the hydrostatic case \epsilon_{nh}=0, allowing us to find {\eta}^{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 (:filelink:`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:
- The RHS term of equation :eq:`eta-n+1-CrankNick` corresponds the
contribution of fresh water flux ({mathcal{P-E}}) to the free-surface variations
(\epsilon_{fw}=1, :varlink:`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_{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_{max}^2 (\gamma - 1/2)(\beta - 1/2) + 1 \geq 0 with c_{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_{nh}
(0 < \gamma_{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_{nh}) being evaluated at
time step n (forward in time). The run-time parameter
:varlink:`implicitNHPress` corresponding to the implicit fraction
\gamma_{nh} of the non-hydrostatic pressure is set by default
to the implicit fraction \gamma of surface pressure
(:varlink:`implicSurfPress`), but can also be specified independently (in
main parameter file
data, namelistPARM01).