 \def\EPS{\local/FIGURES}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le transport grande échelle\label{sec:volumesfinis}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\der#1#2{\frac{\partial #1}{\partial #2}}
\def\dep#1{\left(#1\right)}
\def\depb#1{\left[#1\right]}
\def\dem{1/2}
\def\dq{{\dep{\delta \q}}}
\def\dt{\delta t}


\subsection{Les différentes approches}

De nombreuses méthodes ont été développées au cours des dernières décennies
pour représenter l'advection, c'est à dire les parties non turbulentes
des Eqs~\ref{eq:direct12}
ou \ref{eq:direct11} formellement équivalentes
aux Eqs~\ref{eq:adv1} et \ref{eq:adv2}.

On distingue notamment~:
\begin{itemize}
\item{\bf Les méthodes lagrangiennes:} elles consistent à ``ensemencer" le
domaine d'étude (où la région source) avec des particules et à les
suivre individuellement.
Ces méthodes peuvent être très précises mais il faut généralement
utiliser un grand nombre de particules, notamment pour éviter que des
trous se créent dans la distribution spatiale de ces particules.
\item{\bf Les méthodes semi-Lagrangiennes:} (ou méthodes des caractéristiques
dans le monde des sciences de l'ingénieur)
on suppose qu'à un pas de temps
$t$ donné, $q$ est connu sur une grille. On prédit la valeur sur le même
maillage au temps $t+\delta t$ en remontant à $t$ ``le long du vent":
$q(\X,t+\delta t)=q(\X-\delta t \V, t)$. Le problème est que, si le
point $\X$ se trouve sur le maillage, le point
$\X-\delta t\V$ a peu de chances de s'y trouver.
On interpole donc cette valeur à partir des valeurs sur le maillage au
temps $t$. Cette méthode évite le problème principal des méthodes
lagrangiennes. Sa précision est déterminée par la précision des routines
d'interpolation.
Les méthodes semi-Lagrangiennes ont été très largement utilisées dans la
communauté atmosphérique. Elles présentent un défaut important: il est
difficile de garantir la conservation de la quantité totale de traceur.
\item{\bf L'advection des contours:} elle consiste à advecter des particules
uniquement sur des ``iso-niveaux" du champs considéré.
Cette dernière méthode, potentiellement très précise, se couple
difficilement à d'autres processus (il faut alors calculer l'effet des
autres processus sur la topologie des ``iso-niveaux").
\item{\bf Les différences finies~:}
il existe également tout un tas de méthodes en
différences finies \cite[cf. par exemple][]{Rood:87}.
Les différences finies permettent de dériver des schémas à des ordres élevés
mais ne garantissent pas à priori la conservation du traceur ou la positivité.
\item{\bf Les volumes finis~:} ils sont de plus en plus utilisés dans la 
communauté. Basés sur l'écriture de bilans sur des volumes de contrôle,
ils sont a priori appropriés pour la modélisation de processus
physiques couplés et complexes pour laquelle le respect des lois de 
conservation s'avère souvent essentiel. C'est cette approche qui
a été suivie pour introduire les traceurs dans LMDZ.
\item{\bf Les éléments finis~:} basés sur une décomposition sur des
fonctions de bases analytiques, ils sont encore relativement peu
utilisés dans la communauté.
\end{itemize}
\medskip

\subsection{Les schémas en volumes finis}
\footnote{Ceux qui rechercheraient
une présentation plus systématique et mathématique peuvent se reporter à
\cite{Roux:02}.}
\def\UU{{\cal U}}
\def\FF{{\cal F}}
L'intégration de l'équation de transport (\ref{eq:adv2})
et de l'équation de conservation de la masse (\ref{eq:masse})
sur un polyèdre (volume de contrôle) à $N$ faces conduit aux formulations:
\begin{equation}
\der{m}{t}=\sum_{n=1}^N {\UU}_n  \label{eq:m}
\end{equation}
et
\begin{equation}
\der{\q m}{t}=\sum_{n=1}^N {\FF}_n \label{eq:q}
\end{equation}
où $m$ est la masse totale d'air dans le polyèdre, $\q$ est la
concentration massique moyenne du traceur dans le volume,
${\UU}_n$ et ${\FF}_n$ sont les flux de masse et flux de traceur vers
l'intérieur du polyèdre pour la face $n$.

Les formulations en ``volumes finis" présentent l'avantage très important
d'être conservatives par nature, dès lors que le même flux est utilisé
pour les volumes en amont et en aval de l'interface considérée.

Pour le plus connu de ces schémas, introduit par \cite{Godu:59},
${\FF}$ est simplement estimé comme le produit de $\UU$ par la valeur de
$\q$ dans le volume amont (dans la maille d'où l'air provient).
Ce schéma simple, souvent appelé schéma amont, 
garantit
\begin{enumerate}
\item la conservation de la quantité totale de traceur,
\item la positivité
(un traceur positif partout reste positif),
\item
la monotonie (pour l'advection en dimension 1, une distribution monotone reste monotone),
\item la ``décroissance de la variation totale" du champ de traceur,
c'est à dire la décroissance de la somme des écarts absolus entre deux
concentrations successives.
Cette propriété garantit la stabilité numérique du schéma.
\item l'invariance par addition d'une constante au champ de traceur
(le modèle de Godunov est en fait linéaire par rapport au champ de traceurs, ce qui
garantit en plus la symétrie temporelle du transport comme on l'explique
au \ch{retro}).
\end{enumerate}
Cependant, ces propriétés physiques fondamentales ne sont obtenues
dans le schéma amont qu'au 
prix d'une très grande diffusion numérique.

\cite{VanL:77} a proposé d'utiliser, pour la valeur amont de $\q$, non pas
la valeur moyenne dans la maille amont mais une valeur extrapolée
à la frontière de celle-ci, en utilisant une approximation
polynomiale de la distribution sous-maille du traceur dans la maille
amont.
Deux des schémas proposés par Van Leer en 1977 ont été introduits
apparemment indépendamment dans la littérature météorologique par
\cite{Russ:81} -- le schéma III dans la classification de Van Leer --
et \cite{Prat:86} -- le schéma VI --.
Alors que dans son article original, Van Leer conclut que la complexité
du schéma VI est trop grande par rapport au gain en précision, ce schéma
est devenu une référence dans la communauté météorologique.

Les principes de dérivation de ces différents schémas sont exposés
ci-dessous.


\subsection{Description des schémas en dimension un\label{sec:general}}

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


En dimension 1 et après intégration sur un pas de temps, les 
\eq{m} et \ref{eq:q} s'écrivent simplement
\begin{equation}
\label{eq:dm}
\dt (m)_i=U_{i-\dem}-U_{i+\dem}
\end{equation}
et
\begin{equation}
\label{eq:dqm}
\dt(\q m)_i=F_{i-\dem}-F_{i+\dem}
\end{equation}
où les indices correspondent à la position sur l'axe des $x$,
 $\dt$ représente la différence finie en temps et
$U_{i+\dem}$ et $F_{i+\dem}$ sont les transferts de masse d'air
et de traceur à travers l'interface
$i+\dem$ durant un pas de temps.

Si on connaît complètement la distribution spatio-temporelle du vent, de la
densité de l'air et de la concentration de traceur, on a
\begin{equation}
m_i= \int_{\Delta x} \rho(x,t) dx
\end{equation}
\begin{equation}
\q_i= \int_{\Delta x} \q\rho(x,t) dx/\int_{\Delta x} \rho(x,t) dx
\end{equation}
\begin{equation}
U_{i+1/2}= \int_{t}^{t+\Delta t} \rho(x_{i+1/2},t) u(x_{i+1/2},t) dt
\end{equation}
et
\begin{equation}
F_{i+1/2}= \int_{t}^{t+\Delta t} \rho(x_{i+1/2},t) \q(x_{i+1/2},t) u(x_{i+1/2},t) dt
\end{equation}

Dans la présentation faite ici, on considère que les transferts
de masse $U_{i+1/2}$ sont connus à chaque instant et sur tout le maillage.
L'évolution temporelle de la masse d'air
$m_i$
est donc également complètement déterminée.
La seule chose qui reste à estimer est le flux de traceur.

De façon générale, ce flux $F_{i+\dem}$ peut être écrit comme le produit
du transfert de masse $U_{i+\dem}$ par la valeur moyenne de la concentration
du traceur $\bq_{i+\dem}$ dans l'air qui traverse l'interface au cours du pas
de temps.

\begin{figure}
\centerline{\includegraphics[width=16cm]{\EPS/schema2.eps}}
\caption{Principe du schéma de Van Leer et notations. On montre le cas
où la distribution sous-maille est représentée au moyen
d'un polynôme du premier degré. L'axe vertical correspond à la concentration
massique du traceur. Les deux axes horizontaux correspondent à l'indiçage
des variables et à la masse d'air comptée à partir du centre de la maille
$i$ et normalisée par la masse totale de cette même maille, $m_i$.
La surface grisée correspond à la quantité de traceur qui est transféré
au travers de l'interface durant un pas de temps.
\label{fg:schema2b}}
\end{figure}

La méthode proposée par Van Leer consiste à approximer la distribution
sous-maille par un polynôme pour lequel le calcul de  $\bq_{i+\dem}$ -- et
donc du transfert de traceur -- peut être fait exactement.
Le principe et les notations sont illustrées sur la \fig{schema2b}
pour le cas d'un polynôme du premier degré.
Avec ces notations, les valeurs moyennes de $\q$ et $m$ après un pas de temps
d'advection s'écrivent simplement
\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}
et
\begin{equation}\label{eq:m1d}
m^*_i=m_i+U_{i-\dem}-U_{i+\dem}
\end{equation}

Pour des valeurs positives de $U_{i-\dem}$ et $U_{i+\dem}$, 
l'\eq{q1d} peut se récrire
\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}
où $\bq_i$ est la concentration massique moyenne dans l'air
qui reste dans le volume $i$ pendant le pas de temps
(cf. \fig{schema2b}).
Avec ces notations, le schéma de Van Leer peut être décrit comme la 
séquence d'opérations suivante~:
on commence tout d'abord par définir 
une nouvelle distribution $\bq$ (les croix sur la \fig{schema2b})
à partir de la distribution initiale $\q$ (les points) en utilisant
une approximation polynomiale~: on advecte ensuite cette nouvelle
distribution approchée à travers les mailles.

Un des points essentiels de l'article original de \cite{VanL:77} est
de remarquer qu'on peut garantir la monotonie du schéma d'advection
-- en dimension 1, une distribution monotone reste monotone après advection --
en imposant que {\bf $\bq$ prenne une valeur intermédiaire entre les
deux valeurs voisines de $\q$ dans les régions où $\q$ est monotone}.


Supposons pour fixer les idées que $\q_{i-1}\leq \q_i \leq \q_{i+1}$ et $U\geq 0$.
La condition ci-dessus s'écrit simplement
\begin{equation}
\label{eq:monotone1}
\q_j \leq \bq_{j+\dem} \leq \q_{j+1} \mbox{ pour } j=i-1 \mbox{ et } j=i
\end{equation}
et
\begin{equation}
\label{eq:monotone2}
\q_{i-1} \leq \bq_i \leq \q_i.
\end{equation}
En introduisant la condition~\ref{eq:monotone1} dans l'\eq{q1d},
on assure que $\q^*_i \leq \q_{i}$ tandis que l'\eq{q1db} et les conditions
\ref{eq:monotone1} et \ref{eq:monotone2} garantissent que $\q_{i-1}\leq \q^*_i$.
On voit finalement que la condition suffisante énoncée ci-dessus garantit que
{\bf $\q^*_i$ est compris entre
$\q_{i-1}$ et $\q_i$ (resp. $\q_i$ et $\q_{i+1}$) pour $U>0$ (resp. $U<0$)}.
Or, comme le remarque \cite{VanL:77}, cette proposition implique la monotonie.

Si on rajoute en plus la condition que 
$\bq_i=\q_i$ (distribution sous-maille constante)
quand $\q_i$ est un extremum, on interdit la croissance des extrema.
Ceci implique de fait la {\bf positivité} du schéma et interdit la création
d'oscillations provenant du caractère dispersif du schéma numérique
\cite[provenant de l'advection avec des vitesses différentes des
différentes composantes de fourier de la distribution, cf. \eg\ ][]{Rood:87}.

Noter aussi qu'avec les définitions ci-dessus, une distribution uniforme
de traceur sera inchangée par l'advection, même avec des champs de vents
divergents (on s'en convainc en remplaçant  $\bq_{i-\dem}$ et
 $\bq_{i+\dem}$ par une valeur constante dans l'\eq{q1d}).

Dans les dérivations ci-dessus, on a supposé implicitement que le pas
de temps était suffisamment petit pour éviter de transférer plus que la
maille d'une cellule en un pas de temps (nombre de Courant $U/m$ inférieur
à 1). Si le flux de masse est exactement égal à la masse d'air dans le
volume amont, on calculera la distribution de traceur exactement, quelque
soit l'approximation choisie pour la distribution sous-maille.


\subsection{Le schéma de Godunov\label{sec:Godunov}}

La première approximation consiste à supposer que le traceur est constant
dans chaque maille (polynôme de degré zéro).
$\bq$ est alors simplement la valeur de $\q$ dans la maille amont
($\bq_{i+\dem}=\q_i$ si $U_{i+\dem}>0$ et $\q_{i+1}$ sinon).
Ce schéma proposé à l'origine par \cite{Godu:59} est bon marché, positif
et monotone mais au prix d'une diffusion numérique très forte.
La diffusivité du schéma peut d'ailleurs être quantifiée.
Dans le cas d'un champ de vent non divergent
($u=U\delta x/(m\delta t)=$cste à une dimension),
le schéma \ref{eq:q1d} s'écrit
\begin{equation}
\q^*_i=\q_i +\frac{U}{m}\q_{i-1}-\frac{U}{m}\q_{i}
\end{equation}
et peut se récrire
\begin{eqnarray}
\frac{\q^*_i-\q_i}{\delta t} &=&
u\frac{\q_{i-1}-\q_{i}}{\delta x}\\
&=& u \frac{\q_{i-1}-\q_{i+1}}{2\delta x}+\frac{u\delta x}{2}
\frac{\q_{i-1}-2\q_i+\q_{i+1}}{{\delta x}^2}
\end{eqnarray}
comme la somme d'un schéma numérique centré d'ordre 2
et d'un terme de diffusion numérique avec une diffusivité $u\delta x/2$.
\footnote{Le dernier terme de l'équation est une approximation numérique
de $u\delta x/2\ \partial^2\q/\partial x^2$.}


\subsection{Schémas du second ordre}

Pour passer à un ordre supérieur, on
suppose que la distribution sous-maille est linéaire avec une pente
$\dq_i$ donnant les valeurs aux bords gauche et droit de la maille $i$
comme $\q_i^\pm=\q_i\pm\dq_i/2$ 
(voir l'illustration de la \fig{schema2b}).

Dans ce cas, $\bq_{i+\dem}$ est donné par
\begin{eqnarray}
\bq_{i+1/2} & = & \q_i + \frac{1}{2}\dep{1-\frac{U_{i+\dem}}{m_{i}}} \dq_i \mbox{\ si\ } U_{i+\dem}>0\\
& = & \q_{i+1} - \frac{1}{2}\dep{1+\frac{U_{i+\dem}}{m_{i+1}}} \dq_{i+1} \mbox{ , sinon.}
\end{eqnarray}

Différents schémas ont été proposés \cite[]{VanL:77,VanL:79}
correspondant à différentes estimations de  $\dq$, dont deux schémas
particulièrement intéressants, l'un pour son faible coût et l'autre pour
sa précision.

\begin{figure}
\centerline{\includegraphics[width=9.cm]{\EPS/schema3.eps}
\includegraphics[width=7.cm]{\EPS/1dn.eps}}
\caption{Schémas du 2nd ordre. 
Illustration de l'estimation de la pente par différences
finies pour le schéma I de Van Leer ({\bf a}) et par calcul
par moindres carrés à partir de la distribution en ligne brisée résultant
de l'advection au pas précédent ({\bf b}). Cette seconde estimation
correspond  au schéma des pentes de
\cite{Russ:81}
ou au schéma III de Van Leer
(se reporter au texte pour plus de détails).
{\bf c :} 
Exemples de calculs
d'advection unidimensionnelle sur un domaine périodique de 70 points
(axe horizontal), avec une
vitesse constante $u$ et trois distributions initiales.
Les concentrations de traceur ($c$, unité arbitraire)
sont montrées
après une révolution complète, au pas de temps 350 pour un nombre de Courant
$U/m=0,2$.
\label{fg:1}}
\end{figure}

Dans le schéma I de Van Leer, la pente est simplement estimée à chaque pas
de temps par différences finies
($\delta \q_i=\dep{\q_{i+1}-\q_{i-1}}/2$) comme l'illustre
le graphique {\bf a} de la \fig{1}.
Ce schéma peut être considéré comme une version volumes-finis du schéma de
\cite{From:68}.

Dans le second schéma, la pente au pas de temps $t+\delta t$ dans une maille
donnée est calculée à partir de la distribution en ligne brisée résultant de
l'advection au temps $t$
(illustration sur le graphique {\bf b} de la \fig{1}).
La nouvelle distribution sous-maille 
minimise la distance quadratique par rapport à cette distribution.
Ce second schéma (le schéma III dans l'article original de Van Leer)
a en fait été redécouvert quelques années plus tard par
\cite{Russ:81}\footnote{La mise à jour de la pente, Eq.~(11) de \cite{VanL:77},
est la restriction exacte au cas unidimensionnel et non divergent de
l'Eq~(23) de \cite{Russ:81}.} 
dans un contexte météorologique tridimensionnel.
Ce schéma a été popularisé dans le modèle de circulation générale
du NASA/GISS. On l'appelle souvent
schéma des pentes (``slopes scheme" en anglais). On retiendra ce nom dans ce
qui suit.

On présente, sur le graphique {\bf c} de la
\fig{1}, un exemple typique de calcul d'advection unidimensionnel avec
le schéma de Godunov, Van Leer I et le schéma des pentes.
Le schéma des pentes est bien sûr le plus précis mais la différence principale
se situe entre le schéma de Godunov, extrêmement diffusif, et les deux
autres. On voit également apparaître pour les deux
schémas d'ordre 2 des oscillations provenant de la non monotonie du schéma et la
création de valeurs négatives.

\subsection{Limiteurs de pentes\label{sec:limiteurs}}

L'un des intérêts principaux du travail de \cite{VanL:77} réside
dans la possibilité qu'il offre d'assurer facilement la monotonie des
schémas. Une première façon consiste à appliquer directement la contrainte
explicitée à la fin de la \sec{general}.

Pour les schémas du second ordre,
\cite{VanL:77} propose aussi une condition suffisante à  la fois plus
brutale et plus simple. Il suffit, pour que le schéma soit monotone, d'imposer
que la distribution dans une maille $i$ soit entièrement comprise entre 
les valeurs moyennes des deux mailles adjacentes et que la pente soit
du même signe que dans ces deux mailles (à noter que ce dernier critère
est automatiquement satisfait par le schéma I).
Cette condition suffisante peut s'exprimer facilement comme un
limiteur de pente.
Pour le schéma I, la formulation complète du calcul de la pente avec limiteur
s'écrit simplement
\begin{eqnarray} \label{eq:slope0}
\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}
si $\q_i$ est compris entre $\q_{i-1}$ et $\q_{i+1}$ et 0 sinon.

\begin{figure}
\centerline{\includegraphics[width=15.cm]{\EPS/schema4b.eps}}
\caption{{\bf a)} Illustration de l'application d'un limiteur de pente.
A droite, impact des limiteurs de pentes sur le schéma des
pentes ({\bf b}) et sur le schéma I de Van Leer ({\bf c}).
Les cas sans limiteurs (les mêmes que sur la \fig{1}) et les limiteurs
faibles et forts (se reporter au texte pour plus de détails) sont
montrés pour les distributions carrée et gaussienne.
\label{fg:limit}}
\end{figure}

L'effet de l'application d'un limiteur fort (décrit ici) ou d'un limiteur
faible (application directe des \eq{monotone1} et \eq{monotone2}),
qui correspondent respectivement aux Eq. (66) et (74) données
par \cite{VanL:77} est illustré sur la \fig{limit}{\bf a}.
Un test numérique de l'impact de ces limiteurs est montré sur les
\fig{limit}{\bf c} et {\bf d}.
Le limiteur fort dégrade de façon significative la précision du schéma
des pentes tandis que le limiteur faible corrige ce schéma de façon
sélective au niveau des oscillations, sans émousser par exemple le sommet
de la gaussienne. Pour le schéma I, la différence entre les deux limiteurs
est nettement plus marginale.


Il existe une alternative élégante à l'\eq{slope0} qui consiste à 
utiliser la moyenne géométrique des deux pentes voisines
\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}
quand $\dep{\q_{i+1}-\q_i}\dep{\q_i-\q_{i-1}}>0$ et 0 sinon.
Cette formulation satisfait automatiquement la limitation forte.
Ce schéma est seulement légèrement plus diffusif que le schéma I mais le gain
en temps est également relativement faible, le calcul de la pente
ne représentant qu'une partie du coût du schéma.

\subsection{Schémas du troisième ordre}

La description de la distribution sous-maille du traceur peut
être encore améliorée en utilisant un polynôme du second ordre.
L'équivalent du schéma I de Van Leer (le schéma IV dans l'article original)
consiste à évaluer les coefficients de ce polynôme par différences
finies à partir des valeurs moyennes dans les mailles.
Van Leer a attaché peu d'importance à ce schéma qui présente par exemple
comme défaut de n'être pas plus précis que le schéma I pour un nombre
de Courant de 0,5. Cependant, pour des valeurs plus petites, il peut
devenir nettement plus précis.
Comme le schéma I, ce schéma présente l'intérêt de n'utiliser qu'un 
traceur par point de grille.
\cite{Wood:81} ont développé une alternative au schéma IV en utilisant
un filtre non-linéaire assez élaboré permettant de renforcer les pentes
en cas de choc \cite[se reporter également à][]{Wood:84,Cole:84}.
La méthode qui en résulte a été baptisée Piecewise Parabolic Method.
Elle a été utilisée pour les applications atmosphériques par exemple par
\cite{Carp:90}, \cite{Lin:96}  ou \cite{Vaut:01}.
Dans ce manuscrit on ne présente pas les résultats avec PPM.
Mais des tests ont été réalisés récemment avec ce schéma qui se comporte
de façon assez similaire au schéma des pentes pour un coût numérique
intermédiaire entre les schémas I et III (pentes) en termes de temps de
calcul. Le coût en stockage est le même que pour le schéma I.

L'équivalent à l'ordre 3 du schéma III de Van Leer est le schéma VI 
connu dans la littérature météorologique sous le nom de schéma
de \cite{Prat:86}.
Le schéma de Prather est évidemment beaucoup plus précis que tous les
schémas présentés jusque là comme le montrent les illustrations de la 
section suivante.
Cependant, il nécessite la conservation de 10 traceurs indépendants~:
la moyenne, les 3 pentes dans les 3 directions d'espace $\q_x$, $\q_y$,
$\q_z$ et les 6 moments du second ordre $\q_{xx}$, $\q_{yy}$, $\q_{zz}$,
$\q_{xy}$, $\q_{yz}$ et $\q_{xz}$.
Il est également nécessaire de lui adjoindre des algorithmes du type limiteurs
de pente pour éviter complètement les oscillations numériques ou les valeurs
négatives.


\subsection{Introduction des schémas dans LMDZ}

L'introduction des schémas en volumes finis est relativement facile
dans le modèle de circulation générale LMDZ car les flux de masse
sont déjà définis sur une grille décalée.
Ces flux de masse sont ceux utilisés dans la partie météorologique pour
intégrer l'équation de continuité
pour l'air, à savoir l'\eq{dm} ou l'\eq{cont} donnée plus loin.

\subsubsection*{Flux alternés}

\begin{figure}
\centerline{\includegraphics[height=8cm]{\lmdztdir/schemsplit.eps}
            \includegraphics[height=9cm]{\lmdztdir/split.eps}}
\caption{Advection uniforme d'un pic gaussien
de concentration le long de la diagonale d'un
maillage bidimensionnel régulier.
En haut, le calcul de la divergence des flux est effectué
à partir de flux estimés indépendamment en $x$ et en $y$.
La ligne du bas montre un calcul alterné avec d'abord une advection
en $x$ puis une advection en $y$.
Ces illustrations numériques ont été réalisées avec le schéma I de Van Leer,
sur un maillage de 60 points dans chaque direction horizontale.
\label{fg:split}}
\end{figure}

Les présentations ci-dessus ont été faites avec une seule dimension d'espace.
Une solution naturelle pour passer en 3 dimensions consiste à calculer
d'abord les flux dans chaque dimension puis à calculer la divergence
des flux.

Une alternative classique
consiste à faire trois calculs d'advection, successivement
dans les trois directions.
On parlera de flux alternés ou de directions alternées.

Cette seconde méthode est curieusement plus précise dans le cas d'une advection
en diagonale par rapport au maillage.
Ceci est illustré sur la \fig{split} pour une configuration bidimensionnelle. 
Si on considère un pas de temps d'advection, le calcul direct ne tient
pas compte de l'air passant directement dans la maille située en
diagonale par rapport à la maille du milieu (le petit carré sombre
sur les figures). En revanche, cet air est compté deux fois, une fois dans
l'advection en $x$ et une fois dans l'advection en $y$.
Ceci produit une très forte diffusion latéral en réduisant la propagation
dans la direction du mouvement. En partant d'un panache gaussien, la
distribution obtenue s'étire dans la direction transverse à l'écoulement
(en haut à droite sur la figure).

Au contraire, pour l'advection en flux alternés, le carré sombre commence
par passer dans la maille située à droite de la maille d'origine puis
en haut et arrive bien finalement dans la maille située en diagonale
par rapport à la maille d'origine.
Si on pouvait à chaque instant connaître exactement la distribution sous-maille
du traceur, on aurait donc un calcul exact.

Dans la méthode des flux alternés, on intègre successivement non seulement
l'équation d'advection du traceur (\eq{q1d}) mais aussi l'advection de transport
de l'air (\eq{m1d})\footnote{
Une certaine confusion sur ce point est entretenue dans la littérature
\cite[cf. \eg][]{Carp:90} dans laquelle ont présente souvent
les flux alternés en ne découpant que l'équation de transport
des traceurs puis en introduisant comme des astuces numériques les corrections
nécessaires qui, si elle sont bien choisies, consistent simplement à découper
simultanément l'équation de conservation pour l'air.
}.
Avec cette approche, une distribution uniforme de traceur est inchangée
par l'advection indépendamment du caractère divergent ou non du
champ de vent \cite[][montrent que ce n'est pas garanti par toutes
les formulations en flux alternés]{Lin:96}.

