\documentclass{article}
\usepackage{a4wide}
\usepackage{graphicx}
\usepackage{epsf}

    \usepackage[francais]{babel}
         \usepackage[latin1]{inputenc}




\begin{document}
\def\abs#1{\left| #1 \right|}
\def\dep#1{\left( #1 \right) }
\def\depb#1{\left[ #1 \right] }
\def\der#1#2{\frac{\partial #1}{\partial #2}}
\def\eq#1{Eq.~\ref{eq:#1}}
\def\localdir{.}


\title{Développement d'un modèle physico-numérique de la couche limite convective
\author{
Cours M2 de modélisation numérique, Fr\'ed\'eric Hourdin\\
hourdin@lmd.ens.fr
}\date{\today}}

\def\grad{\vec{\mbox{grad}}\ }
\def\afaire#1{{\bf A FAIRE : #1}}
\maketitle

\def\xmax{x_{\mbox{max}}}
\def\zmax{z_{\mbox{max}}}
\def\mod{\mbox{mod}}

Préambule :

Dans la première heure : réaliser les questions 1 et 2 et réfléchir au
choix du sujet (1 à 4) à réaliser par la suite. 
Pour cela, prenez le temps de lire l'ensemble du sujet.

Il faudra ensuite pour chacun coder des composantes complémentaires
(typiquement dans la fin de la séance d'aujourd'hui et lors de la prochaine) et
effectuer des tests indépendants de ces nouvelles composantes.

Dans une dernière phase (dernière séance), on  couplera tout ou partie
de ces développements pour effectuer des expériences numériques fascinantes.

Veiller à bien coder vos développement dans des sous-programmes indépendants.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le cas ARM cumulus de convection peu profonde continentale}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

On va s'intéresser ici au développement d'un modèle physico-numérique
de la couche limite convective, en s'appuyant sur un cas ARM (Atmospheric
Radiation Measurement, site d'observation en Oklaheoma)
de cycle diurne de couche limite convective avec cumulus l'après-midi.


Pour ce mini-projet, on va travailler avec la version uni-colonne (SCM pour Single Column Model en anglais) ou 1D du modèle LMDZ.

