\section{Implementation in the GCM}

\subsection{The LMD-Z GCM}

The model used for this study is derived from the LMD GCM first described by
\cite{Sado:84}. In addition to several significant modifications in the
"physical package" through time, the model has been recently recoded in a more
organized and flexible manner. It includes in particular a zooming capability~:
longitudes and latitudes are fixed arbitrarily at the beginning of a
simulation. The GCM is briefly described below.

\subsubsection*{Dynamics}


The dynamical part of the code is based on a finite-difference
formulation of the primitive equations of meteorology
developed by
R.~Sadourny \cite[see \eg ][]{Sado:84} and coded by P.~Le Van.
The dynamical equations are discretized on a staggered latitude-longitude
Arakawa C-grid \cite[see \eg ][]{Kasa:77}~.
The C-grid is in fact a node-centered finite volume grid~: the scalar variables
are defined at the center of meshes whereas wind components are defined at
mesh interfaces.

Note that this dynamical code has been widely used not only for Earth but
also for the numerical simulation of the general circulation of
other planetary atmospheres, in particular for Mars \cite[]{Hour:93,Hour:95}
and Titan \cite[]{Hour:95b}.

\subsubsection*{Physical parameterizations}

The physical package used for the present study is very close to that
described by \cite{Letr:94}.
The radiation scheme is the same as used in the model of European Centre for
Medium-Range Weather Forecasts (ECMWF)~: the solar part is a refined version of
the scheme developed by \cite{Fouq:80} and the thermal infra-red part is
due to \cite{Morc:86}. Condensation is parameterized separately for convective and non-convective clouds. Convection is parameterized using in sequence 
a moist adiabatic adjustment and a modified version of the \cite{Kuo:1965}
algorithm. A prognostic equation for cloud water is included
\cite[]{Letr:91}. The large scale transport of vapor and liquid water
is taken into account
with Godunov scheme on the horizontal and a second order centered scheme
with a special treatment of negative concentrations on the vertical.
The planetary boundary layer scheme is based on a second order closure
model.
In this version of the model, the surface conditions over land are
treated using a simple bucket parameterization for moisture and a
single layer for heat and moisture. Note also that diurnal cycle
is not included.

The realism of the climate simulated by this revised version of the model
has been carefully checked with respect to that of the previous version
presented for instance by \cite{Harzallah:1994}.

\subsection{Implementation}

The numerical implementation of finite volume schemes is particularly
easy in the LMD GCM since mass fluxes are already defined at the mesh
interfaces (C-grid) and used to solve the continuity equations
which guaranties a complete consistency between tracer and mass
advection.

\subsubsection*{Time splitting}

The adaptation to 3-Dimensional grid is based on time splitting~: advection is
done sequentially in the 3 directions.
This approach has already been used successfully for  finite volume schemes
\cite[]{VanL:79,Alle:91,Russ:81,Prat:86}.

Here we adopt the sequence proposed by \cite{Russ:81}~:

\centerline{
\begin{tabular}{ll}
Direction & time step \\\hline
X & $\delta t$/2 \\
Y & $\delta t$/2 \\
Z & $\delta t$   \\
Y & $\delta t$/2 \\
X & $\delta t$/2 
\end{tabular}
}

For each step, the values of both mass and tracer mixing ratio are
updated following \eq{m1d} and \eq{q1d}.

Note that as for 1-Dimensional advection, a uniform tracer field is unchanged
by this sequence of operations, independently of the divergent
character of the wind field. This important property is apparently
not satisfied for all kind of time-splitting implementations
 \cite[see the discussion by][]{Lin:96}.

\subsubsection*{Advection over more than one mesh in one time-step}

The refinement of the longitudinal
discretization near the poles for longitude-latitude grids requires
smaller and smaller time steps if one wants to satisfy the criterion
that $U<m$.
One possibility is to use an additional time splitting in the $x$ direction
near the pole \cite[see \eg][]{Alle:91}.