Cette approche en flux alternés est utilisée assez systématiquement
dans les calculs d'advection en volumes finis
\cite[]{VanL:79,Alle:91,Russ:81,Prat:86}.

On utilise ici la séquence proposée par \cite{Russ:81}~:

\centerline{
\begin{tabular}{ll}
Direction & pas de temps \\\hline
X & $\delta t$/2 \\
Y & $\delta t$/2 \\
Z & $\delta t$   \\
Y & $\delta t$/2 \\
X & $\delta t$/2 
\end{tabular}
}

\subsubsection*{Sauter des mailles}

Si on veut ne travailler qu'avec des nombres de Courant
($\delta t\ u/\delta x$ ou $U/m$) inférieurs à 1,
le raffinement de la discrétisation près des pôles, inhérent aux
grilles globales longitudes-latitudes, nécessite l'utilisation de pas de temps
extrêmement petits dans la direction longitudinale.
Une possibilité pour résoudre ce problème consiste à découper encore
davantage le pas de temps pour l'advection longitudinale.

Il existe une alternative utilisant le fait que le calcul est
exact pour un nombre de courant de 1. Dans ce cas en effet, le flux
est simplement le produit du flux de masse par la concentration
moyenne du traceur dans la maille amont.
Si le nombre de courant est plus grand que 1, 
(pour fixer les idées, si on prend $U_{i+\dem}>0$ et
$m_i<U_{i+\dem}< m_i+m_{i-1}$)
on peut calculer le transfert de traceur comme la somme de la quantité
de traceur dans la maille amont ($m_i \q_i$) et de la quantité
de traceur transférée depuis la maille $i-1$ avec un flux
de masse $U_{i+\dem}-m_i$.