%Pour activer ce cas, il faut d'abord installer le modèle. Vous pouvez partir de la version de LMDZ installée à la séance précédente.
%Puis se déplacer sur le répertoire \texttt+LMDZ+ pour éventuellement mettre à jour le modèle grâce à la commande 'svn update'.
%Ensuite, installer le modèle 1D standard :
%\begin{verbatim}
%cd
%wget http://www.lmd.jussieu.fr/~lmdz/Distrib/install_lmdz.sh
%chmod +x install_lmdz.sh
%./install_lmdz.sh -name LMDZ -SCM
%\end{verbatim}
%(mais c'est djà fait normalement).
%Si le modèle  a déjà été installé, se rendre sur le répertoire LMDZ (ou un
%autre nom) sur lequel on trouve le répertoire :
%\begin{verbatim}
%cd ~/LMDZ/1D
%\end{verbatim}

\def\LMDZmain{\verb+LMDZmini+}
\def\LMDZname{\texttt{LMDZ20250209}}
\def\casscm{\texttt{arm\_cu}}
\def\pathcasscm{\LMDZname\texttt{/1D/EXEC/6AL79lmdworkdt1080/}\casscm}
\def\libf{\LMDZname\texttt{/modipsl/modeles/LMDZ/libf}}
\def\workdir{\LMDZname\texttt{/modipsl/modeles/LMDZ/libf/phylmdwork/work}}

\subsection{Installation du modèle}

Placez vous de préférence dans le répertoire dans le quel vous avez lancer l'installation du modèle lors de la première séance (mais ça marchera aussi dans un autre répertoire).
Récupérer un fichier tar contenant les script d'installation de LMDZ, de la version 1D, et des modification nécessaires pour le mini-projet.
Ces modifications consistent à remplacer le jeu de paramétrisations physiques de LMDZ par une K-diffusion (ou diffusion turbulente) avec un coefficient de mélange K constant.

\begin{verbatim}
wget http://web.lmd.jussieu.fr/~hourdin/COURS/M2/MINI/ATM/ParamDev/param_dev.tar
tar xvf param_dev.tar
./param_dev_install.sh
\end{verbatim}
Cette commande va successivement :
\begin{enumerate}
\item installer une nouvelle version de LMDZ, ici \LMDZname.
\item installer les cas 1D sur \LMDZname\texttt{/1D}
\item lancer la commande \texttt{run.sh} pour le cas \casscm. Les résultats se trouveront dans le répertoire \LMDZname\texttt{/1D/EXEC/6AL79/}\casscm.
\item dupliquer le répertoire contetenant les paramétrisations physiques du modèle, \texttt{phylmd}, se trouvant sur \libf\ en \texttt{phylmdwork} et effectuer sur ce dernier les modifications nécessaires pour le mini projet.
\item lancer à nouveau la commande \texttt{run.sh} avec des options différentes pour réaliser une première simulations avec la physique réduite à la K-diffusion dont les résultats se trouveront sur \pathcasscm.
\end{enumerate}

\subsection{Familiarisation avec le cas ARM cumulus}


Sur le répertoire \LMDZname\texttt{/1D/EXEC/6AL79/}\casscm,
on pourra regarder les résultats contenus dans le fichier \texttt{hourly.nc} et les comparer avec les résultats d'une simulations LES contenus dans \texttt{LES.nc}.
Le fichier \texttt{hourly.nc} contient les résultats d'une première simulation 1D effectuée avec la physique du modèle LMDZ.
On pourra regarder par exemple la nénbulosité (\texttt{rneb}), la température potentielle
(\texttt{theta}), l'humidité spécifique (\texttt{ovap}) ou relative (\texttt{rhum}).
\begin{verbatim}
ferret
use hourly.nc
reg/k=50:79
set v ul
fill rneb
set v ur
fill theta
set v ll
plot/l=1 theta
plot/l=5/o theta
plot/l=10/o theta
set v lr
plot/l=1/hlim=0:20 1000*ovap
plot/l=5/o 1000*ovap
plot/l=10/o 1000*ovap
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Le cadre de travail pour la physique simplifiée}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Pour le mini-projet, on partira d'une maquette qui tourne déjà.
%On installera ensuite la maquette du mini-projet. Pour ce faire,
%\begin{verbatim}
%cd ../../modipsl/modeles/LMDZ5/libf
%cp -r phylmd phylmdMON
%cd phylmdMON
%wget http://www.lmd.jussieu.fr/~hourdin/COURS/UTILS/MINIORIG.tar
%tar xvf MINIORIG.tar
%ln -sf MINI/* .
%\end{verbatim}

Le script \texttt{param\_dev\_install.sh} a donc installé d'une part une nouvelle version de la physique et fait tourné un premier cas 1D.

Vous pouvez regarder les résultats du cas 1D sur \pathcasscm,  en lançant sous \texttt{ferret} les commandes
\begin{verbatim}
use phys
fill/k=59:79 temp
\end{verbatim}

Le travail de codage s'effectuera dans le répertoire\\
\workdir.
Le programme central est \texttt{phys\_param\_dev.F90}. Commencez à regarder son contenu.
Ce programme contient pour le moment essentiellement le calcul d'une diffusion verticale (contenu dans le sous programme \texttt{diffkz}). 
Vous vous trouvez (virtuellement) dans le répertoire où sont placés les fichiers fortran de la maquette du mini-projet.
Quand vous voudrez modifier développer ou modifier des codes, vous le ferez dans ce répertoire.
Si vous créez un nouveau fichier contenant un sous-programme, par exemple avec l'éditeur gedit (et c'est fortement conseillé de créer des sous-profgrammes dans de fichier différents), vous devez faire une opération supplémentaire~:
\begin{verbatim}
gedit new_file.f90
cd ..
ln -s work/new_file.f90 .
cd -
\end{verbatim}

Pour compiler et faire tourner ce programme, vous devez vous placer dans le répertoire\\
\pathcasscm.
Commencez par sauvegarder le fichier phys.nc créé initialement :
\begin{verbatim}
mv phys.nc phys_original.nc
\end{verbatim}

Puis vous executez la commande
\begin{verbatim}
./compile.sh
\end{verbatim}
qui va mettre à jour l'executable \texttt{lmdz1d.e} que vous pourrez ensuite executer à noveau pour créer un nouveau fichier \texttt{phys.nc}.

Vous pourrez comparer avec les résultats de la simulation de référence

\texttt{ln -s ../../6AL79/arm\_cu/hourly.nc .}
\begin{verbatim}
ferret
use hourly.nc
use phys.nc
reg/k=50:79
set v ul
fill/d=1 rhum
set v ur
fill/d=2 rh
set v ll
fill/d=1 temp
set v lr
fill/d=2 temp
\end{verbatim}

Si vous tracez le flux sensible et le flux latent en surface 'plot sens' ou 'plot flat' (même nom de variable dan les deux fichiers, mais l'axe des temps est différent),
vous verrez que les flux croissent au cours de la journée.
Vous pouvez également relancer la simulation en jouant avec les forçages.
Vous pouvez en particulier imposer un flux constant à la surface en ajoutant en haut
du fichier config.def par exemple, les lignes
\begin{verbatim}
flux_sensible=120.
flux_latent=700.
flag_forcing=1
\end{verbatim}
Les deux  premières lignes indiquent qu'on impose les flux sensible et
latent à des valeurs constantes.
La troisième clé fait qu'on supprime le forçage "grande échelle" de la
simulation 1D.
Avec ces réglages, on voit qu'on arrive à avoir une humidité relative
qui dépasse 1.

Avec 
\begin{verbatim}
flux_sensible=0.
flux_latent=0.
flag_forcing=1
flag_init=2
\end{verbatim}
on se retrouve avec des flux en surface nuls et une atmosphère sèche
et neutre ($\theta=300K$).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Diffusion turbulente et couche d'Eckman}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{\bf SUJET 1}

On commencera par coder les équations de la couche d'Eckman en
rajoutant au modèle le terme de Coriolis :
\begin{eqnarray}
\partial_t u&=& f ( v-v_g)\\
\partial_t v&=& -f ( u-u_g)\\
\end{eqnarray}
On prendra sur toute la hauteur de l'atmosphère $(u_g,v_g)=(10,0)$ m/s, ce
qui correspond également aux conditions initiales pour le vent.
En reprenant les conditions de cas neutre ci-dessus, faire un calcul de 
couche d'Eckman. Comparer éventuellement au résultat théorique.

Dans un second temps, on modifiera le calcul du coefficient de diffusion turbulente.
On se propose ici de calculer la diffusion turbulente avec un coefficient
de mélange qui dépend du nombre de Richardson. On pourra procéder en deux temps,
en codant d'abord le coefficient de diffusion pour un Richardson nul (atmosphère
neutre) puis en prenant en compte la stratification.

Les équations du modèle de la diffusion turbulente sont déjà codées (dans 'diffkz.F90').
Le flux turbulent d'une espèce conservée de concentration massique $q$ s'écrit :
\begin{equation}
F_q(z)=\rho\overline{w'q'}=-\rho K_z \frac{\partial q}{\partial z}
\end{equation}
(déjà codée).
En surface, on suppose que le flux de température potentielle
$\theta$ est imposé, $F_\theta(0)=H / (\rho C_p)$.
\def\dt#1{\frac{\partial #1}{\partial t}}
\def\cisail{\frac{\partial\bf{v}}{\partial z}}
\def\emin{e_{\mbox{min}}}
La tendance pour 
\begin{equation}
\dt{q}=-\frac{1}{\rho}\frac{\partial q}{\partial z}
\end{equation}

On commencera par coder ce schéma sur la base de différences finies
centrées en $z$ et d'un schéma explicite temporel.

On testera notamment un schéma dans lequel le coefficient de diffusion
turbulente dépend du cisaillement de vent, d'une longueur de mélange
\def\lmix{l_{\mbox{mix}}}
\def\lref{l_0}
\begin{equation}
\lmix=\frac{\lref z}{\lref+z}
\end{equation}
et du nombre de Richardson 
\def\ri{Ri}
\def\ric{Ri_c}
\begin{equation}
\ri=\frac{N^2}{M^2}
\end{equation}
avec 
\begin{equation}
N^2=\frac{g}{\theta}\frac{\partial \theta}{\partial z}
\end{equation}
et
\begin{equation}
M^2=||\cisail||^2
\end{equation}
\begin{equation}
K_z=\lmix \sqrt{\max\depb{\lmix M^2 (1- \ri/\ric),\emin}}
\end{equation}

\def\lm{klev}

\paragraph{Convention d'indexage pour le codage}
Dans la maquette préétablie pour ce mini projet, on considère
$\lm$ couches au milieu desquelles on détermine les variables
d'état du système. Ces variables peuvent en fait être considérées
comme la moyenne de la variable sur la couche en question.
Les flux ou les gradients verticaux des variables sont
calculés aux interfaces entre deux couches.
On supposera que la surface (interface au bas de la première
couche) correspond à l'indice 0. Les interfaces en bas et en
haut de la couche $k$ sont donc indexées par $k$ et $k+1$.

{\bf
Tester le schéma avec des pas de temps de 1 à 100 s pour voir la
limite de stabilité. Relier cette limite aux valeurs de $K_z$ explorées
pendant la simulation.

Effectuer des simulations en conditions neutres ('iflag\_init=1' et 'flux\_sensible=0') avec 
rappel géostrophique pour regarder l'effet de différentes hypothèses sur les
profils de vent dans la couche limite. Comparer les cas $K_z=$cste et $K_z=l M$. Refaire
les calculs avec un flux latent en surface non nul, puis avec un calcul réaliste du
coefficient de traîné en surface, récupéré du sujet 2. On pourra enfin refaire les calculs
en branchant le modèle de panache thermique du sujet 3.}

\subsection{Diffusion turbulente schéma implicite en temps}

On montrera que ce schéma se ramène à une matrice tridiagonale.
On utilisera cette propriété pour développer une version implicite
en temps de ce schéma.

On peut pour ce faire déterminer les coefficients de la matrice 
et utiliser ou coder une méthode classique d'inversion.

On peut aussi inverser cette matrice avec une méthode particulière,
utilisée en pratique dans les modèles de climat car elle permet
un couplage avec par exemple un modèle de surface.

Dans ce second cas on commencera par introduire les coefficents
$A_k$ et $B_k$ tels que 
$q_{k}=A_k q_{k-1} + B_k$ et on établiera à partir de l'équation
de la diffusion une relation récurrente permettant de calculer
$(A_k,B_k)$ à partir de $(A_{k+1},B_{k+1})$.
La condition de flux turbulent nul au sommet permet de déterminer
les coefficients pour $klev$.

On calucle ensuite tous les coefficients par un schéma récursif descendant pour calculer finalement les flux au sol.

En imposant le flux en surface, on peut ensuite déterminer
la valeur de la grandeur au pas de temps $t+\delta t$ dans la
couche 1 puis remonter sur la colonne atmosphérique pour calculer
les variables d'état à partir des coefficients $A_k$ et $B_k$.

A partir de simulations avec différents pas de temps, tester la 
stabilité et la précision du schéma.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Calcul interactif de la température de surface et des flux}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{\bf SUJET 2}

Dans la maquette et plus généralement dans le cas ARM, les
flux latent et sensible sont imposés.

Il s'agira dans ce second sujet de développer un modèle de sol.

\def\wm{W/m$^2$}
\def\swdn{SW^\downarrow}
\def\swmx{F_0}
\def\lwdn{LW^\downarrow}
\def\lwup{LW^\uparrow}
On supposera que le flux radiatif descendant à la surface est la somme
d'un flux infra-rouge constant $LW^\downarrow$=400~\wm et d'un 
flux solaire avec un cycle diurne
$\swdn=\swmx\mbox{max}(\cos(2\pi t/T+\pi),0)$ où $t$ est le temps et
$T$ la durée du jour, avec $\swmx$=900~\wm.

\def\Cd{C_D}
\def\Cdh{C_{D,h}}
\def\modv{||\vec{V_1}||}
On calculera le flux montant radiatif à la surface $\lwup=\sigma T_s^4$, ainsi que les flux sensible $H$, latent $LE$ et les tensions de vent $\tau_x=\rho \overline{w'u'}$ et 
$\tau_y$ à partir des formules
\begin{eqnarray}
\tau_x &=& \Cd \modv (u_1 -u)\\
\tau_y &=& \Cd \modv (v_1 -v)\\
H     &=& C_p \Cdh \modv ( T-s-T_1)\\
LE    &=& \beta L \Cdh \modv (qsat(T_s)-q_1)\\
\Cd   &=& \left({\frac{\kappa}{\log {z/z_0}}}\right)^2\\
\Cdh  &=& 0,7 \Cd
\end{eqnarray}
où $\modv$ est le module du vent horizontal dans la première couche du modèle.
Le coefficient $\beta$ vaut 1 pour une surface d'eau libre et 0 pour un désert sans
eau dans le sol.

On commencera par coder ces calculs de flux en supposant que la température
de la surface $T_s$ est 1~K plus chaude que la température de la première couche du modèle
$T_1$.

Dans un second temps, on assimilera la surface à une couche unique de capacité
calorifique surfacique $S_s$ de telle sorte que la température de surface évolue
suivant l'équation
\begin{equation}
\partial_t T_s= \Sigma F_s / C_s
\end{equation}
où $\Sigma F_s$ représente la somme des flux en surface.

On pourra ensuite de façon optionnelle coder un modèle de conduction thermique dans le sol.

Pour ce faire, on calculera explicitement la conduction 
thermique dans le sol. Pour ce faire, on considérera les équations
donnant le flux thermique $F_s$ et l'évolution de la température du sol
\begin{eqnarray}
F_s&=&- K_s \partial_z T_s\\
C_t \partial_t T_s &=& - \partial_z F_s
\end{eqnarray}
où $K_s$ et $C_t$ sont respectivement la conductivité et la capacité thermique du
sol. On prendra comme valeur 2~W/mK $C_s=2\ 10^6$~J~m$^{-3}$~K$^{-1}$.

On discrétisera le premier mètre de sol avec des couches de 1cm.
On supposera que le centre de la première couche se trouve à la surface
et que ça température correspond à la température d'émission de la surface
dans l'infra-rouge.

On écrira en surface la continuité des flux
\begin{equation}
\lwdn+\swdn -\sigma T_s^4 -H-LE = F_s
\end{equation}
(avec cette convention de signe, $F_s$ est positif vers le bas, donc on suppose
qu'on a orienté les profondeurs vers le bas aussi).

On supposera que la température de la dernière couche a une valeur constante
(300~K) par exemple, ou alors que le flux sous le dernier niveau de température
est nul.

{\bf On prolongera les simulations sur 3 jours pour analyser la représentation
du cycle diurne. 
On branchera le calcul de la diffusion verticale issu du sujet 2.
On regardera l'influence de la capacité calorifique (en multipliant
ou divisant la valeur nominale par 2). On regardera également l'influence du calcul
du coefficient d'aridité $\beta$ en testant des valeurs entre 0 et 1. Enfin
on activera le modèle du panache thermique pour regarder l'influence du transport
convectif dans l'atmosphère sur la représentation de ce cycle diurne.}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Développement d'un modèle de panache thermique}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{\bf SUJET 3}

On va développer ici un schéma en flux de masse des structures
convectives de la couche limite.
Pour ce faire, on va séparer la maille horizontalement en deux
sous-partie, l'une associée aux "ascendances thermiques" et
l'autre aux subsidences.
On supposera pour simplifier le modèle que la partie ascendante couvre
une fraction imposée $\alpha=0.1$ de la maille.


\def\th{\mbox{th}}
\def\wth{w}
\def\fth{f}
\def\eth{e}
\def\dth{d}
\def\qth{q_{\th}}
\def\thetath{\theta_{\th}}

La partie ascendante sera caractériée par une vitesse verticale
$\wth$ et un flux de masse $\fth= \alpha \rho \wth$.

Le flux de masse échange de l'air avec l'extérieur au travers d'un entrainement latéral $\eth$ et d'un détraimenent $\dth$, de sorte
que le panache, supposé stationaire obéit à l'équation de continuité
pour l'air
\begin{equation}\label{eq:conser}
\der{f}{z}=e-d
\end{equation}

On pourra avec cette équation calculer la concentration d'une variable
$q$ dans l'ascendance $\qth$, en écrivant l'équation de continuité
pour le traceur :
\begin{equation}\label{eq:conserq}
\der{f\qth}{z}=e q-d \qth
\end{equation}
où on suppose que l'air entrainé dans le panache à comme concentration
la concentration moyenne dans la maille, $q$.

Sur la base de cette équation, on pourra calculer en particulier
la température potentielle dans le panache ascendant, $\thetath$, qui
permettra de calculer la flotabilité du panache par rapport à
l'environnement
\begin{equation}
\gamma=g\frac{\thetath-\theta}{\theta}
\end{equation}
et d'en déduire la vitesse verticale au travers de l'équation de
conservation du moment sur la verticale
\begin{equation}
\der{f\wth}{z}= -d \wth + \rho \gamma
\end{equation}
Qu'on simplifiera en supposant $\rho \alpha=cste$ et $\dth=0$.
On écrira les équations de conservation avec des schémas amont pour
le calcul des termes de transport sur la verticale.

Pour calculer la flotabilité dans la première couche, on supposera
que l'air dans la partie montante de la maille a comme température
potentielle une extrapolation linéaire au sol de la température
potentielle des deux premières couches.

Commencer par écrire le schéma sur le papier.
Il faudra intégrer dans un premier temps 
depuis la première couche vers le sommet de
l'atmosphère le calcul de $\wth$, de $\gamma$ et de $\thetath$.
On déterminera indirectement $e$ et $d$ à partir de l'équation
\eq{conser} comme
\begin{eqnarray}
e&=&\max\dep{\der{f}{z},0}+\epsilon\\
d&=& \max\dep{-\der{f}{z},0}+\epsilon
\end{eqnarray}

Une fois les propriétés du thermique calculées, la contribution 
aux flux s'écrira sous la forme
\begin{equation}
F_q= \fth (\qth -q)
\end{equation}
où la valeur de $\qth$ dans l'ascendance aura été calculée au 
travers de l'\eq{conserq}.

{\bf On branchera dans le modèle le calcul du coefficient $K_z$ réalisé dans
le sujet 1. On comparera sur un cas simple avec flux de surface constants
le transport vertical dans la couche limite, avec et sans modèle de panache
thermique. On regardera en particulier la relation entre le signe du gradient vertical
de température potentielle et du flux vertical de chaleur $w'\theta'$ (en prenant
en compte la somme du flux turbulent et du flux dans le panache thermique).
On comparera également dans les deux cas le transport d'eau (on mettra pour
ce faire une vapeur d'eau nulle initialement et un flux latent non nul en surface).
On pourra réitérer ce calcul avec les un cycle diurne en branchant le calcul 
des bilans de surface du sujet 1}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Développement d'un panache humide et schéma de nuage}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{\bf SUJET 4}

On commencera par calculer la vitesse verticale d'une ascendance
en supposant qu'une particule d'air par de la première couche de surface
avec une vitesse ascendante nulle, et qu'elle accélère en fonction de la
flotabilité :
\begin{eqnarray}
w \partial_z w &=& B \\
B&=&g\frac{\theta_{asc} -\theta}{\theta}
\end{eqnarray}
On coupera le calcul quand la vitesse verticale $w$ devient négative.
Sortir les valeurs de w.

On supposera ensuite que le coefficient de diffusion turbulente
est $K_z=l w$ avec l=10~m.
$w$ étant nul en haut de la première couche a priori, on prendra ${K_z}_3={K_z}_2/2.$.
Refaire le calcul de diffusion turbulente.

On prendra par exemple pour ces question comme forçage 
\begin{verbatim}
flux_sensible=50.
flux_latent=300.
flag_forcing=1
\end{verbatim}
Il faut éviter notamment que les ascendances aillent jusqu'au sommet
du modèle.

On introduira ensuite une largeur de distribution sous maille de l'eau
dans la maille de largeur $\sigma$. On considérera par exemple
une distribution "chapeau", uniforme entre $q-\sigma/2$ et $q+\sigma/2$.

On utilisera cette distribution pour calculer la fraction nuageuse $f$,
c'est à dire la fraction de la maille couverte par des humidités plus grandes
que $q_{\mbox{sat}}$ (la fonction 'f\_qsat' de la maquette).

On contrôlera la largeur de la distribution par le rapport $r=q/\sigma$.
On considérera  successivement trois hypothèses : dans la première $r=0.2$ partout, dans
la seconde $r=0.2*z$ et dans la troisième, on supposera que ce paramètre
est proportionnel à la différence d'humidité entre l'humidité moyenne 
de la couche et l'humidité dans le panache ascendant, qu'on supposera
elle-même égale à l'humidité dans la première couche du modèle.

{\bf
On pourra ensuite brancher le schéma de condensation  dans la version avec
le modèle des thermiques
(SUJET 3) en supposant que l'humidité dans le panache est celle calculée
par le panache thermique.

On branchera dans ce modèle le modèle de surface du sujet 2 pour regarder le
cycle diurne des cumulus de couche limite. On ajustera le coefficient $\beta$ 
d'aridité de la surface pour produire plus ou moins de nuages. 

On pourra également coder une advection verticale vers le bas de l'eau totale,
$\partial_t q=-w\partial z q$ pour singer l'effet d'une subsidence à grande
échelle comme on en rencontre sur les régions de strato-cumulus.

On pourra également prendre en compte le changement de phase de l'eau
dans le calcul de la flotabilité dans le panache thermique.}



\section{Transfert radiatif}

\def\ir{{_{IR}}}
\def\vis{{_{Vis}}}
\def\upw{^{\uparrow}}
\def\dnw{^{\downarrow}}
%\newcommand{\trans}[2]{ \epsilon \left( #1 , #2 \right) }
\newcommand{\trans}[2]{ \tau_{  #1 , #2 } }
\def\dz{\delta_z}

On se placera pour le calcul du transfert radiatif dans l'hypothèse dite "grise", qui consiste à ne pas tenir compte du spectre d'absorption des gaz, c'est à dire à supposer que l'absorption est aussi efficace quelque soit la longueur d'onde.
On suppose en revanche que le coefficient d'absorption peut être différent dans le visible et dans l'infra-rouge.
On commencera à coder la partie infra-rouge dans la mesure où on peut traiter simplement la partie solaire en pensant qu'une fraction du rayonnement est réfléchit vers l'espace avec l'albédo dit planétaire, et que le reste arrive directement à la surface où il est entièrement absorbé. Dans ce cas, il n'y a pas de chauffage direct de l'atmosphère par le rayonnement solaire.

L'évolution temporele de la température
($T_{k+\frac{1}{2}}$)
au milieu d'une couche due à l'absorption et à l'émission de rayonnement thermique infra-rouge s'écrit alors :
\begin{equation}
\left( \dt{ T_{k+\frac{1}{2}} } \right)\ir =
\frac{g}{C_p}
\frac{ F_k^{\upw}-F_k^{\dnw}  - F_{k+1}^{\upw}+F_{k+1}^{\dnw} }
{ p_k -p_{k+1} }
\end{equation}
où les flux montant $F_k^{\upw}$ et descendant $F_k^{\dnw}$ sont calculés en utilisant une transmissivité intégrée s
$\trans{}{}$ in the following form
\begin{equation}
F_k^{\upw}  =  \epsilon_{s} \sigma {T_s}^4  \trans{0}{k}
+ \sum_{l=0}^k { \sigma T_{l+\frac{1}{2}}^4
\left[ \trans{l+1}{k} -\trans{l}{k} \right]  }
\end{equation}
\begin{equation}
F_k^{\dnw} = \sum_{l=k}^K 
{ \sigma T_{l+\frac{1}{2}}^4
\left[ \trans{k}{l} -\trans{k}{l+1} \right]  }
\end{equation}
(où $\epsilon_s$ est l'émissivité de surface).


Différentes formes simples peuvent être choisies pour la transmissivité.  L'utilisation d'une approximation d'absorption grise à proprement parlé, conduit à de fortes discontinuités de température à la surface (entre la température du sol et celle de l'air au-dessus de la surface), comme l'avaient déjà constaté d'autres auteurs au début de la modélisation numérique du climat~\cite{Mana:61}.

Elle s'écrit :
\def\up{\overline{up}}
\def\u{\overline{u}}
\begin{equation}
\trans{k}{l}=\exp{ \left( - k_\ir\u \right) }
\end{equation}
%
où $\up$ est la masse d'air traversée
%
\begin{equation}
\u = \int { \rho dz}=\frac{\left( p_k -p_l \right)}{g}
\end{equation}
%
et $k_\ir$ est une constante rendant compte de l'efficacité d'absorption du gaz.


Les transmissivités peuvent également être calculées dans le cadre de l'approximation de l'absorption forte, mieux adaptée aux conditions de la troposphère où la plupart des bandes moléculaires efficaces sont fortement saturées. Dans ce cas, en supposant que les raies sont distribuées aléatoirement sur le spectre, et en négligeant les variations de température des coefficients d'absorption le long du chemin optique, la transmissivité de la bande entière, entre deux couches $k$ et $l$, peut être calculée, en supposant un rayonnement complètement isotrope, sous la forme (pour ceux qui ne sont pas familiers avec de tels modèles, voir Goody, 1964) :

\def\up{\overline{up}}
\def\u{\overline{u}}
\begin{equation}
\trans{k}{l}=\exp{ \left( - k_\ir\up \right) }
\end{equation}
%
où $\up$ est une quantité d'absorbant pondérée par la pression
%
\begin{equation}
\up = \int { p \rho dz}=\left( p_k -p_l \right) 
\frac{\left( p_k +p_l \right) }{2g}
\end{equation}
%
et $k_\ir$ est à nouveau une constante rendant compte de l'efficacité d'absorption du gaz.

Cette formule produit à nouveau une certaine discontinuité de température à la surface, mais inférieure à 10~K.  La transmissivité dépend d'un seul paramètre libre, par exemple la transmission optique totale de l'atmosphère
$\tau_\ir=\exp\dep{-k_\ir p_s^2/2}$
(fraction du rayonnement thermique de surface directement émise vers l'espace). Le modèle (que ce soit en approximation grise ou d'abosrption forte) prend en compte l'effet fondamental des variations de température liées à la fonction de Planck, essentielles pour l'équilibre thermique et le forçage radiatif de la circulation, mais pas les dépendances beaucoup plus complexes de pression et de température liées aux coefficients d'absorption moléculaire.

Une approximation classique d'absorption grise est utilisée également pour l'absorption des radiations solaires.  Pour des raisons de simplicité, on ne tentera pas de prendre en compte la diffusion qui complique énormément le calcul du transfert radiatif.
Dans ces conditions, le rayonnement solaire descendant consiste en un faisceau collimaté unique et la transmission grise prend la forme simple suivante
\begin{equation}
\trans{k}{l}=\exp\dep{-\frac{k_\vis}{\cos\theta}\left| p_k - p_l \right| }
\end{equation}
où $\theta$ est l'angle zénithal.
Comme pour le transfert radiatif dans l'infrarouge, le modèle dépend d'un paramètre libre, par exemple la transmission atmosphérique totale
$\tau_\vis$, qu'on pourra calculer pour une incidence moyenne
$\cos\theta=1/2$
(fraction du rayonnement solaire entrant global atteignant la surface), qui est la fraction du flux solaire atteignant la surface de la planète. Elle influence grandement la stabilité verticale de l'atmosphère, un paramètre dynamique important.

Le rayonnement réfléchi, qui dépend de l'albédo de la surface $A_s$, est supposé totalement isotrope dans la direction ascendante et est calculé pour un angle zénithal moyen $\theta_0$ avec $1/\cos\theta_0=1.66$.

Le forçage radiatif de l'atmosphère dépend également du flux solaire entrant au sommet de l'atmosphère $F_0$ et de ceux des paramètres orbitaux qui déterminent les variations spatiales et temporelles de l'insolation (durée du jour, obliquité et durée de l'année).

Sur le site de ARM, on pourra tester des formules comme
\begin{verbatim}
    swdntoa=fsw0*max(cos(2*rpi*(timesec-43200.)/(1.15*86400.)),0.)
\end{verbatim}
qu'on calera sur des simulations avec le modèle 1D standard en réactivant le rayonnement dans la simulation de référence, tournée avec la physique \texttt{phylmd}. Pour cela, on changera \texttt{iflag\_radia=1} dans le fichier \texttt{lmdz1d.def}.


\end{document}

