Un point sur LMDZ en prévision des simulations CMIP du GIEC

Equipe développement LMDZ

9 juillet 2009

Table des matières

Les configurations

Deux versions de "la physique"

La stratégie de la participation au simulations CMIP5 (pour le rapport 5 du Giec) avec le modèle couplé de l'IPSL est basée sur l'utilisation de deux versions très différentes de la "physique" de LMDZ :

une version proche de celle utilisée pour l'AR4 :: cette version doit être figée courant juillet 2009 pour lancer en août les simulations de contrôle.
la version dite "nouvelle physique" :: cette version est encore en cours de rodage et sera figée en décembre 2009.

Deux grilles horizontales

Le modèle sera décliné dans deux configurations horizontales :

96x95: points en longitudes x latitudes (ou 3,75x1,875 en degrés)
144x142: points en longitudes x latitudes (ou 2,5x1,27 en degrés)

On peut se reporter pour regarder des tests de résolution horizonale à une page de texte simulations réalisées par Marie-Alice Foujols

Une nouvelle discrétisation verticale

Le modèle sera principalement décliné sous une forme L39 avec 39 niveaux sur la verticale. Par rapport à l'ancienne configuration à 19 niveaux, les 20 niveaux additionnels sont répartis en gros entre 6 niveaux "au-dessus" de l'ancien modèle, un niveau plus près du sol (33m) et une résolution plus fine d'un facteur environ 1.7 dans la masse de l'atmosphère.

Ce passage à 39 niveaux a nécessité un certain travail d'ajustement sur le modèle qui avait tendance à devenir instable. L'activation de certaines paramétrisations de la version stratosphérique (Lott et al., 2005) avec certaines adaptations a permis de résoudre ces problèmes.

Les otuils pour stabiliser le modèle consistent à

activer la couche éponge qui absorbe les ondes dans les couches les plus hautes du modèle,
augmenter la dissipation horizontale dans les couches les plus hautes,
activer la paramétrisation de Hines pour le freinage de la stratopshère par les ondes de gravité ("ok_hines=y" dans "physiq.def").

Les coefficients de relaxation de la couche éponge et la dépendance verticale des coefficients de la dissipation horizontale ont été revus à cette occasion comme on l'explique plus loin.

Divers changements

A noter un certain nombre d'évolutions depuis la version AR4. En particulier,

La parallélisation MIP/OpenMP.
L'automatisation de la gestion des IO avec des flags pour chaque variable et pour les différents fichiers.
Le passage au calendrier vrai, ou 365 jours qui sera utilisé pour les simulations CMIP5
La réécriture de l'interface avec les surfaces
La prise en compte du forçage par l'ozone et les aérosols.
La récriture de l'interface traceurs.

Ancienne physique : LES DERNIERS REGLAGES

Des séries de tests sont en cours pour régler le modèle avec l'ancienne physique. C'est le cas notamment de la série de simulations réalisées par Sébastien Denvil pour laquelle on trouvera un tableau d'acces aux atlas et différents diagnostics

Les diagnostics, effectués avec les atlas de Patrick ont été enrichis récemment avec de nouvelles climatologies et nouveaux diagnostics par Ionela Musat et Abderrahmane Idelkadi.

Laurent Fairhead a également lancé un test en couplé

Calcul des niveaux verticaux du modèle

Calcul des épaisseurs en niveaux $\sigma$

La discrétisation verticale des version standard du modèle est préscrite en calculant d'abord des niveaux de pression normalisée, $\sigma_l$ où $\sigma=p_l/p_s$ . On commence pour cela par calculer des épaisseurs de couches à partir d'une fonction analytique :

$\begin{displaymath} \delta {\sigma_l}^*=\left[ \alpha + 7 \left( \sin Z \right)^2 \right] \end{displaymath}$

(1)

avec

$\begin{displaymath} Z=\pi \frac{l-1/2}{lm+1} \end{displaymath}$

(2)

qu'on normalise pour déduire $\sigma_l$ :

$\displaystyle \delta\sigma_l$	$\textstyle =$	$\displaystyle {\delta {\sigma_l}^}/{\sum_{l=1}^{l=lm}\delta {\sigma_l}^}$	(3)
$\displaystyle \sigma_L$	$\textstyle =$	$\displaystyle 1-\sum_{l=1}^{l=L}\delta\sigma_l$	(4)

Ainsi définie, l'épaisseur des couches en pression varie d'un facteur $7/\alpha$ entre les bords (en haut et à la surface) et le milieu de l'atmosphère, $Z=\pi/2$ . Dans la version initiale, $\alpha =1$ .

Passage aux coordonnées hybrides

A partir de là, on calcule les coefficients et des coordonnées hybrides qui donnent la pression à un point de pression de surface et au niveau sous la forme . Les coefficients sont calculés selon les formules :