De façon générale, pour
\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}
on prend
\begin{equation}
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{equation}
avec
\begin{equation}
\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{equation}

Cette façon de voir les schémas en volumes finis présente de fortes analogies
avec les approches semi-Lagrangiennes pour lesquelles ont prédit
la valeur ponctuelle en un point du maillage en remontant en amont
la trajectoire de la particule.
Cette analogie a conduit \cite{Lin:96} à introduire la notion
de flux-form Semi-Lagrangian method.
A noter qu'en parallèle, les
adeptes des schémas semi-Lagrangiens classiques modifient leurs
formulations pour garantir la conservation et les font ressembler de plus
en plus à des schémas en volumes finis \cite[\eg][]{Yabe:01}.



On pourrait de façon générale et pour un pas de temps aussi grand qu'on veut,
estimer les contours du maillage transportés à rebours sur un
pas de temps. A partir d'une estimation polynomiale des distributions sous
mailles, comme celles proposées par Van Leer, on pourrait calculer
alors la distribution de traceur sur ce maillage déformé par le transport
à rebours.
Mais on paie vite en complexité du schéma et en coût informatique (introduction
de branchements supplémentaires et gestion de la mémoire,
particulièrement pénalisant 
sur des ordinateurs vectoriels ou parallèles) ce qu'on gagne avec l'utilisation 
d'un pas de temps plus long. De plus, on ne tient pas forcément à utiliser
des pas de temps plus longs que les constantes caractéristiques des
autres processus (transport turbulent dans la couche limite, chimie, etc...).
On restreint donc ici le traitement des nombres de Courant plus grands
que 1 à la direction longitudinale pour traiter spécifiquement le problème
de raffinement du maillage en longitude à l'approche des pôles.
Pour une résolution typique du modèle LMDZT, avec une grille
d'environ 100~000 points et un pas de temps de 15 minutes,
la condition $U/m>1$  n'est rencontrée typiquement que quelques dizaines
de fois par pas de temps ce qui rend le coût marginal raisonnable.

