s_algorithm/text/time_stepping.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/time_stepping.tex,v 1.3 2001/08/09 20:06:18 jmc Exp $
% $Name:  $

The convention used in this section is as follows:
Time is "discretize" using a time step $\Delta t$   
and $\Phi^n$ refers to the variable $\Phi$ 
at time $t = n \Delta t$ . We used the notation $\Phi^{(n)}$
when time interpolation is required to estimate the value of $\phi$
at the time $n \Delta t$.

\section{Time integration}

The discretization in time of the model equations (cf section I )
does not depend of the discretization in space of each
term, so that this section can be read independently.
Also for this purpose, we will refers to the continuous 
space-derivative form of model equations, and focus on 
the time discretization.
 
The continuous form of the model equations is:

\begin{eqnarray}
\partial_t \theta & = & G_\theta
\label{eq-tCsC-theta}
\\
\partial_t S & = & G_s
\label{eq-tCsC-salt}
\\
b' & = & b'(\theta,S,r)
\\
\partial_r \phi'_{hyd} & = & -b'
\label{eq-tCsC-hyd}
\\
\partial_t \vec{\bf v}
+ {\bf \nabla}_h b_s \eta
+ \epsilon_{nh} {\bf \nabla}_h \phi'_{nh}
& = & \vec{\bf G}_{\vec{\bf v}} 
- {\bf \nabla}_h \phi'_{hyd}
\label{eq-tCsC-Hmom}
\\
\epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}
+ \epsilon_{nh} \partial_r \phi'_{nh}
& = & \epsilon_{nh} G_{\dot{r}} 
\label{eq-tCsC-Vmom}
\\
{\bf \nabla}_h \cdot \vec{\bf v} + \partial_r \dot{r}
& = & 0
\label{eq-tCsC-cont}
\end{eqnarray}
where
\begin{eqnarray*}
G_\theta & = &
- \vec{\bf v} \cdot {\bf \nabla} \theta + {\cal Q}_\theta
\\
G_S & = &
- \vec{\bf v} \cdot {\bf \nabla} S + {\cal Q}_S
\\
\vec{\bf G}_{\vec{\bf v}}
& = &
- \vec{\bf v} \cdot {\bf \nabla} \vec{\bf v}
- f \hat{\bf k} \wedge \vec{\bf v}
+ \vec{\cal F}_{\vec{\bf v}}
\\
G_{\dot{r}}
& = &
- \vec{\bf v} \cdot {\bf \nabla} \dot{r}
+ {\cal F}_{\dot{r}}
\end{eqnarray*}
The exact form of all the "{\it G}"s terms is described in the next
section (). Here its sufficient to mention that they contains
all the RHS terms except the pressure / geo- potential terms.

The switch $\epsilon_{nh}$ allows to activate the non hydrostatic
mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise, 
in the hydrostatic limit $\epsilon_{nh} = 0$ 
and equation \ref{eq-tCsC-Vmom} vanishes.

The equation for $\eta$ is obtained by integrating the 
continuity equation over the entire depth of the fluid, 
from $R_{min}(x,y)$ up to $R_o(x,y)$ 
(Linear free surface):
\begin{eqnarray}
\epsilon_{fs} \partial_t \eta =
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
- {\bf \nabla} \cdot \int_{R_{min}}^{R_o} \vec{\bf v} dr
+ \epsilon_{fw} (P-E)
\label{eq-tCsC-eta}
\end{eqnarray}

Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to 
distinguish between a free-surface equation ($\epsilon_{fs}=1$) 
or the rigid-lid approximation ($\epsilon_{fs}=0$);  
and to distinguish when exchange of Fresh-Water is included 
at the ocean surface (natural BC) ($\epsilon_{fw} = 1$) 
or not ($\epsilon_{fw} = 0$).

