s_algorithm/text/tracer.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/tracer.tex,v 1.2 2001/08/09 20:45:27 adcroft Exp $
% $Name:  $

\section{Tracer equations}

The basic discretization used for the tracer equations is the second
order piece-wise constant finite volume form of the forced
advection-diussion equations. There are many alternatives to second
order method for advection and alternative parameterizations for the
sub-grid scale processes. The Gent-McWilliams eddy parameterization,
KPP mixing scheme and PV flux parameterization are all dealt with in
separate sections. The basic discretization of the advection-diffusion
part of the tracer equations and the various advection schemes will be
described here.

\subsection{Time-stepping of tracers: ABII}

The default advection scheme is the centered second order method which
requires a second order or quasi-second order time-stepping scheme to
be stable. Historically this has been the quasi-second order
Adams-Bashforth method (ABII) and applied to all terms. For an
arbitrary tracer, $\tau$, the forced advection-diffusion equation
reads:
\begin{equation}
\partial_t \tau + G_{adv}^\tau = G_{diff}^\tau + G_{forc}^\tau
\end{equation}
where $G_{adv}^\tau$, $G_{diff}^\tau$ and $G_{forc}^\tau$ are the
tendencies due to advection, diffusion and forcing, respectively,
namely:
\begin{eqnarray}
G_{adv}^\tau & = & \partial_x u \tau + \partial_y v \tau + \partial_r w \tau
- \tau \nabla \cdot {\bf v} \\
G_{diff}^\tau & = & \nabla \cdot {\bf K} \nabla \tau
\end{eqnarray}
and the forcing can be some arbitrary function of state, time and
space.

The term, $\tau \nabla \cdot {\bf v}$, is required to retain local
conservation in conjunction with the linear implicit free-surface. It
only affects the surface layer since the flow is non-divergent
everywhere else. This term is therefore referred to as the surface
correction term. Global conservation is not possible using the
flux-form (as here) and a linearized free-surface
(\cite{Griffies00,Campin02}).

The continuity equation can be recovered by setting
$G_{diff}=G_{forc}=0$ and $\tau=1$.

The driver routine that calls the routines to calculate tendancies are
{\em S/R CALC\_GT} and {\em S/R CALC\_GS} for temperature and salt
(moisture), respectively. These in turn call a generic advection
diffusion routine {\em S/R GAD\_CALC\_RHS} that is called with the
flow field and relevent tracer as arguments and returns the collective
tendancy due to advection and diffusion. Forcing is add subsequently
in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendancy
array.

\fbox{ \begin{minipage}{4.75in}
{\em S/R GAD\_CALC\_RHS} ({\em pkg/generic\_advdiff/gad\_calc\_rhs.F})

$\tau$: {\bf tracer} (argument)

$G^{(n)}$: {\bf gTracer} (argument)

$F_r$: {\bf fVerT} (argument)

\end{minipage} }


The space and time discretizations are treated seperately (method of
lines). The Adams-Bashforth time discretization reads:
\marginpar{$\epsilon$: {\bf AB\_eps}}
\marginpar{$\Delta t$: {\bf deltaTtracer}}
\begin{equation}
\tau^{(n+1)} = \tau^{(n)} + \Delta t \left(
(\frac{3}{2} + \epsilon) G^{(n)} - (\frac{1}{2} + \epsilon) G^{(n-1)}
\right)
\end{equation}
where $G^{(n)} = G_{adv}^\tau + G_{diff}^\tau + G_{src}^\tau$ at time
step $n$.

Strictly speaking the ABII scheme should be applied only to the
advection terms. However, this scheme is only used in conjuction with
the standard second, third and fourth order advection
schemes. Selection of any other advection scheme disables
Adams-Bashforth for tracers so that explicit diffusion and forcing use
the forward method.

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

$\tau$: {\bf tracer} (argument)

$G^{(n)}$: {\bf gTracer} (argument)

$G^{(n-1)}$: {\bf gTrNm1} (argument)

$\Delta t$: {\bf deltaTtracer} (PARAMS.h)

\end{minipage} }

\begin{figure}
\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}
\caption{
Comparison of 1-D advection schemes. Courant number is 0.05 with 60
points and solutions are shown for T=1 (one complete period).
a) Shows the upwind biased schemes; first order upwind, DST3,
third order upwind and second order upwind.
b) Shows the centered schemes; Lax-Wendroff, DST4, centered second order,
centered fourth order and finite volume fourth order.
c) Shows the second order flux limiters: minmod, Superbee,
MC limiter and the van Leer limiter.
d) Shows the DST3 method with flux limiters due to Sweby with
$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter,
$\mu=c/(1-c)$.
\label{fig:advect-1d-lo}
}
\end{figure}