Here, we use the fact that advection can be computed exactly if the mass
transfer $U$ is exactly equal to the mass of the upstream mesh.
For $U_{i+\dem}>0$ to fix idea, and if $m_i<U_{i+\dem}< m_i+m_{i-1}$
(which means that all the air of mesh $i$ passes through interface
$i+\dem$ in one time-step)
then the flux $F_i$ is the sum of the tracer amount in the upstream box
$m_i q_i$ and of the tracer advected form the mesh $i-1$ with a mass flux
$U_{i+\dem}-m_i$.
More generally, for
\begin{equation}\label{eq:shift}
\sum_{j=i-n+1}^i m_j < U_{i+\dem}\leq \sum_{j=i-n}^i m_j
\end{equation}
we use
\begin{equation}
\begin{eqalign}
F_{i+\dem}&=\sum_{j=i-n+1}^i m_j q_j\\
&+ \dep{U_{i+\dem}-\sum_{j=i-n+1}^i m_j}
\bq_{i-n+\dem}
\end{eqalign}
\end{equation}
with
\begin{equation}
\begin{eqalign}
\bq_{i-n+\dem}&=q_{i-n}\\
&+\frac{1}{2}\dep{1-\frac{U_{i+\dem}-\sum_{j=i-n+1}^i m_j}{m_{i-n}}} \dq_i
\end{eqalign}
\end{equation}

This approach is analogous to the flux-form Semi-Lagrangian method
proposed by \cite{Lin:96} except that we conserve
the semi-Lagrangian approach in the $x$ direction only
to treat the problem of longitudinal grid refinement near the poles.
For the basic resolution of the LMD GCM, with a grid of about 
100~000 points and a time step of 15~min for advection, the condition
$n>1$ is encountered about several tens of times per time step which
makes the cost negligible.

\subsubsection*{Polar discretization}

In the LMD GCM, the poles correspond to a mesh center, the 
horizontal mesh being twice smaller than current meshes in the
meridional direction.
The model uses a unique value of the mass and other scalar
fields at the pole
and the time evolution of the mass is computed from the integrated
meridional mass fluxes over all the meshes around the pole.

The same scheme is retained for the tracer advection. However, no way
of insuring strict monotonicity at the pole was found except by using
zero slopes. Thus the scheme retained for polar points 
is equivalent to Godunov's.

A numerical test of the solid rotation advection of a cosine distribution 
through the poles is shown in \fg{pole}.
The lower panels show the results of the advection after one full revolution
along a meridian
(for the exact solution it also corresponds to the initial distribution).
The shape of the advection with Van Leer scheme I is elongated in the
$y$ direction because of the numerical diffusion along the
flow. At the same resolution, Prather's scheme is of course less diffusive
but the shape of the iso-lines is slightly affected by the transpolar
advection. 
The last panels on the right show an application of Van Leer scheme I
with a much larger time-step for which the shifting number (\eq{shift})
reaches a value of $n=8$ at some points (case b).
 Note that the results are somewhat better
than on the other Van Leer test (case a) because of the smaller number of time
steps (160 against 16000 for the other cases).

\def\rayon{R}
\def\ri{r}



\begin{figure}
% \centerline{\includegraphics[width=16cm,viewport=0 80 560 390]{\EPS/pole.eps}}
\centerline{\includegraphics[width=16cm]{\EPS/pole.eps}}
\caption{Numerical test of transpolar advection.
Following Williamson and Rasch (1989) and Allen et al. (1991)
the initial tracer distribution is given by
$q(\lambda,\phi)=(1+\cos[\mbox{min}(\ri[\lambda,\phi]/\rayon,1)])/2$
with
$\ri=\arccos(\cos\lambda\cos\phi)$ and $\rayon=7\times(2\pi)/128$.
The discretization is based on a 128$\times$64 grid.
We show from the left to right, the exact advection result, a Van
Leer I test with 16000 time-steps for a full rotation, a Prather's test with
the same characteristics and a Van Leer I simulation with only 160 time-steps.
The lower panels correspond to a complete rotation (it is also the initial
condition for the exact advection). The upper panels show the tracer
distribution just after the first trans-polar crossing. 
\label{fg:pole}}
\end{figure}