%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le modèle LMDZ}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Avant d'introduire la composante traceurs, on
présente rapidement le modèle LMDZ et notamment la discrétisation
des équations de grande échelle sur laquelle s'appuie le transport
des traceurs.
On en profite pour donner un bref aperçu du contenu des
paramétrisations physiques du modèle et des évolutions apportées récemment.
\footnote{
Cette section, relativement technique, ne permettra pas à un non initié de
se familiariser avec le monde des modèles de circulation générale.}

LMDZ correspond à la seconde génération d'un modèle de climat développé
depuis une trentaine d'années au LMD
et décrit initialement par \cite{Sado:84}.
Ce modèle est plus modulaire et flexible\footnote{La flexibilité est ici une
notion positive !}
que son prédécesseur.
Le ``Z" de son nom se réfère à la capacité de raffinement de la grille
(Zoom) obtenue grâce à une écriture généralisée de la formulation numérique
avec un maillage dont les facteurs d'élongation dans les deux directions
horizontales peuvent être choisis arbitrairement.
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Description générale}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Comme la plupart des modèles de climat, le modèle LMDZ intègre sur la sphère
et dans le temps les ``équations primitives de la météorologie".
Ces équations sont une version des équations de Navier Stokes simplifiées
en supposant que l'atmosphère est à tout moment en équilibre hydrostatique
sur la verticale et en négligeant les variations verticales de la géométrie
horizontale (hypothèse de couche mince).
Le moment cinétique par rapport à l'axe des pôles est
par exemple calculé en utilisant comme distance à l'axe
$a \cos\phi$ -- où $a$ est le rayon de la planète
et $\phi$ la latitude -- plutôt que $(a+z) \cos\phi$, qui tiendrait compte de
l'altitude $z$ d'une particule d'air au dessus de la surface.

Le modèle du LMD est bâti sur une discrétisation en différences finies de ces
équations avec certaines propriétés intéressantes du schéma numérique comme
la conservation de la masse, la conservation du moment cinétique 
par la composante axi-symétrique de l'écoulement et de
la vorticité barotrope. On revient sur ces aspects en toute fin de chapitre.

Ces équations ne peuvent pas être intégrées jusqu'à l'échelle visqueuse.
En pratique, dans les modèles globaux, on utilise des mailles de quelques
dizaines à quelques centaines de kilomètres suivant les applications.
L'impact des échelles sous-mailles sur la grande échelle doit donc être
représenté au travers de paramétrisations.
Il faut également représenter des processus fondamentaux comme le transfert
de rayonnement visible et infrarouge dans l'atmosphère, les processus
nuageux ou les interactions avec la surface.
Dans le jargon des modélisateurs, ces paramétrisations sont regroupées
dans la ``physique" du modèle,
par opposition
avec le code représentant la ``dynamique" de grande échelle.
Dans le cas de l'utilisation du modèle de circulation pour d'autres planètes,
c'est cette partie physique qui doit être largement modifiée, et notamment
le calcul du transfert radiatif.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Discrétisation des équations primitives}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\def\m{m}
\def\pk{\Pi}

\def\rcp{\kappa}
\def\cp{C_p}
\def\uabs{\tilde{u}_{abs}}
\def\err{\epsilon}
\def\s{ {\dep{\frac{p}{\p0}}}^{\rcp} }
\def\dsig{\dz \sigma}
\def\ps{p_s}
\def\psk{{\ps}^\kappa}
\def\ucov{\tilde{u}}
\def\vcov{\tilde{v}}
\def\ucont{\tilde{\ucov}}
\def\vcont{\tilde{\vcov}}
\def\cu{c_u}
\def\cv{c_v}
\def\h{\theta}
\def\pext{\tilde{p}_s}
\def\fext{f}
\def\dx{\delta_x}
\def\dy{\delta_y}
\def\dz{\delta_z}
\def\K{\frac{1}{2}
\left( \av{\ucov \ucont}{X} + \av{\vcov \vcont}{Y} \right)}
\def\Z{\frac{\dx \vcov - \dy \ucov + \fext}{\av{m}{X,Y}}}
\def\Zm{\frac{- \dy \ucov + \fext}{\av{\pext}{Y}}}