\subsubsection*{Le point du pôle}

Dans le modèle du LMD, les pôles correspondent à des centres de mailles
triangulaires.
La dimension de ces mailles est 2 fois plus petite dans la direction
méridienne que celle des mailles normales.
Toutes ces mailles ne sont en fait pas indépendantes et on impose
que les valeurs scalaires (pression de surface, température, concentration
de traceurs) soient toutes identiques.
L'évolution de la masse d'air ou de la quantité d'un traceur au pôle est
estimée à partir de la convergence totale des flux méridiens.
En fait, on peut considérer le pôle comme un volume de contrôle consistant
en un polygone avec autant d'arêtes que le nombre de points longitudinaux 
de la grille.

Il semble impossible de garantir strictement la monotonie du schéma
en conservant des pentes non nulles pour ces mailles polaires.
On retient donc aux pôles un schéma de Godunov.

\subsubsection*{Advection transpolaire}

La \fig{pole} montre un exemple d'advection transpolaire d'une distribution
sinusoïdale avec un écoulement en rotation solide le long du méridien
de Greenwich.
Les figures du bas montrent le résultat de l'advection après une révolution
complète (la solution exacte, à gauche, est identique à la distribution
initiale). La distribution obtenue avec le schéma I~(cas a) est allongée dans la
direction méridienne à cause de la diffusion numérique, dans la
direction de l'écoulement.
Pour la même résolution, le schéma de Prather est beaucoup moins diffusif.
La forme des iso-lignes est cependant légèrement altérée par le passage au pôle.
Pour ces calculs, on s'est arrangé pour que le nombre de Courant soit toujours
inférieur à 1. 
Le dernier exemple sur la droite (cas b)
correspond à un calcul avec le schéma I
utilisant un pas de temps beaucoup plus grand, avec un nombre de Courant
qui atteint 8 dans la direction longitudinale près du pôle.
On remarque que les résultats sont plutôt meilleurs que ceux du cas {\bf a}
à cause du plus petit nombre de pas de temps nécessaire (160 au lieu
de 16000).

