\documentclass{article}

\usepackage{graphicx}

\begin{document}

\input{def0}
\input{def}


\begin{center}
{\large
Proposition d'une formulation en volumes finis des équations primitives
}
\end{center}


\begin{center}
{\large
Laboratoire
de Météorologie
Dynamique du CNRS \\
Ecole Normale Supérieure \\
24 rue Lhomond\\
75231 PARIS cedex 05\\
FRANCE\\
}
\end{center}

\def\localdir{.}

\def\h{\theta}

\subsection*{Modèle original}

La grille vectorielle, avec des points de vent zonal et méridien décalés
est commune aux deux versions et est définie comme suit.

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}.

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:
%
%
\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}
%
$\fext=2\Omega \sin{\phi} \cu \cv$\\
où $\Omega$ est la vitesse de rotation de la planète.
%
\begin{equation}
Z=\Z
\end{equation}
%
\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:

\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}
%
%
\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}
%
\begin{equation}
\dz \Phi=-\ps^\rcp  \av{\h}{z} \dz s
\end{equation}
%
%
\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.

On peut, à partir des équations sur $u$ et $v$, déduire une équation
sur la vorticité:
\begin{equation}
\end{equation}


\section{Nouvelle formulation}

La 

L'analyse précédente dégage deux erreurs numériques
indépendantes sur la conservation du moment cinétique,
l'une liée à la définition du facteur de Coriolis et l'autre
à la mauvaise représentation de l'advection méridienne
de $u$.
La première erreur est éliminée en utilisant pour $f$
sa forme rotationnelle (équation~\ref{eq:frot}).
La seconde (plus fondamentale)
peut être éliminée également
en tirant partie du fait que la conservation numérique
de l'enstrophie ne
concerne que les termes d'advection horizontale.
Il est en fait possible de transformer le terme d'advection
verticale afin d'éliminer exactement les erreurs sur l'advection
méridienne du moment cinétique présentées dans la section
précédente.
On peut pour ce faire utiliser les deux formulations suivantes:
\begin{equation} \label{eq:newu1}
\dt{\ucov} -
\av{Z}{Y} \av{V}{X,Y} 
+ \dx \left( \Phi + K \right)
+s \av{\h}{X} \dx \psk
- \frac{\av{\ucov}{Y,Y} \dz \av{W}{X} }
{\av{\pext}{X} \dsig }
+ \frac{\dz \left( \av{W}{X} \av{\ucov}{Z} \right) }
{\av{\pext}{X} \dsig}
=S_{\ucov}
\end{equation}
ou
\begin{equation} \label{eq:newu2}
\dt{\ucov} -
\av{Z}{Y} \av{V}{X,Y} 
+ \dx \left( \Phi + K \right)
+s \av{\h}{X} \dx \psk
- \frac{\av{\uabs}{Y,Y} \dz \av{W}{X} }
{\av{\pext}{X} \dsig }
+ \frac{\dz \left( \av{W}{X} \av{\uabs}{Z} \right) }
{\av{\pext}{X} \dsig}
=S_{\ucov}
\end{equation}
o\'u $\uabs$ est la composante zonale covariante
du vecteur vent absolu:
\begin{equation}
\uabs=\ucov+\cu a \Omega \cos{\phi} =\dep{\frac{d\lambda}{dX}}M
\end{equation}
La deuxième formulation est plus proche
d'une équation d'évolution
du moment cinétique.
On vérifie facilement que ces deux formulations conservent
le moment cinétique dans le cas simple étudié précédemment
d'un écoulement axi-symétrique avec $\pext=cte$;
cette conservation est obtenue grâce à l'introduction du terme
d'advection verticale de $\av{\uabs}{YY}$.
En contre-partie, les nouvelles équations ne conservent plus
l'énergie globale, mais celle-ci n'était conservée, dans
la formulation originale, que pour des écoulements irrotationnels.

On utilise pour la deuxième composante de l'équation du mouvement
une formulation équivalente:
\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{\vcov}{X,X} \dz \av{W}{Y} }
{\av{\pext}{Y} \dsig}
+ \frac{\dz \left( \av{W}{Y} \av{\vcov}{Z} \right) }
{\av{\pext}{X} \dsig}
=S_{\vcov}
\end{equation}

\subsection*{Formulation en volumes finis}

On se rend compte en fait que le schéma ci-dessus est écrit
de sorte qu'on a en fait un schéma en volumes fini sur la vorticité.

En effet, si on écr


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

\end{document}
