\section{Formulation originale et propriétés de conservation}

\subsection{Le schéma numérique}


Le modèle climatique du LMD est bâti, comme tous les
modèles de circulation générale atmosphérique,
sur la résolution numérique des {équations primitives
de la météorologie} décrites dans de nombreux
ouvrages~\cite{Holt:79}.
L'analyse présentée ici a été menée sur la nouvelle
version de la dynamique du LMD écrite par Phu Le Van~\cite{LeVa:89}
sur une formulation de Robert Sadourny.
Cette formulation diffère de l'ancienne essentiellement
par deux points:
dans la nouvelle formulation, la répartition des points en
longitude et en latitude peut être changée arbitrairement.
L'autre modification porte sur la répartition des points
aux pôles\footnote{Aux pôles sont calculés:
le vent méridien dans l'ancienne formulation et les variables
scalaires dans la nouvelle.}.

L'analyse présentée ici peut être transposée sans efforts à
l'ancienne formulation. En fait, l'ancienne configuration de la grille,
dans laquelle les mailles ont toutes la même surface\footnote{
Dans l'ancienne formulation
de la dynamique, les points étaient équidistants
en sinus de la latitude et régulièrement répartis en longitude.
Toutes les mailles avaient par conséquent la même aire.},
est même d'une certaine façon
optimale pour la conservation du moment
cinétique.\\

La coordonnée  verticale du modèle est la pression normalisée
par sa valeur à la surface: $\sigma=p/\ps$.
On utilise en fait $\sigma$ aux niveaux inter-couches
et $s=\sigma^\kappa$ au milieu des couches.
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
(température potentielle $\h = c_p T/\psk$, géopotentiel $\Phi$
et pression de surface $\ps$) sont évaluées aux points
correspondant à des couples de valeurs entières $(X,Y)=(i,j)$.
Les variables dynamiques sont décalées par rapport aux variables
scalaires en utilisant une grille $C$ dans la définition de
Arakawa~\cite{Arak:77}: 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
Figure~\ref{fg:grille}.