\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{Test numérique d'advection transpolaire.
Suivant \cite{Will:89} et \cite{Alle:91}, la distribution
initiale de traceur est donnée par
$\q(\lambda,\phi)=(1+\cos[\mbox{min}(\ri[\lambda,\phi]/\rayon,1)])/2$
avec
$\ri=\arccos(\cos\lambda\cos\phi)$ et $\rayon=7\times(2\pi)/128$.
La grille utilisée comprend 128 points en longitude et 64 en latitude.
On montre de gauche à droite, la solution exacte, un test du schéma I
de Van Leer avec 16000 pas de temps (pour avoir un nombre de Courant plus
petit que 1) pour une révolution complète,
un test du schéma de Prather avec le même pas de temps et enfin une
simulation avec le schéma I mais seulement 160 pas de temps.
Les graphiques du haut montrent la distribution de traceur
juste après le premier passage par un pôle.
Les graphiques du bas montrent le résultat obtenu après une révolution complète
autour du globe.
\label{fg:pole}}
\end{figure}

\subsection{Tests bidimensionnels}

Dans cette section, nous présentons des tests bidimensionnels
des schémas d'advection tels qu'ils sont codés dans le modèle
LMDZT. On teste le schéma de Godunov, le schéma I de Van Leer, le
schéma des pentes et le schéma de Prather.
Pour le schéma I, on utilise le limiteur fort. Pour les schémas de pentes
et de Prather, on assure seulement la positivité en utilisant un
limiteur de flux
\cite[suivant ][on se contente d'imposer de ne pas sortir plus de traceur d'une maille que ce qu'elle contient]{Prat:86}.
Les schémas de Pentes et de Prather que nous utilisons ont été optimisés
et interfacés avec le modèle ancien du LMD par Pascal Simon
et Christophe Genthon.

On insiste dans les tests présentés ci-dessous sur le rapport entre
précision et coût numérique. On montre en effet que l'arbitrage entre
une résolution plus fine, qui rend tous les schémas plus précis, et 
l'utilisation d'un schéma intrinsèquement plus précis, comme Prather,
n'est pas évident et peut dépendre du type de machine utilisé
ou  d'autres considérations relatives à d'autres composantes du modèle.



\subsubsection*{Coût numérique des différents schémas}

De façon générale, on peut gagner en précision en utilisant une résolution
spatiale plus fine, ce qui se fait évidemment au prix d'un coût numérique
plus grand.
Une question pratique importante en termes d'efficacité des schémas est de
savoir ce qu'on gagne en changeant de schéma d'advection à coût numérique
inchangé.
Puisque notre but est l'advection tridimensionnelle dans un modèle
de circulation générale, il faut d'abord se faire une idée du coût
relatif des différents schémas dans une telle configuration.

La comparaison la plus facile est celle de l'occupation en mémoire.
Les schémas de Godunov, Van Leer I et PPM ne nécessitent de conserver
qu'une variable indépendante par maille et par champ de traceur.

Les schémas des pentes et de Prather sont nettement plus coûteux
avec respectivement 4 et 10 variables  d'état pour décrire un traceur
physique.

La comparaison en termes de rapidité est moins évidente et peut dépendre
du type de machine utilisé.
Des test sur machines scalaires (stations SUN) et vectorielles montrent
qu'il y a typiquement un facteur 2 entre le schéma I et les pentes
et entre les pentes et Prather.

Les machines vectorielles favorisent de façon générale l'utilisation
de schémas plus grossier sur des grilles plus fines.
Pour les test présentés ci-dessous, effectués sur un CRAY-90,
les pentes ont exactement le même coût numérique pour une grille
horizontale de 60 par 43 points que le schéma I de Van Leer pour
une grille de 120 par 85 points.

On voit que le fait de doubler la résolution horizontale dans chaque direction
rend le schéma I équivalent à celui des pentes à la fois en termes
de stockage et en termes de rapidité sur une machine vectorielle.
De même, les pentes en résolution double sont comparables à Prather
sur une grille deux fois plus grossière.

De façon générale, les machines scalaires sont plus favorables
aux schémas plus précis.


\subsubsection*{Test avec solution exacte connue}

Nous présentons ici des tests numériques bidimensionnels effectués
avec un champ de vent analytique présentant une rotation différentielle
et pour lequel l'advection peut être calculée exactement.
Ce test est effectué sur la sphère en utilisant la discrétisation du
modèle de circulation, pour permettre de valider directement les
codes utilisés dans le modèle de circulation générale.

Pour un champ de vent horizontal non divergent dans le plan 
longitude-latitude ($\lambda$, $\phi$), les composantes zonale ($u$) et
méridienne ($v$) du vent satisfont la relation suivante~:
\begin{equation}
\frac{\partial u}{\partial \lambda}+\frac{\partial v\cos\phi}{\partial \phi}=0
\end{equation}
\def\ur{U_0}
Le champ de vent que nous utilisons dérive du potentiel suivant~:
\begin{equation}
\label{eq:psi}
\Psi= a \ur \cos^2\frac{\lambda}{2}\cos^2\phi
\end{equation}
(où $\ur$ est une vitesse caractéristique et $a$ est le rayon planétaire)
de sorte qu'on a
\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}}
En suivant une trajectoire (valeur constante de $\Psi$), la vitesse méridienne
$v=a \ddd\phi/\ddd t$ peut être récrite en combinant les équations
\ref{eq:psi} et \ref{eq:vpsi}~:
\begin{equation}
a \ddd\phi/\ddd t= -\sqrt{\Psi \ur/a}\mbox{ sign}\lambda\sqrt{1-\cos^2(\lambda/2)}
\end{equation}
Si on introduit $\alpha=\sqrt{\Psi\ur/a^3}$ et qu'on remplace
 $\cos^2(\lambda/2)$ selon  l'\eq{psi}, cette équation devient
