1 |
\documentclass[a4paper]{article} |
2 |
|
3 |
% Packages |
4 |
|
5 |
\usepackage[latin1]{inputenc} |
6 |
\usepackage{ifpdf} |
7 |
|
8 |
\ifpdf |
9 |
\usepackage{aeguill} |
10 |
\usepackage[pdftex]{color,graphicx} |
11 |
\else |
12 |
\usepackage[T1]{fontenc} |
13 |
\usepackage[dvips]{color,graphicx} |
14 |
\usepackage[active]{srcltx} |
15 |
\fi |
16 |
|
17 |
\usepackage[frenchb]{babel} |
18 |
\usepackage{mathrsfs} |
19 |
\usepackage{amsmath} |
20 |
\usepackage{amsfonts} |
21 |
\usepackage{mathabx} |
22 |
\usepackage{algorithmic} |
23 |
\usepackage[all]{xy} |
24 |
|
25 |
|
26 |
% Load last: |
27 |
\ifpdf |
28 |
\usepackage{hyperref} |
29 |
\hypersetup{pdftex, pdftitle={Manuel pour LMDZE}, pdfauthor={Lionel |
30 |
Guez}, pdfstartview=FitBH} |
31 |
\else |
32 |
\usepackage[hypertex]{hyperref} |
33 |
\fi |
34 |
|
35 |
%--------------- |
36 |
|
37 |
\newcommand{\ud}{\mathrm{d}} |
38 |
\DeclareMathOperator{\diverg}{div} |
39 |
\DeclareMathOperator{\tgh}{th} |
40 |
|
41 |
\newcommand{\Eng}[1]{\textit{\foreignlanguage{english}{#1}}} |
42 |
% (English language) |
43 |
|
44 |
\newcommand{\pseudov}[1]{\overset{\curvearrowbotright}{#1}} |
45 |
|
46 |
\renewcommand{\algorithmicdo}{\textbf{faire}} |
47 |
\renewcommand{\algorithmicend}{\textbf{fin}} |
48 |
\renewcommand{\algorithmicfor}{\textbf{pour}} |
49 |
|
50 |
\author{Lionel GUEZ} |
51 |
\title{Manuel pour LMDZE} |
52 |
|
53 |
%--------------- |
54 |
|
55 |
\begin{document} |
56 |
|
57 |
\maketitle |
58 |
\ifpdf |
59 |
\else |
60 |
\tableofcontents |
61 |
\fi |
62 |
\listoffigures |
63 |
\newpage |
64 |
|
65 |
\part{Installation} |
66 |
|
67 |
LMDZE utilise quatre bibliothèques. Cf. figure (\ref{fig:libraries}). |
68 |
\begin{figure}[htbp] |
69 |
\centering |
70 |
\includegraphics{Graphiques/libraries} |
71 |
\caption{Bibliothèques utilisées par LMDZE.} |
72 |
\label{fig:libraries} |
73 |
\end{figure} |
74 |
|
75 |
Pour télécharger le programme : |
76 |
\begin{verbatim} |
77 |
svn checkout \ |
78 |
svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/lmdze/svn/trunk \ |
79 |
LMDZE_program |
80 |
\end{verbatim} |
81 |
|
82 |
Pour compiler \verb+gcm+ : |
83 |
\begin{enumerate} |
84 |
\item Choisir la résolution spatiale dans \verb+dimens_m.f+. Choisir |
85 |
la valeur de \verb+kdlon+ dans \verb+phylmd/raddim.f+. Pour cela, |
86 |
consulter le tableau (\ref{tab:kdlon}), les |
87 |
\href{http://www.lmd.jussieu.fr/~lmdz/LMDZ4/choisir_une_resolution.html}{conseils} |
88 |
et le programme |
89 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Tests/choixdim.f}{\texttt{choixdim}}. |
90 |
(Il faudra adapter les pas de temps à l'exécution.) |
91 |
\item Vérifier la cible du lien \verb+compiler.mk+. Vérifier dans |
92 |
\verb+make.sh+ la valeur de \verb+dest_dir+. |
93 |
\item \verb+make.sh+. (Si on a modifié les dimensions alors il faudra |
94 |
créer un nouvel état initial adapté à ces dimensions. Il faut donc |
95 |
recompiler \verb+etat0_lim+.) |
96 |
\end{enumerate} |
97 |
\begin{table}[htbp] |
98 |
\centering |
99 |
\begin{tabular}{lll} |
100 |
\texttt{iim} & \texttt{jjm} & \texttt{kdlon} \\ |
101 |
\hline |
102 |
32 & 24 & 41 \\ |
103 |
48 & 32 & 149 \\ |
104 |
72 & 46 & 1621 \\ |
105 |
96 & 72 & 487 \\ |
106 |
160 & 98 & 78 \\ |
107 |
192 & 43 & 10 |
108 |
\end{tabular} |
109 |
\caption{Choix de \texttt{kdlon} dans \texttt{phylmd/raddim.f}.} |
110 |
\label{tab:kdlon} |
111 |
\end{table} |
112 |
|
113 |
La compilation de \verb+clmain.F+ avec le compilateur IBM XL Fortran |
114 |
Enterprise Edition V9.1 Version: 09.01.0000.0005 et les options |
115 |
\verb+-O5 -qnostrict -qessl+ n'est pas finie en 10 h. Avec les options |
116 |
\verb+-O4 -qnostrict -qessl+, elle n'est pas finie en 1 h. |
117 |
|
118 |
\part{Utilisation comme boîte noire} |
119 |
|
120 |
Cf. liste des fichiers d'entrée nécessaires sur le |
121 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
122 |
IOIPSL LMDZ/inout LMDZE.odg}{schéma des entrées-sorties}. Les |
123 |
fichiers \verb+ECDYN.nc+ et \verb+ECPHY.nc+ contiennent des données du |
124 |
centre européen, pour un jour particulier, avec une résolution de 1°. |
125 |
|
126 |
Le fichier \verb+amipbc_sic_1x1.nc+ donne la fraction de |
127 |
la surface de la mer couverte par de la glace. (Cette définition |
128 |
suppose que la fraction en question est indépendante de l'échelle |
129 |
spatiale considérée.) |
130 |
|
131 |
Cf. la |
132 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
133 |
IOIPSL LMDZ/namelists.ods}{liste des paramètres d'entrée}. Avant de |
134 |
lancer l'exécution de \verb+gcm+, il faut en particulier choisir |
135 |
\verb+nday+ (par exemple 1 ou 30). Je peux aussi modifier |
136 |
\verb+iecri+, qui contrôle la fréquence d'écriture dans le fichier |
137 |
\verb+histday.nc+, et \verb+periodav+, qui contrôle la fréquence |
138 |
d'écriture dans le fichier \verb+histmth.nc+ (?). |
139 |
|
140 |
Pour ne pas utiliser Orchidée, choisir \verb+VEGET = n+ |
141 |
dans \verb+physiq.def+. Si on intègre sur 1 jour, choisir |
142 |
\verb+OK_journe = y+, \verb+OK_mensuel = n+. |
143 |
|
144 |
Si on a changé la résolution spatiale, il peut être nécessaire de |
145 |
modifier les variables \verb+nfilun+, \verb+nfilus+, \verb+nfilvn+ et |
146 |
\verb+nfilvs+, définies dans \verb+filtrez/parafilt.h+. |
147 |
\verb+etat0_lim+ se plante et indique alors les bonnes valeurs. |
148 |
|
149 |
\verb+idissip+ dans \verb+gcm.def+ n'est pas pris en compte. Il est |
150 |
automatiquement calculé par le programme. Phu L. V. conseille de |
151 |
diminuer à 3600 les valeurs de \verb+tetagdiv+ et \verb+tetagrot+. |
152 |
|
153 |
\verb+grossismx+ et \verb+grossismy+ doivent être inférieurs à 4 pour |
154 |
que le modèle converge. |
155 |
|
156 |
Phu L. V. me conseille de choisir \verb+ysinus=n+ dans \verb+gcm.def+. |
157 |
|
158 |
\section{Paramètres liés au pas de temps} |
159 |
|
160 |
Les paramètres à choisir au moment de l'exécution et liés au pas de |
161 |
temps sont \verb+day_step+, \verb+iphysiq+, \verb+iperiod+ et |
162 |
\verb+iconser+. Il faut en particulier penser à ajuster ces paramètres |
163 |
lorsque l'on a changé la résolution spatiale. Le principe de base est |
164 |
que le pas de temps doit être suffisamment petit pour que le modèle |
165 |
soit stable. Il doit être d'autant plus petit que la résolution |
166 |
spatiale horizontale est fine. Si les milieux de couches ne montent |
167 |
pas trop haut dans la stratosphère alors la résolution verticale n'a |
168 |
pas à être prise en compte. La référence est un pas de temps maximal |
169 |
de 6 mn environ pour une grille de 64 longitudes par 50 latitudes. |
170 |
Lorsqu'on multiplie par un même facteur le nombre de longitudes et de |
171 |
latitudes, on divise par ce facteur le pas de temps maximal. Pour des |
172 |
milieux de couches haut dans la stratosphère, cf. |
173 |
\href{file:///user/guez/Documents/Utilisation_LMDZ/Utilisation_LMDZ.texfol/Utilisation-LMDZ.dvi}{notes |
174 |
sur l'utilisation de LMDZ}. |
175 |
|
176 |
\verb+day_step+ est le nombre de pas de temps par jour. \verb+iphys+ |
177 |
est le nombre de pas de temps entre deux appels de la physique. |
178 |
\verb+iperiod+ est le nombre de pas de temps entre deux pas Matsuno. |
179 |
\verb+iconser+ est le nombre de pas de temps entre deux écritures des |
180 |
variables de contrôle. Ces quatre paramètres doivent recevoir des |
181 |
valeurs entières. Apparemment, on laisse toujours \verb+iperiod+ à 5 |
182 |
(quelle que soit la résolution). |
183 |
|
184 |
On a en outre la contrainte que \verb+day_step+ soit un multiple de |
185 |
\verb+iperiod+ (cf. l'unité de programme principale \verb+gcm+), c'est-à-dire |
186 |
normalement un multiple de 5. |
187 |
|
188 |
Enfin, on a la contrainte suivante sur \verb+iphysiq+. Dans le fichier |
189 |
\verb+fisrtilp.F+, ligne 133, on trouve : |
190 |
\begin{verbatim} |
191 |
IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN |
192 |
PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime |
193 |
PRINT*, 'Je prefere un sous-intervalle de 6 minutes' |
194 |
ENDIF |
195 |
\end{verbatim} |
196 |
Donc, comme \verb+ninter+ vaut 5, il faut que \verb+dtime+, le pas de |
197 |
temps de la physique, vaille 1800 s à 5.10$^{-3}$ s près. Notons $d$ |
198 |
le pas de temps de la dynamique, en mn. Il faut donc que |
199 |
$\mathtt{iphysiq} \times d = 30$. |
200 |
|
201 |
En résumé, en notant $d_M$ le pas de temps maximal imposé par la |
202 |
résolution, nous avons le système de contraintes suivant : |
203 |
\begin{displaymath} |
204 |
\left\{ |
205 |
\begin{array}{l} |
206 |
\mathtt{iphysiq} \times d = 30 \\ |
207 |
d \le d_M \\ |
208 |
\mathtt{day\_step} \times d = 1440 \\ |
209 |
\mathtt{day\_step} \equiv 0[5] |
210 |
\end{array} |
211 |
\right. |
212 |
\end{displaymath} |
213 |
D'où : |
214 |
\begin{displaymath} |
215 |
\left\{ |
216 |
\begin{array}{l} |
217 |
\mathtt{iphysiq} \ge = 30 / d_M \\ |
218 |
48 \times \mathtt{iphysiq} \equiv 0[5] |
219 |
\end{array} |
220 |
\right. |
221 |
\end{displaymath} |
222 |
C'est-à-dire : |
223 |
\begin{displaymath} |
224 |
\left\{ |
225 |
\begin{array}{l} |
226 |
\mathtt{iphysiq} \ge = 30 / d_M \\ |
227 |
\mathtt{iphysiq} \equiv 0[5] |
228 |
\end{array} |
229 |
\right. |
230 |
\end{displaymath} |
231 |
On peut donc adopter l'algorithme suivant : |
232 |
\begin{verbatim} |
233 |
entrer(iim, jjm) |
234 |
n := ceiling(max(iim / 64., jjm / 50.)) |
235 |
iphysiq := 5 * n |
236 |
day_step := 240 * n |
237 |
\end{verbatim} |
238 |
|
239 |
\section{Ressources nécessaires} |
240 |
|
241 |
Espace disque par utilisateur : |
242 |
\begin{itemize} |
243 |
\item fichiers sources LMDZ (sans IOIPSL ni NetCDF) 3 Moctets ; |
244 |
\item fichiers objets, exécutables (sans IOIPSL ni NetCDF) (résolution |
245 |
$32 \times 24 \times 9$) 4 Moctets ; |
246 |
\item fichiers sortie \verb+etat0_lim+ 9 Moctets ; |
247 |
\item fichiers sortie \verb+gcm+ (principalement \verb+histday.nc+) 26 |
248 |
Moctets. |
249 |
\end{itemize} |
250 |
Par ailleurs, on peut centraliser : |
251 |
\begin{itemize} |
252 |
\item fichiers NetCDF entrée \verb+etat0_lim+ 40 Moctets ; |
253 |
\item bibliothèque IOIPSL 600 koctets ; |
254 |
\item ensemble NetCDF (bibliothèque, utilitaires, manuels, modules, |
255 |
etc.) 3 Moctets. |
256 |
\end{itemize} |
257 |
|
258 |
Mémoire principale (résolution $32 \times 24 \times 9$) : 30 Moctets. |
259 |
|
260 |
Durées. Temps de compilation avec \verb+g95 -O3+ : 5 mn. Temps |
261 |
d'exécution de \verb+gcm+ avec \verb+g95 -O3+ (résolution $32 \times |
262 |
24 \times 9$) : environ 10 mn par mois simulé. |
263 |
|
264 |
\part{Fonctionnement : derrière les coulisses} |
265 |
|
266 |
\begin{verbatim} |
267 |
wc -l `find . -name "*.[fh]" -o -name "*.f90"` |
268 |
\end{verbatim} |
269 |
donne environ $7 \cdot 10^{4}$ lignes. |
270 |
|
271 |
\section{Intérêt de placer un makefile dans chaque répertoire} |
272 |
|
273 |
J'admets l'intérêt de séparer les fichiers sources dans les |
274 |
répertoires \verb+dyn3d+, \verb+bibio+ etc. Avec un aussi grand nombre |
275 |
de fichiers, ce classement apporte de la clarté. À partir de là se |
276 |
pose la question de l'intérêt de placer un \verb+makefile+ dans chaque |
277 |
répertoire. Cette organisation suppose déjà l'absence de dépendance |
278 |
circulaire entre les répertoires pour la compilation. Admettons qu'il |
279 |
n'y en ait pas. Il y aurait logiquement un \verb+makefile+ au dessus |
280 |
de ces répertoires pour les éditions de liens. Je suis arrivé à la |
281 |
conclusion qu'en bon style de \verb+makefile+ des archives |
282 |
n'apparaissent que dans les options d'édition de liens et ne sont pas |
283 |
connues directement du \verb+makefile+ du programme utilisateur. |
284 |
L'utilisation des archives ne permettrait donc pas la mise à jour des |
285 |
deux programmes, \verb+gcm+ et \verb+etat0_lim+, lorsqu'un fichier |
286 |
d'un répertoire est modifié. J'écarte donc la création d'une archive |
287 |
pour chaque répertoire. Par ailleurs, chacun des deux programmes |
288 |
utilise des fichiers répartis dans les différents répertoires. Pour |
289 |
chacun des deux programmes, la liste des objets pertinents doit donc |
290 |
apparaître en prérequis. On ne pourrait mettre simplement en prérequis |
291 |
des répertoires. Mais alors comment le \verb+makefile+ supérieur |
292 |
sait-il mettre à jour ces objets ? En conclusion, je crois qu'il ne |
293 |
faut pas placer un \verb+makefile+ dans chaque répertoire. Un unique |
294 |
\verb+makefile+ au dessus des répertoires est la bonne solution. |
295 |
|
296 |
\section{Description des programmes} |
297 |
|
298 |
Sur le suffixe des fichiers. Pour pouvoir utiliser les outils NAG, |
299 |
tant que le programme est formé d'un mélange de fichiers aux formats |
300 |
fixe et libre, je suis obligé d'associer le format libre au suffixe |
301 |
\verb+f90+ ou \verb+f95+. \verb+sxf90+ m'impose le suffixe \verb+f90+. |
302 |
Pour le format fixe, \verb+sxf90+ m'impose \verb+f+. Cf. |
303 |
\hyperref{file:///user/guez/Documents/Informatique_fonctionnement/Informatique, |
304 |
fonctionnement.texfol/fonctionnement.dvi}{text}{suffixes}{contraintes |
305 |
sur les suffixes de fichiers de Fortran}. |
306 |
|
307 |
Deux unités de programme principales en Fortran, dans les fichiers |
308 |
\verb+etat0_lim.f+ et \verb+gcm.f+. |
309 |
|
310 |
Au plus haut niveau, parties du programme \verb+gcm+ : |
311 |
calcul des tendances (\verb+caldyn+), intégration (\verb+integr+), |
312 |
physique, dissipation (\verb+dissip+), advection. Variables |
313 |
principales : \verb+ucov+, \verb+vcov+, température, |
314 |
pression. \verb+w+ se déduit par quelque chose comme une |
315 |
relation de fermeture. Dans la physique, il n'y a pas d'échange |
316 |
horizontal entre les mailles. Un problème fondamental est l'utilisation |
317 |
d'une grille rectangulaire en latitudes et longitudes pour couvrir |
318 |
la sphère. Les mailles de la grille n'ont pas la même surface. |
319 |
La solution choisie pour résoudre ce problème est le filtrage. |
320 |
Principe du filtrage : décomposition d'un champ en harmoniques |
321 |
sphériques et conservation d'une partie seulement du spectre. |
322 |
Le filtre est très coûteux : il faut connaître le champ |
323 |
en entier pour le calculer en un point particulier. Par ailleurs, |
324 |
le programme calcule à de nombreuses occasions des valeurs moyennes |
325 |
aux pôles, à partir des mailles pôlaires (à chaque pôle, autant |
326 |
de mailles pôlaires que de longitudes dans la grille) |
327 |
|
328 |
\verb+etat0_lim+ crée l'état initial. \verb+etat0_lim+ appelle |
329 |
\verb+etat0+ et \verb+limit+. \verb+etat0+ crée les fichiers |
330 |
\verb+start.nc+ et \verb+startphy.nc+. \verb+limit+ crée |
331 |
\verb+limit.nc+. |
332 |
|
333 |
Le répertoire \verb+phylmd+ contient ce qui concerne la |
334 |
physique. Le répertoire \verb+dyn3d+ contient ce qui concerne |
335 |
la dynamique. |
336 |
|
337 |
\verb+histday.nc+ contient des moyennes journalières. |
338 |
\verb+ini_histday+ et \verb+write_histday+ sont responsables de |
339 |
l'écriture de \verb+histday.nc+. \verb+ini_histday+ initialise le |
340 |
fichier \verb+histday.nc+, puis \verb+write_histday+ écrit les |
341 |
contenus des variables dans ce fichier. Rôles analogues pour |
342 |
\verb+ini_histmth+ et \verb+write_histmth+ avec le fichier |
343 |
\verb+histmth.nc+. L'initialisation d'un fichier inclut |
344 |
l'initialisation de variables avec le sous-programme \verb+histdef+. |
345 |
L'écriture du contenu d'une variable se fait avec le sous-programme |
346 |
\verb+histwrite+. |
347 |
|
348 |
La lecture des fichiers d'entrée \verb+*.def+ se fait entre autres par |
349 |
l'intermédiaire du sous-programme \verb+getin+ de IOIPSL. \verb+getin+ |
350 |
n'est appelé que dans les fichiers \verb+phylmd/conf_phys.f+, |
351 |
\verb+dyn3d/conf_gcm.F+ et \verb+dyn3d/getparam.f+. |
352 |
|
353 |
La résolution des fichiers de conditions initiales et de conditions |
354 |
aux limites étant a priori différente de celle de LMDZ, LMDZ re-maille |
355 |
ces données par interpolation ou moyenne. |
356 |
|
357 |
Dans le module \verb+start_init_orog_m+, la variable \verb+phis+ est |
358 |
allouée et modifiée par \verb+start_init_orog+. Dans |
359 |
\verb+start_init_orog+, \verb+phis+ est passée en argument à |
360 |
\verb+grid_noro+ et remplie par \verb+grid_noro+. |
361 |
|
362 |
Dans \verb+grid_noro+, \verb+mask+ est rempli avec \verb+num_lan+ |
363 |
\verb+/+ \verb+num_tot+. On peut voir les valeurs obtenues en |
364 |
regardant la variable \verb+masque+ dans \verb+startphy.nc+. Le nom de |
365 |
masque devrait être réservé à des variables logiques. \verb+mask+ |
366 |
calculé par \verb+grid_noro+ n'est pas une variable logique et ne |
367 |
contient pas que les valeurs 0 ou 1. Ce doit être la fraction de terre |
368 |
dans la maille considérée. \verb+zdata+ sert à calculer \verb+zusn+. |
369 |
|
370 |
La variable \verb+iperiod+, qui est importante pour la structure de la |
371 |
procédure \verb+leapfrog+, n'est pas modifiée au cours de l'exécution |
372 |
de cette procédure. Cf. figure (\ref{fig:leapfrog}). |
373 |
\begin{figure}[htbp] |
374 |
\centering |
375 |
\includegraphics{Graphiques/leapfrog} |
376 |
\caption[Principe de l'intégration temporelle selon le schéma |
377 |
``Matsuno-leapfrog'']{Principe de l'intégration temporelle selon le |
378 |
schéma ``Matsuno-leapfrog'', dans \texttt{leapfrog}. Cas où |
379 |
\texttt{iphysiq} = 10. ``FM'' : \Eng{forward} Matsuno, ``BM'' : |
380 |
\Eng{backward} Matsuno. Les flèches courbent représentent des pas |
381 |
``leapfrog''. Les numéros sur les flèches indiquent l'ordre des |
382 |
intégrations. Six intégrations sont nécessaires pour avancer de |
383 |
cinq pas de temps.} |
384 |
\label{fig:leapfrog} |
385 |
\end{figure} |
386 |
|
387 |
\subsection{Maillage horizontal} |
388 |
|
389 |
Trois grilles à trois dimensions. Une grille de base pour, entre |
390 |
autres, la température et l'humidité relative. François L. l'appelle |
391 |
la grille physique. Une grille décalée en longitude par rapport |
392 |
à la grille de base, d'un demi-pas de longitude, pour la vitesse |
393 |
zonale \textit{u}. Une grille décalée en latitude par rapport à la |
394 |
grille de base, d'un demi-pas de latitude, pour la vitesse \textit{v}. |
395 |
François L. appelle ces deux dernières grilles les grilles dynamiques. |
396 |
Les points de la grille de base de dernier indice de longitude |
397 |
et ceux de premier indice de longitude ont la même longitude. |
398 |
De même, les points de la grille de base de dernier indice de |
399 |
latitude et ceux de premier indice de latitude ont la même latitude. |
400 |
On peut dire que la grille de base est « repliée sur elle |
401 |
même » en longitude et latitude. Les deux grilles décalées |
402 |
ont autant de points que la grille de base. Elles sont aussi |
403 |
repliées sur elles-mêmes en longitude et latitude. \textit{Cf}. : |
404 |
\begin{itemize} |
405 |
\item Arakawa et Lamb [1977 \#737] ; |
406 |
\item description du modèle par Phu L. V. (juin 1989) ; |
407 |
\item |
408 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
409 |
IOIPSL LMDZ/Manuel_LMDZE.texfol/Discrétisation équations |
410 |
LMDZ.ps}{Discrétisation des équations de la dynamique dans LMDZ} |
411 |
(16 novembre 1995) ; |
412 |
\item manuel de LMDZ pour Mars (5 avril 2004) ; |
413 |
\item commentaires dans la procédure \verb+inigeom +; |
414 |
\item |
415 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
416 |
IOIPSL LMDZ/Grilles horizontales LMDZ.odg}{schéma des grilles |
417 |
horizontales}. |
418 |
\end{itemize} |
419 |
\verb+iim+ est le nombre de longitudes distinctes dans la |
420 |
grille de base (scalaire). Les valeurs de ces longitudes sont |
421 |
dans le tableau \verb+rlonv+. La taille de rlonv est \verb-iim+1- |
422 |
et \verb-rlonv(iim+1)=rlonv(1)-. \verb+jjm+ est le |
423 |
nombre d'intervalles de latitude dans la grille de base. Les |
424 |
valeurs de ces longitudes sont dans le tableau \verb+rlatu+. |
425 |
La taille de rlatu est \verb-jjm+1-. Pour la partie physique |
426 |
du modèle, les points distincts de la grille scalaire sont classés |
427 |
par ordre de latitude décroissante et, pour une même latitude, |
428 |
par ordre de longitude croissante (cf. figure 2.2 du manuel de |
429 |
LMDZ pour Mars). Chacun des points distincts de la grille scalaire |
430 |
peut donc être identifié par un indice dans ce classement : |
431 |
un seul indice au lieu de 2 indices faisant référence à la latitude |
432 |
et à la longitude. Les tableaux \verb+latfi+ et \verb+lonfi+ |
433 |
contiennent les coordonnées associées à cet indice simple. |
434 |
|
435 |
Pour le remplissage des tableaux \verb+latfi+, \verb+lonfi+, |
436 |
\verb+zcufi+ et \verb+zcvfi+ dans \verb+gcm+, cf. |
437 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
438 |
IOIPSL LMDZ/Grilles horizontales LMDZ.odg}{schéma des grilles |
439 |
horizontales}. Dans cette figure, dans le tableau \verb+mask_v+, le |
440 |
point d'indice \verb-(1, jjm + 1)- est marqué |
441 |
{\textquotedbl}T(2){\textquotedbl} pour rappeler qu'il est recopié |
442 |
dans deux éléments de \verb+zcvfi+. zcufi et zcvfi sont les valeurs de |
443 |
cu et cv affectées aux points de la grille « physique ». Par exemple, |
444 |
regardons le |
445 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
446 |
IOIPSL LMDZ/Grilles horizontales LMDZ.odg}{schéma} et considérons le |
447 |
cas de \verb+cv+ et \verb+zcvfi+. \verb+cv+ est bien défini sur les |
448 |
points rouges. Il s'agit de donner des valeurs de \verb+cv+ aux points |
449 |
(noirs) de la grille physique. Pour tous les points de la grille |
450 |
physique sauf le pôle sud, on choisit la valeur au point rouge juste |
451 |
en dessous. Pour le pôle sud, on prend le point rouge juste au dessus. |
452 |
On pourrait imaginer d'interpoler les valeurs en plusieurs points |
453 |
rouges pour chaque point de la grille physique. Et on pourrait de même |
454 |
imaginer des interpolations pour les valeurs de \verb+zcufi+. Cela |
455 |
aurait plus de sens. Mais Phu L. V. semble dire que la faible |
456 |
importance des tableaux \verb+zcufi+ et \verb+zcvfi+ ne justifie pas |
457 |
de compliquer ainsi. Cf. son message du 27/10/6. |
458 |
|
459 |
Dans la version que j'ai téléchargée, la résolution (modifiable) est |
460 |
de 96 longitudes \ensuremath{\times} 71 latitudes \ensuremath{\times} |
461 |
19 niveaux. François L. préfère 72 latitudes pour avoir une grille |
462 |
sur des valeurs rondes. Laurent Li utilise LMDZ avec la résolution 72 |
463 |
longitudes \ensuremath{\times} 45 latitudes \ensuremath{\times} 19 |
464 |
niveaux. Soit une résolution de 5° en longitude et 4° en latitude. |
465 |
LMDZ ne converge plus lorsque la résolution devient trop fine. |
466 |
\textit{Cf.} aussi la |
467 |
\href{http://www.lmd.jussieu.fr/~lmdz/LMDZ4/choisir_une_resolution.html}{documentation |
468 |
officielle de LMDZ4}. |
469 |
|
470 |
\subsection{Maillage vertical} |
471 |
|
472 |
L'indice de niveau vertical est croissant de la surface au sommet de |
473 |
l'atmosphère. Cf. |
474 |
\href{http://www.lmd.jussieu.fr/~lmdz/manuelGCM/main.ps}{manuel de |
475 |
LMDZ pour Mars}, figure 2.4. |
476 |
|
477 |
D'après François L., le temps passé dans les calculs de transfert |
478 |
radiatif est proportionnel au cube du nombre de niveaux verticaux. |
479 |
|
480 |
Dans \verb+disvert+, $s$ peut être calculé de quatre façons |
481 |
selon la valeur de \verb+s_sampling+. D'après Phu L. V., le cas |
482 |
\verb+param+ n'est plus utilisé (il n'est pas utilisé non plus par |
483 |
François L. pour son extension à la haute atmosphère). Les quatre |
484 |
échantillonnages de $s$ sont comparés pour quelques valeurs |
485 |
de \verb+llm+ dans les figures (\ref{fig:compare_sampl_9}), |
486 |
(\ref{fig:compare_sampl_19}), (\ref{fig:compare_sampl_30}) et |
487 |
(\ref{fig:compare_sampl_41}). |
488 |
\begin{figure}[htbp] |
489 |
\centering |
490 |
\input{Graphiques/compare_sampl_9.tex} |
491 |
\caption[Maillage vertical, échantillonnages de $s$ pour : |
492 |
\texttt{llm} = 9]{Maillage vertical. Comparaison des quatre |
493 |
échantillonnages de $s$ pour : \texttt{llm} = 9 ; \texttt{h} = 7 ; |
494 |
\texttt{deltaz} = \nombre{0,04} ; \texttt{beta} = \nombre{1,3} ; |
495 |
\texttt{k0} = 20 ; \texttt{k1} = \nombre{1,2}. La pression à la |
496 |
surface choisie est de 101325 Pa.} |
497 |
\label{fig:compare_sampl_9} |
498 |
\end{figure} |
499 |
\begin{figure}[htbp] |
500 |
\centering |
501 |
\input{Graphiques/compare_sampl_19.tex} |
502 |
\caption[Maillage vertical, échantillonnages de $s$ pour : |
503 |
\texttt{llm} = 19]{Maillage vertical. Comparaison des quatre |
504 |
échantillonnages de $s$ pour : \texttt{llm} = 19 ; \texttt{h} = 7 |
505 |
; \texttt{deltaz} = \nombre{0,04} ; \texttt{beta} = \nombre{1,3} ; |
506 |
\texttt{k0} = 20 ; \texttt{k1} = \nombre{1,2}. La pression à la |
507 |
surface choisie est de 101325 Pa.} |
508 |
\label{fig:compare_sampl_19} |
509 |
\end{figure} |
510 |
\begin{figure}[htbp] |
511 |
\centering |
512 |
\input{Graphiques/compare_sampl_30.tex} |
513 |
\caption[Maillage vertical, échantillonnages de $s$ pour : |
514 |
\texttt{llm} = 30]{Maillage vertical. Comparaison des quatre |
515 |
échantillonnages de $s$ pour : \texttt{llm} = 30 ; \texttt{h} = 7 |
516 |
; \texttt{deltaz} = \nombre{0,04} ; \texttt{beta} = \nombre{1,3} ; |
517 |
\texttt{k0} = 20 ; \texttt{k1} = \nombre{1,2}. La pression à la |
518 |
surface choisie est de 101325 Pa.} |
519 |
\label{fig:compare_sampl_30} |
520 |
\end{figure} |
521 |
\begin{figure}[htbp] |
522 |
\centering |
523 |
\input{Graphiques/compare_sampl_41.tex} |
524 |
\caption[Maillage vertical, échantillonnages de $s$ pour : |
525 |
\texttt{llm} = 41]{Maillage vertical. Comparaison des quatre |
526 |
échantillonnages de $s$ pour : \texttt{llm} = 41 ; \texttt{h} = 7 |
527 |
; \texttt{deltaz} = \nombre{0,04} ; \texttt{beta} = \nombre{1,3} ; |
528 |
\texttt{k0} = 20 ; \texttt{k1} = \nombre{1,2}. La pression à la |
529 |
surface choisie est de 101325 Pa.} |
530 |
\label{fig:compare_sampl_41} |
531 |
\end{figure} |
532 |
|
533 |
On pose : |
534 |
\begin{align*} |
535 |
& b(0) = 0 \\ |
536 |
& b(s) = \exp(1 - 1 / s^2) |
537 |
\qquad \mathrm{pour}\ s \in \mathbb{R}^* |
538 |
\end{align*} |
539 |
$b$ est $\mathscr{C}^\infty$ sur $\mathbb{R}$, paire, strictement croissante |
540 |
sur $\mathbb{R}^+$. On a : |
541 |
\begin{align*} |
542 |
& \exists \lim_{s \to +\infty} b(s) = e \\ |
543 |
& b'(0) = 0 \\ |
544 |
& b'(1) = 2 |
545 |
\end{align*} |
546 |
et $b$ admet un point d'inflexion en $\sqrt{6} / 3$ ($\approx |
547 |
\nombre{0,8}$). Cf. figure (\ref{fig:disvert}). |
548 |
\begin{figure}[htbp] |
549 |
\centering |
550 |
\input{Graphiques/disvert.tex} |
551 |
\caption[Maillage vertical : $a$ et $b$]{Maillage vertical. Pour |
552 |
$p_a = 5 \cdot 10^4$ Pa et $p_s = \nombre{101325}$ Pa.} |
553 |
\label{fig:disvert} |
554 |
\end{figure} |
555 |
Puis, étant donné $p_a \in \mathbb{R}^{+*}$, on pose, pour tout |
556 |
$s$ : |
557 |
\begin{displaymath} |
558 |
a(s) = p_a (s - b(s)) |
559 |
\end{displaymath} |
560 |
Pour $s > 0$ : |
561 |
\begin{displaymath} |
562 |
a'(s) \ge 0 \Leftrightarrow b(s) \le s^3 / 2 |
563 |
\end{displaymath} |
564 |
L'équation $b(s) = s^3 / 2$ admet deux solutions strictement |
565 |
positives. Notons-les $s_1$ et $s_2$, $s_1 < 1 < |
566 |
s_2$. $a$ est strictement croissante sur $[0, s_1]$, |
567 |
strictement décroissante sur $[s_1, s_2]$, strictement |
568 |
croissante sur $[s_2, + \infty[$. Cf. figure (\ref{fig:disvert}). |
569 |
Enfin, étant donné $p_s \in \mathbb{R}^{+*}$, on pose, pour tout |
570 |
$s$ : |
571 |
\begin{align*} |
572 |
p(s) & = a(s) + b(s) p_s \\ |
573 |
& = p_a s + b(s) (p_s - p_a) |
574 |
\end{align*} |
575 |
Cf. figure (\ref{fig:disvert}). On a : |
576 |
\begin{displaymath} |
577 |
s \le \nombre{0,3} \Rightarrow b(s) \ll s |
578 |
\end{displaymath} |
579 |
Cf. figure (\ref{fig:b_s}). |
580 |
\begin{figure}[htbp] |
581 |
\centering |
582 |
\input{Graphiques/b_s.tex} |
583 |
\caption[Maillage vertical : $b(s) / s$]{Maillage vertical. Les |
584 |
niveaux de $s$ sont des niveaux de pression lorsque $b(s) / s$ est |
585 |
très petit.} |
586 |
\label{fig:b_s} |
587 |
\end{figure} |
588 |
En supposant que $p_s - p_a$ soit de l'ordre de $p_a$, on a donc, aux |
589 |
plus hauts niveaux de l'atmosphère, tels que $s(l) \le |
590 |
\nombre{0,3}$ : |
591 |
\begin{displaymath} |
592 |
p(s(l)) \approx a(s(l)) |
593 |
\approx p_a s(l) |
594 |
\end{displaymath} |
595 |
En haut de l'atmosphère, les niveaux de $s$ sont des niveaux de |
596 |
pression. Ces niveaux de pression sont les valeurs de \verb+ap+. |
597 |
(Après une exécution de \verb+etat0_lim+, \verb+ap+ est dans le |
598 |
fichier \verb+start.nc+.) |
599 |
|
600 |
\subsubsection{Cas \texttt{LMD5}} |
601 |
|
602 |
Pour $l \in \{1, \dots, |
603 |
\mathtt{llm}\}$, posons : |
604 |
\begin{displaymath} |
605 |
x_l = \frac{\pi (l - 1/2)}{\mathtt{llm} + 1} |
606 |
\end{displaymath} |
607 |
Pour $x \in \mathbb{R}$, posons : |
608 |
\begin{displaymath} |
609 |
f(x) = 1 + 7 \sin^2 x |
610 |
\end{displaymath} |
611 |
$f$ est périodique de période $\pi$. Cf. figure (\ref{fig:f_x}). |
612 |
\begin{figure}[htbp] |
613 |
\centering |
614 |
\input{Graphiques/f_x.tex} |
615 |
\caption{Maillage vertical, fonction utilisée dans le cas |
616 |
\texttt{LMD5}.} |
617 |
\label{fig:f_x} |
618 |
\end{figure} |
619 |
$(x_l)$ est un échantillonnage de $[0, \pi]$ d'autant plus fin que |
620 |
\verb+llm+ est grand : |
621 |
\begin{align*} |
622 |
& \forall l \in \{1, \dots, \mathtt{llm}\}, x_l \in ]0, \pi[ \\ |
623 |
& \exists \lim_{\mathtt{llm} \to +\infty} x_1 = 0 \\ |
624 |
& \exists \lim_{\mathtt{llm} \to +\infty} x_\mathtt{llm} = \pi \\ |
625 |
& x_{l+1} - x_l = \frac{\pi}{\mathtt{llm} + 1} |
626 |
\underset{\mathtt{llm} \to + \infty}{\longrightarrow} 0 |
627 |
\end{align*} |
628 |
Le programme calcule : |
629 |
\begin{displaymath} |
630 |
\mathtt{ds(l)} = \frac{f(x_l)}{\sum_{l'=1}^\mathtt{llm} f(x_{l'})} |
631 |
\end{displaymath} |
632 |
Or : |
633 |
\begin{displaymath} |
634 |
\sum_{l=1}^\mathtt{llm} f(x_{l}) \frac{\pi}{\mathtt{llm} + 1} |
635 |
\underset{\mathtt{llm} \to + \infty}{\longrightarrow} \int_0 ^\pi f |
636 |
\end{displaymath} |
637 |
Notons : |
638 |
\begin{displaymath} |
639 |
I := \int_0 ^\pi f = 9 \pi / 2 |
640 |
\end{displaymath} |
641 |
On a donc : |
642 |
\begin{displaymath} |
643 |
\sum_{l=1}^\mathtt{llm} f(x_{l}) |
644 |
\underset{\mathtt{llm} \to +\infty}{\sim} \frac{I}{\pi} \mathtt{llm} |
645 |
\end{displaymath} |
646 |
(L'erreur est de 2 \% environ pour \verb+llm+ = 41.) Le niveau |
647 |
$s$ le plus petit est 0, le deuxième plus petit est : |
648 |
\begin{displaymath} |
649 |
s(\mathtt{llm}) = \mathtt{ds(llm)} |
650 |
\underset{\mathtt{llm} \to +\infty}{\sim} \frac{\pi f(\pi)}{I\ \mathtt{llm}} |
651 |
\end{displaymath} |
652 |
|
653 |
Donc : |
654 |
\begin{displaymath} |
655 |
p(s(\mathtt{llm})) |
656 |
\underset{\mathtt{llm} \to +\infty}{\sim} |
657 |
\frac{\pi f(\pi) p_a}{I\ \mathtt{llm}} |
658 |
= \frac{2 p_a}{9\ \mathtt{llm}} |
659 |
\end{displaymath} |
660 |
On a : |
661 |
\begin{align*} |
662 |
& \mathtt{llm} = 9 \Rightarrow \mathtt{ap(llm)} |
663 |
\approx 3 \cdot 10\ \mathrm{hPa} \\ |
664 |
& \mathtt{llm} = 19 \Rightarrow \mathtt{ap(llm)} \approx 8\ \mathrm{hPa} \\ |
665 |
& \mathtt{llm} = 50 \Rightarrow \mathtt{ap(llm)} \approx 2\ \mathrm{hPa} |
666 |
\end{align*} |
667 |
Il faudrait prendre \verb+llm+ de l'ordre de $10^3$ pour résoudre la |
668 |
région décrite par la paramétrisation de la chimie de l'ozone. |
669 |
|
670 |
\subsubsection{Cas \texttt{strato1}} |
671 |
|
672 |
On pose, pour $(x, y) \in |
673 |
\mathbb{R}^2$ : |
674 |
\begin{displaymath} |
675 |
f(x, y) = \tgh x + \frac{1}{2} \tgh \frac{x - y}{2} |
676 |
\end{displaymath} |
677 |
Relevons quelques propriétés de $f$, considérée comme une fonction de |
678 |
$x$, à $y$ fixé. $f$ est une fonction strictement croissante sur |
679 |
$\mathbb{R}$, de $- 3/2$ à $3/2$. Par rapport à la fonction $\tgh$, le |
680 |
second terme de $f$ est 2 fois plus large, translaté de $y$ et 2 fois |
681 |
moins haut. Puisque $\tgh 3$ est déjà très proche de 1, les domaines |
682 |
de variation des deux termes de $f$ sont clairement séparés dès que $y |
683 |
\ge 9$. Cf. figure (\ref{fig:strato1}). |
684 |
\begin{figure}[htbp] |
685 |
\centering |
686 |
\input{Graphiques/strato1.tex} |
687 |
\caption[Maillage vertical, cas \texttt{strato1}]{Maillage vertical, |
688 |
cas \texttt{strato1}. Fonction utilisée pour \texttt{llm} = 70.} |
689 |
\label{fig:strato1} |
690 |
\end{figure} |
691 |
|
692 |
Le programme calcule : |
693 |
\begin{displaymath} |
694 |
\mathtt{dz(l)} |
695 |
= \nombre{1,56} |
696 |
+ f\left(\frac{\mathtt{l} - 12}{5}, \frac{\mathtt{llm} - 12}{5}\right) |
697 |
\end{displaymath} |
698 |
D'après l'étude de $f$, $\mathtt{dz(l)} \in ]\nombre{0,06} ; |
699 |
\nombre{3,06}[$, \verb+dz(l)+ croît avec $l$. Pour \verb+llm+ |
700 |
suffisamment grand : |
701 |
\begin{align*} |
702 |
& \mathtt{dz(1)} \approx \nombre{1,06} - \tgh\nombre{2,2} |
703 |
\approx \nombre{8,43} \cdot 10^{-2} \\ |
704 |
& \mathtt{dz(llm)} \approx \nombre{2,56} |
705 |
\end{align*} |
706 |
(l'erreur relative sur \verb+dz(1)+ devient inférieure à 1 \% si |
707 |
\verb+llm+ $\ge 37$ ; l'erreur relative sur \verb+dz(llm)+ devient |
708 |
inférieure à 1 \% si \verb+llm+ $\ge 23$). Les domaines de variation |
709 |
des deux termes $\tgh$ de \verb+dz+ sont clairement séparés dès que |
710 |
\verb+llm+ $\ge 57$. |
711 |
|
712 |
Posons, pour $z \in \mathbb{R}$ : |
713 |
\begin{displaymath} |
714 |
s(z) |
715 |
= \frac{\exp(- z / h) - \exp(- \mathtt{zz(llm + 1)} / h)} |
716 |
{1 - \exp(- \mathtt{zz(llm + 1)} / h)} |
717 |
\end{displaymath} |
718 |
Le programme calcule $s(\mathtt{zz})$. On a, pour $h = 7$ km : |
719 |
\begin{align*} |
720 |
& \mathtt{llm} = 40 |
721 |
\Rightarrow \mathtt{ap(llm)} \approx \nombre{2,70}\ \mathrm{Pa} \\ |
722 |
& \mathtt{llm} = 41 |
723 |
\Rightarrow \mathtt{ap(llm)} \approx \nombre{2,01}\ \mathrm{Pa} \\ |
724 |
& \mathtt{llm} = 50 |
725 |
\Rightarrow \mathtt{ap(llm)} \approx \nombre{0,14}\ \mathrm{Pa} |
726 |
\end{align*} |
727 |
|
728 |
Si $\exp[- \mathtt{zz(llm + 1)} / h] \ll 1$, pour des valeurs de $z$ |
729 |
telles que : |
730 |
\begin{equation} |
731 |
\label{eq:z_interm} |
732 |
\left\{ |
733 |
\begin{array}{l} |
734 |
\exp(- z / h) \gg \exp[- \mathtt{zz(llm + 1)} / h] \\ |
735 |
\exp(- z / h) \le \nombre{0,3} |
736 |
\end{array} |
737 |
\right. |
738 |
\end{equation} |
739 |
on a : |
740 |
\begin{align*} |
741 |
& s(z) \approx \exp(- z / h) \\ |
742 |
& b[s(z)] \ll s(z) |
743 |
\end{align*} |
744 |
Donc, en supposant que $p_s - p_a$ soit de l'ordre de $p_a$ : |
745 |
\begin{displaymath} |
746 |
p[s(z)] \approx a[s(z)] \approx p_a \exp(- z / h) |
747 |
\end{displaymath} |
748 |
Si $h$ est proche de la hauteur d'échelle réelle dans la région |
749 |
considérée alors les variations d'altitude entre deux niveaux |
750 |
$s$ sont proches de \verb+dz+. Par exemple, pour \verb+llm+ = |
751 |
70, \verb=zz(llm + 1)= $\approx 125$ km. Si $h = 7$ km, la condition |
752 |
(\ref{eq:z_interm}) devient $z \in [9\ \mathrm{km}, 92\ \mathrm{km}]$, |
753 |
qui correspond aux niveaux 17 à 56. |
754 |
|
755 |
\subsubsection{Cas \texttt{strato2}} |
756 |
|
757 |
C'est le cas conseillé par François L. pour un domaine incluant la |
758 |
stratosphère. Ce cas est analogue au cas \verb+LMD5+. Seule la |
759 |
fonction $f$ utilisée change. Pour $x \in \mathbb{R}$, posons : |
760 |
\begin{displaymath} |
761 |
g(x) |
762 |
= \frac{1}{2} \left[1 - \tgh\left(\frac{x - \pi / 2}{\pi / 2}\right)\right] |
763 |
\end{displaymath} |
764 |
Par rapport à la fonction $\tgh$, $g$ est translatée de $\pi / 2$ vers |
765 |
la droite, dilatée horizontalement d'un facteur $\pi / 2$, inversée |
766 |
verticalement, translatée verticalement de 1 et contractée |
767 |
verticalement d'un facteur 2. Elle décroît donc de 1 à $-1$. La |
768 |
fonction $f$ utilisée est maintenant : |
769 |
\begin{displaymath} |
770 |
f(x) = (1 + 7 \sin^2 x) g(x)^2 |
771 |
\end{displaymath} |
772 |
\begin{figure}[htbp] |
773 |
\centering |
774 |
\input{Graphiques/strato2.tex} |
775 |
\caption[Maillage vertical, cas \texttt{strato2}]{Maillage vertical, |
776 |
cas \texttt{strato2}. $f(\pi) \approx \nombre{0,014}$.} |
777 |
\label{fig:strato2} |
778 |
\end{figure} |
779 |
On a : |
780 |
\begin{displaymath} |
781 |
I := \int_0 ^\pi f \approx \nombre{4,02} |
782 |
\end{displaymath} |
783 |
donc : |
784 |
\begin{displaymath} |
785 |
p(s(\mathtt{llm})) |
786 |
\underset{\mathtt{llm} \to +\infty}{\sim} |
787 |
\frac{\pi f(\pi) p_a}{I\ \mathtt{llm}} |
788 |
= \frac{\nombre{0,011} p_a}{\mathtt{llm}} |
789 |
\end{displaymath} |
790 |
|
791 |
\subsection{Module \texttt{inter\_barxy\_m}} |
792 |
|
793 |
Cf. |
794 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Documentation |
795 |
IOIPSL LMDZ/Mouvement de données/inter_barxy.odg}{graphique de |
796 |
transmission de données}. |
797 |
|
798 |
Le seul cas où \verb+size(champint, 2)+ vaut \verb+jjm+ dans |
799 |
\verb+inter_barxy+ se produit suite à l'appel : |
800 |
\begin{verbatim} |
801 |
vvent = start_inter_3d("V", ...) |
802 |
\end{verbatim} |
803 |
C'est aussi le seul cas où \verb+rlatimod+ reçoit \verb+rlatu(:jjm)+. |
804 |
Dans tous les autres cas, \verb+rlatimod+ dans \verb+inter_barxy+ |
805 |
reçoit \verb+rlatv+. |
806 |
|
807 |
Dans le cas le plus fréquent, dans \verb+inter_barxy+ : |
808 |
\begin{verbatim} |
809 |
size(champint, 2) == jjm + 1 |
810 |
\end{verbatim} |
811 |
et \verb+rlatimod+ reçoit \verb+rlatv+. On veut interpoler sur les |
812 |
\verb+jjm+ + 1 latitudes \verb+rlatu+ mais on reçoit \verb+rlatv+. |
813 |
Bizarre. |
814 |
|
815 |
Je crois que \verb+inter_bary+ réalise une moyenne et non une |
816 |
interpolation. (Par exemple, pour l'initialisation de l'abondance de |
817 |
vapeur d'eau, on veut calculer \verb+q3d+ aux points \verb+rlatu+.) |
818 |
Cf. |
819 |
\href{file:///user/guez/Documents/Utilisation_LMDZ/Utilisation_LMDZ.texfol/Utilisation-LMDZ.pdf}{analyse}. |
820 |
Je crois que le calcul suppose une fonction en escalier, dont les |
821 |
marches sont limitées par \verb+yjdat+. Pour $i \ge 2$, \verb+fdat(i)+ |
822 |
est interprété comme la hauteur de la marche $[\mathtt{yjdat(i-1)}, |
823 |
\mathtt{yjdat(i)}]$. Pour $i \ge 2$, \verb+inter_bary(i)+ est la |
824 |
moyenne sur l'intervalle $[\mathtt{yjmod(i-1)}, \mathtt{yjmod(i)}]$. |
825 |
|
826 |
Cf. le remaillage \verb+@ave+ dans Ferret. On pourrait imaginer une |
827 |
procédure qui reçoit les tableaux \verb+x1edge+, \verb+x2edge+ et |
828 |
\verb+y1+, qui retourne un tableau \verb+y2+, où \verb+x1edge+ et |
829 |
\verb+x2edge+ sont deux listes de frontières d'intervalles, \verb+y1+ |
830 |
et \verb+y2+ ont un élément de moins que \verb+x1edge+ et |
831 |
\verb+x2edge+, \verb+y1+ donne les valeurs d'une fonction en escalier |
832 |
dont les marches sont limitées par \verb+x1edge+, et \verb+y2+ est la |
833 |
moyenne de cette fonction sur les intervalles limités par |
834 |
\verb+x2edge+. |
835 |
|
836 |
\subsection{Les traceurs en général} |
837 |
|
838 |
Procédures qui concernent les traceurs dans le programme |
839 |
\verb+etat0_lim+ : \verb+iniadvtrac+, \verb+etat0+. Procédures qui |
840 |
concernent les traceurs dans le programme \verb+gcm+ : |
841 |
\verb+iniadvtrac+, \verb+dynetat0+, \verb+caladvtrac+, \verb+advtrac+, |
842 |
\verb+physiq+, \verb+phytrac+. Laurent L. connaît peut-être bien le |
843 |
traitement des traceurs. Il a écrit une grande partie de la |
844 |
``physique'' du modèle. |
845 |
|
846 |
Dans le fichier \verb+traceur.def+, le premier nombre est le nombre de |
847 |
traceurs. Chaque ligne suivante correspond à un traceur. Pour chaque |
848 |
traceur, on indique le numéro du schéma |
849 |
d'advection à utiliser. Les numéros sont définis dans |
850 |
\verb+iniadvtrac+. |
851 |
|
852 |
\verb+traceur.def+ est lu uniquement dans le sous-programme |
853 |
\verb+iniadvtrac+. \verb+nq+ est lu dans \verb+traceur.def+. |
854 |
\verb+iniadvtrac+ impose \verb+nqmx+ $\ge$ \verb+nq+. \verb+nq+ est |
855 |
passé à \verb+calfis+ et à \verb+physiq+. \verb+physiq+ impose |
856 |
\verb+nq+ $\ge 2$. Donc \verb+nqmx+ doit être supérieur à 2. |
857 |
\verb+nbtr+ (variable de \verb+dimphy+) vaut 1 pour \verb+nqmx+ |
858 |
\verb+=+ \verb+2+ et \verb+nqmx+ -- \verb+2+ pour \verb+nqmx+ $\ge 3$. |
859 |
\verb+nq-2+ est passé à \verb+phytrac+ et devient \verb+nqmax+ dans |
860 |
\verb+phytrac+. À cause de l'appel de \verb+dynredem1+ par |
861 |
\verb+etat0+, \verb+nq+ doit être égal à \verb+nqmx+. \verb+initrrnpb+ |
862 |
impose \verb+nbtr+ $\ge 2$, mais il n'est appelé que si \verb+rnpb+ |
863 |
est vrai. Cf. aussi message de Ionela M. du 29/11/6. |
864 |
|
865 |
\verb+phytrac+ peut effectuer un traitement spécial pour les traceurs |
866 |
plomb et radon (décroissance radioactive etc.). \verb+phytrac+ |
867 |
suppose que ce sont les traceurs numérotés 3 et 4. Donc si on veut |
868 |
laisser les traceurs plomb et radon, avec le traitement spécial qui |
869 |
leur est réservé, alors il faut les laisser en troisième et quatrième |
870 |
position dans \verb+traceur.def+ et laisser \verb+rnpb+ à |
871 |
\verb+.true.+ dans \verb+physiq+. Cette variable est définie dans |
872 |
\verb+physiq+ et \verb+phytrac+ uniquement. Cf. le message de |
873 |
Frédéric H. du 26/1/7. |
874 |
|
875 |
Dans le programme \verb+etat0_lim+, la variable qui contient la |
876 |
distribution initiale des traceurs est \verb+q3d+. Cette variable est |
877 |
définie dans \verb+etat0+. Pour la vapeur d'eau (traceur numéro 1), la |
878 |
distribution initiale est calculée à partir de la variable \verb+R+ du |
879 |
fichier \verb+ECDYN.nc+. La distribution initiale des traceurs est |
880 |
écrite dans \verb+start.nc+. |
881 |
|
882 |
Dans le programme \verb+gcm+, la procédure \verb+iniadvtrac+ lit les |
883 |
noms des traceurs. La procédure \verb+dynetat0+ lit les distributions |
884 |
initiales à partir du fichier \verb+start.nc+ vers la variable |
885 |
\verb+q+. Cette variable se retrouve intacte dans \verb+leapfrog+. Par |
886 |
ailleurs, la concentration de radon dans le sol est lue dans |
887 |
\verb+starttrac+ par \verb+phytrac+. Dans \verb+phytrac+, on calcule |
888 |
la fraction massique à $t + \ud t$ : |
889 |
\begin{displaymath} |
890 |
q(\vec r, t + \ud t) = q(\vec r, t) + (\partial_t q)(\vec r, t)\ud t |
891 |
\end{displaymath} |
892 |
Cf. |
893 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Guide_IPSL_climate_models/Traceurs.odg}{graphique |
894 |
de flux de données}. Les distributions de traceurs calculées sont |
895 |
écrites dans \verb+dyn_hist_ave.nc+ et \verb+histrac.nc+. Les |
896 |
distributions finales sont écrites dans \verb+restart.nc+ par |
897 |
\verb+dynredem1+. Par ailleurs, la concentration de radon dans le sol |
898 |
est écrite dans \verb+restarttrac+ par \verb+phytrac+. |
899 |
|
900 |
\subsection{L'ozone en tant que traceur} |
901 |
|
902 |
Résumé. On part d'une certaine paramétrisation de la chimie |
903 |
stratosphérique de l'ozone, préparée par Daniel Cariolle et Hubert |
904 |
Teyssedre. La paramétrisation calcule en chaque point une production |
905 |
chimique d'ozone. On suit l'ozone ainsi produit comme un traceur |
906 |
passif, c'est-à-dire qu'il n'a aucun effet sur le reste du modèle. |
907 |
Les paramètres chimiques sont remaillés pour LMDZ, en latitude et en |
908 |
temps, dans \verb+etat0_lim+. Le remaillage vertical est fait une fois |
909 |
par jour dans \verb+gcm+. On néglige donc l'effet sur le remaillage de |
910 |
la variation temporelle de pression à la surface pendant une journée. |
911 |
Le champ initial d'ozone spécifié dans \verb+etat0_lim+ est le champ |
912 |
de référence de D. Cariolle (un des paramètres). Pour les passages de |
913 |
fraction molaire d'ozone à fraction massique, on néglige la présence |
914 |
de vapeur d'eau dans l'air. L'évolution de la fraction massique |
915 |
d'ozone pendant un pas de temps de la physique est calculée dans |
916 |
\verb+phytrac+ : les contributions du transport et de la chimie sont |
917 |
intégrées séparément. La température, la masse volumique de l'air et |
918 |
la direction du Soleil sont considérées comme constantes pendant un |
919 |
pas de temps de la physique. Un des paramètres de D. Cariolle est la |
920 |
densité-colonne d'ozone au dessus d'un point donné. Cette |
921 |
densité-colonne est calculée dans l'approximation ``plans |
922 |
parallèles''. La production chimique est intégrée pendant un pas de |
923 |
temps de la physique par la méthode du point milieu. Une fois par jour |
924 |
de simulation, dans \verb+gcm+, une procédure lit dans un fichier tous |
925 |
les paramètres chimiques au bon jour, remaillés en latitude pour LMDZ. |
926 |
|
927 |
Cf. Cariolle (1986 \#761), |
928 |
\href{file:/user/guez/Documents/Commentaires_lectures/Autres_publications/-1991/Cariolle_1986_761.texfol/Cariolle_1986_761.dvi}{commentaires |
929 |
sur Cariolle (1986 \# 761)}, texte des courriels (François L., |
930 |
29/11/6 ; D. Cariolle, 29/11/6 ; François L., 12/1/7), fichiers |
931 |
\verb+O3_readme.txt+ de D. Cariolle et \verb+coefoz_v2_3.nc+. |
932 |
Attention : les symboles ne désignent apparemment pas les mêmes |
933 |
grandeurs dans Cariolle (1986, \#761) et dans la paramétrisation de |
934 |
2005. $r$ passe d'une concentration massique à une fraction molaire. |
935 |
$P$ et $L$ passent de la dimension masse$^{-1}$ temps$^{-1}$ à la |
936 |
dimension temps$^{-1}$. |
937 |
|
938 |
Phu L. V. conseille d'utiliser le schéma de transport numéro 10 (choix |
939 |
dans \verb+traceur.def+). |
940 |
|
941 |
Nous notons avec un indice ``Mob'' les valeurs des paramètres tirés de |
942 |
Mobidic. |
943 |
|
944 |
\subsubsection{La paramétrisation de la chimie} |
945 |
|
946 |
Le taux de production photochimique par unité de volume (ou de masse |
947 |
d'air) est le nombre de molécules créées par unité de volume (ou de |
948 |
masse d'air) et par unité de temps. On peut aussi définir un taux de |
949 |
destruction chimique d'une espèce, de dimension temps$^{-1}$, sachant |
950 |
que le nombre de molécules de cette espèce détruites par unité de |
951 |
volume et par unité de temps est proportionnel à la densité de |
952 |
l'espèce. Comment $P$ est-il donc défini dans la paramétrisation de |
953 |
2005 ? En divisant par $N$, la densité totale de la phase gazeuse ? |
954 |
Quelle définition pour $L$ ? Par ailleurs, quelle grandeur la |
955 |
paramétrisation de 2005 permet-elle de calculer ? Certainement pas la |
956 |
dérivée partielle comme noté dans l'équation (5) de Cariolle (1986 |
957 |
\#761), car la dérivée partielle doit dépendre du mouvement de |
958 |
l'atmosphère. Cf. |
959 |
\href{file:/user/guez/Documents/Commentaires_lectures/Autres_publications/-1991/Cariolle_1986_761.texfol/Cariolle_1986_761.dvi}{commentaires |
960 |
sur Cariolle (1986 \# 761)}. La paramétrisation donne-t-elle la |
961 |
dérivée lagrangienne de la fraction molaire, comme noté dans le |
962 |
fichier \verb+O3_readme.txt+ ? Notons : |
963 |
\begin{description} |
964 |
\item[$r$] : fraction molaire d'ozone |
965 |
\item[$q$] : fraction massique d'ozone |
966 |
\item[$s_{O_3}$] : masse nette d'ozone créée par la chimie par unité de |
967 |
volume et par unité de temps |
968 |
\item[$S$] : nombre net de molécules d'ozone créées par la chimie par |
969 |
unité de volume et par unité de temps |
970 |
\item[$N$] : densité totale de la phase gazeuse |
971 |
\item[$\rho$] : masse volumique totale de la phase gazeuse |
972 |
\end{description} |
973 |
On a : |
974 |
\begin{equation} |
975 |
\label{eq:o3} |
976 |
\left\{ |
977 |
\begin{array}{l} |
978 |
\partial_t \rho_{\mathrm{O}_3} + \diverg(\rho_{\mathrm{O}_3} \vec v) |
979 |
= s_{O_3} \\ |
980 |
\partial_t \rho + \diverg(\rho \vec v) = 0 |
981 |
\end{array} |
982 |
\right. |
983 |
\end{equation} |
984 |
Remarquer que l'équation bilan sur le nombre total de molécules en |
985 |
phase gazeuse est plus compliquée que celle sur la masse totale : le |
986 |
second membre n'est a priori pas nul. D'après (\ref{eq:o3}) : |
987 |
\begin{equation} |
988 |
\label{eq:dqdt} |
989 |
\frac{\ud q}{\ud t} = s_{O_3} / \rho |
990 |
\end{equation} |
991 |
D'où : |
992 |
\begin{displaymath} |
993 |
\frac{\ud r}{\ud t} |
994 |
= |
995 |
\frac{1}{\mu_\mathrm{moy}} \frac{\ud \mu_\mathrm{moy}}{\ud t} r + S / N |
996 |
\end{displaymath} |
997 |
Donc il serait bizarre que la paramétrisation donne $\ud r/\ud t$. Je |
998 |
vais supposer qu'elle donne $S/N$. $P$ et $L$ sont alors les taux de |
999 |
production et de destruction chimique par unité de volume divisés par |
1000 |
$N$. Je note : $P_\mathrm{net} := S / N$, taux de production chimique |
1001 |
d'ozone, par molécule d'air. |
1002 |
|
1003 |
Désignons par $R_\mathrm{het}$ le taux de production chimique |
1004 |
hétérogène par molécule d'ozone, interpolé journalièrement, remaillé |
1005 |
horizontalement et verticalement pour LMDZ. Notons $\phi$ la latitude |
1006 |
et $\lambda$ la longitude. Notons $d(t)$ le jour contenant une date |
1007 |
$t$, $s$ la coordonnée verticale hybride et $\alpha$ l'angle entre la |
1008 |
direction du Soleil et la verticale (\Eng{solar zenith angle?}). Nous |
1009 |
pouvons définir un taux de destruction par chimie hétérogène corrigé : |
1010 |
\begin{multline*} |
1011 |
R_\mathrm{het,corr}(\lambda, \phi, s, t) |
1012 |
= |
1013 |
R_\mathrm{het}(\lambda, \phi, s, d(t)) |
1014 |
1_{T(\lambda, \phi, s, t) > 195\ \mathrm{K}} |
1015 |
1_{\phi \notin [- 45°, 45°]} \\ |
1016 |
\times 1_{\alpha(\lambda, \phi, t) < 87°} |
1017 |
\left(\frac{\mathrm{Clx}}{\nombre{3,8} \cdot 10^{-9}}\right)^2 |
1018 |
\end{multline*} |
1019 |
(La définition de Clx n'est pas claire pour moi.) Clx est supposé être |
1020 |
une constante. On calcule alors : |
1021 |
\begin{multline*} |
1022 |
P_\mathrm{net}(\lambda, \phi, s, t) |
1023 |
= |
1024 |
P_\mathrm{net,Mob}(\lambda, \phi, s, d(t)) \\ |
1025 |
+ \frac{\partial P_\mathrm{net}}{\partial r}(\lambda, \phi, s, d(t)) |
1026 |
[r(\lambda, \phi, s, t) - r_\mathrm{Mob}(\lambda, \phi, s, d(t))] \\ |
1027 |
+ \frac{\partial P_\mathrm{net}}{\partial T}(\lambda, \phi, s, d(t)) |
1028 |
[T(\lambda, \phi, s, t) - T_\mathrm{Mob}(\lambda, \phi, s, d(t))] \\ |
1029 |
+ \frac{\partial P_\mathrm{net}}{\partial \Sigma} |
1030 |
(\lambda, \phi, s, d(t)) |
1031 |
[\Sigma(\lambda, \phi, s, t) - \Sigma_\mathrm{Mob} |
1032 |
(\lambda, \phi, s, d(t))] \\ |
1033 |
+ R_\mathrm{het,corr}(\lambda, \phi, s, t) r(\lambda, \phi, s, t) |
1034 |
\end{multline*} |
1035 |
Ou encore : |
1036 |
\begin{displaymath} |
1037 |
S = P_\mathrm{net,Mob} N + \ldots |
1038 |
+ R_\mathrm{het,corr} n_{\mathrm{O}_3} |
1039 |
\end{displaymath} |
1040 |
Donc $R_\mathrm{het,corr}$ est le taux de production (algébrique) par |
1041 |
chimie hétérogène par unité de volume et par molécule d'ozone. |
1042 |
|
1043 |
Soit une position dans l'atmosphère. Soit $\pseudov{\ud^2 S} = \ud^2 S |
1044 |
\vec u_r$ une surface élémentaire à cette position. Notons $\ud^2 |
1045 |
\tau$ le volume compris dans l'angle solide de sommet le centre de la |
1046 |
Terre, limitant $\ud^2 S$, et aux rayons $r'$ supérieurs au rayon $r$ |
1047 |
de la position considérée. J'interprète $\Sigma$ et |
1048 |
$\Sigma_\mathrm{Mob}$ comme les nombres de molécules d'ozone dans |
1049 |
$\ud^2 \tau$, divisés par $\ud^2 S$. D'où : |
1050 |
\begin{displaymath} |
1051 |
\mu_{\mathrm{O}_3} \Sigma |
1052 |
= \frac{1}{r^2} \int_r ^{+ \infty} q \rho r'^2 \ud r' |
1053 |
\end{displaymath} |
1054 |
Avec l'équilibre hydrostatique vertical : |
1055 |
\begin{displaymath} |
1056 |
\mu_{\mathrm{O}_3} \Sigma |
1057 |
= \frac{1}{r^2} \int_0 ^p q r'^2 \frac{\ud p'}{g} |
1058 |
\end{displaymath} |
1059 |
|
1060 |
$a_2$, $a_4$ et $a_6$ sont constants sur les quatre niveaux de |
1061 |
pression les plus proches de la surface, sur toutes les latitudes et |
1062 |
tous les mois : $a_2 \approx - 2 \cdot 10^{-6}$, $a_4 \equiv 0$, $a_6 |
1063 |
\equiv 0$. L'intervalle de pression correspond à une couche de 400 ou |
1064 |
500 m au dessus de la surface. Pour $a_2$, il y a une discontinuité |
1065 |
entre la valeur dans cette couche proche de la surface et les valeurs |
1066 |
au-dessus. Cf. |
1067 |
\href{file:///user/guez/Documents/Utilisation_LMDZ/Input_etat0_lim/Ozone/a2.ps}{exemple |
1068 |
de variation de $a_2$}. |
1069 |
|
1070 |
Ordres de grandeur dans \verb+coefoz_v2_3.nc+ : |
1071 |
\begin{displaymath} |
1072 |
\begin{array}{rcl} |
1073 |
- 10^{-11} \lesssim & P_\mathrm{net,Mob} & \lesssim 10^{-12} \\ |
1074 |
- 10^{-5} \lesssim & a_2 & \lesssim 10^{-11} \\ |
1075 |
10^{-8} \lesssim & r_\mathrm{Mob} & \lesssim 10^{-5} \\ |
1076 |
- 10^{-12} \lesssim & a_4 & \lesssim 10^{-13} \\ |
1077 |
184 \le & T_\mathrm{Mob} & \le 320 \\ |
1078 |
- 10^{-27} \lesssim & a_6 & \lesssim 10^{-28} \\ |
1079 |
10^{12} \lesssim & \Sigma_\mathrm{Mob} & \lesssim 10^{19} \\ |
1080 |
10^{-13} \lesssim & - R_\mathrm{het} & \lesssim 10^{-5} |
1081 |
\end{array} |
1082 |
\end{displaymath} |
1083 |
Deux valeurs de $|a_6|$ sont non nulles inférieures à $1,2 \cdot |
1084 |
10^{-38}$ : environ $7 \cdot 10^{-39}$ et $2 \cdot 10^{-39}$. $1,2 |
1085 |
\cdot 10^{-38}$ est le plus petit nombre strictement positif normalisé |
1086 |
dans l'ensemble IEEE des réels simple précision. Les deux petites |
1087 |
valeurs de $|a_6|$ sont donc dénormalisées dans le fichier NetCDF. La |
1088 |
perte de précision n'est pas importante. Ce problème ne se produit pas |
1089 |
pour les autres coefficients. |
1090 |
|
1091 |
\subsubsection{Passage de la fraction molaire à la fraction massique d'ozone} |
1092 |
|
1093 |
Nous avons besoin de la masse moléculaire moyenne de l'air. Nous |
1094 |
pourrions considérer l'air comme un mélange de trois espèces : une |
1095 |
espèce fictive ``air sec'', l'eau (vapeur) et l'ozone. Nous supposons |
1096 |
connue la fraction massique d'eau (elle est calculée par |
1097 |
\verb+etat0+). Notons ici $q_i$ la fraction massique de l'espèce $i$ |
1098 |
et $x_i$ sa fraction molaire. Nous avons donc le système de quatre |
1099 |
équations : |
1100 |
\begin{displaymath} |
1101 |
\left\{ |
1102 |
\begin{array}{ll} |
1103 |
q_{\mathrm{H}_2\mathrm{O}} |
1104 |
& = \frac{\mu_{\mathrm{H}_2\mathrm{O}}}{\mu_\mathrm{moy}} |
1105 |
x_{\mathrm{H}_2\mathrm{O}} \\ |
1106 |
\mu_\mathrm{moy} |
1107 |
& = \mu_{\mathrm{H}_2\mathrm{O}} x_{\mathrm{H}_2\mathrm{O}} |
1108 |
+ \mu_{\mathrm{O}_3} x_{\mathrm{O}_3} |
1109 |
+ \mu_\textrm{air sec} x_\textrm{air sec} \\ |
1110 |
1 & = x_{\mathrm{H}_2\mathrm{O}} + x_{\mathrm{O}_3} |
1111 |
+ x_\textrm{air sec} \\ |
1112 |
q_{\mathrm{O}_3} & = \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1113 |
x_{\mathrm{O}_3} |
1114 |
\end{array} |
1115 |
\right. |
1116 |
\end{displaymath} |
1117 |
pour les quatre inconnues : $\mu_\mathrm{moy}$, |
1118 |
$x_{\mathrm{H}_2\mathrm{O}}$, $x_\textrm{air sec}$, $q_{\mathrm{O}_3}$. La |
1119 |
solution de ce système est : |
1120 |
\begin{equation} |
1121 |
\label{eq:qO3} |
1122 |
\left\{ |
1123 |
\begin{array}{ll} |
1124 |
\mu_\mathrm{moy} |
1125 |
& = |
1126 |
\frac{[(\mu_{\mathrm{O}_3} - \mu_\textrm{air sec}) x_{\mathrm{O}_3} |
1127 |
+ \mu_\textrm{air sec}] \mu_{\mathrm{H}_2\mathrm{O}}} |
1128 |
{\mu_{\mathrm{H}_2\mathrm{O}} |
1129 |
+ (\mu_\textrm{air sec} - \mu_{\mathrm{H}_2\mathrm{O}}) |
1130 |
q_{\mathrm{H}_2\mathrm{O}}} \\ |
1131 |
q_{\mathrm{O}_3} |
1132 |
& = |
1133 |
\frac{\mu_{\mathrm{O}_3} |
1134 |
[(\mu_\textrm{air sec} - \mu_{\mathrm{H}_2\mathrm{O}}) |
1135 |
q_{\mathrm{H}_2\mathrm{O}} |
1136 |
+ \mu_{\mathrm{H}_2\mathrm{O}}]} |
1137 |
{\mu_{\mathrm{H}_2\mathrm{O}} |
1138 |
\left[ |
1139 |
\mu_{\mathrm{O}_3} |
1140 |
+ \mu_\textrm{air sec} |
1141 |
\left( |
1142 |
\frac{1}{x_{\mathrm{O}_3}} - 1 |
1143 |
\right) |
1144 |
\right]} |
1145 |
\end{array} |
1146 |
\right. |
1147 |
\end{equation} |
1148 |
Les valeurs typiques semblent être : |
1149 |
\begin{align*} |
1150 |
\max q_{\mathrm{H}_2\mathrm{O}} \# 10^{-2} \\ |
1151 |
\max x_{\mathrm{O}_3} \# 10^{-5} |
1152 |
\end{align*} |
1153 |
L'écart entre $\mu_\mathrm{moy}$ et $\mu_\textrm{air sec}$ est alors |
1154 |
$\lesssim 10^{-2}$, et essentiellement dû à la présence de l'eau. Une |
1155 |
approximation intermédiaire entre : |
1156 |
\begin{displaymath} |
1157 |
\left\{ |
1158 |
\begin{array}{ll} |
1159 |
\mu_\mathrm{moy} & = \mu_\textrm{air sec} \\ |
1160 |
q_{\mathrm{O}_3} & = \frac{\mu_{\mathrm{O}_3}}{\mu_\textrm{air sec}} |
1161 |
x_{\mathrm{O}_3} |
1162 |
\end{array} |
1163 |
\right. |
1164 |
\end{displaymath} |
1165 |
et le système (\ref{eq:qO3}) serait de calculer $\mu_\mathrm{moy}$ en |
1166 |
ne tenant compte que de la présence de l'eau : |
1167 |
\begin{displaymath} |
1168 |
\mu_\mathrm{moy} |
1169 |
= |
1170 |
\frac{\mu_\textrm{air sec}} |
1171 |
{1 + |
1172 |
\left( |
1173 |
\frac{\mu_\textrm{air sec}}{\mu_{\mathrm{H}_2\mathrm{O}}} - 1 |
1174 |
\right) |
1175 |
q_{\mathrm{H}_2\mathrm{O}}} |
1176 |
\end{displaymath} |
1177 |
|
1178 |
\subsubsection{Choix du type de remaillage} |
1179 |
|
1180 |
Copiant ce qui était déjà fait pour la vapeur d'eau, j'ai d'abord |
1181 |
choisi de remailler verticalement les paramètres de Cariolle en |
1182 |
interpolant une spline cubique de la pression. (Noter que ce n'est pas |
1183 |
une interpolation en logarithme de pression.) Si on interpole |
1184 |
simplement les données de Cariolle et que la pression à laquelle on |
1185 |
interpole est à l'extérieur du tableau de pressions de Cariolle alors |
1186 |
la formule d'interpolation devient sans crier gare une formule |
1187 |
d'extrapolation. C'est la formule pour l'intervalle extrême qui est |
1188 |
utilisée. Pour mieux contrôler cette extrapolation, j'ai ajouté aux |
1189 |
données de Cariolle un point correspondant à une pression nulle. Le |
1190 |
but était d'avoir approximativement une fraction molaire constante |
1191 |
au-dessus du point le plus haut de Cariolle et une production qui tend |
1192 |
vers 0. |
1193 |
|
1194 |
Mais l'interpolation avec une spline cubique produit de fortes |
1195 |
oscillations. (Elle donne d'ailleurs souvent des fractions molaires |
1196 |
négatives.) Cf. exécutions |
1197 |
\href{file:/user/guez/Documents/Utilisation_LMDZ/Utilisation_LMDZ.texfol/Utilisation-LMDZ.dvi}{28} |
1198 |
et |
1199 |
\href{file:///user/guez/Documents/Utilisation_LMDZ/Results_etat0_lim/check_coefoz_38.ps}{38} |
1200 |
de \verb+etat0_lim+. Le remaillage vertical est donc maintenant fait |
1201 |
par moyenne pour la plupart des paramètres chimiques, comme le |
1202 |
remaillage horizontal. |
1203 |
|
1204 |
Considérons le remaillage de la fraction molaire $r_\mathrm{Mob}$ (le |
1205 |
raisonnement est le même pour $P_\mathrm{net,Mob}$). Nous avons choisi |
1206 |
un remaillage par moyenne donc nous voulons associer à une cellule de |
1207 |
la grille tri-dimensionnelle de LMDZ une moyenne de $r_\mathrm{Mob}$ |
1208 |
sur cette cellule. Plus généralement, considérons un volume quelconque |
1209 |
de l'espace réel humain. Le plus sensé est d'associer à ce volume une |
1210 |
moyenne de $r_\mathrm{Mob}$ pondérée par la densité de l'air : |
1211 |
\begin{align*} |
1212 |
\langle r_\mathrm{Mob} \rangle |
1213 |
& = \frac{\iiint r_\mathrm{Mob} N \ud^3 \tau}{\iiint N \ud^3 \tau} \\ |
1214 |
& = \frac{\iiint (a + z)^2 r_\mathrm{Mob} N \cos\phi\,\ud z\, \ud \phi\, \ud \lambda} |
1215 |
{\iiint (a + z)^2 N \cos\phi\, \ud z\, \ud \phi\, \ud \lambda} |
1216 |
\end{align*} |
1217 |
Opérons le changement de variable $(z, \phi, \lambda) \to (p, \phi, |
1218 |
\lambda)$, avec l'hypothèse d'équilibre hydrostatique vertical. Notons |
1219 |
$\ud(z, \phi, \lambda) / \ud(p, \phi, \lambda)$ le jacobien de ce |
1220 |
changement de variable : |
1221 |
\begin{align*} |
1222 |
\ud z\, \ud \phi\, \ud \lambda |
1223 |
& = | \frac{\ud(z, \phi, \lambda)}{\ud(p, \phi, \lambda)} | |
1224 |
\ud p\, \ud \phi\, \ud \lambda \\ |
1225 |
& = | \partial_p z | \ud p\, \ud \phi\, \ud \lambda \\ |
1226 |
& = \frac{\ud p}{\rho g} \ud \phi\, \ud \lambda |
1227 |
\end{align*} |
1228 |
D'où : |
1229 |
\begin{displaymath} |
1230 |
\langle r_\mathrm{Mob} \rangle = \frac |
1231 |
{ |
1232 |
\iiint (a + z)^2 r_\mathrm{Mob} \cos\phi \frac{\ud p}{\mu_\mathrm{moy} g} |
1233 |
\ud \phi\, \ud \lambda |
1234 |
} |
1235 |
{ |
1236 |
\iiint (a + z)^2 \cos\phi \frac{\ud p}{\mu_\mathrm{moy} g} |
1237 |
\ud \phi\, \ud \lambda |
1238 |
} |
1239 |
\end{displaymath} |
1240 |
Si $z \ll a$ et $\mu_\mathrm{moy}$ et $g$ sont approximativement |
1241 |
constants sur le domaine considéré alors : |
1242 |
\begin{displaymath} |
1243 |
\langle r_\mathrm{Mob} \rangle \approx \frac |
1244 |
{\iiint r_\mathrm{Mob} \cos\phi\, \ud p\, \ud \phi\, \ud \lambda} |
1245 |
{\iiint \cos\phi \, \ud p\, \ud \phi\, \ud \lambda} |
1246 |
\end{displaymath} |
1247 |
(Attention : même si $r_\mathrm{Mob}$ ne dépend pas de $\lambda$, l'intégrale sur |
1248 |
$\lambda$ ne se simplifie que si le domaine en $(p, \phi)$ est |
1249 |
indépendant de $\lambda$.) Considérons maintenant le cas particulier |
1250 |
où le domaine est une cellule de LMDZ. Nous considérons cette cellule |
1251 |
comme un parallélépipède rectangle dans l'espace $(p, \phi, \lambda)$, |
1252 |
d'où : |
1253 |
\begin{displaymath} |
1254 |
\langle r_\mathrm{Mob} \rangle = \frac{1}{\Delta p \Delta(\sin\phi)} |
1255 |
\iint r_\mathrm{Mob} \cos\phi\, \ud p\, \ud \phi |
1256 |
\end{displaymath} |
1257 |
Cf. figure (\ref{fig:cell}). |
1258 |
\begin{figure}[htbp] |
1259 |
\centering |
1260 |
\includegraphics{Graphiques/Cellule} |
1261 |
\caption[Remaillage de la paramétrisation de la chimie de |
1262 |
l'ozone]{En noir, cellule sur laquelle on moyenne les données |
1263 |
concernant l'ozone. En rouge, cellules de Mobidic.} |
1264 |
\label{fig:cell} |
1265 |
\end{figure} |
1266 |
Nous n'avons pas a priori $r_\mathrm{Mob}$ en tout point $(p, \phi)$, |
1267 |
nous avons un nombre fini de valeurs de $r_\mathrm{Mob}$. Il nous |
1268 |
reste donc à supposer une forme de variation de $r_\mathrm{Mob}$ en |
1269 |
fonction de $(p, \phi)$. Le plus simple est de supposer que |
1270 |
$r_\mathrm{Mob}$ est constant sur des cellules rectangulaires du plan |
1271 |
$(p, \phi)$. Nous devons choisir les limites de ces cellules. En |
1272 |
latitude, nous prenons les milieux des valeurs de Mobidic. En |
1273 |
pression, nous prenons les moyennes géométriques des valeurs de |
1274 |
Mobidic (si bien que les limites de cellule en pression sont à |
1275 |
mi-chemin des valeurs de Mobidic sur un axe logarithmique). |
1276 |
|
1277 |
Pour les coefficients $a_2$ et $R_\mathrm{het}$, la démarche est moins évidente. Supposons |
1278 |
connu $a_2$ en tout point de ``l'espace humain'', quelle valeur |
1279 |
moyenne de $a_2$, $\overline{a_2}$, associer à un volume $\tau$ ? |
1280 |
Considérons la paramétrisation simplifiée : |
1281 |
\begin{displaymath} |
1282 |
P_\mathrm{net} = P_\mathrm{net,Mob} + a_2 (r - r_\mathrm{Mob}) |
1283 |
\end{displaymath} |
1284 |
Nous voudrions que : |
1285 |
\begin{displaymath} |
1286 |
\overline{P_\mathrm{net}} |
1287 |
= \overline{P_\mathrm{net,Mob}} |
1288 |
+ \overline{a_2} (\bar r - \overline{r_\mathrm{Mob}}) |
1289 |
\end{displaymath} |
1290 |
Ce qui donne : |
1291 |
\begin{displaymath} |
1292 |
\overline{a_2} |
1293 |
= \frac{\iiint a_2 (r - r_\mathrm{Mob}) N \ud^3 \tau} |
1294 |
{\iiint (r - r_\mathrm{Mob})N \ud^3 \tau} |
1295 |
\end{displaymath} |
1296 |
On peut adopter : |
1297 |
\begin{displaymath} |
1298 |
\overline{a_2} = \frac{\iiint a_2 N \ud^3 \tau}{\iiint N \ud^3 \tau} |
1299 |
\end{displaymath} |
1300 |
et garder ainsi une certaine compatibilité avec l'hypothèse selon |
1301 |
laquelle $r_\mathrm{Mob}$ est constant sur une cellule de |
1302 |
Mobidic. Idem pour $R_\mathrm{het}$. |
1303 |
|
1304 |
Pour $a_4$ et $T_\mathrm{Mob}$, nous calculons aussi la moyenne |
1305 |
pondérée par $N$. |
1306 |
|
1307 |
Pour $\Sigma_\mathrm{Mob}$, verticalement, une moyenne ne me |
1308 |
semblerait pas avoir de sens. Je choisis de faire une interpolation au |
1309 |
milieu de couche de LMDZ (pression \verb+pls+). Quelle interpolation ? |
1310 |
On a : |
1311 |
\begin{equation} |
1312 |
\label{eq:sigma} |
1313 |
\frac{\partial \Sigma_\mathrm{Mob}}{\partial z} = - r_\mathrm{Mob} N \approx - \frac{r_\mathrm{Mob} p}{k_B T} |
1314 |
\end{equation} |
1315 |
En remaillant $r_\mathrm{Mob}$ et $T$, nous avons supposé qu'ils |
1316 |
étaient constants sur des intervalles de pression centrés sur les |
1317 |
niveaux de pression de Mobidic. Considérons donc un tel intervalle, |
1318 |
centré sur la pression $p_1$, correspondant à une altitude $z_1$. |
1319 |
Supposons aussi $\mu_\mathrm{moy}$ et $g$ constants. Sur l'intervalle, |
1320 |
l'échelle $H$ de gradient de pression est constante : |
1321 |
\begin{displaymath} |
1322 |
p(z) = p(z_1) \exp\left(- \frac{z - z_1}{H}\right) |
1323 |
\end{displaymath} |
1324 |
Donc, en intégrant (\ref{eq:sigma}) sur $z$ : |
1325 |
\begin{displaymath} |
1326 |
\Sigma_\mathrm{Mob}(p) - \Sigma_\mathrm{Mob}(p_1) = \frac{r_\mathrm{Mob}}{\mu_\mathrm{moy} g} (p - p_1) |
1327 |
\end{displaymath} |
1328 |
$\Sigma_\mathrm{Mob}$ varie linéairement avec la pression. Nous connaissons les |
1329 |
valeurs de $\Sigma_\mathrm{Mob}$ aux niveaux de Mobidic, que nous supposons être |
1330 |
les centres des intervalles de fraction molaire constante, et non les |
1331 |
limites des intervalles. Nous pourrions commencer par calculer |
1332 |
$\Sigma_\mathrm{Mob}$ aux limites des intervalles. Notons $p_1$ et $p_2$ deux |
1333 |
niveaux de pression successifs de Mobidic, $r_{\mathrm{Mob},1}$ et $r_{\mathrm{Mob},2}$ les |
1334 |
fractions molaires correspondantes, $p_l = \sqrt{p_1 p_2}$ la |
1335 |
limite entre les deux intervalles de centres $p_1$ et $p_2$. Posons : |
1336 |
\begin{displaymath} |
1337 |
\alpha = \frac{r_{\mathrm{Mob},1}}{r_{\mathrm{Mob},2}} \frac{p_l - p_1}{p_2 - p_l} |
1338 |
\end{displaymath} |
1339 |
On a : |
1340 |
\begin{align*} |
1341 |
& \Sigma_\mathrm{Mob}(p_l) - \Sigma_{\mathrm{Mob},1} |
1342 |
= \frac{r_{\mathrm{Mob},1}}{\mu_\mathrm{moy} g} (p_l - p_1) \\ |
1343 |
& \Sigma_{\mathrm{Mob},2} - \Sigma_\mathrm{Mob}(p_l) |
1344 |
= \frac{r_{\mathrm{Mob},2}}{\mu_\mathrm{moy} g} (p_2 - p_l) |
1345 |
\end{align*} |
1346 |
donc : |
1347 |
\begin{displaymath} |
1348 |
\Sigma_\mathrm{Mob}(p_l) = \frac{\Sigma_{\mathrm{Mob},1} + \alpha \Sigma_{\mathrm{Mob},2}}{1 + \alpha} |
1349 |
\end{displaymath} |
1350 |
Plus simplement, je choisis d'interpoler linéairement en pression |
1351 |
entre les valeurs de Mobidic. |
1352 |
|
1353 |
Pour $\Sigma_\mathrm{Mob}$ horizontalement, nous moyennons simplement pour |
1354 |
conserver le nombre de molécules d'ozone au-dessus d'une surface : |
1355 |
\begin{align*} |
1356 |
\bar \Sigma_\mathrm{Mob} & = \frac{\iint \Sigma_\mathrm{Mob} \ud^2 S}{S} \\ |
1357 |
& = |
1358 |
\frac{\iint \Sigma_\mathrm{Mob} (a + z^2) \cos\phi \ud \phi \ud \lambda} |
1359 |
{\iint (a + z^2) \cos\phi \ud \phi \ud \lambda} \\ |
1360 |
& \approx \frac{\iint \Sigma_\mathrm{Mob} \cos\phi \ud \phi \ud \lambda} |
1361 |
{\iint \cos\phi \ud \phi \ud \lambda} |
1362 |
\end{align*} |
1363 |
|
1364 |
Pour $a_6$, une moyenne pondérée par $N$ semble censée. $a_6(p) (\Sigma(p) - |
1365 |
\Sigma_\mathrm{Mob}(p))$ est une contribution à |
1366 |
$P_\mathrm{net}(p)$. LMDZ calcule la moyenne |
1367 |
$\overline{P_\mathrm{net}}$ pondérée par $N$ sur une cellule. Une |
1368 |
contribution à cette moyenne est : |
1369 |
\begin{displaymath} |
1370 |
\frac{\int a_6 N (\Sigma - \Sigma_\mathrm{Mob}) \ud \tau}{\int N \ud \tau} |
1371 |
\end{displaymath} |
1372 |
$\overline{a_6} (\Sigma - \Sigma_\mathrm{Mob})$, avec $\Sigma$ la |
1373 |
valeur en milieu de couche dans LMDZ et $\Sigma_\mathrm{Mob}$ la |
1374 |
valeur interpolée, se rapproche de cette contribution. |
1375 |
|
1376 |
Cf. figures (\ref{fig:regr_latit}), (\ref{fig:regr_p_stepav}) et |
1377 |
(\ref{fig:regr_p_lint}). |
1378 |
\begin{figure}[htbp] |
1379 |
\centering |
1380 |
\includegraphics{Graphiques/regr_latit} |
1381 |
\caption[Remaillage en latitude de la paramétrisation de la chimie |
1382 |
de l'ozone]{Remaillage en latitude. En bleu, les valeurs d'origine, |
1383 |
en rouge les valeurs utilisées.} |
1384 |
\label{fig:regr_latit} |
1385 |
\end{figure} |
1386 |
\begin{figure}[htbp] |
1387 |
\centering |
1388 |
\includegraphics{Graphiques/regr_p_stepav} |
1389 |
\caption[Remaillage de la paramétrisation de la chimie de l'ozone, |
1390 |
en pression par moyenne de fonction en escalier] {Remaillage en |
1391 |
pression par moyenne de fonction en escalier. En bleu, les |
1392 |
valeurs d'origine, en rouge les valeurs utilisées. La position de |
1393 |
\texttt{p3d(i, j, 1)} par rapport à \texttt{plev(n\_plev)} est |
1394 |
quelconque.} |
1395 |
\label{fig:regr_p_stepav} |
1396 |
\end{figure} |
1397 |
\begin{figure}[htbp] |
1398 |
\centering |
1399 |
\includegraphics{Graphiques/regr_p_lint} |
1400 |
\caption[Remaillage de la paramétrisation de la chimie de l'ozone, |
1401 |
en pression par interpolation linéaire]{Remaillage en pression par |
1402 |
interpolation linéaire. Les calculs aux valeurs élevées de |
1403 |
\texttt{pls} (premiers indices de niveau vertical) peuvent être |
1404 |
des extrapolations. La valeur supplémentaire à pression nulle du |
1405 |
maillage source contrôle l'extrapolation aux basses pressions.} |
1406 |
\label{fig:regr_p_lint} |
1407 |
\end{figure} |
1408 |
|
1409 |
\subsubsection{Remaillage : dans quel programme, à quel moment ? Quelles |
1410 |
entrées-sorties ?} |
1411 |
|
1412 |
Les paramètres chimiques doivent être remaillés en latitude pour LMDZ. |
1413 |
Ce remaillage en latitude peut être fait une fois pour toutes, dans |
1414 |
\verb+etat0_lim+. Par ailleurs, au cours d'une simulation, la pression |
1415 |
à une longitude, une latitude et un niveau vertical donnés change |
1416 |
parce que la pression à la surface change. Donc il faut en principe |
1417 |
faire le remaillage en pression des données de Cariolle à chaque pas |
1418 |
de temps physique. Enfin, nous voulons lisser les variations |
1419 |
temporelles des paramètres chimiques en les interpolant |
1420 |
journalièrement. Ce remaillage temporel peut être fait une fois pour |
1421 |
toutes, dans \verb+etat0_lim+. |
1422 |
|
1423 |
En pratique, le remaillage en pression à chaque pas de temps n'est pas |
1424 |
indispensable parce que l'essentiel de la chimie de l'ozone, avec |
1425 |
notre paramétrisation, se déroule à haute altitude, où les niveaux |
1426 |
verticaux de LMDZ sont pratiquement des niveaux de pression fixés. |
1427 |
Nous pouvons donc nous demander si, pour simplifier le programme, nous |
1428 |
pouvons faire les trois remaillages (latitude, pression, temps) dans |
1429 |
\verb+etat0_lim+. Nous négligerions alors pour le remaillage vertical |
1430 |
l'évolution des pressions aux interfaces des couches. Mais la |
1431 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Guide_IPSL_climate_models/taille.ods}{taille |
1432 |
du fichier contenant les paramètres remaillés} serait alors trop |
1433 |
grande.\footnote{Le remaillage en pression dans \texttt{etat0\_lim} |
1434 |
serait une bonne solution si nous utilisions les paramètres mensuels |
1435 |
sans interpolation journalière.} |
1436 |
|
1437 |
Nous devons donc faire un des trois remaillages dans \verb+gcm+. |
1438 |
Lequel ? On ne peut faire le remaillage vertical sans avoir fait le |
1439 |
remaillage en latitude. Le choix dans \verb+gcm+ se restreint donc aux |
1440 |
remaillages vertical et temporel. Pour le calcul de la production |
1441 |
chimique d'ozone, nous avons besoin à chaque pas de temps de cinq |
1442 |
coefficients tri-dimensionnels (taille totale : $(\mathtt{iim} + 1) |
1443 |
\times (\mathtt{jjm} + 1) \times \mathtt{llm} \times 5$). Si nous |
1444 |
partons des coefficients remaillés en pression dans \verb+etat0_lim+, |
1445 |
nous devons calculer ces cinq coefficients par interpolation |
1446 |
temporelle pour chaque jour. Si nous partons des coefficients |
1447 |
interpolés journalièrement dans \verb+etat0_lim+, nous pouvons aussi |
1448 |
nous contenter de calculer ces cinq coefficients par remaillage |
1449 |
vertical pour chaque jour. Les deux transformations : |
1450 |
\begin{displaymath} |
1451 |
\xymatrix{ |
1452 |
(\mathtt{jjm} + 1) \times 45 \ar[d]_{\textrm{remaillage vertical}} |
1453 |
& (\mathtt{iim} + 1) \times (\mathtt{jjm} + 1) \times \mathtt{llm} |
1454 |
\times 12 \ar[ld]^{\textrm{remaillage temporel}} \\ |
1455 |
(\mathtt{iim} + 1) \times (\mathtt{jjm} + 1) \times \mathtt{llm} |
1456 |
} |
1457 |
\end{displaymath} |
1458 |
demandent des nombres d'opérations du même ordre de grandeur : |
1459 |
$(\mathtt{iim} + 1) \times (\mathtt{jjm} + 1) \times \mathtt{llm}$. Le |
1460 |
plus sensé est donc de faire le remaillage en pression dans |
1461 |
\verb+gcm+. Pour que l'arrêt et le redémarrage du programme \verb+gcm+ |
1462 |
soient sans conséquence, nous ne pouvons faire le remaillage en |
1463 |
pression moins d'une fois par jour. (En particulier, nous ne pouvons |
1464 |
le faire une fois par exéuction de \verb+gcm+, au début de \verb+gcm+, |
1465 |
à partir de la valeur courante du champ de pression.) Nous faisons |
1466 |
donc le remaillage vertical pour chaque jour, négligeant l'évolution |
1467 |
des pressions aux interfaces des couches pendant un jour. |
1468 |
|
1469 |
Nous devons encore nous demander s'il vaut mieux lire tout le contenu |
1470 |
de \verb+coefoz_LMDZ.nc+ au début de l'exécution de \verb+gcm+ ou lire |
1471 |
dans le fichier une fois par jour les paramètres au bon jour. Si nous |
1472 |
lisons tout en une seule fois, nous gagnons en temps d'exécution (nous |
1473 |
limitons les accès au disque) mais nous augmentons la |
1474 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Guide_IPSL_climate_models/taille.ods}{mémoire |
1475 |
principale consommée}. La lecture une fois par jour simulé est-elle |
1476 |
trop fréquente ? On peut se référer à une durée de 40 s de temps |
1477 |
écoulé par jour simulé en résolution $96 \times 72 \times 50$, sur |
1478 |
deux processus MPI. Pour le moment, la programmation de la lecture une |
1479 |
fois par jour me paraît plus simple. |
1480 |
|
1481 |
\verb+etat0_lim+ fait le remaillage en latitude et en temps des huit |
1482 |
coefficients et le remaillage vertical du coefficient $r$ seulement, |
1483 |
au jour initial seulement, pour servir d'état initial de l'abondance |
1484 |
d'ozone. Nous ne pouvons donc pas faire le remaillage en latitude et |
1485 |
temps après \verb+etat0+. |
1486 |
|
1487 |
Dans quel fichier \verb+etat0_lim+ doit-il écrire les coefficients |
1488 |
remaillés ? Le fichier \verb+limit.nc+ n'est pas très adapté car les |
1489 |
coefficients ne sont pas sur la grille de la ``physique''. De plus, |
1490 |
les paramètres chimiques ne sont pas des conditions à la limite. Par |
1491 |
ailleurs, les paramètres chimiques sont les mêmes dans plusieurs |
1492 |
exécutions successives de \verb+gcm+, donc ils ne peuvent pas être |
1493 |
dans \verb+start+ ou \verb+startphy+. Nous créons donc un nouveau |
1494 |
fichier \verb+coefoz_LMDZ.nc+. Cf. |
1495 |
\href{file:///user/guez/Documents/Informatique_fonctionnement/Programs/Guide_IPSL_climate_models/inout |
1496 |
LMDZE.odg}{schéma des entrées-sorties de \texttt{etat0\_lim} et |
1497 |
\texttt{gcm}}. \verb+coefoz_LMDZ.nc+ contient donc les données de |
1498 |
Cariolle remaillées en latitude pour LMDZ et interpolées |
1499 |
journalièrement. |
1500 |
|
1501 |
Principe de la procédure \verb+regr_lat_time_coefoz+ : |
1502 |
\begin{algorithmic} |
1503 |
\FOR{chaque coefficient à remailler} |
1504 |
\STATE remaillage en latitude |
1505 |
\STATE remaillage en temps |
1506 |
\ENDFOR |
1507 |
\end{algorithmic} |
1508 |
Cf. figure (\ref{fig:regr_lat_time_coefoz}). |
1509 |
\begin{figure}[htbp] |
1510 |
\centering |
1511 |
\includegraphics{Graphiques/regr_lat_time_coefoz} |
1512 |
\caption{Schéma des entrées-sorties de |
1513 |
\texttt{regr\_lat\_time\_coefoz}.} |
1514 |
\label{fig:regr_lat_time_coefoz} |
1515 |
\end{figure} |
1516 |
|
1517 |
\subsubsection{Calcul de l'évolution de la distribution d'ozone dans \texttt{phytrac}} |
1518 |
|
1519 |
D'après l'équation (\ref{eq:dqdt}) : |
1520 |
\begin{displaymath} |
1521 |
\partial_t q |
1522 |
= - \vec U \cdot \vec\nabla q |
1523 |
+ P_\mathrm{net} \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1524 |
\end{displaymath} |
1525 |
Notons $\vec r$ le vecteur position. \verb+phytrac+ calcule : |
1526 |
\begin{displaymath} |
1527 |
q(\vec r, t + \mathtt{pdtphys}) = q(\vec r, t) + \textrm{termes de transport} |
1528 |
+ \int_t ^{t + \mathtt{pdtphys}} P_\mathrm{net}(\vec r, t') |
1529 |
\frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} \ud t' |
1530 |
\end{displaymath} |
1531 |
Le dernier terme est la contribution de la chime de l'ozone. Nous |
1532 |
avons besoin de $P_\mathrm{net}$ à chaque point de la grille |
1533 |
horizontale de LMDZ et à chaque niveau vertical de LMDZ. Le fichier |
1534 |
\verb+coefoz_LMDZ.nc+ contenant les grandeurs $P_\mathrm{net,Mob}$, |
1535 |
$a_2$, $r_\mathrm{Mob}$, $a_4$, $T _\mathrm{Mob}$, $a_6$, |
1536 |
$\Sigma_\mathrm{Mob}$ et $R_\mathrm{het}$, \verb+phytrac+ doit |
1537 |
intégrer : |
1538 |
\begin{multline*} |
1539 |
P_\mathrm{net} \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1540 |
= P_\mathrm{net,Mob} \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1541 |
+ a_2 \left(q - r_\mathrm{Mob} \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}}\right) |
1542 |
+ a_4 \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1543 |
(T - T_\mathrm{Mob}) \\ |
1544 |
+ \frac{a_6}{\mu_\mathrm{moy}} |
1545 |
\left( |
1546 |
\frac{1}{r^2} \int_r ^{+ \infty} q \rho r'^2 \ud r' |
1547 |
- \mu_{\mathrm{O}_3} \Sigma_\mathrm{Mob} |
1548 |
\right) |
1549 |
+ R_\mathrm{het,corr} q |
1550 |
\end{multline*} |
1551 |
|
1552 |
Dans \verb+phytrac+, nous n'avons accès qu'à la température et à la |
1553 |
masse volumique à l'instant courant, et à aucun autre instant. |
1554 |
Calculons donc l'intégrale comme si $T$ et $\rho$ étaient des |
1555 |
constantes pendant \verb+pdtphys+. La direction du Soleil, qui |
1556 |
intervient dans le terme de chimie hétérogène, pourrait être |
1557 |
considérée comme une fonction du temps. Je choisis pour simplifier de |
1558 |
considérer aussi la direction du Soleil comme une constante pendant un |
1559 |
pas de temps de la physique. Récrivons alors $P_\mathrm{net} |
1560 |
\mu_{\mathrm{O}_3} / \mu_\mathrm{moy}$ en mettant en évidence les |
1561 |
constantes. Posons : |
1562 |
\begin{align*} |
1563 |
& c := [P_\mathrm{net,Mob} - a_2 r_\mathrm{Mob} |
1564 |
+ a_4 (T - T_\mathrm{Mob}) - a_6 \Sigma_\mathrm{Mob}] |
1565 |
\frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} \\ |
1566 |
& b := a_2 + R_\mathrm{het,corr} \\ |
1567 |
& a_{6m} := \frac{a_6}{\mu_\mathrm{moy}} |
1568 |
\end{align*} |
1569 |
$c$, $b$ et $a_{6m}$ sont des champs constants connus pendant |
1570 |
\verb+pdtphys+. On a : |
1571 |
\begin{equation} |
1572 |
\label{eq:dq_dt} |
1573 |
P_\mathrm{net} \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1574 |
= c + b q + a_{6m} \frac{1}{r^2} \int_r ^{+ \infty} q \rho r'^2 \ud r' |
1575 |
\end{equation} |
1576 |
Si on calcule l'évolution de $q$ en ne tenant compte |
1577 |
que de la chimie (c'est-à-dire sans tenir compte des termes de |
1578 |
transport) alors on obtient : |
1579 |
\begin{displaymath} |
1580 |
\partial_t q = P_\mathrm{net} \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1581 |
\end{displaymath} |
1582 |
|
1583 |
Limitons-nous d'abord aux deux premiers termes dans le membre de |
1584 |
droite de (\ref{eq:dq_dt}) : |
1585 |
\begin{equation} |
1586 |
\label{eq:dq_dt_c_a2} |
1587 |
\partial_t q = c + b q |
1588 |
\end{equation} |
1589 |
La solution exacte de cette équation différentielle sur $q$ donne : |
1590 |
\begin{equation} |
1591 |
\label{eq:delta_q_exact} |
1592 |
\begin{array}{|ll} |
1593 |
\delta q = (c + b q) \frac{e^{b \delta t} - 1}{b} |
1594 |
& \textrm{si } b \ne 0 \\ |
1595 |
\delta q = c \delta t & \textrm{si } b = 0 |
1596 |
\end{array} |
1597 |
\end{equation} |
1598 |
pour un intervalle de temps $\delta t$ quelconque. \`A l'ordre 2 en |
1599 |
$b \delta t$, nous avons donc, que $b$ soit nul ou non : |
1600 |
\begin{equation} |
1601 |
\label{eq:delta_q} |
1602 |
\delta q = (c + b q) (1 + b \delta t / 2) \delta t |
1603 |
\end{equation} |
1604 |
(Ce qui revient à avoir intégré l'équation différentielle par la |
1605 |
méthode du point milieu. Cf. Press et al. [1992 \# 318, équation |
1606 |
(16.1.2)].) Pour $\delta t = 1800$ s, on a : |
1607 |
\begin{displaymath} |
1608 |
4 \cdot 10^{-11} \le |a_2| \delta t \le \nombre{0,02} |
1609 |
\end{displaymath} |
1610 |
Entre la valeur calculée de $e^x -1$ et la valeur calculée de $x + x^2 |
1611 |
/ 2$, laquelle est la plus proche de la valeur exacte de $e^x -1$ ? |
1612 |
Pour $x \ge 10^{-2}$, c'est la valeur calculée de $e^x -1$. Pour $0 |
1613 |
\le x \le 10^{-3}$, c'est la valeur calculée de $x + x^2 / 2$. En |
1614 |
conséquence, l'utilisation de (\ref{eq:delta_q}) à la place de |
1615 |
(\ref{eq:delta_q_exact}) est un bon choix. |
1616 |
|
1617 |
Si on tient maintenant compte du terme de densité-colonne dans |
1618 |
(\ref{eq:dq_dt}), nous avons une équation intégro-différentielle aux |
1619 |
dérivées partielles. Nous commençons par une discrétisation spatiale |
1620 |
verticale, à une position horizontale fixée. Nous avons \verb+llm+ |
1621 |
couches d'atmosphère et nous supposons que $q$ est uniforme dans |
1622 |
chaque couche. La fonction inconnue $q(\vec r, t)$ devient alors un |
1623 |
vecteur de fonctions du temps : $[q_1(t), \dots, q_\mathtt{llm}(t)]$. |
1624 |
La fonction connue $c(\vec r)$ devient un vecteur de constantes |
1625 |
connues : $[c_1, \dots, c_\mathtt{llm}]$, idem pour $b$ et $a_{6m}$. |
1626 |
Notons $r_k$ la position de la limite entre les couches $k-1$ |
1627 |
et $k$. Pour chaque couche $k$, nous avons : |
1628 |
\begin{displaymath} |
1629 |
\frac{\ud q_k}{\ud t} |
1630 |
= c_k + b_k q_k |
1631 |
+ a_{6m,k} \frac{1}{r_k^2} |
1632 |
\sum_{k'=k} ^{\mathtt{llm}} q_{k'} \int_{r_{k'}} ^{r_{k'+1}} \rho r^2 \ud r |
1633 |
\end{displaymath} |
1634 |
Supposons que la différence entre $r_\mathtt{llm}$ et |
1635 |
le rayon de la Terre soit négligeable. Alors pour $k \le k' \le \mathtt{llm} |
1636 |
- 1$ : |
1637 |
\begin{equation} |
1638 |
\label{eq:approx_column} |
1639 |
\frac{1}{r_k^2} \int_{r_{k'}} ^{r_{k'+1}}\rho r'^2 \ud r' |
1640 |
\approx \int_{z_{k'}} ^{z_{k'+1}} \rho \ud z |
1641 |
\end{equation} |
1642 |
Mais le cas de la couche supérieure est plus problématique : |
1643 |
$r_{\mathtt{llm}+1} = + \infty$. Nous pouvons faire l'approximation |
1644 |
(\ref{eq:approx_column}) pour $k' = \mathtt{llm}$ si $\rho$ décroît |
1645 |
suffisamment rapidement au dessus de $r_\mathtt{llm}$. C'est |
1646 |
l'approximation qui est faite dans \verb+phytrac+ pour le calcul de |
1647 |
\verb+zmasse+, et que je garde. Alors : |
1648 |
\begin{equation} |
1649 |
\label{eq:dqk_dt} |
1650 |
\frac{\ud q_k}{\ud t} |
1651 |
= c_k + b_k q_k |
1652 |
+ a_{6m,k} \sum_{k'=k} ^{\mathtt{llm}} q_{k'} \mathtt{zmasse}_{k'} |
1653 |
\end{equation} |
1654 |
Nous nous sommes ramenés à un système d'équations différentielles du |
1655 |
premier ordre. |
1656 |
|
1657 |
La somme partielle dans (\ref{eq:dqk_dt}) est la densité-colonne |
1658 |
d'ozone au dessus de la base de la couche $k$. Nous pourrions raffiner |
1659 |
en calculant à partir de là des densités-colonnes au dessus des |
1660 |
milieux de couches (par interpolation linéaire en pression). Je reste |
1661 |
pour l'instant à l'équation (\ref{eq:dqk_dt}). |
1662 |
|
1663 |
Les équations différentielles du premier ordre du système |
1664 |
(\ref{eq:dqk_dt}) sont linéaires à coefficients constants. Nous |
1665 |
pourrions chercher une solution exacte. Considérant l'étude que j'ai |
1666 |
faite de l'équation (\ref{eq:dq_dt_c_a2}), je choisis d'appliquer |
1667 |
simplement la méthode du point milieu. |
1668 |
|
1669 |
$c$ et $b$ doivent être recalculés à chaque pas de temps de la |
1670 |
physique mais, parmi les termes composant $c$, les termes provenant de |
1671 |
Mobidic ne varient que pour le remaillage en pression: leur |
1672 |
combinaison peut être calculée une fois par jour. Distinguons donc : |
1673 |
\begin{displaymath} |
1674 |
c_\mathrm{Mob} |
1675 |
:= |
1676 |
(P_\mathrm{net,Mob} - a_2 r_\mathrm{Mob} |
1677 |
- a_6 \Sigma_\mathrm{Mob} - a_4 T_\mathrm{Mob}) |
1678 |
\frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} |
1679 |
\end{displaymath} |
1680 |
$c_\mathrm{Mob}$ regroupe les quatre paramètres ``climatologiques''. |
1681 |
On a : |
1682 |
\begin{displaymath} |
1683 |
c = c_\mathrm{Mob} + a_4 \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} T |
1684 |
\end{displaymath} |
1685 |
De même, pour le calcul de $b$, le filtrage de la chimie hétérogène |
1686 |
selon la latitude et la prise en compte du coefficient Clx peuvent |
1687 |
être faits une fois par jour. En résumé, nous avons distingué : |
1688 |
\begin{itemize} |
1689 |
\item la variable $q$ à intégrer pendant un pas de temps de la |
1690 |
physique, qui est donc considérée comme dépendant du temps pendant |
1691 |
\verb+pdtphys+ ; |
1692 |
\item les variables considérées comme constantes pendant |
1693 |
\verb+pdtphys+ mais qui changent d'un pas de temps de la physique à |
1694 |
l'autre ; |
1695 |
\item les variables qui ne changent qu'une fois par jour. |
1696 |
\end{itemize} |
1697 |
|
1698 |
\subsubsection{Les procédures \texttt{regr\_pr\_comb\_coefoz} et |
1699 |
\texttt{regr\_pr\_(av|int)\_coefoz}} |
1700 |
|
1701 |
Les procédures du module \verb+regr_pr+ (cf. figure |
1702 |
(\ref{fig:regridding_proced})) fournissent des variables ayant un |
1703 |
double indice horizontal alors que \verb+tr_seri+ dans \verb+phytrac+ |
1704 |
a un indice horizontal simple. |
1705 |
\begin{figure}[htbp] |
1706 |
\centering |
1707 |
\includegraphics{Graphiques/regridding_proced} |
1708 |
\caption{Arbre des appels des procédures de remaillage. Les cadres |
1709 |
verts symbolisent des modules. Les procédures du plus bas niveau |
1710 |
sont celles qui effectuent réellement le remaillage. Elles sont |
1711 |
les plus générales et ne font pas référence aux notions de |
1712 |
pression, latitude ou temps.} |
1713 |
\label{fig:regridding_proced} |
1714 |
\end{figure} |
1715 |
Il faut lire dans le fichier \verb+coefoz_LMDZ.nc+, remailler en |
1716 |
pression et faire la transformation d'indiçage une seule fois par jour |
1717 |
de simulation. De même, $c_\mathrm{Mob}$, $a_4 \mu_{\mathrm{O}_3} / |
1718 |
\mu_\mathrm{moy}$ et $a_{6m}$ doivent être calculés une seule fois par |
1719 |
jour. Toutes ces opérations sont faites dans la procédure |
1720 |
\verb+regr_pr_comb_coefoz+. Quand pouvons-nous appeler |
1721 |
\verb+regr_pr_comb_coefoz+ ? \verb+calfis+ appelle forcément une seule |
1722 |
fois \verb+physiq+, qui appelle forcément une seule fois |
1723 |
\verb+phytrac+. Nous pourrions donc faire l'appel dans n'importe |
1724 |
laquelle de ces trois procédures. Nous choisissons \verb+phytrac+ pour |
1725 |
une plus grande modularité. |
1726 |
|
1727 |
La procédure \verb+regr_pr_comb_coefoz+ lit en entier les huit |
1728 |
coefficients chimiques au jour courant et met à jour les cinq |
1729 |
variables de module : \verb+a2+, \verb+a4_mass+, \verb+c_mob+, |
1730 |
\verb+a6_mass+, \verb+r_het_interm+. Nous regroupons ici les |
1731 |
définitions des quatre dernières : |
1732 |
\begin{align*} |
1733 |
& a_{4m} := a_4 \frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} \\ |
1734 |
& c_\mathrm{Mob} |
1735 |
:= (P_\mathrm{net,Mob} - a_2 r_\mathrm{Mob} - a_6 \Sigma_\mathrm{Mob}) |
1736 |
\frac{\mu_{\mathrm{O}_3}}{\mu_\mathrm{moy}} - a_{4m} T_\mathrm{Mob} \\ |
1737 |
& a_{6m} := \frac{a_6}{\mu_\mathrm{moy}} \\ |
1738 |
& R_\mathrm{het,interm} = R_\mathrm{het} 1_{\phi \notin [- 45°, 45°]} |
1739 |
\left(\frac{\mathrm{Clx}}{\nombre{3,8} \cdot 10^{-9}}\right)^2 |
1740 |
\end{align*} |
1741 |
Chacune des cinq variables de module est de profil \verb+(klon, llm)+. |
1742 |
|
1743 |
La procédure \verb+regr_pr_(av|int)_coefoz+ est appelée une fois par |
1744 |
jour, pour chaque coefficient chimimque. Cf. figure (\ref{fig:regr_pr_coefoz}). |
1745 |
\begin{figure}[htbp] |
1746 |
\centering |
1747 |
\includegraphics{Graphiques/regr_pr_coefoz} |
1748 |
\caption{Entrées-sorties de la procédure |
1749 |
\texttt{regr\_pr\_(av|int)\_coefoz}. En rouge : les données, en |
1750 |
noir : les traitements.} |
1751 |
\label{fig:regr_pr_coefoz} |
1752 |
\end{figure} |
1753 |
|
1754 |
|
1755 |
\subsubsection{Vue d'ensemble informatique} |
1756 |
|
1757 |
Cf. figure (\ref{fig:ozone}). |
1758 |
\begin{figure}[htbp] |
1759 |
\centering |
1760 |
\includegraphics{Graphiques/ozone} |
1761 |
\caption{Vue d'ensemble de ce qui concerne le traceur ozone.} |
1762 |
\label{fig:ozone} |
1763 |
\end{figure} |
1764 |
Les résultats du programme \texttt{gcm} pour l'ozone apparaissent dans |
1765 |
cinq fichiers (cf. figure (\ref{fig:ozone})). Dans \verb+histrac.nc+, |
1766 |
ce sont les variables : |
1767 |
\begin{itemize} |
1768 |
\item \verb+float O3(time_counter, presnivs, lat, lon)+ |
1769 |
\item \verb+flO3+ flux |
1770 |
\item \verb+d_tr_th_O3+ tendance thermique |
1771 |
\item \verb+d_tr_cv_O3+ tendance convection |
1772 |
\item \verb+d_tr_cl_O3+ tendance couche limite |
1773 |
\end{itemize} |
1774 |
La différence entre \verb+O3+ dans \verb+histrac.nc+ et : |
1775 |
\begin{verbatim} |
1776 |
float O3VL1(time_counter, sigss, lat, lon) |
1777 |
\end{verbatim} |
1778 |
dans \verb+dyn_hist_ave.nc+ n'est pas claire pour moi. Une longitude |
1779 |
en plus, 180°, dans \verb+dyn_hist_ave.nc+. Les latitudes sont les |
1780 |
mêmes. Cf. |
1781 |
\hyperref{file:/user/guez/Documents/Utilisation_LMDZ/Utilisation_LMDZ.texfol/Utilisation-LMDZ.dvi}{section}{ozone}{résultats}. |
1782 |
Faire varier les paramètres dans \verb+physiq.def+ : |
1783 |
\begin{verbatim} |
1784 |
OK_mensuel=y |
1785 |
OK_instan=y |
1786 |
lev_histhf=4 |
1787 |
lev_histmth=2 |
1788 |
\end{verbatim} |
1789 |
n'apporte aucune sortie supplémentaire sur l'ozone. Cf. exécutions 58 |
1790 |
à 63 de \verb+gcm+. |
1791 |
|
1792 |
Diverses procédures s'occupent de remaillage. Elles sont de |
1793 |
généralités différentes. Cf. figure (\ref{fig:regridding_proced}). |
1794 |
\begin{description} |
1795 |
\item[\texttt{regr\_pr\_comb\_coefoz}] : lecture, remaillage en pression, |
1796 |
empaquetage, combinaison |
1797 |
\item[\texttt{regr\_pr\_coefoz}] : lecture, remaillage en pression, empaquetage |
1798 |
\item[\texttt{regr\_lat\_time\_coefoz}] : lecture, remaillage en temps et |
1799 |
latitude, écriture |
1800 |
\item[\texttt{regr\_pr\_o3}] : lecture, remaillage en pression |
1801 |
\item[\texttt{regr\_pr}] : remaillage en pression |
1802 |
\item[\texttt{regr1\_lint}, \texttt{regr1\_step\_av}, |
1803 |
\texttt{regr3\_lint}] : remaillage |
1804 |
\end{description} |
1805 |
|
1806 |
|
1807 |
\section{Versions de \texttt{etat0\_lim}} |
1808 |
|
1809 |
Problème de l'appel à \verb+startget_dyn+ pour \verb+vvent+. Cf. |
1810 |
Flyspray, tâche numéro 29. Dans mon rapport, j'aurais dû ajouter que |
1811 |
le même problème que pour \verb+pls+ se posait pour \verb+workvar+. Si |
1812 |
on corrige en sélectionnant les indices de latitude de 1 à \verb+jjm+ |
1813 |
dans les arguments effectifs \verb+pls+ et \verb+workvar+ alors on |
1814 |
obtient des modifications dans la sortie standard et dans le fichier |
1815 |
\verb+start.nc+ uniquement. Dans la sortie standard, des petites |
1816 |
différences seulement, entre autres dans les valeurs de \verb+V+ |
1817 |
\verb+min+ \verb+max+ affichées après l'appel incriminé de |
1818 |
\verb+startget_dyn+. Dans le fichier \verb+start.nc+, les seules |
1819 |
variables différentes sont \verb+controle+ et \verb+vcov+. |
1820 |
|
1821 |
Modifications provisoires : |
1822 |
\begin{description} |
1823 |
\item[M1] Mise en commentaire de l'appel à \verb+startget+ pour |
1824 |
\verb+phis+ dans \verb+etat0+. |
1825 |
\item[M2] \verb+q3d(:, :, :, 5) = 3.+ dans \verb+etat0+. |
1826 |
\end{description} |
1827 |
|
1828 |
Modifications cumulatives : |
1829 |
\begin{description} |
1830 |
\item[L1] Passage à \verb+netcdf95+ dans \verb+limit.f+. Correction |
1831 |
lecture \verb+Rugos+. Correction dans \verb+etat0+, appel à |
1832 |
\verb+startget_dyn+ pour \verb+vvent+. Version 2.10a de |
1833 |
\verb+spline+ et \verb+splint+. |
1834 |
\item[L2] Initialisation d'un traceur ozone (sans rétroaction sur |
1835 |
l'atmosphère) à partir des données de Mobidic. Modification de |
1836 |
\verb+nqmx+, 5 au lieu de 4 dans \verb+dimens_m.f+. |
1837 |
\item[L3] Correction dans \verb+etat0+ : \verb+q3d+ doit contenir des |
1838 |
fractions massiques et non molaires. |
1839 |
\item[L4] Correction : attribut \verb+save+ dans les modules |
1840 |
\verb+comvert+, \verb+serre+, \verb+clesphys+, \verb+iniprint+, |
1841 |
\verb+ener+, \verb+comdissnew+, \verb+dimphy+, \verb+logic+, \verb+comgeom+. |
1842 |
\item[L5] \'Ecriture du fichier \verb+coefoz_LMDZ.nc+ (ajout de la |
1843 |
procédure \verb+aver_interp_coefoz+), lecture de la distribution |
1844 |
initiale d'ozone dans ce fichier (suppression de la procédure |
1845 |
\verb+initial_o3+). |
1846 |
\item[L6] Pour l'écriture de \verb+coefoz_LMDZ.nc+, modification de |
1847 |
l'extrapolation aux basses pressions. Nouveau cas \verb+strato1+ dans |
1848 |
\verb+disvert+. |
1849 |
\item[L7] Création des modules \verb+regr_coefoz_m+, \verb+pressure_m+ |
1850 |
et \verb+regrid+. Clarification de \verb+inter_barxy+. Modification |
1851 |
du remaillage en latitude des coefficients pour l'ozone : en |
1852 |
supposant une fonction en escalier et non une fonction linéaire par |
1853 |
morceaux. |
1854 |
\item[L8] Modification du remaillage en latitude des coefficients pour |
1855 |
l'ozone : pondération avec le cosinus de la latitude. |
1856 |
\item[L9] Modification du remaillage en pression des coefficients pour |
1857 |
l'ozone : moyenne au lieu d'une interpolation par spline cubique. |
1858 |
\item[L10] Nouveau cas \verb+strato2+ dans \verb+disvert+. |
1859 |
Correction de bogue dans \verb+inter_bary+ : \verb+y0+ initialisé à |
1860 |
$-90$ au lieu de 0. |
1861 |
\item[L11] Remaillage de la variable \verb+a2+ dans \verb+regr_coefoz+. |
1862 |
\item[L12] Changement de nom : \verb+sigma+ en \verb+s+ dans |
1863 |
\verb+disvert+. Changement de noms : \verb+r+ et \verb+P_net+ en |
1864 |
\verb+r_Mob+ et \verb+P_net_Mob+ dans \verb+coefoz_LMDZ.nc+. |
1865 |
\item[L13] Remaillage de \verb+a4+ et \verb+temperature+ dans |
1866 |
\verb+regr_coefoz+. |
1867 |
\item[L14] Remaillage de \verb+a6+ et \verb+Sigma+ dans |
1868 |
\verb+regr_coefoz+. |
1869 |
\item[L15] Remaillage de \verb+R_Het+ dans \verb+regr_coefoz+. |
1870 |
Suppression de la variable \verb+newlmt+ dans \verb+limit+. |
1871 |
\end{description} |
1872 |
|
1873 |
\section{Versions de \texttt{gcm}} |
1874 |
|
1875 |
Modifications provisoires : |
1876 |
\begin{description} |
1877 |
\item[M1] Modifications de François L. dans \verb+ini_histday.h+, |
1878 |
\verb+physiq.F+, \verb+phytrac.F+, \verb+write_histday+\verb+.h+. |
1879 |
Avec ces modifications, les fichiers \verb+ini_histhf.h+, |
1880 |
\verb+write_histhf.h+, \verb+ini_histmth.h+, \verb+write_histmth.h+, |
1881 |
\verb+ini_histmthNMC.h+ et \verb+write_histmthNMC.h+ ne sont plus |
1882 |
utilisés. Les variables \verb+lev_histhf+, \verb+lev_histmth+ et |
1883 |
\verb+ok_mensuel+ ne sont donc pas utilisées. |
1884 |
\item[M2] Dans \verb+gcm+, après l'appel à \verb+dynetat0+ : |
1885 |
\verb+teta+ moyenné horizontalement, \verb+ucov=vcov=0+, |
1886 |
\verb+q=0+. |
1887 |
\item[M3] Dans \verb+gcm+, après l'appel à \verb+dynetat0+ : |
1888 |
\verb+ps+ = 101325. |
1889 |
\item[M4] Dans \verb+leapfrog.F+, mise en commentaire des appels |
1890 |
à \verb+writedynav+ et \verb+bilan_dyn+ (deux appels |
1891 |
de chaque sous-programme). Pour que le fichier \verb+dyn_hist_ave.nc+ |
1892 |
soit presque vide et que le fichier \verb+dynzon.nc+ ne |
1893 |
soit pas créé. |
1894 |
\end{description} |
1895 |
|
1896 |
Modifications cumulatives : |
1897 |
\begin{description} |
1898 |
\item[L1] Sur la base de la version téléchargée le 2/3/6, corrections |
1899 |
de Phu L. V. dans \verb+physiq.F+ et \verb+write_histmthNMC.h+, \textit{cf.} |
1900 |
messages des 13/3/6 et 14/3/6. |
1901 |
\item[L2] Décroissance radioactive de l'ozone (modification dans |
1902 |
\verb+phytrac+). |
1903 |
\item[L3] Correction : attribut \verb+save+ dans les modules |
1904 |
\verb+comvert+, \verb+conema3_m+, \verb+serre+, \verb+clesphys+, |
1905 |
\verb+iniprint+, \verb+ener+, \verb+comgeomphy+, \verb+yoethf+, |
1906 |
\verb+comdissnew+, \verb+dimphy+, \verb+logic+, \verb+comgeom+. |
1907 |
\item[L4] \'Evolution du traceur ozone selon la paramétrisation de |
1908 |
Cariolle, à l'ordre 0 (terme $P_\mathrm{net}$ uniquement) et sans la |
1909 |
chimie hétérogène. |
1910 |
\item[L5] Correction dans \verb+phytrac+, pour empêcher la fraction |
1911 |
massique d'ozone de devenir négative. Nouveaux cas \verb+strato1+ et |
1912 |
\verb+strato2+ dans \verb+disvert+. |
1913 |
\item[L6] Augmentation de \verb+lmx+ de 20 à 100 dans \verb+gradsdef+. |
1914 |
\item[L7] Ajout du terme en $a_2$ dans le calcul de la production |
1915 |
chimique d'ozone, dans \verb+phytrac+. |
1916 |
\item[L8] Ajout du terme en $a_4$ dans le calcul de la production |
1917 |
chimique d'ozone, dans \verb+phytrac+. |
1918 |
\item[L9] Ajout du terme en $a_6$ dans le calcul de la production |
1919 |
chimique d'ozone, dans \verb+phytrac+. |
1920 |
\item[L10] Ajout du terme de chimie hétérogène dans le calcul de la |
1921 |
production chimique d'ozone, dans \verb+o3_chem_m+. |
1922 |
\item[L11] Correction d'un bogue dans \verb+o3_chem+ : le filtrage |
1923 |
selon la direction du Soleil doit affecter $b$ et non \verb+r_het+. |
1924 |
\end{description} |
1925 |
|
1926 |
\section{Utilisation de module à la place de l'inclusion de fichier} |
1927 |
|
1928 |
Pour remplacer les \verb+include+ de fichiers particuliers, |
1929 |
par exemple \verb+temp.h+ et \verb+temps.inc+, par |
1930 |
des \verb+use temps+ dans tous les fichiers, dans chacun |
1931 |
des répertoires \verb+bibio+, \verb+dyn3d+, \verb+filtrez+ |
1932 |
et \verb+phylmd+ : |
1933 |
\begin{verbatim} |
1934 |
%M |
1935 |
.+\.f\(90\)?$ |
1936 |
Q |
1937 |
include +\(['"]\)temps\.\(h\|inc\)\1 |
1938 |
use temps |
1939 |
\end{verbatim} |
1940 |
Pour chaque fichier, taper ensuite ``!''. Enregistrer tous |
1941 |
les fichiers modifiés avec la commande \verb+Ctrl-x s !+. |
1942 |
Tous les fichiers modifiés étant encore ouverts, dans chacun |
1943 |
d'entre eux, déplacer la ligne \verb+use+ au bon endroit, |
1944 |
enregistrer et fermer le \Eng{buffer}. |
1945 |
|
1946 |
Effacer le fichier \verb+.h+. Ajouter le fichier \verb+.f+ dans la |
1947 |
liste des objets pour l'un ou l'autre exécutable (selon qu'ils |
1948 |
l'utilisent ou non). |
1949 |
\begin{verbatim} |
1950 |
make depend |
1951 |
make.sh |
1952 |
\end{verbatim} |
1953 |
Corriger les mauvaises positions restantes des instructions \verb+use+. |
1954 |
|
1955 |
\end{document} |