\begin{figure}
\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-hi.eps}}
\caption{
Comparison of 1-D advection schemes. Courant number is 0.89 with 60
points and solutions are shown for T=1 (one complete period).
a) Shows the upwind biased schemes; first order upwind and DST3.
Third order upwind and second order upwind are unstable at this Courant number.
b) Shows the centered schemes; Lax-Wendroff, DST4. Centered second order,
centered fourth order and finite volume fourth order and unstable at this
Courant number.
c) Shows the second order flux limiters: minmod, Superbee,
MC limiter and the van Leer limiter.
d) Shows the DST3 method with flux limiters due to Sweby with
$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter,
$\mu=c/(1-c)$.
\label{fig:advect-1d-hi}
}
\end{figure}

\section{Linear advection schemes}

The advection schemes known as centered second order, centered fourth
order, first order upwind and upwind biased third order are known as
linear advection schemes because the coefficient for interpolation of
the advected tracer are linear and a function only of the flow, not
the tracer field it self. We discuss these first since they are most
commonly used in the field and most familiar.

\subsection{Centered second order advection-diffusion}

The basic discretization, centered second order, is the default. It is
designed to be consistant with the continuity equation to facilitate
conservation properties analogous to the continuum. However, centered
second order advection is notoriously noisey and must be used in
conjuction with some finite amount of diffusion to produce a sensible
solution.

The advection operator is discretized:
\begin{equation}
{\cal A}_c \Delta r_f h_c G_{adv}^\tau = 
\delta_i F_x + \delta_j F_y + \delta_k F_r
\end{equation}
where the area integrated fluxes are given by:
\begin{eqnarray}
F_x & = & U \overline{ \tau }^i \\
F_y & = & V \overline{ \tau }^j \\
F_r & = & W \overline{ \tau }^k
\end{eqnarray}
The quantities $U$, $V$ and $W$ are volume fluxes defined:
\marginpar{$U$: {\bf uTrans} }
\marginpar{$V$: {\bf vTrans} }
\marginpar{$W$: {\bf rTrans} }
\begin{eqnarray}
U & = & \Delta y_g \Delta r_f h_w u \\
V & = & \Delta x_g \Delta r_f h_s v \\
W & = & {\cal A}_c w
\end{eqnarray}

For non-divergent flow, this discretization can be shown to conserve
the tracer both locally and globally and to globally conserve tracer
variance, $\tau^2$. The proof is given in \cite{Adcroft95,Adcroft97}.

\fbox{ \begin{minipage}{4.75in}
{\em S/R GAD\_C2\_ADV\_X} ({\em gad\_c2\_adv\_x.F})

$F_x$: {\bf uT} (argument)

$U$: {\bf uTrans} (argument)

$\tau$: {\bf tracer} (argument)

{\em S/R GAD\_C2\_ADV\_Y} ({\em gad\_c2\_adv\_y.F})

$F_y$: {\bf vT} (argument)

$V$: {\bf vTrans} (argument)

$\tau$: {\bf tracer} (argument)

{\em S/R GAD\_C2\_ADV\_R} ({\em gad\_c2\_adv\_r.F})

$F_r$: {\bf wT} (argument)

$W$: {\bf rTrans} (argument)

$\tau$: {\bf tracer} (argument)

\end{minipage} }


\subsection{Third order upwind bias advection}

Upwind biased third order advection offers a relatively good
compromise between accuracy and smoothness. It is not a ``positive''
scheme meaning false extrema are permitted but the amplitude of such
are significantly reduced over the centered second order method.

The third order upwind fluxes are discretized:
\begin{eqnarray}
F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i
         + \frac{1}{2} |U| \delta_i \frac{1}{6} \delta_{ii} \tau \\
F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j
         + \frac{1}{2} |V| \delta_j \frac{1}{6} \delta_{jj} \tau \\
F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k
         + \frac{1}{2} |W| \delta_k \frac{1}{6} \delta_{kk} \tau 
\end{eqnarray}

At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing
$\delta_{nn}$ to be evaluated. We are currently examing the accuracy
of this boundary condition and the effect on the solution.


\subsection{Centered fourth order advection}

Centered fourth order advection is formally the most accurate scheme
we have implemented and can be used to great effect in high resolution
simultation where dynamical scales are well resolved. However, the
scheme is noisey like the centered second order method and so must be
used with some finite amount of diffusion. Bi-harmonic is recommended
since it is more scale selective and less likely to diffuse away the
well resolved gradient the fourth order scheme worked so hard to
create.