\begin{equation}
-\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{equation}
Les trajectoires, iso-valeurs de $\Psi$, sont parcourues avec la loi
horaire donnée ci-dessus. Le temps mis par une particule pour revenir
à sa position de départ est donc
$2\pi/\alpha=2\pi a/(\ur\cos\phi\cos\lambda/2)$.
Ce temps est plus court au centre du domaine et infini sur les bords.

\subsubsection*{Résultats numériques}

\begin{figure}
% Viewport original de uvref 44 80 586 239
\centerline{\includegraphics[width=16cm,trim=0 10 0 0]{\EPS/uvn.eps}}
\caption{Tests d'advection bidimensionnelle avec un champ de vent analytique.
La distribution initiale de traceur est une fonction gaussienne de la longitude.
Les graphiques du haut montrent l'état initial (à gauche)
et la solution exacte aux instants $T/2$ et $T$.
Le temps $T$ correspond à une révolution complète au centre du domaine.
En dessous, on montre les résultats numériques au temps $T$ pour différents
schémas d'advection et 3 résolutions~: résolution pleine, 1/2 et 1/3,
correspondant respectivement à des grilles longitude-latitude de
$120\times 60$, $60 \times 30$ et $40\times 20$ points.
\label{fg:uv}}
\end{figure}

Sur le graphique en haut à gauche de la \fig{uv},
on montre à la fois la distribution
initiale du traceur, une gaussienne ne dépendant que de la longitude,
et le champ de vent (dans une unité arbitraire).
Toujours en haut, les figures à droite montrent l'évolution exacte
du traceur sous l'effet de l'advection.
La rotation est plus rapide au centre qu'aux bords du domaine,
ce qui produit cette forme 
en spirale et une filamentation.