$\displaystyle B_l$	$\textstyle =$	$\displaystyle \exp\left[ 1 -\frac{1}{{\sigma_l}^2} \right]$	(5)
$\displaystyle A_l$	$\textstyle =$	$\displaystyle P_{\mbox{hyb}}\left( \sigma_l - B_l \right)$	(6)

Le lien entre les niveaux

et $\sigma_l$ est, avec ce choix, $P_l=B_l (P_s -P_{\mbox{hyb}}) + P_{\mbox{hyb}}\sigma_l$ .

On choisit pour la terre $P_{\mbox{hyb}}=500$ hPa=.

Le cas de la discrétisation avec stratosphère

Pour l'extension à la stratosphère, l'Eq. 1 est modifiée selon

$\begin{displaymath} \delta {\sigma_l}^*=\left[ \alpha + 7 \left( \sin Z \right)... ...1-\tanh\left( \left( Z-\pi \right)/{\pi} \right)}{2} \right]^2 \end{displaymath}$

(7)

Configurations standards

**Figure 1:** Epaisseur en altitude (à gauche, pour une hauteur d'échelle de 7 km) et en pression (à droite) des couches pour 3 discrétisations (L=19, 39 et 50), soit avec "ok_strato=n" (en utilisant l'Eq. 1), soit avec "ok_strato=y" (Eq. 7), et avec $\alpha =1$ ou 0,3 pour la version à 39 niveaux.
$\includegraphics[width=6cm]{vertz.eps}$ $\includegraphics[width=6cm]{vertp.eps}$

Le facteur $\alpha$ a été introduit (et sorti dans les fichiers .def) pour contrôler l'épaisseur de la première couche quand on est passé de la version non troposphérique à 19 couches utilisée pour l'exercice CMIP4 pour le rapport AR4 du GIEC à la version à 39 niveaux avec stratosphère utilisée pour CMIP5, pour la quelle on a fixée $\alpha=0,3$ .

On compare sur la Fig. 1 ces deux configurations de référence, la version stratosphérique développé par Lott et al. (2005) ainsi qu'une version avec stratosphère à 39 niveaux pour la quelle on aurait fixé $\alpha =1$ .

Dans le code, le calcul de la discrétisation verticale est effectué dans

dyn3d/disvert.F

Le passage de l'Eq. 1 à l'Eq. 7 se fait en passant la variable "ok_strato" de la valeur "n" à "y". En revanche pour le moment (juillet 2009), et en attendant un choix définitif sur la façon de le faire, le paramètre $\alpha$ n'est pas contrôlable par les fichiers ".def". En pratique, $\alpha =1$ pour tous les cas sauf pour "llm=39" et "ok_strato=y" où on impose $\alpha$ =0.3".

Les paramètres de contrôle de la dissipation temporelle dans le modèle

On a applique trois opérateurs de dissipation différents, l'un sur la température (assimilé à la lettre pour enthalpie), l'un pour la divergence et le dernier pour le rotationnel du vent. On peut fixer dans "gcm.def" d'abord part les nombre d'itération avec les choix conseillés :

## nombre d'iterations de l'operateur de dissipation   gradiv
nitergdiv=1
## nombre d'iterations de l'operateur de dissipation  nxgradrot
nitergrot=2
## nombre d'iterations de l'operateur de dissipation  divgrad
niterh=2

Les constantes de temps sont également contrôlées dans "gcm.def"

## temps (s) de dissipation des plus petites long.d ondes pour u,v (gradiv)
tetagdiv=7200.
## temps (s) de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
tetagrot=7200.
## temps (s) de dissipation des plus petites long.d ondes pour  h ( divgrad)
tetatemp=7200.

Ces constantes sont toujours de l'ordre de quelques dizaines de minutes à quelques heures dans les configurations classiques du modèle, mais dépendent davantage de la configuration (en général, on prend des constantes de temps un peu plus courte pour des grilles plus fines).

En fait, ces constantes de temps sont ensuite affectées d'un facteur qui dépend de la verticale afin de rendre la diffusion plus efficace dans la haute atmosphère.

Dans les versions d'origine du modèle, ce facteur était calculé comme

$\begin{displaymath} f_l=\alpha -\frac{\alpha -1}{1+Z^2} \end{displaymath}$

(8)

avec

$\begin{displaymath} Z=1-\frac{P_0}{P_l} \end{displaymath}$

(9)

où

est la pression approximative au milieu des couches.

Dans la version originale, $\alpha =2$ .

Dans la nouvelle version, la formulation est changée pour

$\begin{displaymath} f_l=\frac{1}{2}\left[ \tanh{\frac{{z^*}_l-Z_d}{\delta z}}+1 \right] \times (\alpha -1) \end{displaymath}$