The centered fourth order fluxes are discretized:
\begin{eqnarray}
F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i \\
F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j \\
F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k
\end{eqnarray}

As for the third order scheme, the best discretization near boundaries
is under investigation but currenlty $\delta_i \tau=0$ on a boundary.

\subsection{First order upwind advection}

Although the upwind scheme is the underlying scheme for the robust or
non-linear methods given later, we haven't actually supplied this
method for general use. It would be very diffusive and it is unlikely
that it could ever produce more useful results than the positive
higher order schemes.

Upwind bias is introduced into many schemes using the {\em abs}
function and is allows the first order upwind flux to be written:
\begin{eqnarray}
F_x & = & U \overline{ \tau }^i - \frac{1}{2} |U| \delta_i \tau \\
F_y & = & V \overline{ \tau }^j - \frac{1}{2} |V| \delta_j \tau \\
F_r & = & W \overline{ \tau }^k - \frac{1}{2} |W| \delta_k \tau
\end{eqnarray}

If for some reason, the above method is required, then the second
order flux limiter scheme described later reduces to the above scheme
if the limiter is set to zero.


\section{Non-linear advection schemes}

Non-linear advection schemes invoke non-linear interpolation and are
widely used in computational fluid dynamics (non-linear does not refer
to the non-linearity of the advection operator). The flux limited
advection schemes belong to the class of finite volume methods which
neatly ties into the spatial discretization of the model.

When employing the flux limited schemes, first order upwind or
direct-space-time method the time-stepping is switched to forward in
time.

\subsection{Second order flux limiters}

The second order flux limiter method can be cast in several ways but
is generally expressed in terms of other flux approximations. For
example, in terms of a first order upwind flux and second order
Lax-Wendroff flux, the limited flux is given as:
\begin{equation}
F = F_1 + \psi(r) F_{LW}
\end{equation}
where $\psi(r)$ is the limiter function,
\begin{equation}
F_1 = u \overline{\tau}^i - \frac{1}{2} |u| \delta_i \tau
\end{equation}
is the upwind flux,
\begin{equation}
F_{LW} = F_1 + \frac{|u|}{2} (1-c) \delta_i \tau
\end{equation}
is the Lax-Wendroff flux and $c = \frac{u \Delta t}{\Delta x}$ is the
Courant (CFL) number.

The limiter function, $\psi(r)$, takes the slope ratio
\begin{eqnarray}
r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \tau_{i} - \tau_{i-1} } & \forall & u > 0
\\
r = \frac{ \tau_{i+1} - \tau_{i} }{ \tau_{i} - \tau_{i-1} } & \forall & u < 0
\end{eqnarray}
as it's argument. There are many choices of limiter function but we
only provide the Superbee limiter \cite{Roe85}:
\begin{equation}
\psi(r) = \max[0,\min[1,2r],\min[2,r]]
\end{equation}


\subsection{Third order direct space time}

The direct-space-time method deals with space and time discretization
together (other methods that treat space and time seperately are known
collectively as the ``Method of Lines''). The Lax-Wendroff scheme
falls into this category; it adds sufficient diffusion to a second
order flux that the forward-in-time method is stable. The upwind
biased third order DST scheme is:
\begin{eqnarray}
F = u \left( \tau_{i-1}
        + d_0 (\tau_{i}-\tau_{i-1}) + d_1 (\tau_{i-1}-\tau_{i-2}) \right)
& \forall & u > 0 \\
F = u \left( \tau_{i}
        - d_0 (\tau_{i}-\tau_{i-1}) - d_1 (\tau_{i+1}-\tau_{i}) \right)
& \forall & u < 0
\end{eqnarray}
where
\begin{eqnarray}
d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\
d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| )
\end{eqnarray}
The coefficients $d_0$ and $d_1$ approach $1/3$ and $1/6$ repectively
as the Courant number, $c$, vanishes. In this limit, the conventional
third order upwind method is recovered. For finite Courant number, the
deviations from the linear method are analogous to the diffusion added
to centered second order advection in the Lax-Wendroff scheme.

