\section{Numerical tests}

This section presents a series of numerical tests performed with
the Godunov, Van Leer I, Slopes and Prather's schemes.
For the Van Leer scheme, we use the strong limit to insure monotonicity
and positivity. For the Russell and Prather's scheme, only positivity
is insured by limiting fluxes following \cite{Prat:86}.
The Prather's and Slopes schemes were originally optimized and interfaced 
for the LMD model by Pascal Simon and Christophe Genthon.

We first derive an estimation of the numerical cost of the various
schemes in order to compare their accuracy at a constant numerical cost.
We then present a series of numerical experiments including
2D advection with an analytical wind field, tests of advection of an
idealized tracer with the GCM winds, and a more realistic test on the Radon
cycle.


\subsection{Numerical cost of the various schemes}

For practical purposes, the interesting question about accuracy is :
"What is gained or lost at a constant numerical cost when
changing the advection algorithm?".
Since our goal is 3D advection in a GCM, we must first have an idea of
how the different algorithms compare in terms of 
numerical cost in this particular context.

The easiest and most rigorous comparison is that of storage requirements
(note that for chemical
applications for instance, this issue may dominate the problem of CPU,
at least for short term simulations).
Godunov, Van Leer I and Van Leer IV schemes only need the storage
of one value (the mean mixing ratio) per mesh.
The Slopes scheme requires the storage of the first order moment
in each direction which means 4 times more.
Finally, Prather's scheme requires the additional storage of the 3
second order moments
in each direction and the 3 cross derivatives which makes a total of 10
variables for each grid point.

The comparison in terms of CPU is less straightforward.
Simulations on SUN Ultra workstations and on a CRAY-90 show
that there is typically a factor 2 between Slopes and Van Leer I
and between Prather and Slopes.

The main difference on vector machines is that you gain in running
models on a finer grid. So there is a larger advantage for using
lower order schemes on a finer grid.
For the tests presented hereafter, run on a CRAY-90, Slopes had
exactly the same cost for a 60x43 horizontal grid as Van Leer I for 120x85.

Finally, if playing with the horizontal resolution only, a doubling of
resolution in each horizontal direction for Van
Leer scheme I makes it equivalent to Slopes schemes in terms of storage
and similar in terms of CPU on vector machines.
The Slopes scheme in double resolution is also similar to
Prather's for both CPU and storage.
The doubling in each of the three directions makes Van Leer Scheme
more similar to Prather's.
Scalar machines are more favorable to higher order codes in terms of CPU.


\subsection{2-Dimensional tests}

We first present a test of horizontal advection
with a simple analytical wind field for which the advection can also be computed analytically.
For a non diverging horizontal flow in the longitude-latitude
($\lambda$, $\phi$) reference frame, the zonal ($u$) and meridional ($v$)
wind components satisfy the following relation:
\begin{equation}
\frac{\partial u}{\partial \lambda}+\frac{\partial v\cos\phi}{\partial \phi}=0
\end{equation}
\def\ur{U_0}
The analytical wind field used for advection is derived from the following potential
\begin{equation}
\label{eq:psi}
\Psi= a \ur \cos^2\frac{\lambda}{2}\cos^2\phi
\end{equation}
(where $\ur$ is a normalizing velocity and $a$ the planetary radius) as
\begin{eqnarray} \label{eq:upsi}
u = \frac{1}{a} \frac{\partial \Psi}{\partial \phi}  =- 2 \ur \cos\phi\sin\phi \cos^2\frac{\lambda}{2}\\ \label{eq:vpsi}
v= \frac{1}{a\cos\phi}\frac{\partial \Psi}{\partial \lambda} =
-\ur\cos\phi\cos\frac{\lambda}{2}\sin\frac{\lambda}{2}
\end{eqnarray}

