s_algorithm/text/time_stepping.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/time_stepping.tex,v 1.5 2001/09/24 19:30:40 jmc Exp $
% $Name:  $

The convention used in this section is as follows:
Time is ``discretized'' using a time step $\Delta t$   
and $\Phi^n$ refers to the variable $\Phi$ 
at time $t = n \Delta t$ . We use 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 and so  this section can be read independently.
 
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 \ref{sect:discrete}. Here its sufficient to mention that they contains
all the RHS terms except the pressure/geo-potential terms.

The switch $\epsilon_{nh}$ allows one 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}
is not used.

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

Here, $\epsilon_{fs}$ is a flag to switch between the free-surface,
$\epsilon_{fs}=1$, and a rigid-lid, $\epsilon_{fs}=0$. The flag
$\epsilon_{fw}$ determines whether an 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 (equation
\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 now presented with explicit
references to the Fortran code. We often refer to the discretized
equations of the model that are detailed in the following sections.

The general algorithm consist of a ``predictor step'' that computes
the forward tendencies ("G" terms") comprising of explicit-in-time
terms and 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}.

Terms that are integrated implicitly in time are handled at the end of
the {\it THERMODYNAMICS} and {\it DYNAMICS} routines. Then the
surface pressure and non hydrostatic pressure are solved for in ({\it
SOLVE\_FOR\_PRESSURE}).

Finally, in the ``corrector step'', (routine {\it
THE\_CORRECTION\_STEP}) the new values of $u,v,w,\theta,S$ are
determined (see details in \ref{sect:II.1.3}).

At this point, the regular time stepping process is complete. However,
after the correction step there are optional adjustments such as
convective adjustment or filters (zonal FFT filter, shapiro filter)
that can be applied on both momentum and tracer fields, just prior to
incrementing the time level to $n+1$.

Since the pressure solver precision is of the order of the ``target
residual'' and can 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 in {\it CALC\_EXACT\_ETA}. This ensures exact volume
conservation. Note that there is no need for an equivalent
non-hydrostatic ``exact conservation'' step, since by default $w$ is
already computed after the filters are applied.

Optional forcing terms (usually part of a physics ``package''), that
account for a specific source or sink process (e.g. condensation as a
sink of water vapor Q) are generally incorporated in the main
algorithm as follows: at the the beginning of the time step, the
additional tendencies are computed as a 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 (an implicit method).  Equations \ref{eq-tCsC-theta} to
\ref{eq-tCsC-cont} are 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_{fixed}}^{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 viscosity and diffusivity terms are not
considered as part of the ``{\it G}'' terms, but are written
separately here.

The default time-stepping method is the Adams-Bashforth quasi-second
order scheme in which the ``G'' terms are extrapolated forward in time
from time-levels $n-1$ and $n$ to time-level $n+1/2$ and provides a
simple but stable algorithm:
\begin{equation}
G^{(n+1/2)} = G^n + (1/2+\epsilon_{AB}) (G^n - G^{n-1})
\end{equation}
where $\epsilon_{AB}$ is a small number used to stabilize the time
stepping.

In the standard non-staggered formulation, the Adams-Bashforth time
stepping is also applied to the hydrostatic pressure/geo-potential
gradient, $\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}).
\marginpar{JMC: Clarify this term?}


\subsection{Stagger baroclinic time stepping}

An alternative to synchronous time-stepping is to stagger the momentum
and tracer fields in time. This allows the evaluation and gradient of
$\phi'_{hyd}$ to be centered in time with out needing to use the
Adams-Bashforth extrapoltion. This option is known as staggered
baroclinic time stepping because tracer and momentum are stepped
forward-in-time one after the other.  It can be activated by turning
on a run-time parameter {\bf staggerTimeStep} in namelist ``{\it
PARM01}''.

The main advantage of staggered time-stepping compared to synchronous,
is improved stability to internal gravity wave motions and a very
natural implementation of a 2nd order in time hydrostatic
pressure/geo-potential gradient term. However, synchronous
time-stepping might be better for convection problems, time dependent
forcing and diagnosing the model are also easier and it allows a more
efficient parallel decomposition.

The staggered baroclinic time-stepping scheme is equations
\ref{eq-tDsC-theta} to \ref{eq-tDsC-cont} except that \ref{eq-tDsC-hyd} is replaced with:
\begin{equation}
{\phi'_{hyd}}^{n+1/2} = \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r)
dr
\end{equation}
and the time-level for tracers has been shifted back by half:
\begin{eqnarray*}
\theta^* & = &
\theta ^{(n-1/2)} + \Delta t G_{\theta} ^{(n)}
\\
S^* & = &
S ^{(n-1/2)} + \Delta t G_{S} ^{(n)}
\\
\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*}


\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_{fixed})
{\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_{fixed}}^{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} yields $\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 the convenience of the restart
procedure, the vertical integration of the continuity equation has
been moved to the beginning of the time step (instead of at the end),
without any consequence on the solution.

Regarding the implementation of the surface pressure solver, all
computation are done within the routine {\it SOLVE\_FOR\_PRESSURE} and
its dependent calls.  The standard method to solve the 2D elliptic
problem (\ref{eq-solve2D}) uses the conjugate gradient method (routine
{\it CG2D}); 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
at the same point in the code.


\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 potential energy.  The modification presented now allows one
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_{fixed}}^{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_{fixed}}^{R_o} 
[ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
\end{eqnarray*}
\\
In the hydrostatic case ($\epsilon_{nh}=0$), allowing us to find
${\eta}^{n+1}$, thus:
$$
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
{\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. 

Note that:
\begin{enumerate}
\item The non-hydrostatic part of the code has not yet been 
updated, so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
\item The stability criteria with Crank-Nickelson time stepping
for the pure linear gravity wave problem in cartesian coordinates is:
\begin{itemize}
\item $\beta + \gamma < 1$ : unstable
\item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
\item $\beta + \gamma \geq 1$ : stable if
$$ 
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} }
$$
\end{itemize}
\end{enumerate}

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