The DST3 method described above must be used in a forward-in-time
manner and is stable for $0 \le |c| \le 1$. Although the scheme
appears to be forward-in-time, it is in fact second order in time and
the accuracy increases with the Courant number! For low Courant
number, DST3 produces very similar results (indistinguishable in
Fig.~\ref{fig:advect-1d-lo}) to the linear third order method but for
large Courant number, where the linear upwind third order method is
unstable, the scheme is extremely accurate
(Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots.

\subsection{Third order direct space time with flux limiting}

The overshoots in the DST3 method can be controlled with a flux limiter.
The limited flux is written:
\begin{equation}
F =
\frac{1}{2}(u+|u|)\left( \tau_{i-1} + \psi(r^+)(\tau_{i} - \tau_{i-1} )\right)
+
\frac{1}{2}(u-|u|)\left( \tau_{i-1} + \psi(r^-)(\tau_{i} - \tau_{i-1} )\right)
\end{equation}
where
\begin{eqnarray}
r^+ & = & \frac{\tau_{i-1} - \tau_{i-2}}{\tau_{i} - \tau_{i-1}} \\
r^- & = & \frac{\tau_{i+1} - \tau_{i}}{\tau_{i} - \tau_{i-1}}
\end{eqnarray}
and the limiter is the Sweby limiter:
\begin{equation}
\psi(r) = \max[0, \min[\min(1,d_0+d_1r],\frac{1-c}{c}r ]]
\end{equation}

\subsection{Multi-dimensional advection}

In many of the aforementioned advection schemes the behaviour in
multiple dimensions is not necessarily as good as the one dimensional
behaviour. For instance, a shape preserving monotonic scheme in one
dimension can have severe shape distortion in two dimensions if the
two components of horizontal fluxes are treated independently. There
is a large body of literature on the subject dealing with this problem
and among the fixes are operator and flux splitting methods, corner
flux methods and more. We have adopted a variant on the standard
splitting methods that allows the flux calculations to be implemented
as if in one dimension:
\begin{eqnarray}
\tau^{n+1/3} & = & \tau^{n}
- \Delta t \left( \frac{1}{\Delta x} \delta_i F^x(\tau^{n})
           + \tau^{n} \frac{1}{\Delta x} \delta_i u \right) \\
\tau^{n+2/3} & = & \tau^{n}
- \Delta t \left( \frac{1}{\Delta y} \delta_j F^y(\tau^{n+1/3})
           + \tau^{n} \frac{1}{\Delta y} \delta_i v \right) \\
\tau^{n+3/3} & = & \tau^{n}
- \Delta t \left( \frac{1}{\Delta r} \delta_k F^x(\tau^{n+2/3})
           + \tau^{n} \frac{1}{\Delta r} \delta_i w \right)
\end{eqnarray}

In order to incorporate this method into the general model algorithm,
we compute the effective tendancy rather than update the tracer so
that other terms such as diffusion are using the $n$ time-level and
not the updated $n+3/3$ quantities:
\begin{equation}
G^{n+1/2}_{adv} = \frac{1}{\Delta t} ( \tau^{n+3/3} - \tau^{n} )
\end{equation}
So that the over all time-stepping looks likes:
\begin{equation}
\tau^{n+1} = \tau^{n} + \Delta t \left( G^{n+1/2}_{adv} + G_{diff}(\tau^{n}) + G^{n}_{forcing} \right)
\end{equation}
1	adcroft	1.3	% $Header: /u/gcmpack/mitgcmdoc/part2/tracer.tex,v 1.2 2001/08/09 20:45:27 adcroft Exp $
2	adcroft	1.2	% $Name: $
3	adcroft	1.1
4			\section{Tracer equations}
5
6	adcroft	1.2	The basic discretization used for the tracer equations is the second
7			order piece-wise constant finite volume form of the forced
8			advection-diussion equations. There are many alternatives to second
9			order method for advection and alternative parameterizations for the
10			sub-grid scale processes. The Gent-McWilliams eddy parameterization,
11			KPP mixing scheme and PV flux parameterization are all dealt with in
12			separate sections. The basic discretization of the advection-diffusion
13			part of the tracer equations and the various advection schemes will be
14			described here.
15
16	adcroft	1.3	\subsection{Time-stepping of tracers: ABII}
17
18			The default advection scheme is the centered second order method which
19			requires a second order or quasi-second order time-stepping scheme to
20			be stable. Historically this has been the quasi-second order
21			Adams-Bashforth method (ABII) and applied to all terms. For an
22			arbitrary tracer, $\tau$, the forced advection-diffusion equation
23			reads:
24			\begin{equation}
25			\partial_t \tau + G_{adv}^\tau = G_{diff}^\tau + G_{forc}^\tau
26			\end{equation}
27			where $G_{adv}^\tau$, $G_{diff}^\tau$ and $G_{forc}^\tau$ are the
28			tendencies due to advection, diffusion and forcing, respectively,
29			namely:
30			\begin{eqnarray}
31			G_{adv}^\tau & = & \partial_x u \tau + \partial_y v \tau + \partial_r w \tau
32			- \tau \nabla \cdot {\bf v} \\
33			G_{diff}^\tau & = & \nabla \cdot {\bf K} \nabla \tau
34			\end{eqnarray}
35			and the forcing can be some arbitrary function of state, time and
36			space.
37
38			The term, $\tau \nabla \cdot {\bf v}$, is required to retain local
39			conservation in conjunction with the linear implicit free-surface. It
40			only affects the surface layer since the flow is non-divergent
41			everywhere else. This term is therefore referred to as the surface
42			correction term. Global conservation is not possible using the
43			flux-form (as here) and a linearized free-surface
44			(\cite{Griffies00,Campin02}).
45
46			The continuity equation can be recovered by setting
47			$G_{diff}=G_{forc}=0$ and $\tau=1$.
48
49			The driver routine that calls the routines to calculate tendancies are
50			{\em S/R CALC\_GT} and {\em S/R CALC\_GS} for temperature and salt
51			(moisture), respectively. These in turn call a generic advection
52			diffusion routine {\em S/R GAD\_CALC\_RHS} that is called with the
53			flow field and relevent tracer as arguments and returns the collective
54			tendancy due to advection and diffusion. Forcing is add subsequently
55			in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendancy
56			array.
57
58			\fbox{ \begin{minipage}{4.75in}
59			{\em S/R GAD\_CALC\_RHS} ({\em pkg/generic\_advdiff/gad\_calc\_rhs.F})
60
61			$\tau$: {\bf tracer} (argument)
62
63			$G^{(n)}$: {\bf gTracer} (argument)
64
65			$F_r$: {\bf fVerT} (argument)
66
67			\end{minipage} }
68
69
70			The space and time discretizations are treated seperately (method of
71			lines). The Adams-Bashforth time discretization reads:
72			\marginpar{$\epsilon$: {\bf AB\_eps}}
73			\marginpar{$\Delta t$: {\bf deltaTtracer}}
74			\begin{equation}
75			\tau^{(n+1)} = \tau^{(n)} + \Delta t \left(
76			(\frac{3}{2} + \epsilon) G^{(n)} - (\frac{1}{2} + \epsilon) G^{(n-1)}
77			\right)
78			\end{equation}
79			where $G^{(n)} = G_{adv}^\tau + G_{diff}^\tau + G_{src}^\tau$ at time
80			step $n$.
81
82			Strictly speaking the ABII scheme should be applied only to the
83			advection terms. However, this scheme is only used in conjuction with
84			the standard second, third and fourth order advection
85			schemes. Selection of any other advection scheme disables
86			Adams-Bashforth for tracers so that explicit diffusion and forcing use
87			the forward method.
88
89			\fbox{ \begin{minipage}{4.75in}
90			{\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F})
91
92			$\tau$: {\bf tracer} (argument)
93
94			$G^{(n)}$: {\bf gTracer} (argument)
95
96			$G^{(n-1)}$: {\bf gTrNm1} (argument)
97
98			$\Delta t$: {\bf deltaTtracer} (PARAMS.h)
99
100			\end{minipage} }
101
102			\begin{figure}
103			\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}
104			\caption{
105			Comparison of 1-D advection schemes. Courant number is 0.05 with 60
106			points and solutions are shown for T=1 (one complete period).
107			a) Shows the upwind biased schemes; first order upwind, DST3,
108			third order upwind and second order upwind.
109			b) Shows the centered schemes; Lax-Wendroff, DST4, centered second order,
110			centered fourth order and finite volume fourth order.
111			c) Shows the second order flux limiters: minmod, Superbee,
112			MC limiter and the van Leer limiter.
113			d) Shows the DST3 method with flux limiters due to Sweby with
114			$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter,
115			$\mu=c/(1-c)$.
116			\label{fig:advect-1d-lo}
117			}
118			\end{figure}
119
120			\begin{figure}
121			\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-hi.eps}}
122			\caption{
123			Comparison of 1-D advection schemes. Courant number is 0.89 with 60
124			points and solutions are shown for T=1 (one complete period).
125			a) Shows the upwind biased schemes; first order upwind and DST3.
126			Third order upwind and second order upwind are unstable at this Courant number.
127			b) Shows the centered schemes; Lax-Wendroff, DST4. Centered second order,
128			centered fourth order and finite volume fourth order and unstable at this
129			Courant number.
130			c) Shows the second order flux limiters: minmod, Superbee,
131			MC limiter and the van Leer limiter.
132			d) Shows the DST3 method with flux limiters due to Sweby with
133			$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter,
134			$\mu=c/(1-c)$.
135			\label{fig:advect-1d-hi}
136			}
137			\end{figure}
138
139			\section{Linear advection schemes}
140
141			The advection schemes known as centered second order, centered fourth
142			order, first order upwind and upwind biased third order are known as
143			linear advection schemes because the coefficient for interpolation of
144			the advected tracer are linear and a function only of the flow, not
145			the tracer field it self. We discuss these first since they are most
146			commonly used in the field and most familiar.
147
148	adcroft	1.2	\subsection{Centered second order advection-diffusion}
149
150			The basic discretization, centered second order, is the default. It is
151			designed to be consistant with the continuity equation to facilitate
152	adcroft	1.3	conservation properties analogous to the continuum. However, centered
153			second order advection is notoriously noisey and must be used in
154			conjuction with some finite amount of diffusion to produce a sensible
155			solution.
156
157			The advection operator is discretized:
158	adcroft	1.1	\begin{equation}
159	adcroft	1.3	{\cal A}_c \Delta r_f h_c G_{adv}^\tau =
160			\delta_i F_x + \delta_j F_y + \delta_k F_r
161	adcroft	1.1	\end{equation}
162	adcroft	1.2	where the area integrated fluxes are given by:
163			\begin{eqnarray}
164	adcroft	1.3	F_x & = & U \overline{ \tau }^i \\
165			F_y & = & V \overline{ \tau }^j \\
166			F_r & = & W \overline{ \tau }^k
167	adcroft	1.2	\end{eqnarray}
168	adcroft	1.1	The quantities $U$, $V$ and $W$ are volume fluxes defined:
169			\marginpar{$U$: {\bf uTrans} }
170			\marginpar{$V$: {\bf vTrans} }
171			\marginpar{$W$: {\bf rTrans} }
172			\begin{eqnarray}
173			U & = & \Delta y_g \Delta r_f h_w u \\
174			V & = & \Delta x_g \Delta r_f h_s v \\
175			W & = & {\cal A}_c w
176			\end{eqnarray}
177
178	adcroft	1.3	For non-divergent flow, this discretization can be shown to conserve
179			the tracer both locally and globally and to globally conserve tracer
180			variance, $\tau^2$. The proof is given in \cite{Adcroft95,Adcroft97}.
181
182			\fbox{ \begin{minipage}{4.75in}
183			{\em S/R GAD\_C2\_ADV\_X} ({\em gad\_c2\_adv\_x.F})
184
185			$F_x$: {\bf uT} (argument)
186
187			$U$: {\bf uTrans} (argument)
188
189			$\tau$: {\bf tracer} (argument)
190
191			{\em S/R GAD\_C2\_ADV\_Y} ({\em gad\_c2\_adv\_y.F})
192
193			$F_y$: {\bf vT} (argument)
194
195			$V$: {\bf vTrans} (argument)
196
197			$\tau$: {\bf tracer} (argument)
198
199			{\em S/R GAD\_C2\_ADV\_R} ({\em gad\_c2\_adv\_r.F})
200
201			$F_r$: {\bf wT} (argument)
202
203			$W$: {\bf rTrans} (argument)
204
205			$\tau$: {\bf tracer} (argument)
206
207			\end{minipage} }
208
209
210			\subsection{Third order upwind bias advection}
211
212			Upwind biased third order advection offers a relatively good
213			compromise between accuracy and smoothness. It is not a ``positive''
214			scheme meaning false extrema are permitted but the amplitude of such
215			are significantly reduced over the centered second order method.
216
217			The third order upwind fluxes are discretized:
218			\begin{eqnarray}
219			F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i
220			+ \frac{1}{2} \|U\| \delta_i \frac{1}{6} \delta_{ii} \tau \\
221			F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j
222			+ \frac{1}{2} \|V\| \delta_j \frac{1}{6} \delta_{jj} \tau \\
223			F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k
224			+ \frac{1}{2} \|W\| \delta_k \frac{1}{6} \delta_{kk} \tau
225			\end{eqnarray}
226
227			At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing
228			$\delta_{nn}$ to be evaluated. We are currently examing the accuracy
229			of this boundary condition and the effect on the solution.
230
231
232			\subsection{Centered fourth order advection}
233
234			Centered fourth order advection is formally the most accurate scheme
235			we have implemented and can be used to great effect in high resolution
236			simultation where dynamical scales are well resolved. However, the
237			scheme is noisey like the centered second order method and so must be
238			used with some finite amount of diffusion. Bi-harmonic is recommended
239			since it is more scale selective and less likely to diffuse away the
240			well resolved gradient the fourth order scheme worked so hard to
241			create.
242
243			The centered fourth order fluxes are discretized:
244			\begin{eqnarray}
245			F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i \\
246			F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j \\
247			F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k
248			\end{eqnarray}
249
250			As for the third order scheme, the best discretization near boundaries
251			is under investigation but currenlty $\delta_i \tau=0$ on a boundary.
252
253			\subsection{First order upwind advection}
254
255			Although the upwind scheme is the underlying scheme for the robust or
256			non-linear methods given later, we haven't actually supplied this
257			method for general use. It would be very diffusive and it is unlikely
258			that it could ever produce more useful results than the positive
259			higher order schemes.
260
261			Upwind bias is introduced into many schemes using the {\em abs}
262			function and is allows the first order upwind flux to be written:
263			\begin{eqnarray}
264			F_x & = & U \overline{ \tau }^i - \frac{1}{2} \|U\| \delta_i \tau \\
265			F_y & = & V \overline{ \tau }^j - \frac{1}{2} \|V\| \delta_j \tau \\
266			F_r & = & W \overline{ \tau }^k - \frac{1}{2} \|W\| \delta_k \tau
267			\end{eqnarray}
268
269			If for some reason, the above method is required, then the second
270			order flux limiter scheme described later reduces to the above scheme
271			if the limiter is set to zero.
272
273
274			\section{Non-linear advection schemes}
275
276			Non-linear advection schemes invoke non-linear interpolation and are
277			widely used in computational fluid dynamics (non-linear does not refer
278			to the non-linearity of the advection operator). The flux limited
279			advection schemes belong to the class of finite volume methods which
280			neatly ties into the spatial discretization of the model.
281
282			When employing the flux limited schemes, first order upwind or
283			direct-space-time method the time-stepping is switched to forward in
284			time.
285
286			\subsection{Second order flux limiters}
287
288			The second order flux limiter method can be cast in several ways but
289			is generally expressed in terms of other flux approximations. For
290			example, in terms of a first order upwind flux and second order
291			Lax-Wendroff flux, the limited flux is given as:
292			\begin{equation}
293			F = F_1 + \psi(r) F_{LW}
294			\end{equation}
295			where $\psi(r)$ is the limiter function,
296			\begin{equation}
297			F_1 = u \overline{\tau}^i - \frac{1}{2} \|u\| \delta_i \tau
298			\end{equation}
299			is the upwind flux,
300			\begin{equation}
301			F_{LW} = F_1 + \frac{\|u\|}{2} (1-c) \delta_i \tau
302			\end{equation}
303			is the Lax-Wendroff flux and $c = \frac{u \Delta t}{\Delta x}$ is the
304			Courant (CFL) number.
305
306			The limiter function, $\psi(r)$, takes the slope ratio
307			\begin{eqnarray}
308			r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \tau_{i} - \tau_{i-1} } & \forall & u > 0
309			\\
310			r = \frac{ \tau_{i+1} - \tau_{i} }{ \tau_{i} - \tau_{i-1} } & \forall & u < 0
311			\end{eqnarray}
312			as it's argument. There are many choices of limiter function but we
313			only provide the Superbee limiter \cite{Roe85}:
314			\begin{equation}
315			\psi(r) = \max[0,\min[1,2r],\min[2,r]]
316			\end{equation}
317
318
319			\subsection{Third order direct space time}
320
321			The direct-space-time method deals with space and time discretization
322			together (other methods that treat space and time seperately are known
323			collectively as the ``Method of Lines''). The Lax-Wendroff scheme
324			falls into this category; it adds sufficient diffusion to a second
325			order flux that the forward-in-time method is stable. The upwind
326			biased third order DST scheme is:
327			\begin{eqnarray}
328			F = u \left( \tau_{i-1}
329			+ d_0 (\tau_{i}-\tau_{i-1}) + d_1 (\tau_{i-1}-\tau_{i-2}) \right)
330			& \forall & u > 0 \\
331			F = u \left( \tau_{i}
332			- d_0 (\tau_{i}-\tau_{i-1}) - d_1 (\tau_{i+1}-\tau_{i}) \right)
333			& \forall & u < 0
334			\end{eqnarray}
335			where
336			\begin{eqnarray}
337			d_1 & = & \frac{1}{6} ( 2 - \|c\| ) ( 1 - \|c\| ) \\
338			d_2 & = & \frac{1}{6} ( 1 - \|c\| ) ( 1 + \|c\| )
339			\end{eqnarray}
340			The coefficients $d_0$ and $d_1$ approach $1/3$ and $1/6$ repectively
341			as the Courant number, $c$, vanishes. In this limit, the conventional
342			third order upwind method is recovered. For finite Courant number, the
343			deviations from the linear method are analogous to the diffusion added
344			to centered second order advection in the Lax-Wendroff scheme.
345
346			The DST3 method described above must be used in a forward-in-time
347			manner and is stable for $0 \le \|c\| \le 1$. Although the scheme
348			appears to be forward-in-time, it is in fact second order in time and
349			the accuracy increases with the Courant number! For low Courant
350			number, DST3 produces very similar results (indistinguishable in
351			Fig.~\ref{fig:advect-1d-lo}) to the linear third order method but for
352			large Courant number, where the linear upwind third order method is
353			unstable, the scheme is extremely accurate
354			(Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots.
355
356			\subsection{Third order direct space time with flux limiting}
357
358			The overshoots in the DST3 method can be controlled with a flux limiter.
359			The limited flux is written:
360			\begin{equation}
361			F =
362			\frac{1}{2}(u+\|u\|)\left( \tau_{i-1} + \psi(r^+)(\tau_{i} - \tau_{i-1} )\right)
363			+
364			\frac{1}{2}(u-\|u\|)\left( \tau_{i-1} + \psi(r^-)(\tau_{i} - \tau_{i-1} )\right)
365			\end{equation}
366			where
367			\begin{eqnarray}
368			r^+ & = & \frac{\tau_{i-1} - \tau_{i-2}}{\tau_{i} - \tau_{i-1}} \\
369			r^- & = & \frac{\tau_{i+1} - \tau_{i}}{\tau_{i} - \tau_{i-1}}
370			\end{eqnarray}
371			and the limiter is the Sweby limiter:
372			\begin{equation}
373			\psi(r) = \max[0, \min[\min(1,d_0+d_1r],\frac{1-c}{c}r ]]
374			\end{equation}
375	adcroft	1.1
376	adcroft	1.3	\subsection{Multi-dimensional advection}
377	adcroft	1.1
378	adcroft	1.3	In many of the aforementioned advection schemes the behaviour in
379			multiple dimensions is not necessarily as good as the one dimensional
380			behaviour. For instance, a shape preserving monotonic scheme in one
381			dimension can have severe shape distortion in two dimensions if the
382			two components of horizontal fluxes are treated independently. There
383			is a large body of literature on the subject dealing with this problem
384			and among the fixes are operator and flux splitting methods, corner
385			flux methods and more. We have adopted a variant on the standard
386			splitting methods that allows the flux calculations to be implemented
387			as if in one dimension:
388			\begin{eqnarray}
389			\tau^{n+1/3} & = & \tau^{n}
390			- \Delta t \left( \frac{1}{\Delta x} \delta_i F^x(\tau^{n})
391			+ \tau^{n} \frac{1}{\Delta x} \delta_i u \right) \\
392			\tau^{n+2/3} & = & \tau^{n}
393			- \Delta t \left( \frac{1}{\Delta y} \delta_j F^y(\tau^{n+1/3})
394			+ \tau^{n} \frac{1}{\Delta y} \delta_i v \right) \\
395			\tau^{n+3/3} & = & \tau^{n}
396			- \Delta t \left( \frac{1}{\Delta r} \delta_k F^x(\tau^{n+2/3})
397			+ \tau^{n} \frac{1}{\Delta r} \delta_i w \right)
398			\end{eqnarray}
399	adcroft	1.1
400	adcroft	1.3	In order to incorporate this method into the general model algorithm,
401			we compute the effective tendancy rather than update the tracer so
402			that other terms such as diffusion are using the $n$ time-level and
403			not the updated $n+3/3$ quantities:
404			\begin{equation}
405			G^{n+1/2}_{adv} = \frac{1}{\Delta t} ( \tau^{n+3/3} - \tau^{n} )
406			\end{equation}
407			So that the over all time-stepping looks likes:
408			\begin{equation}
409			\tau^{n+1} = \tau^{n} + \Delta t \left( G^{n+1/2}_{adv} + G_{diff}(\tau^{n}) + G^{n}_{forcing} \right)
410			\end{equation}