\def\ddd{\mbox{d}}
Following a trajectory (constant value of $\Psi$) the meridional velocity
$v=a \ddd\phi/\ddd t$ can be rewritten by combining \eq{psi} and \eq{vpsi} as
\begin{equation}
a \ddd\phi/\ddd t= -\sqrt{\Psi \ur/a}\mbox{ sign}\lambda\sqrt{1-\cos^2(\lambda/2)}
\end{equation}
By introducing $\alpha=\sqrt{\Psi\ur/a^3}$ and after replacing $\cos^2(\lambda/2)$ according to \eq{psi}, this equation can be integrated as
\begin{equation}
\begin{eqalign}
-\alpha\int \mbox{ sign}{\lambda} \ddd t & =
\int \frac{\cos\phi \ddd \phi}{\sqrt{\cos^2\phi-\Psi/(a\ur)}} \\ & =
\depb{\arcsin\dep{\frac{\sin\phi}{\sqrt{1-\Psi/(a\ur)}}}}
\end{eqalign}
\end{equation}
The trajectories correspond to iso-values of $\Psi$. The time law is that
given above. The time taken by an air parcel to come back to its
original position is $2\pi/\alpha=2\pi a/(\ur\cos\phi\cos\lambda/2)$.
This period it shorter at the center of the domain and infinite on the boundaries.

For the numerical tests, the time step was chosen small enough to
avoid the problem of near polar advection for the Prather's and Slopes
schemes for which the special treatment in longitude described above is not implemented.
In practice, the length of the simulation was fixed in order to have a full
rotation at the center of the domain. The total number of time-steps is 4000, 6000 and 12000
for increasing resolutions.
Note that a change of time step modifies the accuracy of the results but
much less than a change of resolution.


