\section{1-Dimensional derivation}

\subsection{General description\label{sec:general}}

\def\diradv{\EPS}
\def\dt{\delta_t}
\def\bq{{\hat{q}}} % q au bord de la maille

In one dimension, and after integration over one time step,
\eq{m} and \eq{q} just read
\begin{equation}
\label{eq:dm}
\dt (m)_i=U_{i-\dem}-U_{i+\dem}
\end{equation}
and
\begin{equation}
\label{eq:dqm}
\dt(qm)_i=F_{i-\dem}-F_{i+\dem}
\end{equation}
where indices correspond now to position on the $x$ axis,
 $\dt$ represents the finite difference in time,
 and $U_{i+\dem}$ and $F_{i+\dem}$ are the mass and tracer transfers
 across the interface $i+\dem$
 during one time-step (fluxes integrated in time).
$F_{i+\dem}$ can be rewritten as the product of the mass transfer $U_{i+\dem}$ by the mean tracer mixing ratio $\bq_{i+\dem}$ corresponding to the air crossing the interface.

Van Leer's method consists
in approximating the sub-grid scale distribution by a polynomial
for which $\bq_{i+\dem}$ and hence advection can be computed exactly.
The principle and notations are illustrated in \fg{schema2b} for the
case of a first order polynomial.

\begin{figure}
\centerline{\includegraphics[width=8.cm]{\EPS/schema2b.eps}}
\caption{General principle of Van Leer finite volume schemes and notations for the present paper illustrated in the case where a first order polynomial is used to approximate the sub-grid scale distribution (see text for details).
The vertical axis corresponds to the tracer mixing ratio. The horizontal
axis is shown both for grid index and for the mass counted from the
center of mesh $i$ and normalized by the mass of this mesh ($m_i$).
The shaded areas correspond to the quantity of tracer which is advected
through the mesh interfaces during one time step.
\label{fg:schema2b}}
\end{figure}

With these notations, the new values of $q$ and $m$ after one time step are
simply
\begin{equation}\label{eq:q1d}
q^*_i=\dep{q_im_i+U_{i-\dem}\bq_{i-\dem}-U_{i+\dem}\bq_{i+\dem}}/m^*_i
\end{equation}
and
\begin{equation}\label{eq:m1d}
m^*_i=m_i+U_{i-\dem}-U_{i+\dem}
\end{equation}


For positive values of $U_{i-\dem}$ and $U_{i+\dem}$, \eq{q1d} can be rewritten
as
\begin{equation}\label{eq:q1db}
q^*_i=\depb{U_{i-\dem}\bq_{i-\dem}+\dep{m_i-U_{i+\dem}}\bq_i}/m^*_i
\end{equation}
where $\bq_i$
is the mean mixing ratio
in that part of the air which remains in volume $i$ during the time-step
(see \fg{schema2b}). With these notations, the Van Leer scheme can be described
as the sequence of the two following steps~: first, a new distribution
$\bq$ (crosses on \fg{schema2b}) is defined from distribution $q$
(dots) based on the polynomial approximation; then this new distribution
is rearranged among the grid boxes.

A very simple formulation of the {\bf sufficient requirement for monotonicity}
proposed by
\cite{VanL:77} in the general context of a compressible flow is that
{\bf $\bq$ must take its value between the two neighboring $q$ in regions where
$q$ is monotonic}.

Suppose to fix ideas that $q_{i-1}\leq q_i \leq q_{i+1}$ and $U\geq 0$.
The requirement above is just that
\begin{equation}
\label{eq:monotone1}
q_j \leq \bq_{j+\dem} \leq q_{j+1} \mbox{ for } j=i-1 \mbox{ and } j=i
\end{equation}
and that
\begin{equation}
\label{eq:monotone2}
q_{i-1} \leq \bq_i \leq q_i.
\end{equation}
Introducing condition~\ref{eq:monotone1} in \eq{q1d} insures that
$q^*_i \leq q_{i+1}$ while combination of \eq{q1db} and conditions
\ref{eq:monotone1} and \ref{eq:monotone2} leads to $q_{i-1}\leq q^*_i$.
Finally, we see that the condition given above insures that
{\bf $q^*_i$ lies between
$q_{i-1}$ and $q_i$ (resp. $q_i$ and $q_{i+1}$) for $U>0$ (resp. $U<0$)}.
As pointed out by \cite{VanL:77}, this statement
is a sufficient requirement for monotonicity.

With the additional requirement that ``$\bq_i=q_i$ (constant value in the mesh)
when $q_i$ is an extremum'' insures that maxima will not
increase, nor minima decrease. This in turn insures {\bf positivity}
and prevents from the formation of small scale numerical oscillations
coming from the different speeds of advection of the various Fourier
components of the tracer field \cite[dispersion, see e.g., ][]{Rood:87}.

