\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_i1$ 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é.