\begin{figure}
% Viewport original de uvref 44 80 586 239
\centerline{\includegraphics[width=17cm,trim=0 10 0 0]{\EPS/uvref.eps}}
%viewport original de uv.eps 42 63 586 539
\centerline{\includegraphics[width=17cm]{\EPS/uv.eps}}
\caption{Test of horizontal advection with an analytical wind field.
The initial tracer distribution is a Gaussian function of latitude only.
The upper panel shows the tracer distribution at time zero, and at time
$T/2$ and $T$ as computed with the analytical solution. The lower panel
shows the results at time $T$ the result of advection with the
different transport schemes for various resolutions
(full resolution with $120\times 60$ points, 1/2 with $60 \times 30$
and 1/3 with $40\times 20$ in longitude-latitude.
\label{fg:w}}
\end{figure}

The winds are shown in the upper left corner in \fg{w} with an arbitrary scale as well as the initial distribution chosen for numerical tests : a Gaussian in longitude independent of latitude. The two panels on the right show successive distributions computed analytically from the formula above.
The air rotates faster near the center of the figure than close to the boundary which produces a spiral.

The very accurate reconstruction of the fine central filament (values larger than 0.9) by the Prather's scheme for the full resolution is of course very impressive.
However, the Slopes scheme does as well or even somewhat better for the full resolution than Prather's scheme in half resolution. Same thing between Van Leer scheme I on the full grid and the slope scheme in low resolution.

So, for the impact of changing the resolution in two directions only, Van Leer scheme I is somewhat better than Slopes and Slopes better than Prather's
at least in terms of storage to cost ratio.
When changing the resolution in the three directions, the comparison must be made between the high and intermediate resolution with a significant advantage for Prather's with respect to Slopes or Slopes with respect to Van Leer I.

Note also the large difference in performance between the Godunov scheme and Van Leer I despite of 
their equivalent storage requirement.

\subsection{GCM simulations with an idealized tracer\label{sec:ideal}}

Here we present some experiments performed with an idealized tracer,
assumed to be emitted on continents only, with the same source functions
everywhere and no vertical variation.
In addition to its resemblance to known trace species such as 
radon, this allows to have very different spatial scales for the source depending
on the region on the globe.
The tracer is destroyed exponentially with a time constant of 10 days.
Those simulations are also idealized in the sense that the effect of
sub-grid scale processes (convection, turbulent mixing in the planetary
boundary layer) on the tracer distribution are not accounted for.

The GCM simulation is performed on a regular horizontal grid with
120 points in longitude and 85 in latitude and 20 layers on the vertical.

In order to compare the various schemes with different numerical
resolutions but exactly the same wind field, the
advection schemes were run, either at the full model resolution,
or at resolution 1/2 or at resolution 1/3, but in a unique GCM experiment.

To run one of the schemes at resolution 1/2 for instance,
the meshes are  gathered horizontally four by four
(two by two in each horizontal direction)
as illustrated in \fg{groupe}.
The horizontal mass flux at an interface is just the sum (in the $x$ direction
for $V$ and in the $y$ direction for $U$) of two consecutive GCM mass fluxes
taken from the full grid and the vertical mass fluxes are summed four by four.

\begin{figure}
\centerline{\includegraphics[width=8.cm]{\EPS/groupe.eps}}
\caption{Illustration of the mesh gathering\label{fg:groupe}}
\end{figure}

The simulation has been run on 40 days starting on 1st of july 1988.
Sea Surface temperatures were taken from AMIP.
The model was initialized with the ECMWF analysis.

\begin{figure}
\centerline{\includegraphics[width=18.cm]{/d2/hourdin/tex/ARTICLES/vanleer/EPS/new1.eps}}
\caption{Advection of an idealized tracer with continental sources in
a GCM simulation performed with 120 longitudes by 85 latitudes
and 20 vertical layers. ``Res. 1/2'' means that the meshes have
been gathered horizontally 2 by 2 (resp. 3 by 3) in each horizontal
direction for computing advection.
The tracer field is plotted after 20 days at 500 mbar.\label{fg:new1}}
\end{figure}

\fg{new1} shows the tracer field at 500~mbar, as simulated with the
various schemes and various resolutions at day 33.
Although the tracer field
is dominated by the continent-ocean contrasts in the source term, it
clearly shows regions of advection of marine air on continents and
vice versa. This particular date was chosen because of the very nice
filament of continental air exported from south America over south
atlantic.

Once again, the more expensive schemes are more accurate, but, for
a fixed numerical cost, at least in terms of storage,
the results are very similar for the various schemes.

\begin{figure}
\centerline{\includegraphics[width=15cm]{/d2/hourdin/tex/ARTICLES/vanleer/EPS/evol.eps}}
\caption{Time evolutions on the 40 days of integration of the root mean square error computed with respect to Prather's advection at full GCM resolution
at 500~mbar.
Left panels show errors computed on the full grid. Panel on the right
shows errors on the mean tracer mixing ratio as computed
 on the half resolution grid.
\label{fg:evol}}
\end{figure}

A quantitative comparison is shown in \fg{evol}, which displays
the time evolution of the root mean square error computed with respect to
Prather's results considered as the reference.
The left panel shows the time evolution of the root-mean square error
computed by considering Prather's advection in full resolution
as a reference. The curves for Slopes at half resolution and Van Leer I are
almost superimposed which means that,
 when looking at the results at a 60$\times$43 resolution, the results of the 
slopes advection are identical to the results of the 120$\times$85 advection
with Van Leer I for an equivalent storage requirement.

On the right panel, errors are computed on the finest grid.
The error for Van Leer scheme I are slightly larger than on the left
panel but by less than a factor 2, which means that Van leer scheme I
succeeds in catching part of the structure at the grid scale.
The errors are much larger for the schemes computed on a coarser grid.
Of course, this error could be significantly decreased by reconstructing
the distribution of the finest grid using
the higher order moments for the Slopes and Prather's schemes.
However, this figure gives an idea of what is probably lost by using
Prather's or Slopes
scheme on a coarser grid and not being able to correctly  treat the 
higher order moments, for instance in strongly nonlinear processes such
as chemical reactions.

Finally, note that Van Leer IV does not give results much more accurate
than Van Leer I.

% PARTIE D'ALEXANDRE SUR LE RADON
\input{\local/alex_part}