Diffusive Part -- Forward or Backward Scheme

The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. For a tendancy $D_x$, representing a diffusion term or a restoring term to a tracer climatology (when present, see § 5.6), a forward time differencing scheme is used :

$$x^{t+\rdt} = x^{t-\rdt} + 2 \rdt {D_x}^{t-\rdt}$$ (3.3)

This is diffusive in time and conditionally stable. The conditions for stability of second and fourth order horizontal diffusion schemes are [Griffies, 2004]:

$$A^h < \left\{ \begin{aligned}&\frac{e^2}{ 8 \; \rdt } &&\quad \... ...\rdt } &&\quad \text{bilaplacian diffusion} \end{aligned} \right.$$

where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. The linear constraint (3.4) is a necessary condition, but not sufficient. If it is not satisfied, even mildly, then the model soon becomes wildly unstable. The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient.

For the vertical diffusion terms, a forward time differencing scheme can be used, but usually the numerical stability condition imposes a strong constraint on the time step. Two solutions are available in NEMO to overcome the stability constraint: $(a)$ a forward time differencing scheme using a time splitting technique (ln_zdfexp = true) or $(b)$ a backward (or implicit) time differencing scheme (ln_zdfexp = false). In $(a)$, the master time step $\Delta$t is cut into $N$ fractional time steps so that the stability criterion is reduced by a factor of $N$. The computation is performed as follows:

$$\begin{split}& x_\ast ^{t-\rdt} = x^{t-\rdt} & x_\ast ^{t-... ...or $L=1$ to $N$} & x^{t+\rdt} = x_\ast^{t+\rdt} \end{split}$$ (3.5)

with DF a vertical diffusion term. The number of fractional time steps, $N$, is given by setting nn_zdfexp, (namelist parameter). The scheme $(b)$ is unconditionally stable but diffusive. It can be written as follows:

$$x^{t+\rdt} = x^{t-\rdt} + 2 \rdt _x^{t+\rdt}$$ (3.6)

This scheme is rather time consuming since it requires a matrix inversion, but it becomes attractive since a value of 3 or more is needed for N in the forward time differencing scheme. For example, the finite difference approximation of the temperature equation is:

$$\frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\rdt}\equiv +\frac{1}{e_{3t} }\delta _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} \right]$$ (3.7)

where RHS is the right hand side of the equation except for the vertical diffusion term. We rewrite (3.6) as:

$$-c(k+1)\;T^{t+1}(k+1) + d(k)\;T^{t+1}(k) - \;c(k)\;T^{t+1}(k-1) \equiv b(k)$$ (3.8)

where

$$c(k) = A_w^{vT} (k) / e_{3w} (k)$$

$$d(k) = e_{3t} (k) / (2\rdt) + c_k + c_{k+1}$$

$$b(k) = e_{3t} (k) \; \left( T^{t-1}(k) / (2\rdt) + \text{RHS} \right)$$

(3.8) is a linear system of equations with an associated matrix which is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, therefore a special adaptation of the Gauss elimination procedure is used to find the solution (see for example Richtmyer and Morton [1967]).