With the requirements above, note also
the important property that {\bf a uniform tracer distribution is
unchanged by advection, and not only for non
divergent flows} (replace $q_i$, $\bq_{i-\dem}$ and $\bq_{i+\dem}$ by
a constant value in \eq{q1d}).

Finally, note that in all the derivations above,
the time-step is assumed to be small enough so that
the mass transfer $U$ is smaller or equal to the mass of the upstream volume.
Note also that if $U$ is exactly equal to the mass in the upstream mesh,
advection will be computed exactly independently of the polynomial
approximation.


\subsection{First order (Godunov) scheme}

A first approximation is to consider the tracer distribution inside
the mesh as a constant in which case $\bq$ is just the upstream value of $q$
($\bq_{i+\dem}=q_i$ if $U_{i+\dem}>0$ and $q_{i+1}$ otherwise).
As already stated,
this scheme which was first proposed by \cite{Godu:59} is computationally very
cheap, positive and monotonic,
but at the expense of a very large numerical diffusion.

\subsection{Second order scheme}

The next step is to assume that the tracer varies linearly as a function of mass with a slope $\dq_i=q_{i+\dem}-q_{i-\dem}$ where $q_{i\pm\dem}$ are the tracer mixing ratio extrapolated at the mesh boundary (see \fg{schema2b} for illustration).

Then $\bq_{i+\dem}$ is simply given by
\begin{eqnarray}
\bq_{i+1/2} & = & q_i + \frac{1}{2}\dep{1-\frac{U_{i+\dem}}{m_{i}}} \dq_i \mbox{\ if\ } U_{i+\dem}>0\\
& = & q_{i+1} - \frac{1}{2}\dep{1+\frac{U_{i+\dem}}{m_{i+1}}} \dq_{i+1} \mbox{ , otherwise.}
\end{eqnarray}

Various schemes were proposed by \cite{VanL:77,VanL:79} based on various estimations of $\dq$ two of which appear to be very attractive, the first one for its low numerical cost and the second one for its accuracy.

\begin{figure}
\centerline{\includegraphics[width=12.cm]{\EPS/schema3.eps}}
\caption{2nd order schemes : estimation of the slope by finite
differences for Van Leer I (upper panel) and by least-square fit
of the broken line distribution after a one time step integration
for the Slopes scheme (see text for details).\label{fg:schema3}}
\end{figure}

In Van Leer scheme I, slopes are just estimated at each time step by finite differences ($\delta q_i=\dep{q_{i+1}-q_{i-1}}/2$) as illustrated in the upper scheme in \fg{schema3}.
This scheme can be seen as a finite-volume version of the Fromm's algorithm \cite{From:68}.

In the second scheme, the slope at time step $t+\delta t$ in a given mesh is computed from the broken--line distribution corresponding to the same air mass at time $t$ by conserving the first moment of the distribution.
The new slope also corresponds to the linear distribution which minimizes the square root distance with the corresponding distribution at time $t$.
This is illustrated in the lower panels of \fg{schema3}.