\newcommand{\glob}[1]{ \left< #1 \right> }
\newcommand{\av}[2]{{\overline{#1}}^{ #2 }}
\newcommand{\avg}[1]{\left< #1 \right>}
% \newcommand{\emis}[2]{E \left( #1 , #2 \right) }
\newcommand{\der}[2]{\frac{\partial #1 }{\partial #2} }

\begin{figure}
\includegraphics[width=16cm]{\local/FIGURES/vargrille.eps}
\caption{Disposition des variables pour la grille du modèle LMDZ.\label{fg:grille1}}
\end{figure}


Dans le modèle LMDZ, les équations primitives sont discrétisées
horizontalement sur une grille-C dans la classification d'Arakawa
\cite[cf. \eg][]{Kasa:77}.
On note $X$ et $Y$ les coordonnées horizontales~:
$X$ (resp. $Y$) est une fonction biunivoque de la longitude $\lambda$
(resp. de la latitude $\phi$).
Les variables scalaires
(la température potentielle $\theta$, le géopotentiel $\Phi$
et la pression de surface $\ps$) sont évaluées aux points
correspondant à des couples de valeurs entières $(X,Y)=(i,j)$.
Le vent zonal est calculé aux points $(X,Y)=(i+1/2,j)$ et  le vent
méridien aux points $(X,Y)=(i,j+1/2)$.
La disposition des variables sur la grille est illustrée sur la
\fig{grille1}.
De façon à pouvoir modifier la distribution des longitudes et latitudes
de la grille, on utilise en fait les composantes covariantes
($\ucov$ et $\vcov$) et contravariantes ($\ucont$ et $\vcont$)
du vent définies par
\begin{equation}
\begin{array}{llllllllll}
\ucov = \cu u & \mbox{et} & \ucont = u / \cu & \mbox{avec} & 
\cu = a \cos{\phi} \left( d\lambda/dX \right)  \\
\vcov = \cv v & \mbox{et} & \vcont = v / \cv & \mbox{avec} &
\cv = a \left( d\phi / dY \right)
\end{array}
\end{equation}
%
où $u$ et $v$ sont les composantes physiques du vecteur vent
horizontal.

Comme beaucoup de modèles globaux, le modèle du LMD avait initialement
été codé avec une coordonnée verticale de type $\sigma$. La coordonnée
$\sigma=p/\ps$ (où $p$ est la pression et $p_s$ la pression à la surface
au point considéré) est pratique parce qu'elle varie de 1 à la surface à 0
au sommet de l'atmosphère quelque soit le relief
sous-jacent. Cette coordonnée devient cependant problématique plus haut dans
l'atmosphère où on préfère travailler sur des isobares ou, mieux, sur
des isentropes.
Suivant là aussi beaucoup de modèles de climat, nous avons avec Phu LeVan
récrit le modèle en coordonnée hybride $\sigma-p$.
La coordonnée hybride est définie de façon implicite
comme donnant la pression dans la couche $l$ du modèle sous la forme
\begin{equation}
\label{eq:pl}
p_l=A_l + B_l\ps
\end{equation}
On choisit près de la surface $A\sim 0$ et $B\sim 1$. Vers le haut
du modèle, on a $A\sim 0$ et $B\sim 0$ avec $A\gg B \ps$. On retrouve ainsi
la coordonnée $\sigma$ près de la surface et des
coordonnées pression plus haut dans l'atmosphère.

\def\aire{{\cal A}}

Les niveaux de pression définis par la relation~\ref{eq:pl}
correspondent aux interfaces entre les couches du modèle avec $p_1=\ps$
et $p_{N+1}=0$, $N$ étant le nombre de couches dans le modèle.
La masse d'air contenue dans une maille du modèle comprise entre les
niveaux $p_k$ et $p_{k+1}$ est donnée par
\begin{equation}
m_k=\aire\frac{p_k-p_{k+1}}{g}
\end{equation}
($\aire=\cu\cv$ est l'aire de la maille) ou encore
\begin{equation}
\m_k= \frac{\aire}{g} \depb{A_k-A_{k+1} +\dep{B_k-B_{k+1}} \ps}
\end{equation}


On introduit alors les trois composantes du flux de masse~:
\begin{equation}
U=\av{\m}{X} \ucont ,\  V= \av{\m}{Y} \vcont \  \mbox{et} \  
W
\end{equation}
où la notation $\av{a}{X}$ signifie qu'on prend la moyenne arithmétique
de la quantité $a$ suivant la direction $X$ et 
où le flux de masse vertical $W$ est défini à partir de
l'équation de continuité
\begin{equation} \label{eq:cont}
\dt{\m} + \dx U+ \dy V +\dz W=0
\end{equation}
la notation $\delta_X$ correspondant à la différence entre
deux points consécutifs suivant la direction $X$.\footnote{
Pour résoudre cette équation, on commence
par calculer la convergence de masse $\Omega_k$
cumulée depuis le sommet de l'atmosphère
jusqu'au niveau $k$~:
\begin{equation} 
\Omega_k= \sum_{l=k}^N \dep{\dx U+ \dy V}
\end{equation}
La convergence au sol $\Omega_1$ donne accès à l'évolution de la pression
de surface
\begin{equation}
\frac{\aire}{g}\dt{\ps}=\Omega_1
\end{equation}
 Finalement, en remarquant que
\begin{equation}
\dt{m}=-\frac{\aire}{g}\dz B \dt{\ps}=-\dz B \Omega_1
\end{equation}
on obtient le flux de masse vertical $W$ par intégration
depuis le haut de l'atmosphère
\begin{equation}
- \dz B \Omega_1 -\dz \Omega +\dz W=0
\end{equation}
}

On ne présente pas ici la discrétisation des équations du mouvements
qui n'a pas changé par rapport à la version $\sigma$ du modèle
et dont ont discute en toute fin de chapitre.
En revanche, l'introduction de la coordonnée hybride à entraîné la 
modification du schéma de calcul de la fonction d'Exner
$\pk=p^\kappa$.\footnote{$\kappa=R/C_p$ où $R$ et $C_p$ sont les constantes
thermodynamiques de l'air sec.}
L'équilibre hydrostatique est intégré verticalement suivant le schéma~:
\begin{equation}
\label{eq:hydrostat}
\dz \Phi=-  \av{\h}{z} \dz \pk
\end{equation}
où $\Phi$ est le géopotentiel,
avec $\av{\h}{z}_1=\h_1$ et ${\dz \pk}_1=\ps^\kappa-\pk_1$.
Au lieu de calculer $\pk_k$ au
milieu de la couche $k$ à partir d'une moyenne arithmétique des pressions
$p_k$ et $p_{k+1}$ (ce qui est entre autres coûteux numériquement), on 
utilise une relation numérique différente entre  $p$ et $\pk$~:
\begin{equation}\label{eq:exner}
\frac{\aire}{g} \av{p\dz \pk}{z}=-\kappa \pk \m
\end{equation}
Cette relation garantit un équivalent numérique de la relation
\begin{equation}\label{eq:ener1}
\int_0^\infty \Phi \rho dz =\int_0^\infty RT \rho dz
\end{equation}
entre énergie potentielle et énergie interne, à savoir
\footnote{
A partir de l'équation hydrostatique \ref{eq:hydrostat}, on peut
écrire
\begin{equation}
\Phi_l=\Phi_s-\sum_{k=1}^l \depb{\av{\h}{z} \dz \pk}_k
\end{equation}
avec la convention (utilisée dans le modèle) 
$\depb{\av{\h}{z}}_1=\h_1$.
L'énergie potentielle de gravité totale de la colonne s'écrit
donc
\begin{eqnarray}
\sum_{l=1}^N\Phi_l\m -\Phi_s
&=&-\sum_{l=1}^N\m_l \sum_{1\le k \le l} \depb{\av{\h}{z} \dz \pk}_k \\
&=&-\sum_{k=1}^N\sum_{k\le l \le N} \depb{\av{\h}{z} \dz \pk}_k m_l \\
&=&-\sum_{k=1}^N\depb{\av{\h}{z}\dz \pk}_k \times \frac{\aire}{g} p_k\\
&=&-\frac{\aire}{g}\sum_{k=1}^N\h_k\depb{\av{p \dz \pk}{z}}_k
\end{eqnarray}
La dernière transformation, qui consiste à faire glisser la moyenne d'une
demi maille sur la verticale, est valable si on choisit pour les
conditions aux limites inférieure et supérieure~:
\begin{equation}
\depb{\av{p\dz \pk}{z}}_1=\frac{p_2 \dep{\dz \pk}_2}{2}+\ps \dep{\dz \pk}_1
\end{equation}
et
\begin{equation}
\depb{\av{p\dz \pk}{z}}_N=\frac{p_{N-1} \dep{\dz \pk}_{N-1}}{2}
\end{equation}
On voit alors que l'équation~\ref{eq:ener1} peut être satisfaite
simplement si on choisit les niveaux $\pk$ suivant la relation
\ref{eq:exner}, cqfd.
}
\begin{equation} \label{eq:ener1}
\sum_{l=1}^N\Phi_l\m_l -\Phi_s = \sum_{l=1}^N RT_l\m_l =
\sum_{l=1}^N \kappa \h_l \pk_l \m_l
\end{equation}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Les paramétrisations physiques du modèle de climat terrestre\label{sec:phylmd}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


La première utilisation de ce nouveau code hydrodynamique est martienne.
Le modèle du LMD est le premier à simuler un cycle saisonnier entier sur la planète rouge \cite[]{Hour:93}.
Pendant ce temps le modèle climatique terrestre du LMD poursuit une vie trépidante. La version 4ter, qui donnait des résultats excellents sur la Mousson Indienne (présentez une simulation du modèle à un chercheur authentique du LMD, il commencera par regarder la ``précip'' de juillet sur l'Inde), est progressivement remplacée par les versions 5 et 6. Ces versions sont notamment couplées à SECHIBA pour la végétation et à ORCALIM pour l'océan.

C'est lors d'un séminaire interne du LMD, tenu en la royale abbaye de Fontevraud en 1990, qu'est prise collectivement la décision de faire de ce nouveau code dynamique l'ossature du futur modèle de climat du LMD. Ce choix est motivé par la volonté de disposer d'un outil souple et modulaire, permettant facilement l'échange et le test de modules et de procédures. Cette décision s'inscrit dans l'idée de modélisation communautaire alors dans l'air du temps.

Il faut donc adjoindre au code dynamique un jeu de paramétrisations physiques. Dans un premier temps, la décision est prise de partir de la physique du Centre Européen pour les Prévisions Météorologiques à Moyen Terme. Ce modèle contient en effet le code radiatif Fouquart/Morcrette, développé à l'origine pour le modèle du LMD, et les schémas de nuages d'Hervé LeTreut.
Cette décision ne sera finalement pas tenue. Si Laurent Li s'est bien battu pour effectuer les portages nécessaires, le groupe décidera quelques années plus tard de privilégier la continuité avec les versions précédentes.\footnote{Un lourd exercice de ``convergence'' LMDZ/LMD-5ter sera mené à bien en 1996. Michèle Forichon y laissera pas mal d'énergie. Plus généralement, les acteurs de l'époque, qui tentent de mener à bien les décisions de Fontevraud, s'useront entre décisions, attaques et contre-décisions.}
C'est la combinaison de ce nouveau code dynamique et de la réécriture
des paramétrisations physiques qui donnera naissance au modèle LMDZ d'aujourd'hui.

On décrit ci-dessous dans ses grandes lignes la version du modèle
LMDZ3.3 qui a été utilisée pour les applications "traceurs"
décrites dans ce document ainsi que les évolutions apportées
au court de la mise au point du modèle couplé IPSLCM4 et qui
définissent la version LMDZ4 \cite[]{Hour:06LMDZ}.

L'effet de la turbulence de petite échelle dans la couche limite
est pris en compte au moyen d'une fermeture en super-viscosité avec
un flux turbulent proportionnel au gradient vertical de la quantité
transportée. La viscosité turbulente dépend du cisaillement
vertical du vent et du nombre de Richardson, selon les formules
présentées par \cite{Lava:81}.
Un "contre-gradient" est introduit sur la température potentielle pour
permettre un transport de chaleur "en remontant le gradient" dans les
couches limites marginalement stables.
On revient très largement sur ces aspects relatifs à la couche limite
dans le \ch{therm}.

Le code radiatif est celui du modèle du Centre Européen
pour les Prévisions Météorologiques à Moyen Terme \cite[]{Morc:91}.
La partie concernant
l'absorption et la diffusion du rayonnement solaire est une version raffinée
du modèle de \cite{Fouq:80}. La partie infrarouge thermique a été 
développée par \cite{Morc:86}. La condensation est paramétrée de façon
différente pour la convection profonde et pour les autres nuages.
Pour la convection nuageuse, on utilise
la paramétrisation développée par \cite{Tied:89}.
Cette paramétrisation est basée sur une représentation 
en flux de masse d'une colonne convective nuageuse idéalisée. 
La colonne atmosphérique
est divisée en trois sous-colonnes~: une ascendance concentrée,
une subsidence rapide associée à la pluie et enfin l'environnement, généralement
en subsidence lente.
Pour la partie non convective des nuages, on diagnostique une fraction
nuageuse et un contenu en eau des nuages en se donnant a priori une
fonction de distribution sous-maille de l'eau \cite[]{LeTr:91}.
Les nuages ainsi diagnostiqués sont utilisés à la fois pour le calcul 
radiatif et pour calculer un taux de chauffage et un taux de précipitation.
La pluie calculée ainsi peut se réévaporer dans les couches inférieures 
du modèle.
L'eau vapeur et l'eau condensée sont ensuite advectées indépendamment
dans la partie dynamique\footnote{Cette séparation entre eau vapeur et liquide
est en partie fictive car, juste après l'advection, en entrée
des paramétrisations physiques, l'eau liquide est
réévaporée pour travailler en eau totale. La séparation entre eau vapeur
et eau liquide pour l'advection est donc uniquement numérique.}.
Le transport de grande échelle de la vapeur d'eau et de l'eau liquide,
traité à l'origine avec un schéma amont très diffusif sur l'horizontale
et un schéma centré non positif sur la verticale, est maintenant calculé
avec les schémas en volumes finis décrits un peu plus loin.

Dans les expériences présentées dans ce document, la température de surface
de la mer est imposée. La conduction thermique dans les sols
continentaux est quant à elle calculée au moyen d'un schéma multi-couches
d'un sol homogène
développé à l'origine pour la version martienne du modèle
\cite[]{Hour:93}.
Si on note $\xi$ la conduction thermique du sol et $C$ la capacité thermique
volumique, on peut montrer facilement que le flux conductif à la surface
peut s'écrire sous la forme
$F_c=-I\partial T/\partial z'$
et la conduction dans le sol
 $\partial T/\partial t=\partial^2 T/\partial {z'}^2$.
$I=\sqrt{\xi C}$ est appelé inertie thermique et $z'$ est une pseudo profondeur
définie par $z'=z\sqrt{C/\xi}$.

Suivant \cite{Lava:81}, l'évaporation à la surface est calculée comme
$E=1,35 \beta \rho {\ustar}^2 [q_s(T_s)-q]$
où $\ustar$ est la tension  de vent en surface,
$q$ est l'humidité spécifique de la première couche du modèle atmosphérique
et $q_s(T_s)$ est l'humidité spécifique à saturation pour la température de
surface $T_s$.
Le coefficient d'aridité $\beta$ peut être soit imposé soit dépendre directement
d'un contenu en eau du sol, calculé au fil du temps en fonction du bilan
entre précipitation et évaporation,
$d q_s/dt = P-E$.
Le modèle utilisé par défaut, dit modèle du saut d'eau (bucket en anglais),
fait l'hypothèse que $q_s$ peut
croître jusqu'à une valeur de 150~mm d'eau. Entre 75 et 150 mm, le coefficient
d'aridité $\beta$ est égal à 1 et l'évaporation est égale à l'évaporation
potentielle. Ce coefficient varie linéairement entre 0 et 1 pour des contenus
en eau du sol variant de 0 à 75 mm.


La partie physique du modèle est en fait en évolution constante.
La description qui en est faite ci-dessus correspond à une vue instantanée
de la version qui a servi à développer les aspects liés au transport
des traceurs. Sauf mention explicite, c'est cette version de LMDZ3
qui est utilisée
dans les simulations présentées dans ce document.

\figurecaption{
\centerline{\includegraphics[width=10.7cm]{\local/FIGURES/precip1g.eps}}
}
{Climatologie des précipitations en janvier
dans deux versions du modèle LMDZ.
En haut : la version LMDZ3 avec le schéma de convection de Tiedtke
et le modèle de seau d'eau en surface. Cette version a été beaucoup utilisée
pour le développement des aspects chimie et traceurs.
Au milieu : le modèle LMDZ4 avec
le schéma de convection d'Emanuel, des nuages couplés à la convection et
le schéma de surface SECHIBA du modèle ORCHIDEE. C'est la version utilisée pour le modèle
couplé de l'IPSL.
En bas : Observations (Climatologie de Xie-Arkin).
Les simulations sont effectuées en prescrivant les températures de surface
de l'océan et les cartes montrent des moyennes sur 5 années consécutives.}
{fg:precip1}
{1.}{0.73}{0.26}{-15.5cm}

\figurecaption{
\centerline{\includegraphics[width=10.7cm]{\local/FIGURES/precip7g.eps}}
}
{Climatologie des précipitations en juillet
dans deux versions du modèle LMDZ.
En haut : la version LMDZ3 avec le schéma de convection de Tiedtke
et le modèle de seau d'eau en surface. Cette version a été beaucoup utilisée
pour le développement des aspects chimie et traceurs.
Au milieu : le modèle LMDZ4 avec
le schéma de convection d'Emanuel, des nuages couplés à la convection et
le schéma de surface SECHIBA du modèle ORCHIDEE. C'est la version utilisée pour le modèle
couplé de l'IPSL.
En bas : Observations (Climatologie de Xie-Arkin).
Les simulations sont effectuées en prescrivant les températures de surface
de l'océan et les cartes montrent des moyennes sur 5 années consécutives.}
{fg:precip7}
{1.}{0.73}{0.26}{-15.5cm}

Depuis, notamment en vue du développement du modèle couplé IPSLCM4, la physique
du modèle a évolué de façon significative.
Les améliorations principales concernent~:
\begin{enumerate}
\item le remplacement du schéma de convection nuageuse de Tiedtke  par
le schéma d'\cite{Eman:93}. Comme celui de Tiedtke, le schéma d'Emanuel
est basé sur une séparation de la colonne atmosphérique en trois
compartiments.
Les différences entre les deux schémas concernent principalement la
fermeture (paramétrisation
donnant par exemple le flux de masse à la base de l'ascendance adiabatique),
l'échange d'air entre l'ascendance et l'environnement (décrit
de façon plus sophistiquée dans le schéma d'Emanuel) ainsi que
les descentes précipitantes \cite[beaucoup plus actives dans le schéma
d'Emanuel, cf. ][]{Guic:04,Hour:06LMDZ}.
Un travail important est réalisé au LMD sur cette paramétrisation
d'Emanuel \cite[]{Gran:04}.
\item le schéma statistique de nuages qui a été couplé à la paramétrisation de
la convection suivant l'approche de \cite{Bony:01}.
Un effort important a été consacré au réglage de la composante
nuageuse.
\item le modèle ORCHIDEE qui est maintenant utilisé pour représenter la
thermodynamique des surfaces continentales \cite[]{Rosnay:2002,Krinner:gbc2005}.
Ce modèle distingue pour l'eau du sol un réservoir de surface et un réservoir
profond. Les paramètres thermodynamiques de la surface sont prescrits en
fonction de 10  types de végétation. Le routage de l'eau par les grands
bassins fluviaux, depuis les zones de ruissellement jusqu'aux océans,
est représenté. Différents degrés de complexité peuvent ensuite être activés,
afin d'introduire
par exemple une représentation explicite de la croissance des
feuilles voir, à terme, un calcul évolutif de la végétation.
\item la prise en compte de 4 sous-surfaces avec des mailles fractionnaires
pour les terres, les glaciers, les océans  et la banquise. Ces 4 sous-surfaces
sont en fait couplées à 4 sous-colonnes d'atmosphère sur lesquelles
on effectue indépendamment 4 calculs de diffusion turbulente dans la couche
limite.
\end{enumerate}
Le modèle ainsi modifié
-- version LMDZ4 -- couplé au modèle global d'océan et de banquise
ORCALIM, est utilisé notamment pour effectuer des prévisions du changement climatique.

On ne s'appesantit pas dans ce document sur la climatologie du modèle
d'atmosphère. On montre cependant une comparaison de la version précédente
LMDZ3 du modèle et du nouveau modèle LMDZ4 en ce qui concerne la précipitation
(\fig{precip1} et \fig{precip7}).
Les différences principales entre les deux versions proviennent en fait
du changement de schéma de convection. L'ancien modèle avait tendance
à exagérer l'intensité des précipitations tropicales dans les zones de
convergences, notamment à l'Est de Madagascar ainsi que sur l'océan Pacifique.
Ce biais est nettement diminué dans le nouveau modèle.
La diminution (amélioration) des précipitations sur les continents nord
en été (Sibérie, Canada), est elle davantage due au schéma de surface.
La détérioration la plus notable avec le nouveau modèle est sans doute
la distribution des précipitations de mousson autour de l'Inde en 
juillet.
 Lorsqu'on interprète des résultats relatifs à la simulation
des concentrations atmosphériques d'aérosols ou d'espèces chimiques,
il faut bien sûr garder à l'esprit
ces différences quant à la capacité du modèle à simuler correctement
le climat.





\subsection{Organisation informatique - version unidimensionnelle - mode guidé}

\begin{figure}
\centerline{\includegraphics[width=10cm]{\lmdztdir/img283.epsi}}
\caption{Séparation entre les parties dynamique 3D et physique 1D du code
LMDZ.\label{fg:physdyn}}
\end{figure}

\subsubsection*{Modularité}

L'organisation du modèle LMDZ est très étroitement calquée sur cette 
distinction entre partie "dynamique" -- la seule partie où soient pris en
compte des échanges horizontaux entre des mailles du modèle -- et la
partie ``physique" qui peut être vue comme une juxtaposition de ``colonnes"
d'atmosphère.

Cette spécificité de la partie physique est exploitée en ce sens que
le codage de toutes les paramétrisations est fait avec un indice
interne muet qui représente la grille horizontale.
Cette organisation est illustrée sur la \fig{physdyn}.
Chacune des parties, physique et dynamique, a ses propres fichiers
de conditions initiales et de sorties.

Cette écriture se prête à la fois à la vectorisation et
à la parallélisation des codes.
Cette approche permet également de disposer de façon transparente
d'une version unidimensionnelle du modèle de circulation.
Pour disposer d'un modèle unidimensionnel, il suffit d'écrire un programme
dans lequel on initialise des profils météorologiques sur
un point particulier du globe et dans lequel on appelle ensuite
en boucle le moniteur physique.

Autre atout important, cette conception modulaire permet de gérer
en parallèle des ``physiques" différentes interfacées avec le même code
dynamique. Ce point est essentiel pour les études menées au LMD sur Mars et
sur Titan.
Enfin, on peut noter qu'il existe aussi une version bidimensionnelle,
latitude-altitude, du noyau dynamique
du modèle qui a été abondamment utilisée sur Titan comme on le voit
à la fin de ce document.

\subsubsection*{Mode climatique et mode guidé}

Le modèle LMDZ est avant tout développé pour des études climatiques.
On effectue alors des intégrations longues, dans lesquelles l'état
initial est vite ``oublié" et les résultats ne s'interprètent qu'en termes 
statistiques.
Cependant, notamment pour des aspects de validation, il peut s'avérer
utile de contraindre le modèle à suivre une situation météorologique
observée. 
Dans ce mode ``guidé" (on parle souvent de ``nudging" en anglais), les
champs du modèle sont rappelés avec un terme linéaire vers les champs
des analyses ou des réanalyses
\begin{equation}
\dt{X}=M(X) + \frac{X_a-X}{\tau}
\end{equation}
où $X$ est une variable météorologique ($u$, $v$, $\theta$, $p_s$) et $X_a$
est la variable correspondante issue des analyses et interpolée sur la grille
spatio-temporelle du modèle. $M$ représente le calcul des
dérivées temporelles par le modèle.
Les temps de relaxation $\tau$ peuvent être différents pour les différentes
variables ou dépendre de la région considérée.
Cette approche peut être vue comme une alternative légère à
l'assimilation directe de données météorologiques dans le modèle
telle qu'elle est mise en {\oe}uvre dans les centres de prévision
opérationnelle \cite[]{Jeuk:96}.

Une approche intermédiaire a été mise en {\oe}uvre dans LMDZ
\cite[]{Bona:these}.
Elle consiste à minimiser
une fonction mesurant la distance entre l'état du modèle et les analyses.
Cette approche, qui s'apparente encore davantage à l'assimilation
météorologique opérationnelle,
permet d'introduire, en plus des analyses, des observations
supplémentaires ou des contraintes comme des pénalités sur les modes
de gravité excités par la procédure d'assimilation.
Cette dernière méthode n'est pas utilisée ici.

