\documentclass[11pt]{article}


\usepackage{natbib}
\bibliographystyle{jas}
%\renewcommand{\floatpagefraction}{0.7}

\begin{document}
\cite{Yama:83}. Cette version a été utilisée dans une étude d'intercomparaison
coucou
\end{document}






\def\q{c}
\section{Les fermetures basées sur une équation pronostique de l'énergie cinétique turbulente}


La fermetures turbulente utilisée dans LMDZ repose sur
une version simplifiée des
fermetures d'ordre 2 où on se focalise sur l'équation d'évolution
de l'énergie cinétique turbulente
\begin{equation}
e=\frac{1}{2}\depb{\overline{u'^2}+\overline{v'^2}+\overline{w'^2}}
\end{equation}
cette énergie étant ensuite utilisée pour le calcul du coefficient de 
mélange $K_z$.

En effet,
si on revient à la présentation de la longueur de mélange faite précédemment,
il est naturel de prendre $\sqrt{e}$ comme amplitude des
mouvements turbulents verticaux $\ave{|w'|}$
dans le cas d'une turbulence isotrope en atmosphère neutre.
Dans le cas d'une atmosphère stable, où l'on s'attend à des mouvements
anisotropes plutôt horizontaux, $\sqrt{e}$ fournira plutôt
une surestimation de $\ave{|w'|}$ (respectivement une sous-estimation
pour une atmosphère instable).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\der#1#2{\frac{\partial #1}{\partial #2}}

Sous les hypothèses de couche limite mentionnées plus haut,
on peut montrer
%\cite[voir eg][]{Stul:88}
que l'évolution de
l'énergie cinétique turbulente s'écrit sous la forme
\begin{equation}
\der{e}{t}=P-D
-\frac{1}{\rho}\der{\ave{w'p'}}{z}
-\der{\ave{w'e}}{z}-\epsilon
\end{equation}
Le terme
\begin{equation}
P=-\overline{w'u'}\der{u}{z}-\overline{w'v'}\der{v}{z}
\end{equation}
est généralement positif
(en particulier si $\overline{w'u'}=-K_z\partial u/\partial z$)
et correspond à la production mécanique de turbulence.


Le terme 
\begin{equation}
D=-\frac{g}{\theta}\overline{w'\theta'}
\end{equation}
est en général positif
dans une atmosphère stable. Il correspond alors
à l'inhibition (ou destruction)
de la turbulence par les effets de stratification.
Il peut devenir un terme de production dans les atmosphères instables.
Les termes suivants sont un terme de pression, le transport vertical
turbulent d'énergie cinétique turbulente (le transport par la grande échelle
est négligé car $w$ est supposé petit), et $\epsilon$ est le terme de
dissipation mécanique de la turbulence.
Physiquement, l'énergie contenue principalement dans les plus grosses
structures turbulentes ``cascade" vers les échelles moléculaires où elle
est dissipée.

\subsection{TKE et conservation de l'énergie}

On va partir d'une forme simplifiée de l'équation de la TKE,
\begin{equation}
\der{e}{t}=P-D
-\frac{1}{\rho}\der{\ave{w'p'}}{z}
-\der{\ave{w'e}}{z}-\epsilon
\end{equation}
avec 
\begin{equation}
D=-\frac{g}{\theta_0}\overline{w'\theta'}
\end{equation}

\begin{equation}
\der{}{t} \int_0^\infty \rho e dz =
 \int_0^\infty \rho (P-D-\epsilon) dz
\end{equation}

\begin{eqnarray}
\der{}{t} \int_0^\infty \rho e dz &=&
\depb{ -\rho\overline{w'u'}u -\rho\overline{w'v'}v}_0^\infty + 
\int_0^\infty \dep{\der{\rho\overline{w'u'}}{z}u +\der{\rho\overline{w'v'}}{z}v} dz\\
&+&\depb{\frac{g\rho}{\theta_0}z\overline{w'\theta'}}_0^\infty
-\int_0^\infty \frac{g}{\theta_0} z
\der{\rho\overline{w'\theta'}}{z} dz
-\int_0^\infty \rho \epsilon dz
\end{eqnarray}


Si on introduit la tension de vent en surface $(\overline{w'u'},\overline{w'v'})_s=-C_D (u,v)_s=-C_D {\bf v}$, et si on utilise 
\begin{equation}
\rho{\der{X}{t}}_{|\mbox{diff.}}=-\der{\overline{\rho w'X'}}{z}
\end{equation}

\begin{eqnarray}
\der{}{t} \int_0^\infty \rho e dz &=&
-C_D {\bf v}_s^2 -
\int_0^\infty \rho\der{\bf v}{t}_{|\mbox{diff.}}{\bf . v} dz
+\int_0^\infty \frac{gz}{\theta_0}\rho\der{\theta}{t}_{|\mbox{diff.}} dz
-\int_0^\infty \rho \epsilon dz
\end{eqnarray}

\def\ps{{p_s}}
\begin{eqnarray}
\der{}{t} \int_0^\ps  e dp &=&
-g C_D {\bf v}_s^2 
-
\der{}{t}_{|\mbox{diff.}} \int_0^\ps  \frac{1}{2}{\bf v}^2 dp
-\int_0^\ps \frac{gz}{\theta_0} \der{\theta}{t}_{|\mbox{diff.}}   dp
-\int_0^\ps  \epsilon dp
\end{eqnarray}

Pour le vent, je trouve ce que je cherchais :
l'évolution de l'intégrale verticale de la TKE est opposée à la destruction de l'énergie cinétique grande échelle sous l'effet de la diffusion verticale.

Pour la partie thermique, que peut-on en dire ?

Par rapport aux corrections mises en place par Jean-Louis pour la conservation de l'énergie, si j'ai bien compris, ca revient à diagnostiquer le terme
\begin{equation}
C_D {\bf v}_s^2+ \der{}{t} \int_0^\ps  \frac{1}{2}{\bf v}^2 \frac{dp}{g}
\end{equation}
entre le début et la fin de la physique, et de rajouter ça comme chauffage localement.
Faut-il bien prendre en compte le terme de frottement et l'est il ? ...
Si ce terme est négatif, c'est à dire si on perd plus d'énergie cinétique dans l'atmosphère qu'on en gagne par frotement en surface, alors on réchauffe localement l'atmosphère comme on peut (en ajoutant un terme de chauffage là où l'énergie a beaucoup diminué ?).

Avec l'équation de la TKE, on peut faire beaucoup mieux.
Le chauffage à ajouter est en fait le terme $\epsilon$ de la dissipation.

Avec le modèle de Yamada utilisé dans LMDZ (cf. description détaillée ci-dessous), 
\begin{equation}
\epsilon=\frac{q^3}{l B_1}
\end{equation}
où $l$ est la longueur de mélange ("long" dans le codes yamada4.F90, de l'ordre de 10 à 300 m a priori) et $B_1=16,6$.
La TKE est reliée à $q$ par $e=q^2/2$. Attention, actuellement, c'est $q^2$ qu'on nomme pbltke dans physiq.F et dans les fichiers de sorties ...

Je laisse pour mémoire ci-dessous la description plus précise du modèle de Yamada utilisé dans LMDZ.

\subsection{Le modele de Mellor et Yamada}

Mellor et Yamada, à partir d'une paramétrisation des termes du troisième
ordre (à la base il s'agit donc d'une fermeture à l'ordre 2) proposent
une série de simplifications, avec une hiérarchie de schémas.
Nous présentons ici le schéma de niveau 2.5 (qui malgré cette nomenclature
originale n'est qu'une approximation relativement grossière d'un schéma
d'ordre 2), le plus largement testé dans des
applications météorologiques, et dont certains résultats sont présentés
plus loin.
Dans le modèle 2.5, les termes croisés s'écrivent formellement, comme pour la 
diffusion turbulente,
\begin{equation}
\overline{w'\phi'}=-K_\phi \der{\phi}{z}
\end{equation}
On écrit la diffusivité turbulente $K_\phi=lqS_\phi$ en fonction d'une
longueur de mélange, $l$, d'une vitesse typique $q=\sqrt{2e}$ et d'une
fonction de stabilité, qui peut prendre des valeurs différentes
pour les traceurs ($S_\phi$), la température ($S_h$) ou la quantité de
mouvement ($S_m$).




\def\ca{1,96}
\def\cb{0,1912}
\def\cc{0,2341}
\def\cd{1}
\def\ce{0,2231}
\def\cf{1,318}

\def\ba{0,6588}
\def\bb{0,1776}
\def\bc{0,3221}
\def\be{0,03156}
\def\ce{0,2231}
\def\cf{1,1318}

Un des tours de force de ce travail est de faire sortir directement
de simplifications successives d'une fermeture du second ordre
une formulation analytique pour les fonctions $S_m$ et $S_h$.
Ces grandeurs sont fonctions du seul {\bf nombre de  Richardson de flux},
$\rif=D/P$, mesure de l'importance relative du forçage
mécanique de la turbulence et de l'inhibition par stratification.
On retient ici
les valeurs numériques données par
\cite{Yama:83}. Cette version a été utilisée dans une étude d'intercomparaison
de fermetures turbulentes par \cite{Ayot:96}. Elle est également utilisée
plus loin dans des tests numériques. Pour la quantité de 
mouvement, on obtient~:
\begin{equation}
S_m=\ca\frac{\dep{\cb-\rif}\dep{\cc-\rif}}{\dep{\cd-\rif}\dep{\ce-\rif}}, \rif< 0,16
\end{equation}
et
\begin{equation}
S_m=0,085, \rif\ge 0,16
\end{equation}
Pour la température, on obtient pour l'inverse du nombre de Prandtl turbulent
$\omega=K_h/K_m=S_h/S_m$, l'expression
\begin{equation}
\omega=\cf\frac{\ce-\rif}{\cc-\rif}, \rif< 0,16
\end{equation}
\begin{equation}
\omega=1,12, \rif\ge 0,16
\end{equation}
En remarquant que $\rif=\omega \ri$, où
\def\shear{\left|\left|\der{\V}{z}\right|\right|}
\begin{equation}
\ri=\frac{g}{\theta}\frac{\der{\theta}{z}}{{\shear}^2}
\end{equation}
{\bf est le nombre de Richardson gradient},
on obtient finalement
\begin{equation}
\rif= \ba [\ \ \ri+\bb - \sqrt{\ri^2-\bc\ri+\be}\ \ ]\mbox{,   } \ri<\ric
\end{equation}
\begin{equation}
 \rif=\rifc\mbox{,   } \ri\ge\ric
\end{equation}
où $\rifc=0,191$ et $\ric=0,195$ sont des nombres de Richardson critiques.
On obtient ainsi des expressions explicites de $S_m$ et $\omega$ en fonction
des seules variables de grande échelle.


Dans le modèle 2.5, le terme de pression est négligé et le transport
vertical est traité comme une diffusion turbulente
de sorte que l'évolution de l'énergie
cinétique turbulente ou de $q$ s'écrit finalement
\def\v{{\bf v}}
\begin{equation}
\label{eq:dqorig}
\frac{1}{2}\frac{\partial q^2}{\partial t}=
q l S_m{\shear}^2 \depb{1 - \omega \ri} - \frac{q^3}{lB_1} + \frac{\partial}{\partial z}
\depb{l q S_q \frac{\partial q^2 / 2}{\partial z}}
\end{equation}
avec $B_1=16,6$.
\cite{Yama:83} utilise $S_q=0,2$ à la fois pour l'eau et pour les traceurs
passifs.

La paramétrisation de Mellor et Yamada, même dans ses versions plus
sophistiquées, laisse en fait une grosse zone d'ombre sur la spécification de la
longueur de mélange $l$. Dans certains articles, les auteurs ont
suggéré une équation pour $q^2l$ analogue à l'équation pour $q^2$ mais
cette équation est beaucoup moins fondée que sa grande s{\oe}ur (selon
les auteurs eux-mêmes).
Le plus souvent, la longueur de mélange est spécifiée, par exemple
en utilisant la formule de \cite{Blac:62}
\begin{equation}
l=l_0\frac{\vk z}{\vk z+l_0}
\end{equation}
asymptotique à $\vk z$ près de la surface et à une longueur $l_0$ au-dessus
de la couche de surface.
La longueur $l_0$ peut elle-même être fixée à une 
constante (de l'ordre de 100-200~m la plupart du temps) ou être calculée
en fonction d'autres quantités.
\cite{Yama:83} utilise par exemple pour $l_0$ une altitude moyenne pondérée par
l'intensitée de  la turbulence
\begin{equation}
l_0=0,2\ \frac{\int_0^\infty z q dz}{\int_0^\infty q dz}
\end{equation}
C'est cette formulation qui est retenue pour la version introduite
dans le modèle LMDZ et dont on montre certains résultats plus loin.
\footnote{A noter que l'intégration numérique du modèle de Mellor
et Yamada s'avère souvent délicate.
Une intégration naïve de l'équation d'évolution de l'énergie cinétique
turbulente avec un schéma temporel explicite (on calcule les termes
sources et puits du membre de droite au temps $t$ qu'on ajoute 
à l'énergie cinétique au temps $t$ pour obtenir la nouvelle
valeur à $t+\delta t$) contraint à prendre des pas de temps de quelques
secondes, même avec les discrétisations grossières utilisées dans
le modèle de circulation générale.
Dans la version développée pour LMDZ, on contourne en partie cette
difficulté en récrivant formellement l'équation d'évolution
de l'énergie cinétique turbulente (sans diffusion) sous la forme
\begin{equation}
\frac{1}{2}\frac{\partial q^2}{\partial t}=q^3 \chi
\end{equation}
ou encore
\begin{equation}
\label{eq:dq}
\frac{\partial}{\partial t}\dep{\frac{1}{q}}=- \chi
\end{equation}
avec
\begin{equation}
\chi = \frac{l S_m}{q^2} M^2\dep{1-\rif}-\frac{1}{l B_1}
\end{equation}
Si on suppose que $\chi$ ne varie pas au court d'un pas de temps,
la solution de l'\eq{dq} est
\begin{equation}
q^{(t+\delta t)} =  \frac{q^{(t)}}{1-\chi^{(t)} q^{(t)} \delta t}
\end{equation}
On retient directement cette solution quand $\chi\le 0$.
En revanche, quand $\chi>0$, on utilise une forme approchée
\begin{equation}
q^{(t+\delta t)} =  q^{(t)} \dep{1+\chi^{(t)}q^{(t)}\delta t}
\end{equation}

Cette formulation numérique produit des résultats numériques
presque indiscernables de l'intégration temporelle explicite de
l'équation d'origine mais avec des pas de temps de typiquement
quelques minutes à dizaines de minutes pour les configurations
classiques du modèle de circulation.

La diffusion verticale de l'énergie cinétique turbulente est calculée
a posteriori.
}% Fin de la footnote


\bibliography{/home/hourdin/TEX/BIBLIO/fred,/home/hourdin/TEX/LMDZ/lmdz}

\end{document}

