s_algorithm/text/time_stepping.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/time_stepping.tex,v 1.6 2001/09/26 20:19:52 adcroft 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''?}

\fbox{ \begin{minipage}{4.75in}
{\em S/R TIMESTEP} ({\em timestep.F})

$G_u^n$: {\bf Gu} ({\em DYNVARS.h})

$G_u^{n-1}, u^*$: {\bf GuNm1} ({\em DYNVARS.h})

$G_v^n$: {\bf Gv} ({\em DYNVARS.h})

$G_v^{n-1}, v^*$: {\bf GvNm1} ({\em DYNVARS.h})

$G_u^{(n+1/2)}$: {\bf GuTmp} (local)

$G_v^{(n+1/2)}$: {\bf GvTmp} (local)

\end{minipage} }


\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}

\fbox{ \begin{minipage}{4.75in}
{\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F})

$u^*$: {\bf GuNm1} ({\em DYNVARS.h})

$v^*$: {\bf GvNm1} ({\em DYNVARS.h})

$\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h)

$\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)

\end{minipage} }


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.

\fbox{ \begin{minipage}{4.75in}
{\em S/R CORRECTION\_STEP} ({\em correction\_step.F})

$\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)

$\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em DYNVARS.h)

$u^*$: {\bf GuNm1} ({\em DYNVARS.h})

$v^*$: {\bf GvNm1} ({\em DYNVARS.h})

$u^{n+1}$: {\bf uVel} ({\em DYNVARS.h})

$v^{n+1}$: {\bf vVel} ({\em DYNVARS.h})