Pour les tests numériques on utilise trois résolutions~:
$120\times 60$ (résolution pleine), $60 \times 30$ (résolution 1/2)
et  $40\times 20$ (résolution 1/3).
Le pas de temps
est choisi suffisamment petit pour que le nombre de Courant reste toujours
plus petit que 1 (le traitement spécial en longitude n'étant codé
que pour le schéma I de Van Leer).
L'état final montré sur les figures correspond à une révolution
complète au centre du maillage ce qui nécessite un nombre
de pas de temps de 4000, 6000 et 12000 pour les différentes résolutions
horizontales testées.

La reconstruction très fine des filaments (avec des valeurs maximum
dépassant 0,9 pour une valeur initiale de 1) avec le schéma de Prather
en pleine résolution est très impressionnante.
Cependant, le schéma des pentes fait plutôt mieux à cette résolution que 
Prather dans une résolution deux fois plus grossière dans chaque direction
horizontale.
On observe la même chose entre le schéma I de Van Leer et le schéma des pentes.

Donc, si on change la résolution dans seulement deux directions, le schéma
I se comporte plutôt mieux que les pentes en terme de rapport
qualité/coût. 

Si on change de résolution dans les trois directions, la comparaison
doit être  faite entre la pleine résolution et la résolution intermédiaire
avec un avantage significatif pour les schémas précis.