\begin{figure}
\centerline{\includegraphics[width=7cm,angle=-90]{post/nouveau.epsi}}
\end{figure}

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.
On introduit également:
%
\paragraph{la pression extensive:}
$\pext$ (pression au sol multipliée
par l'aire de la maille).
%
\paragraph{les trois composantes du flux de masse:}
\begin{equation}
U=\av{\pext}{X} \ucont ,\  V= \av{\pext}{Y} \vcont \  \mbox{et} \  
W= \pext \dot{\sigma}
\ \mbox{avec}\ \dot{\sigma}=\frac{d\sigma}{dt}
\end{equation}
%
\paragraph{le facteur de Coriolis multiplié par l'aire de la maille:}
$\fext=2\Omega \sin{\phi} \cu \cv$\\
où $\Omega$ est la vitesse de rotation de la planète.
%
\paragraph{la vorticité potentielle absolue:}
\begin{equation}
Z=\Z
\end{equation}
%
\paragraph{l'énergie cinétique}
\begin{equation}
K=\K
\end{equation}\\
%
La notation $\delta X$ signifie simplement qu'on
effectue la différence entre deux points consécutifs
suivant la direction $X$.
La notation $\av{a}{X}$ signifie qu'on prend la moyenne arithmétique
de la quantité $a$ suivant la direction $X$.
Les équations discrétisées sont écrites sous la forme
suivante:
\paragraph{équations du mouvement:}
\begin{equation} \label{eq:u1}
\dt{\ucov} -
\av{Z}{Y} \av{V}{X,Y} 
+ \dx \left( \Phi + K \right)
+s \av{\h}{X} \dx \psk
+\frac{\av{\av{W}{X} \dz \ucov}{Z}}{\av{\pext}{X}\dz \sigma}
=S_{\ucov}
\end{equation}
\begin{equation} \label{eq:v1}
\dt{\vcov} + \av{Z}{X} \av{U}{X,Y} + \dy \left( \Phi + K \right)
+s \av{\h}{Y} \dy \psk
+\frac{\av{\av{W}{Y} \dz \vcov}{Z}}{\av{\pext}{Y}\dz \sigma}
=S_{\vcov}
\end{equation}
%
\paragraph{équation thermodynamique:}
%
\begin{equation}
\label{eq:thermo}
\dt{\dep{\pext \h}}
+\dx \dep{\av{\h}{X}U} +\dy \dep{\av{\h}{Y}V}
+\frac{\dz \dep{\av{\h}{Z} W}}{\dz \sigma}=S_\h
\end{equation}
%
\paragraph{équation hydrostatique:}
\begin{equation}
\dz \Phi=-\ps^\rcp  \av{\h}{z} \dz s
\end{equation}
%
\paragraph{équations de continuité:}
%
\begin{equation}
\label{eq:cont1}
\dt{\ps}  = \sum_z{\dz \sigma \dep{\dx U+ \dy V}}
\end{equation}
\begin{equation}
\label{eq:cont2}
\dz W   = -\dz \sigma \dep{\dx U+ \dy V + \dt{\ps}}
\end{equation}
%
On a noté $S$ les termes sources dans les différentes équations.

\subsection{Les propriétés de conservation}

On dit qu'une quantité $q$ est conservée globalement par
les équations discrétisées si on peut trouver un équivalent
numérique $[q]$ de $q$ tel que:
\begin{equation}
\dt{\glob{\pext d\sigma [q]}}=0
\end{equation}
où la notation $\glob{}$ indique qu'on effectue une sommation sur
l'ensemble du domaine (on s'intéresse uniquement
à la façon de calculer les dérivées spatiales et pas
à l'intégration temporelle).
Les équations~\ref{eq:cont1} et~\ref{eq:cont2}, par exemple, sont
écrites de façon à ce que la masse totale de l'atmosphère
soit exactement conservée.
De même l'équation thermodynamique~\ref{eq:thermo}
permet de conserver numériquement $\theta$ pour les écoulements
adiabatiques. En fait, la forme particulière de cette équation
est choisie pour conserver à la fois $\theta$ et $\theta^2$
(ce qui se montre en combinant cette équation avec l'équation
de continuité). On montre, plus généralement,
que pour une quantité 
$q$ conservée par les équations physiques
et une fonction $f$ de $q$ donnée,
il existe une formulation numérique permettant
de conserver à la fois $q$ et $f\dep{q}$
(voir par exemple~\cite{Sado:80}).

Le problème de la conservation devient plus complexe lorsqu'on
s'attaque aux quantités faisant intervenir les équations dynamiques.
Il est par exemple impossible de conserver à la fois l'énergie
et l'enstrophie avec les formulations classiques utilisées dans les
GCM~\cite{Sado:75}. Les modèles ont généralement privilégié
la conservation de l'énergie. Robert Sadourny
a souligné l'importance
de la conservation de l'enstrophie et montré qu'elle
conduit à une plus grande stabilité numérique,
au moins dans un modèle des équations de couche
mince~\cite{Sado:75a}.
La conservation de l'enstrophie a donc été privilégiée
dans le modèle du LMD.
Le dilemme porte en fait sur la façon de discrétiser la
partie rotationnelle des termes advectifs dans les équations
horizontales du mouvement.

Dans la formulation originale, l'énergie totale n'est conservée
que pour les écoulements irrotationnels. Cette conservation
elle-même n'est pas immédiate et conduit à l'introduction
d'une contrainte sur la disposition des niveaux
verticaux~\cite{Sado:80}:
les niveaux intercouches (niveaux $\sigma$) et les milieux des
couches (niveaux $s$) sont répartis de façon à satisfaire
la relation suivante:
\begin{equation}
\av{\sigma \dz s}{Z}=\kappa s \dz \sigma
\end{equation}
Cette relation conduit à déplacer le centre des couches par rapport
à la moyenne arithmétique des niveaux inter-couches.
Le décalage peut devenir tel que le centre d'une couche
se trouve à l'extérieur de la couche elle-même.
En pratique, il faut choisir une répartition suffisamment lisse
des niveaux $\sigma$ pour que les niveaux $s$ soient proches du centre
des couches en terme de masse.

Le terme rotationnel est lui conçu pour conserver la vorticité
potentielle dans le cas d'un écoulement barotrope ou dans le
cas d'un modèle à coordonnée verticale
isentropique~\cite{Sado:75a}.

Dans toute cette partie, on a volontairement omis le filtrage
longitudinal introduit
dans les hautes latitudes du modèle pour éviter d'utiliser un pas
de temps trop petit, du fait de la diminution des distances
entre les points de grille près des pôles.
