s_algorithm/text/time_stepping.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/time_stepping.tex,v 1.2 2001/08/09 16:22:09 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
\\
\partial_t S & = & G_s
\\
b' & = & b'(\theta,S,r)
\\
\partial_r \phi'_{hyd} & = & -b'
\label{eq-r-split-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-r-split-hmom}
\\
\epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}
+ \epsilon_{nh} \partial_r \phi'_{nh}
& = & \epsilon_{nh} G_{\dot{r}} 
\label{eq-r-split-rdot}
\\
{\bf \nabla}_h \cdot \vec{\bf v} + \partial_r \dot{r}
& = & 0
\label{eq-r-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-r-split-rdot} 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-cont-2D}
\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-r-split-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}
 
The general algorithm consist in  a "predictor step" that computes
the forward tendencies ("G" terms") and all 
the "first guess" values star notation): 
$\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$
in non-hydrostatic mode). This is done in the routine {\it DYNAMICS}.

Then the implicit terms that appear here on the left hand side (LHS),
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 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-rtd-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.

optionnal forcing terms (package):\\
Regarding optional forcing terms (usually part of a "package"), 
that a 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 additionnal 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.

\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^*
\\
\left[ 1 - \partial_r \kappa_v^S \partial_r \right]
S^{n+1} & = & S^*
\\
%{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-rtd-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-rtd-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-rtd-eta}
\\
\epsilon_{nh} \left( \dot{r} ^{n+1}
+ \Delta t \partial_r {\phi'_{nh}} ^{n+1}
\right)
& = & \epsilon_{nh} \dot{r}^*
\label{eq-rtd-rdot}
\\
{\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
& = & 0
\label{eq-rtd-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 "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-rtd-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-rtd-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-rtd-rdot}
\\
{\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
& = & 0
%\label{eq-rtd-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-rtd-hmom} into \ref{eq-rtd-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{solve_2d}
\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{solve_2d_rhs}
\end{eqnarray}

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