(10)

où ${z^*}_l=8.\ln{P_0/P_l}$ est la pseudo altitude du niveau

calculée pour une hauteur d'échelle de 8 km (ce serait mieux de passer à 7 et surtout de pouvoir contrôler cette valeur. A voir au moment de la convergence avec les versions planeto).

Pour le moment (juillet 2009), pour éviter de trop brouiller les cartes et dans l'attente de généralisation dans le cadre de la convergence avec les versions planétaires, cette formule n'est active que pour "ok_strato=y" et "llm=39" (configuration de référence pour l'exercice CMIP5).

Dans ce cas, les paramètres de la fonctions sont contrôlables dans "gcm.def" grâce aux paramètres $\alpha$ , $\delta z$ et , respectivement

dissip_factz=4.
dissip_deltaz=10.
dissip_zref=30.

avec leurs valeurs par défaut pour CMIP5. Avec ces valeurs, l'idée est qu'on dissipe 4 fois plus efficacement au-dessus de 30 km qu'à la surface, avec une transition sur 10 km.

Pour la stabilité numérique du schéma temporel explicite utilisé pour les opérateurs de dissipation, la fréquence à laquelle est appelée la dissipation est calculée automatiquement et est un multiple de "iperiod", fréquence à laquelle on effectue les pas "matsuno" dans le modèle. En général "iperiod=5". Le pas de temps d'appel à la dissipation ne peut donc descendre en-dessous de "iperiod $\times$ dtvr", où "dtvr" est le pas de temps dynamique du modèle.

D'après la discussion ci-dessus, on voit qu'il faut que la constante de temps de dissipation soit plus grande que ce pas de temps minimum d'appel à la dissipation. En tentant compte du facteur sur la verticale, cette contrainte s'écrit

$\begin{displaymath} \mbox{min(tetatemp,tetagdiv,tetagrot)}/\alpha <\mbox{iperiod}\times \mbox{dtvr} \end{displaymath}$

(11)

La couche en éponge

Description

La "couche en éponge" ou "sponge layer" en anglais, a pour but d'absorber les ondes au sommet du modèle. Les ondes qui se propagent vers le haut voient leur amplitude croître sous l'effet de la diminution de la masse volumique de l'air. Dans un modèle numérique limité en pression (à ), les ondes sont artificiellement réfléchies sur le sommet du modèle.

Dans la couche éponge du LMDZ, appelée "top_bound", on fait décroître la partie non-zonale de l'écoulement avec une constante de temps $\tau$ . Si on note $\left[ X \right]$ la moyenne zonale d'un champ et $X^*=X-\left[ X \right]$ la perturbation associée, on écrit simplement

$\begin{displaymath} \frac{\partial X}{\partial t}_{\vert\mbox{eponge}}=-\nu X^* \end{displaymath}$

(12)

avec $\nu =1/\tau$ . Les constantes de temps sont d'autant plus grande qu'on s'approche du sommet du modèle.

Choix des constantes de temps

Dans la version stratosphérique du modèle à 50 niveaux développée par Lott et al. (2005), la couche éponge est active seulement dans les 4 dernières couches avec $\nu _{lm}=1/\tau_0$ , $\nu _{lm-1}=2/\tau_0$ , $\nu _{lm-2}=4/\nu _0$ , $\nu _{lm-3}=8/\tau_0$ et $\nu _l=0$ pour . La valeur de $\tau_0$ était fixée à 10 s.

Cette formule a été un peu assouplie (des instabilités numériques apparaissaient juste sous la dernière couche affectée pour la couche éponge) en écrivant

$\begin{displaymath} \nu _l= \max(P_{lm}/P_l-0.01,0.)/\tau_0 \end{displaymath}$

(13)

La constante de temps croît ainsi linéairement avec la pression en partant de la dernière couche du modèle vers le bas, juqu'à ce quelle soit 100 fois la valeur dans la couche la plus hautes.

Paramètres de contrôle

La couche éponge est contrôlée par un choix sur le calcul des constantes

# iflag_top_bound : 0 : pas d'action, 1: actif dans les 4 dernieres couches, 2 : nouveau
iflag_tob_bound=2
# Inverse de la constante de temps de relaxation dans la couche la plus haure du modele.
tau_top_bound=1.e-5

Rm : peut-être que ce serait mieux de mettre la constante en secondes ...

À propos de ce document...

Un point sur LMDZ en prévision des simulations CMIP du GIEC

This document was generated using the LaTeX2HTML translator Version 2002-2-1 (1.71)

The command line arguments were:
latex2html -split 0 configurations_cmip

The translation was initiated by hourdin on 2009-07-09

hourdin 2009-07-09