Il est enfin intéressant de noter la différence flagrante de
performances entre les schémas de Godunov et le schéma I en dépit d'un
coût équivalent en termes de stockage.

\subsection{Remarques pour conclure}

Les schémas en volumes finis proposés par \cite{VanL:77}
conduisent facilement à des mises en {\oe}uvre
tridimensionnelles qui satisfont des propriétés essentielles du
transport comme
\begin{itemize}
\item la conservation de la quantité totale de traceur,
\item la non modification d'une distribution constante de traceurs,
\item l'invariance par addition d'une constante au champ de traceur,
\item la monotonie,
\item la non création d'extrema d'origine numérique (d'où la
positivité),
\item la décroissance de la variation totale (TVD en anglais) --
somme des écarts absolus entre deux points consécutifs -- au cours
du temps qui assure la stabilité du schéma. 
On peut montrer que cette propriété (valable uniquement en dimension
1) est respectée à la fois par le schéma I de Van Leer et par le schéma
PPM \cite[cf. ][]{Roux:02}.
\end{itemize}


En pratique, on constate que les schémas plus sophistiqués se comportent
mieux, mais au pris d'un coût numérique additionnel du même ordre que celui
qu'aurait entraîné l'utilisation d'une grille plus fine.
A noter qu'il se peut que les schémas deviennent à partir d'un certain
stade moins diffusifs que l'atmosphère elle-même.
C'est particulièrement vrai pour la basse troposphère, dans laquelle
la turbulence de couche limite induit une très forte diffusion verticale,
qui, couplée à des cisaillements importants du vent horizontal, conduit aussi
à une forte dispersion horizontale effective.
Les tests présentés plus loin illustrent ce point.

Pour finir, il faut noter que nous avons présenté ici les schémas
dérivés à l'origine par Van Leer. Le schéma PPM qui a également
été testé dans LMDZ (résultats non présentés)
semble supérieur à coût numérique
égal (en tous cas en termes de stockage) à ces schémas d'origine.

L'ensemble des schémas décrits ci-dessus a été introduit dans le modèle
LMDZ. Le schéma I de Van Leer a été retenu en standard pour sa robustesse
et sa simplicité mais ce choix ne doit pas être considéré comme définitif
et doit être reconsidéré en fonction du problème abordé.