This second scheme introduced in a 1-Dimensional context as Scheme III by \cite{VanL:77} was published later on by \cite{Russ:81}\footnote{The slope update, eq. (11) in the original Van Leer's paper, is the exact restriction to the 1-Dimensional non-divergent context of eq. (23) from \cite{Russ:81}.} in the more general context of 3-Dimensional advection. It has been used at NASA/GISS for several years in atmospheric transport studies. It is generally referred to as Slopes scheme in the GCM community. This name will be retained hereafter. For a detailed description of the numerical formulation in the 3-Dimensional divergent context, see \cite{Russ:81}.

\fg{1} gives an example of the typical behavior of the Godunov, Van Leer I and Slopes schemes in the 1-Dimensional context. The Slopes scheme is of course the most accurate but the main difference is between the very diffusive first order scheme and the two others. However, it clearly appears in the figure that Van Leer I and Slopes schemes are not monotonic (oscillations leading to negative values
are especially strong in the case of the square distribution).

\begin{figure}
\centerline{\includegraphics[width=6.cm]{\diradv/1d.eps}}
\caption{Comparison of the second order schemes.
This figure gives an example of 1-Dimensional advection on a periodic domain
with a constant velocity $U$ and
70 points of three different tracer distributions. The vertical
axis is in arbitrary units.
Results are shown for time-step 350 and $U/m=0.2$.\label{fg:1}}
\end{figure}

\subsection{Slopes limitation}

One of the main interest of the work by \cite{VanL:77} is to propose
a very simple way of insuring monotonicity. A first way consists in
applying directly the constraint given at the end of \sec{general}.

For the second order schemes,
he also gave the somewhat stronger but simpler requirement that
the entire tracer distribution in a given mesh $i$ does not exceed
the mean value in the two adjacent boxes and that the slope must be
of the same sign of that in the two adjacent boxes (the latter
criterion is automatically satisfied for Van Leer Scheme I).
This requirement can be formulated as a slope limitation.
For Van Leer scheme I, the complete computation of the slope with limitation is simply
\begin{eqnarray} \label{eq:slope}
\dq_i &=& \mbox{sign}\dep{q_{i+1}-q_i}\times \\ \nonumber
&& \mbox{min}\dep{\frac{|q_{i+1}-q_{i-1}|}{2},2|q_{i+1}-q_i|,2|q_i-q_{i-1}|}
\end{eqnarray}
if $q_i$ lies in between $q_{i-1}$ and $q_{i+1}$ and 0 otherwise.

\begin{figure}
\centerline{\includegraphics[width=8.cm]{\EPS/schema4b.epsi}}
\caption{An example of slope limitation.
\label{fg:schema4}}
\end{figure}

The effect of the strong and weak (direct application of \eq{monotone1} and
\eq{monotone2}) limitations,
which correspond respectively to eq. (66) and (74) given by \cite{VanL:77},
is illustrated in \fg{schema4}.

A numerical test of the impact of the weak and strong limits is shown
in \fg{limit}. The strong limit significantly alters the accuracy of the
Slopes scheme whereas the weak constraint corrects this scheme almost only
where oscillations occur, only reducing a little bit the maximum of the
Gaussian distribution.
For Van Leer Scheme I, the difference between the two limits is small 
when compared with the difference between the exact and numerical solutions.


\begin{figure}
\centerline{\includegraphics[width=6.cm]{/d2/hourdin/tex/ARTICLES/vanleer/EPS/limit.eps}}
\caption{Impact of the slope limitation on the Slopes (upper panel)
and Van Leer I (lower panel) schemes.
\label{fg:limit}
The cases with
no limits (same as in \fg{1}) and weak and strong limits
(see text for details) are shown for the square and Gaussian distributions.
}
\end{figure}

A very elegant alternative to \eq{slope} is to compute the slope as a geometric
average
\begin{equation}\label{eq:geom}
\dq_i = \frac{2\dep{q_{i+1}-q_i}\dep{q_i-q_{i-1}}}{q_{i+1}-q_{i-1}}
\end{equation}
which automatically satisfies the strong requirement.
This scheme is slightly more diffusive than direct application of
\eq{slope} but the gain in time is also weak since this equation is only part
of the
 scheme.

In the following tests we use \eq{slope} whereas \cite{Alle:91} retained
\eq{geom} which makes no real difference.

\subsection{Third-order scheme}

Approximation of the tracer distribution can be refined further by using second order 
polynomials approximant of the sub-grid scale tracer distribution.

The equivalent to Van Leer I scheme (Van Leer IV in the original paper), just consists
in evaluating the polynomial coefficients by finite differences computed on the values
of $q$ at the mesh centers. Van Leer did not pay much attention to this scheme which
presents the drawback of being equivalent to Van Leer I for $U/m=0.5$.
For small values of $U/m$ however it can be significantly more accurate. 
As Van Leer I, this scheme presents the advantage of carrying only one variable per
mesh. \cite{Wood:81} derived an alternative to Van Leer IV  with
an elaborate non-linear limiting-steepening filter \cite[See also,]{Wood:84,Cole:84}.
The resulting ``Piecewise Parabolic Method'' has been 
applied for meteorological modeling by \cite{Carp:90} and \cite{Lin:96}.
For the present paper, we only test the scheme originally derived by Van Leer.

The equivalent to Van Leer scheme III 
is a scheme in which first and second order moments are conserved during advection.
This scheme  \cite[scheme VI from][]{VanL:77} was introduced in the GCM
community by \cite{Prat:86}.
Prather's scheme is of  course much more accurate than the
previous schemes (see the next section).
However, limits are still required in order to  completely avoid 
spurious numerical oscillations.

\def\B{{$\bullet$}}
\def\A{{\large $\star$}}
\def\D{{\tiny $S$}}
\def\C{{\tiny $4$}}


\begin{figure}
%{\input{\local/test}}\\

\centerline{Root mean square error}
\caption{1D comparisons of the different schemes.\label{fg:tout1d}}
\B : reference with 70 meshes, regular grid and 1 full revolution.\\
\A : finer grid (140 meshes).\\
\C : 4 full revolutions (with 70 meshes).\\
\D : stretched grid (1 revolution, 70 meshes).\\
The reference case (\B) corresponds to \fg{1} with the same three tracer distributions
(Gaussian, Triangle and Square). As for \fg{1}, the values of the tracer field lie
between 0 and 1. The Courant number, $U/m$, is 0.1 so that one full
revolution corresponds to 700 time-steps.
For the stretched grid,
there is a factor 4 between the smallest and largest mesh.
The largest and smallest meshes are for point ~18 and ~52 respectively with the
initial distribution being centered at point 35.
\end{figure}