\end{minipage} }


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.7	% $Header: /u/gcmpack/mitgcmdoc/part2/time_stepping.tex,v 1.6 2001/09/26 20:19:52 adcroft 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	adcroft	1.7	\marginpar{JMC: Clarify ``this term''?}
241
242			\fbox{ \begin{minipage}{4.75in}
243			{\em S/R TIMESTEP} ({\em timestep.F})
244
245			$G_u^n$: {\bf Gu} ({\em DYNVARS.h})
246
247			$G_u^{n-1}, u^*$: {\bf GuNm1} ({\em DYNVARS.h})
248
249			$G_v^n$: {\bf Gv} ({\em DYNVARS.h})
250
251			$G_v^{n-1}, v^*$: {\bf GvNm1} ({\em DYNVARS.h})
252
253			$G_u^{(n+1/2)}$: {\bf GuTmp} (local)
254
255			$G_v^{(n+1/2)}$: {\bf GvTmp} (local)
256
257			\end{minipage} }
258
259
260
261	adcroft	1.6
262	adcroft	1.1
263	jmc	1.3	\subsection{Stagger baroclinic time stepping}
264	adcroft	1.1
265	adcroft	1.6	An alternative to synchronous time-stepping is to stagger the momentum
266			and tracer fields in time. This allows the evaluation and gradient of
267			$\phi'_{hyd}$ to be centered in time with out needing to use the
268			Adams-Bashforth extrapoltion. This option is known as staggered
269			baroclinic time stepping because tracer and momentum are stepped
270			forward-in-time one after the other. It can be activated by turning
271			on a run-time parameter {\bf staggerTimeStep} in namelist ``{\it
272			PARM01}''.
273
274			The main advantage of staggered time-stepping compared to synchronous,
275			is improved stability to internal gravity wave motions and a very
276			natural implementation of a 2nd order in time hydrostatic
277			pressure/geo-potential gradient term. However, synchronous
278			time-stepping might be better for convection problems, time dependent
279			forcing and diagnosing the model are also easier and it allows a more
280			efficient parallel decomposition.
281	adcroft	1.1
282	adcroft	1.6	The staggered baroclinic time-stepping scheme is equations
283			\ref{eq-tDsC-theta} to \ref{eq-tDsC-cont} except that \ref{eq-tDsC-hyd} is replaced with:
284			\begin{equation}
285			{\phi'_{hyd}}^{n+1/2} = \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r)
286			dr
287			\end{equation}
288			and the time-level for tracers has been shifted back by half:
289	adcroft	1.1	\begin{eqnarray*}
290			\theta^* & = &
291			\theta ^{(n-1/2)} + \Delta t G_{\theta} ^{(n)}
292			\\
293			S^* & = &
294			S ^{(n-1/2)} + \Delta t G_{S} ^{(n)}
295			\\
296	adcroft	1.6	\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
297			\theta^{n+1/2} & = & \theta^*
298	adcroft	1.1	\\
299	adcroft	1.6	\left[ 1 - \partial_r \kappa_v^S \partial_r \right]
300			S^{n+1/2} & = & S^*
301	adcroft	1.1	\end{eqnarray*}
302
303
304	jmc	1.3	\subsection{Surface pressure}
305	adcroft	1.1
306	jmc	1.4	Substituting \ref{eq-tDsC-Hmom} into \ref{eq-tDsC-cont}, assuming
307	adcroft	1.1	$\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$:
308	jmc	1.2	\begin{eqnarray}
309	adcroft	1.1	\epsilon_{fs} {\eta}^{n+1} -
310	jmc	1.5	{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
311	jmc	1.3	{\bf \nabla}_h b_s {\eta}^{n+1}
312	adcroft	1.1	= {\eta}^*
313	jmc	1.4	\label{eq-solve2D}
314	jmc	1.2	\end{eqnarray}
315	adcroft	1.1	where
316	jmc	1.2	\begin{eqnarray}
317	adcroft	1.1	{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
318	jmc	1.5	\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr
319	adcroft	1.1	\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
320	jmc	1.4	\label{eq-solve2D_rhs}
321	jmc	1.2	\end{eqnarray}
322	adcroft	1.1
323	adcroft	1.7	\fbox{ \begin{minipage}{4.75in}
324			{\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F})
325
326			$u^*$: {\bf GuNm1} ({\em DYNVARS.h})
327
328			$v^*$: {\bf GvNm1} ({\em DYNVARS.h})
329
330			$\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h)
331
332			$\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
333
334			\end{minipage} }
335
336
337	adcroft	1.6	Once ${\eta}^{n+1}$ has been found, substituting into
338			\ref{eq-tDsC-Hmom} yields $\vec{\bf v}^{n+1}$ if the model is
339			hydrostatic ($\epsilon_{nh}=0$):
340	adcroft	1.1	$$
341			\vec{\bf v}^{n+1} = \vec{\bf v}^{*}
342	jmc	1.3	- \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
343	adcroft	1.1	$$
344
345			This is known as the correction step. However, when the model is
346			non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
347	adcroft	1.6	additional equation for $\phi'_{nh}$. This is obtained by substituting
348			\ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
349	jmc	1.4	\ref{eq-tDsC-cont}:
350	adcroft	1.1	\begin{equation}
351	jmc	1.3	\left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
352	adcroft	1.1	= \frac{1}{\Delta t} \left(
353	jmc	1.3	{\bf \nabla}_h \cdot \vec{\bf v}^{*} + \partial_r \dot{r}^ \right)
354	adcroft	1.1	\end{equation}
355			where
356			\begin{displaymath}
357	jmc	1.3	\vec{\bf v}^{*} = \vec{\bf v}^ - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
358	adcroft	1.1	\end{displaymath}
359			Note that $\eta^{n+1}$ is also used to update the second RHS term
360			$\partial_r \dot{r}^* $ since
361			the vertical velocity at the surface ($\dot{r}_{surf}$)
362			is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.
363
364			Finally, the horizontal velocities at the new time level are found by:
365			\begin{equation}
366			\vec{\bf v}^{n+1} = \vec{\bf v}^{**}
367	jmc	1.3	- \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
368	adcroft	1.1	\end{equation}
369			and the vertical velocity is found by integrating the continuity
370	adcroft	1.6	equation vertically. Note that, for the convenience of the restart
371			procedure, the vertical integration of the continuity equation has
372			been moved to the beginning of the time step (instead of at the end),
373	adcroft	1.1	without any consequence on the solution.
374
375	adcroft	1.7	\fbox{ \begin{minipage}{4.75in}
376			{\em S/R CORRECTION\_STEP} ({\em correction\_step.F})
377
378			$\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
379
380			$\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em DYNVARS.h)
381
382			$u^*$: {\bf GuNm1} ({\em DYNVARS.h})
383
384			$v^*$: {\bf GvNm1} ({\em DYNVARS.h})
385
386			$u^{n+1}$: {\bf uVel} ({\em DYNVARS.h})
387
388			$v^{n+1}$: {\bf vVel} ({\em DYNVARS.h})
389
390			\end{minipage} }
391
392
393
394	adcroft	1.6	Regarding the implementation of the surface pressure solver, all
395			computation are done within the routine {\it SOLVE\_FOR\_PRESSURE} and
396			its dependent calls. The standard method to solve the 2D elliptic
397			problem (\ref{eq-solve2D}) uses the conjugate gradient method (routine
398			{\it CG2D}); the solver matrix and conjugate gradient operator are
399			only function of the discretized domain and are therefore evaluated
400			separately, before the time iteration loop, within {\it INI\_CG2D}.
401			The computation of the RHS $\eta^*$ is partly done in {\it
402			CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}.
403
404			The same method is applied for the non hydrostatic part, using a
405			conjugate gradient 3D solver ({\it CG3D}) that is initialized in {\it
406			INI\_CG3D}. The RHS terms of 2D and 3D problems are computed together
407			at the same point in the code.
408	adcroft	1.7
409	adcroft	1.6
410	adcroft	1.1
411			\subsection{Crank-Nickelson barotropic time stepping}
412
413	adcroft	1.6	The full implicit time stepping described previously is
414			unconditionally stable but damps the fast gravity waves, resulting in
415			a loss of potential energy. The modification presented now allows one
416			to combine an implicit part ($\beta,\gamma$) and an explicit part
417			($1-\beta,1-\gamma$) for the surface pressure gradient ($\beta$) and
418			for the barotropic flow divergence ($\gamma$).
419	adcroft	1.1	\\
420			For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
421			$\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
422			stable, Crank-Nickelson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$
423			corresponds to the forward - backward scheme that conserves energy but is
424			only stable for small time steps.\\
425			In the code, $\beta,\gamma$ are defined as parameters, respectively
426			{\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
427			the main data file "{\it data}" and are set by default to 1,1.
428
429	jmc	1.4	Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:
430	adcroft	1.1	$$
431			\frac{ \vec{\bf v}^{n+1} }{ \Delta t }
432	jmc	1.3	+ {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
433			+ \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
434	adcroft	1.1	= \frac{ \vec{\bf v}^* }{ \Delta t }
435			$$
436			$$
437			\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
438	jmc	1.5	+ {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
439	adcroft	1.1	[ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
440			= \epsilon_{fw} (P-E)
441			$$
442			where:
443			\begin{eqnarray*}
444			\vec{\bf v}^* & = &
445			\vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
446	jmc	1.3	+ (\beta-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n}
447			+ \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
448	adcroft	1.1	\\
449			{\eta}^* & = &
450			\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
451	jmc	1.5	- \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
452	adcroft	1.1	[ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
453			\end{eqnarray*}
454			\\
455	adcroft	1.6	In the hydrostatic case ($\epsilon_{nh}=0$), allowing us to find
456			${\eta}^{n+1}$, thus:
457	adcroft	1.1	$$
458			\epsilon_{fs} {\eta}^{n+1} -
459	jmc	1.5	{\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
460	jmc	1.3	{\bf \nabla}_h {\eta}^{n+1}
461	adcroft	1.1	= {\eta}^*
462			$$
463			and then to compute (correction step):
464			$$
465			\vec{\bf v}^{n+1} = \vec{\bf v}^{*}
466	jmc	1.3	- \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
467	adcroft	1.1	$$
468
469			The non-hydrostatic part is solved as described previously.
470	adcroft	1.6
471			Note that:
472			\begin{enumerate}
473			\item The non-hydrostatic part of the code has not yet been
474			updated, so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
475			\item The stability criteria with Crank-Nickelson time stepping
476			for the pure linear gravity wave problem in cartesian coordinates is:
477			\begin{itemize}
478			\item $\beta + \gamma < 1$ : unstable
479			\item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
480			\item $\beta + \gamma \geq 1$ : stable if
481	adcroft	1.1	$$
482			c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
483			$$
484			$$
485			\mbox{with }~
486			%c^2 = 2 g H {\Delta t}^2
487			%(\frac{1-cos 2 \pi / k}{\Delta x^2}
488			%+\frac{1-cos 2 \pi / l}{\Delta y^2})
489			%$$
490			%Practically, the most stringent condition is obtained with $k=l=2$ :
491			%$$
492			c_{max} = 2 \Delta t \: \sqrt{g H} \:
493			\sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
494			$$
495	adcroft	1.6	\end{itemize}
496			\end{enumerate}
497