The hydrostatic potential is found by
integrating \ref{eq-tCsC-hyd} with the boundary condition that
$\phi'_{hyd}(r=R_o) = 0$:
\begin{eqnarray*}
& &
\int_{r'}^{R_o} \partial_r \phi'_{hyd} dr =
\left[ \phi'_{hyd} \right]_{r'}^{R_o} =
\int_{r'}^{R_o} - b' dr
\\
\Rightarrow & &
\phi'_{hyd}(x,y,r') = \int_{r'}^{R_o} b' dr
\end{eqnarray*}

\subsection{General method}
 
An overview of the general method is presented hereafter, 
with explicit references to the Fortran code. This part
often refers to the discretized equations of the model
that are detailed in the following sections.

The general algorithm consist in  a "predictor step" that computes
the forward tendencies ("G" terms") and all 
the "first guess" values (star notation): 
$\theta^*, S^*, \vec{\bf v}^*$ (and $\dot{r}^*$
in non-hydrostatic mode). This is done in the two routines 
{\it THERMODYNAMICS} and {\it DYNAMICS}.

Then the implicit terms that appear on the left hand side (LHS)
of equations \ref{eq-tDsC-theta} - \ref{eq-tDsC-cont},
are solved as follows:
Since implicit vertical diffusion and viscosity terms 
are independent from the barotropic flow adjustment,
they are computed first, solving a 3 diagonal Nr x Nr linear system, 
and incorporated at the end of the {\it THERMODYNAMICS} and 
{\it DYNAMICS} routines.
Then the surface pressure and non hydrostatic pressure
are evaluated ({\it SOLVE\_FOR\_PRESSURE}); 

Finally, within a "corrector step', 
(routine {\it THE\_CORRECTION\_STEP})
the new values of $u,v,w,\theta,S$ 
are derived according to the above equations
(see details in II.1.3). 

At this point, the regular time step is over, but  
the correction step contains also other optional
adjustments such as convective adjustment algorithm, or filters 
(zonal FFT filter, shapiro filter)
that applied on both momentum and tracer fields.
just before setting the $n+1$ new time step value.

Since the pressure solver precision is of the order of
the "target residual" that could be lower than the 
the computer truncation error, and also because some filters 
might alter the divergence part of the flow field,
a final evaluation of the surface r anomaly $\eta^{n+1}$
is performed, according to \ref{eq-tDsC-eta} ({\it CALC\_EXACT\_ETA}).
This ensures a perfect volume conservation.
Note that there is no need for an equivalent Non-hydrostatic
"exact conservation" step, since W is already computed after 
the filters applied.

Regarding optional forcing terms (usually part of a "package"), 
that account for a specific source or sink term (e.g.: condensation
as a sink of water vapor Q), they are generally incorporated 
in the main algorithm as follows;
At the the beginning of the time step,
the additional tendencies are computed
as function of the present state (time step $n$) and external forcing ;
Then within the main part of model,
only those new tendencies are added to the model variables.

[more details needed]\\

The atmospheric physics follows this general scheme.

[more about C\_grid, A\_grid conversion \& drag term]\\

\subsection{Standard synchronous time stepping}

In the standard formulation, the surface pressure is 
evaluated at time step n+1 (implicit method).
The above set of equations is then discretized in time 
as follows:
\begin{eqnarray}
\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
\theta^{n+1} & = & \theta^*
\label{eq-tDsC-theta}
\\
\left[ 1 - \partial_r \kappa_v^S \partial_r \right]
S^{n+1} & = & S^*
\label{eq-tDsC-salt}
\\
%{b'}^{n} & = & b'(\theta^{n},S^{n},r)
%\partial_r {\phi'_{hyd}}^{n} & = & {-b'}^{n}
%\\
{\phi'_{hyd}}^{n} & = & \int_{r'}^{R_o} b'(\theta^{n},S^{n},r) dr
\label{eq-tDsC-hyd}
\\
\vec{\bf v} ^{n+1}
+ \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
+ \epsilon_{nh} \Delta t {\bf \nabla} {\phi'_{nh}}^{n+1}
- \partial_r A_v \partial_r \vec{\bf v}^{n+1}
& = &
\vec{\bf v}^*
\label{eq-tDsC-Hmom}
\\
\epsilon_{fs} {\eta}^{n+1} + \Delta t
{\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
& = & 
    \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n} 
\nonumber
\\
% = \epsilon_{fs} {\eta}^{n} & + & \epsilon_{fw} \Delta_t (P-E)^{n} 
\label{eq-tDsC-eta}
\\
\epsilon_{nh} \left( \dot{r} ^{n+1}
+ \Delta t \partial_r {\phi'_{nh}} ^{n+1}
\right)
& = & \epsilon_{nh} \dot{r}^*
\label{eq-tDsC-Vmom}
\\
{\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
& = & 0
\label{eq-tDsC-cont}
\end{eqnarray}
where
\begin{eqnarray}
\theta^* & = &
\theta ^{n} + \Delta t G_{\theta} ^{(n+1/2)}
\\
S^* & = &
S ^{n} + \Delta t G_{S} ^{(n+1/2)}
\\
\vec{\bf v}^* & = &
\vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
+ \Delta t  {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
\\
\dot{r}^* & = &
\dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
\end{eqnarray}

Note that implicit vertical terms (viscosity and diffusivity) are 
not considered as part of the "{\it G}" terms, but are 
written separately here.

To ensure a second order time discretization for both 
momentum and tracer,
The "{\it G}" terms are "extrapolated" forward in time
(Adams Bashforth time stepping)
from the values computed at time step $n$ and $n-1$
to the time $n+1/2$:
$$G^{(n+1/2)} = G^n + (1/2+\epsilon_{AB}) (G^n - G^{n-1})$$
A small number for the parameter $\epsilon_{AB}$ is generally used 
to stabilize this time stepping.

In the standard non-stagger formulation, 
the Adams-Bashforth time stepping is also applied to the 
hydrostatic (pressure / geo-) potential term $\nabla_h \Phi'_{hyd}$. 
Note that presently, this term is in fact incorporated to the 
$\vec{\bf G}_{\vec{\bf v}}$ arrays ({\bf gU,gV}).

\subsection{Stagger baroclinic time stepping}

An alternative is to evaluate $\phi'_{hyd}$ with the 
new tracer fields, and step forward the momentum after.
This option is known as stagger baroclinic time stepping, 
since tracer and momentum are step forward in time one after the other.
It can be activated turning on a running flag parameter 
{\bf staggerTimeStep} in file "{\it data}"). 

The main advantage of this time stepping compared to a synchronous one,
is a better stability, specially regarding internal gravity waves,
and a very natural implementation of a 2nd order in time 
hydrostatic pressure / geo- potential term.
In the other hand, a synchronous time step might be  better
for convection problems; Its also make simpler time dependent forcing
and diagnostic implementation ; and allows a more efficient threading.

Although the stagger time step does not affect deeply the 
structure of the code --- a switch allows to evaluate the 
hydrostatic pressure / geo- potential from new $\theta,S$ 
instead of the Adams-Bashforth estimation ---
this affect the way the time discretization is presented :

\begin{eqnarray*}
\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
\theta^{n+1/2} & = & \theta^*
\\
\left[ 1 - \partial_r \kappa_v^S \partial_r \right]
S^{n+1/2} & = & S^*
\end{eqnarray*}
with
\begin{eqnarray*}
\theta^* & = &
\theta ^{(n-1/2)} + \Delta t G_{\theta} ^{(n)}
\\
S^* & = &
S ^{(n-1/2)} + \Delta t G_{S} ^{(n)}
\end{eqnarray*}
And 
\begin{eqnarray*}
%{b'}^{n+1/2} & = & b'(\theta^{n+1/2},S^{n+1/2},r)
%\\
%\partial_r {\phi'_{hyd}}^{n+1/2} & = & {-b'}^{n+1/2}
{\phi'_{hyd}}^{n+1/2} & = & \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r) dr
%\label{eq-tDsC-hyd}
\\
\vec{\bf v} ^{n+1}
+ \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
+ \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
- \partial_r A_v \partial_r \vec{\bf v}^{n+1}
& = &
\vec{\bf v}^*
%\label{eq-tDsC-Hmom}
\\
\epsilon_{fs} {\eta}^{n+1} + \Delta t
{\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
& = & 
\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n} 
\\
\epsilon_{nh} \left( \dot{r} ^{n+1}
+ \Delta t \partial_r {\phi'_{nh}} ^{n+1}
\right)
& = & \epsilon_{nh} \dot{r}^*
%\label{eq-tDsC-Vmom}
\\
{\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
& = & 0
%\label{eq-tDsC-cont}
\end{eqnarray*}
with
\begin{eqnarray*}
\vec{\bf v}^* & = &
\vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
+ \Delta t  {\bf \nabla}_h {\phi'_{hyd}}^{n+1/2}
\\
\dot{r}^* & = &
\dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
\end{eqnarray*}

%---------------------------------------------------------------------

\subsection{Surface pressure}

Substituting \ref{eq-tDsC-Hmom} into \ref{eq-tDsC-cont}, assuming
$\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$:
\begin{eqnarray}
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{min})
{\bf \nabla}_h b_s {\eta}^{n+1}
= {\eta}^*
\label{eq-solve2D}
\end{eqnarray}
where
\begin{eqnarray}
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
\Delta t {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^* dr
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} 
\label{eq-solve2D_rhs}
\end{eqnarray}

Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-tDsC-Hmom}
would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic
($\epsilon_{nh}=0$):
$$
\vec{\bf v}^{n+1} = \vec{\bf v}^{*}
- \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
$$

This is known as the correction step. However, when the model is
non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
additional equation for $\phi'_{nh}$. This is obtained by
substituting \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
\ref{eq-tDsC-cont}:
\begin{equation}
\left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
= \frac{1}{\Delta t} \left(
{\bf \nabla}_h \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)
\end{equation}
where
\begin{displaymath}
\vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
\end{displaymath}
Note that $\eta^{n+1}$ is also used to update the second RHS term
$\partial_r \dot{r}^* $ since
the vertical velocity at the surface ($\dot{r}_{surf}$) 
is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.

Finally, the horizontal velocities at the new time level are found by:
\begin{equation}
\vec{\bf v}^{n+1} = \vec{\bf v}^{**}
- \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
\end{equation}
and the vertical velocity is found by integrating the continuity
equation vertically.
Note that, for convenience regarding the restart procedure,
the integration of the continuity equation has been 
moved at the beginning of the time step (instead of at the end),
without any consequence on the solution.

Regarding the implementation, all those computation are done
within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent 
{\it CALL}s.
The standard method to solve the 2D elliptic problem (\ref{eq-solve2D})
uses the conjugate gradient method (routine {\it CG2D}); The 
the solver matrix and conjugate gradient operator are only function
of the discretized domain and are therefore evaluated separately,
before the time iteration loop, within {\it INI\_CG2D}. 
The computation of the RHS $\eta^*$ is partly 
done in {\it CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}.

The same method is applied for the non hydrostatic part, using
a conjugate gradient 3D solver ({\it CG3D}) that is initialized 
in {\it INI\_CG3D}. The RHS terms of 2D and 3D problems 
are computed together, within the same part of the code.

\newpage
%-----------------------------------------------------------------------------------
\subsection{Crank-Nickelson barotropic time stepping}

The full implicit time stepping described previously is unconditionally stable
but damps the fast gravity waves, resulting in a loss of 
gravity potential energy.
The modification presented hereafter allows to combine an implicit part
($\beta,\gamma$) and an explicit part ($1-\beta,1-\gamma$) for the surface
pressure gradient ($\beta$) and for the barotropic flow divergence ($\gamma$).
\\
For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
$\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
stable, Crank-Nickelson scheme; $(\beta,\gamma)=(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, $\beta,\gamma$ are defined as parameters, respectively 
{\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
the main data file "{\it data}" and are set by default to 1,1.

Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} 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}^* }{ \Delta t }
$$
$$
\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
+ {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} 
[ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
= \epsilon_{fw} (P-E)
$$
where:
\begin{eqnarray*}
\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_{min}}^{R_o} 
[ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
\end{eqnarray*}
\\
In the hydrostatic case ($\epsilon_{nh}=0$),
this allow to find ${\eta}^{n+1}$, according to:
$$
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{min})
{\bf \nabla}_h {\eta}^{n+1}
= {\eta}^*
$$ 
and then to compute (correction step):
$$
\vec{\bf v}^{n+1} = \vec{\bf v}^{*}
- \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
$$

The non-hydrostatic part is solved as described previously. 
\\ \\
N.B:
\\
 a) The non-hydrostatic part of the code has not yet been 
updated, %since it falls out of the purpose of this test,
so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
\\
b) To remind, the stability criteria with the Crank-Nickelson time stepping
for the pure linear gravity wave problem in cartesian coordinate is:
\\
$\star$~ $\beta + \gamma < 1$ : unstable
\\
$\star$~ $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
\\
$\star$~ $\beta + \gamma \geq 1$ : stable if
%, for all wave length $(k\Delta x,l\Delta y)$
$$ 
c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
$$
$$
\mbox{with }~
%c^2 = 2 g H {\Delta t}^2 
%(\frac{1-cos 2 \pi / k}{\Delta x^2} 
%+\frac{1-cos 2 \pi / l}{\Delta y^2})
%$$
%Practically, the most stringent condition is obtained with $k=l=2$ :
%$$
c_{max} =  2 \Delta t \: \sqrt{g H} \: 
\sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
$$
1	jmc	1.4	% $Header: /u/gcmpack/mitgcmdoc/part2/time_stepping.tex,v 1.3 2001/08/09 20:06:18 jmc Exp $
2	jmc	1.2	% $Name: $
3	adcroft	1.1
4			The convention used in this section is as follows:
5			Time is "discretize" using a time step $\Delta t$
6			and $\Phi^n$ refers to the variable $\Phi$
7			at time $t = n \Delta t$ . We used the notation $\Phi^{(n)}$
8			when time interpolation is required to estimate the value of $\phi$
9			at the time $n \Delta t$.
10
11	jmc	1.3	\section{Time integration}
12	adcroft	1.1
13			The discretization in time of the model equations (cf section I )
14			does not depend of the discretization in space of each
15			term, so that this section can be read independently.
16			Also for this purpose, we will refers to the continuous
17			space-derivative form of model equations, and focus on
18			the time discretization.
19
20			The continuous form of the model equations is:
21
22			\begin{eqnarray}
23			\partial_t \theta & = & G_\theta
24	jmc	1.4	\label{eq-tCsC-theta}
25	adcroft	1.1	\\
26			\partial_t S & = & G_s
27	jmc	1.4	\label{eq-tCsC-salt}
28	adcroft	1.1	\\
29			b' & = & b'(\theta,S,r)
30			\\
31			\partial_r \phi'_{hyd} & = & -b'
32	jmc	1.4	\label{eq-tCsC-hyd}
33	adcroft	1.1	\\
34			\partial_t \vec{\bf v}
35	jmc	1.3	+ {\bf \nabla}_h b_s \eta
36			+ \epsilon_{nh} {\bf \nabla}_h \phi'_{nh}
37	adcroft	1.1	& = & \vec{\bf G}_{\vec{\bf v}}
38	jmc	1.3	- {\bf \nabla}_h \phi'_{hyd}
39	jmc	1.4	\label{eq-tCsC-Hmom}
40	adcroft	1.1	\\
41			\epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}
42			+ \epsilon_{nh} \partial_r \phi'_{nh}
43			& = & \epsilon_{nh} G_{\dot{r}}
44	jmc	1.4	\label{eq-tCsC-Vmom}
45	adcroft	1.1	\\
46	jmc	1.3	{\bf \nabla}_h \cdot \vec{\bf v} + \partial_r \dot{r}
47	adcroft	1.1	& = & 0
48	jmc	1.4	\label{eq-tCsC-cont}
49	adcroft	1.1	\end{eqnarray}
50			where
51			\begin{eqnarray*}
52			G_\theta & = &
53	jmc	1.3	- \vec{\bf v} \cdot {\bf \nabla} \theta + {\cal Q}_\theta
54	adcroft	1.1	\\
55			G_S & = &
56	jmc	1.3	- \vec{\bf v} \cdot {\bf \nabla} S + {\cal Q}_S
57	adcroft	1.1	\\
58			\vec{\bf G}_{\vec{\bf v}}
59			& = &
60	jmc	1.3	- \vec{\bf v} \cdot {\bf \nabla} \vec{\bf v}
61	adcroft	1.1	- f \hat{\bf k} \wedge \vec{\bf v}
62			+ \vec{\cal F}_{\vec{\bf v}}
63			\\
64			G_{\dot{r}}
65			& = &
66	jmc	1.3	- \vec{\bf v} \cdot {\bf \nabla} \dot{r}
67	adcroft	1.1	+ {\cal F}_{\dot{r}}
68			\end{eqnarray*}
69			The exact form of all the "{\it G}"s terms is described in the next
70			section (). Here its sufficient to mention that they contains
71			all the RHS terms except the pressure / geo- potential terms.
72
73			The switch $\epsilon_{nh}$ allows to activate the non hydrostatic
74			mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise,
75	jmc	1.3	in the hydrostatic limit $\epsilon_{nh} = 0$
76	jmc	1.4	and equation \ref{eq-tCsC-Vmom} vanishes.
77	adcroft	1.1
78			The equation for $\eta$ is obtained by integrating the
79			continuity equation over the entire depth of the fluid,
80			from $R_{min}(x,y)$ up to $R_o(x,y)$
81			(Linear free surface):
82	jmc	1.2	\begin{eqnarray}
83	adcroft	1.1	\epsilon_{fs} \partial_t \eta =
84			\left. \dot{r} \right\|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
85	jmc	1.3	- {\bf \nabla} \cdot \int_{R_{min}}^{R_o} \vec{\bf v} dr
86	adcroft	1.1	+ \epsilon_{fw} (P-E)
87	jmc	1.4	\label{eq-tCsC-eta}
88	jmc	1.2	\end{eqnarray}
89	adcroft	1.1
90			Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to
91			distinguish between a free-surface equation ($\epsilon_{fs}=1$)
92			or the rigid-lid approximation ($\epsilon_{fs}=0$);
93			and to distinguish when exchange of Fresh-Water is included
94			at the ocean surface (natural BC) ($\epsilon_{fw} = 1$)
95			or not ($\epsilon_{fw} = 0$).
96
97			The hydrostatic potential is found by
98	jmc	1.4	integrating \ref{eq-tCsC-hyd} with the boundary condition that
99	adcroft	1.1	$\phi'_{hyd}(r=R_o) = 0$:
100			\begin{eqnarray*}
101			& &
102			\int_{r'}^{R_o} \partial_r \phi'_{hyd} dr =
103			\left[ \phi'_{hyd} \right]_{r'}^{R_o} =
104			\int_{r'}^{R_o} - b' dr
105			\\
106			\Rightarrow & &
107			\phi'_{hyd}(x,y,r') = \int_{r'}^{R_o} b' dr
108			\end{eqnarray*}
109
110	jmc	1.3	\subsection{General method}
111
112	jmc	1.4	An overview of the general method is presented hereafter,
113			with explicit references to the Fortran code. This part
114			often refers to the discretized equations of the model
115			that are detailed in the following sections.
116
117	jmc	1.3	The general algorithm consist in a "predictor step" that computes
118			the forward tendencies ("G" terms") and all
119	jmc	1.4	the "first guess" values (star notation):
120			$\theta^, S^, \vec{\bf v}^$ (and $\dot{r}^$
121			in non-hydrostatic mode). This is done in the two routines
122			{\it THERMODYNAMICS} and {\it DYNAMICS}.
123	jmc	1.3
124	jmc	1.4	Then the implicit terms that appear on the left hand side (LHS)
125			of equations \ref{eq-tDsC-theta} - \ref{eq-tDsC-cont},
126	jmc	1.3	are solved as follows:
127			Since implicit vertical diffusion and viscosity terms
128			are independent from the barotropic flow adjustment,
129			they are computed first, solving a 3 diagonal Nr x Nr linear system,
130	jmc	1.4	and incorporated at the end of the {\it THERMODYNAMICS} and
131			{\it DYNAMICS} routines.
132	jmc	1.3	Then the surface pressure and non hydrostatic pressure
133			are evaluated ({\it SOLVE\_FOR\_PRESSURE});
134
135			Finally, within a "corrector step',
136			(routine {\it THE\_CORRECTION\_STEP})
137			the new values of $u,v,w,\theta,S$
138			are derived according to the above equations
139			(see details in II.1.3).
140
141			At this point, the regular time step is over, but
142			the correction step contains also other optional
143			adjustments such as convective adjustment algorithm, or filters
144			(zonal FFT filter, shapiro filter)
145			that applied on both momentum and tracer fields.
146			just before setting the $n+1$ new time step value.
147
148			Since the pressure solver precision is of the order of
149			the "target residual" that could be lower than the
150			the computer truncation error, and also because some filters
151			might alter the divergence part of the flow field,
152			a final evaluation of the surface r anomaly $\eta^{n+1}$
153	jmc	1.4	is performed, according to \ref{eq-tDsC-eta} ({\it CALC\_EXACT\_ETA}).
154	jmc	1.3	This ensures a perfect volume conservation.
155			Note that there is no need for an equivalent Non-hydrostatic
156			"exact conservation" step, since W is already computed after
157			the filters applied.
158
159			Regarding optional forcing terms (usually part of a "package"),
160	jmc	1.4	that account for a specific source or sink term (e.g.: condensation
161	jmc	1.3	as a sink of water vapor Q), they are generally incorporated
162			in the main algorithm as follows;
163			At the the beginning of the time step,
164	jmc	1.4	the additional tendencies are computed
165	jmc	1.3	as function of the present state (time step $n$) and external forcing ;
166			Then within the main part of model,
167			only those new tendencies are added to the model variables.
168
169			[more details needed]\\
170	jmc	1.4
171	jmc	1.3	The atmospheric physics follows this general scheme.
172
173	jmc	1.4	[more about C\_grid, A\_grid conversion \& drag term]\\
174
175	jmc	1.3	\subsection{Standard synchronous time stepping}
176	adcroft	1.1
177			In the standard formulation, the surface pressure is
178			evaluated at time step n+1 (implicit method).
179			The above set of equations is then discretized in time
180			as follows:
181			\begin{eqnarray}
182			\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
183			\theta^{n+1} & = & \theta^*
184	jmc	1.4	\label{eq-tDsC-theta}
185	adcroft	1.1	\\
186			\left[ 1 - \partial_r \kappa_v^S \partial_r \right]
187			S^{n+1} & = & S^*
188	jmc	1.4	\label{eq-tDsC-salt}
189	adcroft	1.1	\\
190			%{b'}^{n} & = & b'(\theta^{n},S^{n},r)
191			%\partial_r {\phi'_{hyd}}^{n} & = & {-b'}^{n}
192			%\\
193			{\phi'_{hyd}}^{n} & = & \int_{r'}^{R_o} b'(\theta^{n},S^{n},r) dr
194	jmc	1.4	\label{eq-tDsC-hyd}
195	adcroft	1.1	\\
196			\vec{\bf v} ^{n+1}
197	jmc	1.3	+ \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
198			+ \epsilon_{nh} \Delta t {\bf \nabla} {\phi'_{nh}}^{n+1}
199	adcroft	1.1	- \partial_r A_v \partial_r \vec{\bf v}^{n+1}
200			& = &
201			\vec{\bf v}^*
202	jmc	1.4	\label{eq-tDsC-Hmom}
203	adcroft	1.1	\\
204			\epsilon_{fs} {\eta}^{n+1} + \Delta t
205	jmc	1.3	{\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
206	adcroft	1.1	& = &
207			\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
208			\nonumber
209			\\
210			% = \epsilon_{fs} {\eta}^{n} & + & \epsilon_{fw} \Delta_t (P-E)^{n}
211	jmc	1.4	\label{eq-tDsC-eta}
212	adcroft	1.1	\\
213			\epsilon_{nh} \left( \dot{r} ^{n+1}
214			+ \Delta t \partial_r {\phi'_{nh}} ^{n+1}
215			\right)
216			& = & \epsilon_{nh} \dot{r}^*
217	jmc	1.4	\label{eq-tDsC-Vmom}
218	adcroft	1.1	\\
219	jmc	1.3	{\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
220	adcroft	1.1	& = & 0
221	jmc	1.4	\label{eq-tDsC-cont}
222	adcroft	1.1	\end{eqnarray}
223			where
224			\begin{eqnarray}
225			\theta^* & = &
226			\theta ^{n} + \Delta t G_{\theta} ^{(n+1/2)}
227			\\
228			S^* & = &
229			S ^{n} + \Delta t G_{S} ^{(n+1/2)}
230			\\
231			\vec{\bf v}^* & = &
232			\vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
233	jmc	1.3	+ \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
234	adcroft	1.1	\\
235			\dot{r}^* & = &
236			\dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
237			\end{eqnarray}
238
239			Note that implicit vertical terms (viscosity and diffusivity) are
240			not considered as part of the "{\it G}" terms, but are
241			written separately here.
242
243			To ensure a second order time discretization for both
244			momentum and tracer,
245	jmc	1.4	The "{\it G}" terms are "extrapolated" forward in time
246	adcroft	1.1	(Adams Bashforth time stepping)
247			from the values computed at time step $n$ and $n-1$
248			to the time $n+1/2$:
249			$$G^{(n+1/2)} = G^n + (1/2+\epsilon_{AB}) (G^n - G^{n-1})$$
250			A small number for the parameter $\epsilon_{AB}$ is generally used
251			to stabilize this time stepping.
252
253			In the standard non-stagger formulation,
254			the Adams-Bashforth time stepping is also applied to the
255			hydrostatic (pressure / geo-) potential term $\nabla_h \Phi'_{hyd}$.
256			Note that presently, this term is in fact incorporated to the
257	jmc	1.3	$\vec{\bf G}_{\vec{\bf v}}$ arrays ({\bf gU,gV}).
258	adcroft	1.1
259	jmc	1.3	\subsection{Stagger baroclinic time stepping}
260	adcroft	1.1
261			An alternative is to evaluate $\phi'_{hyd}$ with the
262			new tracer fields, and step forward the momentum after.
263			This option is known as stagger baroclinic time stepping,
264			since tracer and momentum are step forward in time one after the other.
265			It can be activated turning on a running flag parameter
266	jmc	1.3	{\bf staggerTimeStep} in file "{\it data}").
267	adcroft	1.1
268			The main advantage of this time stepping compared to a synchronous one,
269			is a better stability, specially regarding internal gravity waves,
270			and a very natural implementation of a 2nd order in time
271			hydrostatic pressure / geo- potential term.
272			In the other hand, a synchronous time step might be better
273			for convection problems; Its also make simpler time dependent forcing
274			and diagnostic implementation ; and allows a more efficient threading.
275
276			Although the stagger time step does not affect deeply the
277			structure of the code --- a switch allows to evaluate the
278			hydrostatic pressure / geo- potential from new $\theta,S$
279			instead of the Adams-Bashforth estimation ---
280			this affect the way the time discretization is presented :
281
282			\begin{eqnarray*}
283			\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
284			\theta^{n+1/2} & = & \theta^*
285			\\
286			\left[ 1 - \partial_r \kappa_v^S \partial_r \right]
287			S^{n+1/2} & = & S^*
288			\end{eqnarray*}
289			with
290			\begin{eqnarray*}
291			\theta^* & = &
292			\theta ^{(n-1/2)} + \Delta t G_{\theta} ^{(n)}
293			\\
294			S^* & = &
295			S ^{(n-1/2)} + \Delta t G_{S} ^{(n)}
296			\end{eqnarray*}
297			And
298			\begin{eqnarray*}
299			%{b'}^{n+1/2} & = & b'(\theta^{n+1/2},S^{n+1/2},r)
300			%\\
301			%\partial_r {\phi'_{hyd}}^{n+1/2} & = & {-b'}^{n+1/2}
302			{\phi'_{hyd}}^{n+1/2} & = & \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r) dr
303	jmc	1.4	%\label{eq-tDsC-hyd}
304	adcroft	1.1	\\
305			\vec{\bf v} ^{n+1}
306	jmc	1.3	+ \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
307			+ \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
308	adcroft	1.1	- \partial_r A_v \partial_r \vec{\bf v}^{n+1}
309			& = &
310			\vec{\bf v}^*
311	jmc	1.4	%\label{eq-tDsC-Hmom}
312	adcroft	1.1	\\
313			\epsilon_{fs} {\eta}^{n+1} + \Delta t
314	jmc	1.3	{\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
315	adcroft	1.1	& = &
316			\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
317			\\
318			\epsilon_{nh} \left( \dot{r} ^{n+1}
319			+ \Delta t \partial_r {\phi'_{nh}} ^{n+1}
320			\right)
321			& = & \epsilon_{nh} \dot{r}^*
322	jmc	1.4	%\label{eq-tDsC-Vmom}
323	adcroft	1.1	\\
324	jmc	1.3	{\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
325	adcroft	1.1	& = & 0
326	jmc	1.4	%\label{eq-tDsC-cont}
327	adcroft	1.1	\end{eqnarray*}
328			with
329			\begin{eqnarray*}
330			\vec{\bf v}^* & = &
331			\vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
332	jmc	1.3	+ \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{n+1/2}
333	adcroft	1.1	\\
334			\dot{r}^* & = &
335			\dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
336			\end{eqnarray*}
337
338			%---------------------------------------------------------------------
339
340	jmc	1.3	\subsection{Surface pressure}
341	adcroft	1.1
342	jmc	1.4	Substituting \ref{eq-tDsC-Hmom} into \ref{eq-tDsC-cont}, assuming
343	adcroft	1.1	$\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$:
344	jmc	1.2	\begin{eqnarray}
345	adcroft	1.1	\epsilon_{fs} {\eta}^{n+1} -
346	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{min})
347			{\bf \nabla}_h b_s {\eta}^{n+1}
348	adcroft	1.1	= {\eta}^*
349	jmc	1.4	\label{eq-solve2D}
350	jmc	1.2	\end{eqnarray}
351	adcroft	1.1	where
352	jmc	1.2	\begin{eqnarray}
353	adcroft	1.1	{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
354	jmc	1.3	\Delta t {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^* dr
355	adcroft	1.1	\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
356	jmc	1.4	\label{eq-solve2D_rhs}
357	jmc	1.2	\end{eqnarray}
358	adcroft	1.1
359	jmc	1.4	Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-tDsC-Hmom}
360	adcroft	1.1	would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic
361			($\epsilon_{nh}=0$):
362			$$
363			\vec{\bf v}^{n+1} = \vec{\bf v}^{*}
364	jmc	1.3	- \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
365	adcroft	1.1	$$
366
367			This is known as the correction step. However, when the model is
368			non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
369			additional equation for $\phi'_{nh}$. This is obtained by
370	jmc	1.4	substituting \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
371			\ref{eq-tDsC-cont}:
372	adcroft	1.1	\begin{equation}
373	jmc	1.3	\left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
374	adcroft	1.1	= \frac{1}{\Delta t} \left(
375	jmc	1.3	{\bf \nabla}_h \cdot \vec{\bf v}^{*} + \partial_r \dot{r}^ \right)
376	adcroft	1.1	\end{equation}
377			where
378			\begin{displaymath}
379	jmc	1.3	\vec{\bf v}^{*} = \vec{\bf v}^ - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
380	adcroft	1.1	\end{displaymath}
381			Note that $\eta^{n+1}$ is also used to update the second RHS term
382			$\partial_r \dot{r}^* $ since
383			the vertical velocity at the surface ($\dot{r}_{surf}$)
384			is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.
385
386			Finally, the horizontal velocities at the new time level are found by:
387			\begin{equation}
388			\vec{\bf v}^{n+1} = \vec{\bf v}^{**}
389	jmc	1.3	- \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
390	adcroft	1.1	\end{equation}
391			and the vertical velocity is found by integrating the continuity
392			equation vertically.
393			Note that, for convenience regarding the restart procedure,
394			the integration of the continuity equation has been
395			moved at the beginning of the time step (instead of at the end),
396			without any consequence on the solution.
397
398			Regarding the implementation, all those computation are done
399			within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent
400			{\it CALL}s.
401	jmc	1.4	The standard method to solve the 2D elliptic problem (\ref{eq-solve2D})
402	adcroft	1.1	uses the conjugate gradient method (routine {\it CG2D}); The
403			the solver matrix and conjugate gradient operator are only function
404			of the discretized domain and are therefore evaluated separately,
405			before the time iteration loop, within {\it INI\_CG2D}.
406			The computation of the RHS $\eta^*$ is partly
407			done in {\it CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}.
408
409			The same method is applied for the non hydrostatic part, using
410			a conjugate gradient 3D solver ({\it CG3D}) that is initialized
411			in {\it INI\_CG3D}. The RHS terms of 2D and 3D problems
412			are computed together, within the same part of the code.
413
414			\newpage
415			%-----------------------------------------------------------------------------------
416			\subsection{Crank-Nickelson barotropic time stepping}
417
418			The full implicit time stepping described previously is unconditionally stable
419			but damps the fast gravity waves, resulting in a loss of
420			gravity potential energy.
421			The modification presented hereafter allows to combine an implicit part
422			($\beta,\gamma$) and an explicit part ($1-\beta,1-\gamma$) for the surface
423			pressure gradient ($\beta$) and for the barotropic flow divergence ($\gamma$).
424			\\
425			For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
426			$\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
427			stable, Crank-Nickelson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$
428			corresponds to the forward - backward scheme that conserves energy but is
429			only stable for small time steps.\\
430			In the code, $\beta,\gamma$ are defined as parameters, respectively
431			{\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
432			the main data file "{\it data}" and are set by default to 1,1.
433
434	jmc	1.4	Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:
435	adcroft	1.1	$$
436			\frac{ \vec{\bf v}^{n+1} }{ \Delta t }
437	jmc	1.3	+ {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
438			+ \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
439	adcroft	1.1	= \frac{ \vec{\bf v}^* }{ \Delta t }
440			$$
441			$$
442			\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
443	jmc	1.3	+ {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o}
444	adcroft	1.1	[ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
445			= \epsilon_{fw} (P-E)
446			$$
447			where:
448			\begin{eqnarray*}
449			\vec{\bf v}^* & = &
450			\vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
451	jmc	1.3	+ (\beta-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n}
452			+ \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
453	adcroft	1.1	\\
454			{\eta}^* & = &
455			\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
456	jmc	1.3	- \Delta t {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o}
457	adcroft	1.1	[ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
458			\end{eqnarray*}
459			\\
460			In the hydrostatic case ($\epsilon_{nh}=0$),
461			this allow to find ${\eta}^{n+1}$, according to:
462			$$
463			\epsilon_{fs} {\eta}^{n+1} -
464	jmc	1.3	{\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{min})
465			{\bf \nabla}_h {\eta}^{n+1}
466	adcroft	1.1	= {\eta}^*
467			$$
468			and then to compute (correction step):
469			$$
470			\vec{\bf v}^{n+1} = \vec{\bf v}^{*}
471	jmc	1.3	- \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
472	adcroft	1.1	$$
473
474			The non-hydrostatic part is solved as described previously.
475			\\ \\
476			N.B:
477			\\
478			a) The non-hydrostatic part of the code has not yet been
479			updated, %since it falls out of the purpose of this test,
480			so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
481			\\
482			b) To remind, the stability criteria with the Crank-Nickelson time stepping
483			for the pure linear gravity wave problem in cartesian coordinate is:
484			\\
485			$\star$~ $\beta + \gamma < 1$ : unstable
486			\\
487			$\star$~ $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
488			\\
489			$\star$~ $\beta + \gamma \geq 1$ : stable if
490			%, for all wave length $(k\Delta x,l\Delta y)$
491			$$
492			c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
493			$$
494			$$
495			\mbox{with }~
496			%c^2 = 2 g H {\Delta t}^2
497			%(\frac{1-cos 2 \pi / k}{\Delta x^2}
498			%+\frac{1-cos 2 \pi / l}{\Delta y^2})
499			%$$
500			%Practically, the most stringent condition is obtained with $k=l=2$ :
501			%$$
502			c_{max} = 2 \Delta t \: \sqrt{g H} \:
503			\sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
504			$$