mercredi 30 juillet 2014

Mise en oeuvre des modèles météo WRF-ARW et WRF-NMM - Partie 9 - Aller plus loin

Nous avons fait un long parcours jusqu'ici pour mettre en oeuvre le modèle WRF, que ce soit en version ARW ou NMM. Jusqu'à maintenant, nous n'avons fait qu'utiliser le modèle "tel quel", sans chercher à comprendre davantage son fonctionnement, comment l'optimiser, et améliorer la qualité de la simulation. Cet ultime article, comme promis précédemment, sera une discussion permettant de donner les pistes nécessaires à une mise en exploitation plus sérieuse du modèle.

Les paramétrisations physiques


Nous allons dans un premier temps rechercher comment on peut obtenir la simulation la plus réaliste possible. Pour cela, il faut bien comprendre comment fonctionne un modèle météo, et avoir un aperçu des phénomènes qu'il faut prendre en compte.

Un modèle découpe l'atmosphère en une grille tridimensionnelle plus ou moins fine, et considère que chaque point de grille est un volume d'air dont les propriétés (température, pression, vent...) sont homogènes. L'atmosphère est en quelque sorte discrétisée, et la connaissance que l'on en a, qui provient des observations ou d'un autre modèle comme c'est le cas dans nos articles, est reportée dans cette grille de façon physiquement cohérente. On appelle cette grille initiale une analyse.

Le modèle va calculer l'évolution de cette analyse en appliquant les équations de la dynamique des fluides. Sans rentrer dans les termes barbares des équations, celles-ci décrivent :
  • le bilan des forces physiques : pesanteur, pression,... et leur effet sur l'accélération de l'air et donc la création du vent
  • la prise en compte de la force de coriolis, déviation du vent dûe à la rotation de la terre
  • la conservation de la masse : au cours du mouvement, aucune quantité de matière ne se rajoute ou disparait.
  • le transport d'eau, que ce soit en vapeur, liquide ou solide.
  • la conservation du tourbillon, c'est-à-dire la rotation des volumes d'air sur eux-même (dûe entre autres à la rotation de la terre).
En résolvant ces équations à chaque pas de temps, on en déduit les mouvements de l'atmosphère et donc son évolution. Mais cette description ne prend pas en compte tous les phénomènes. En effet, la taille de maille est généralement trop importante pour résoudre des phénomènes d'échelle aérologique ou de micro-échelle, comme les cumulus de beau temps ou la turbulence. D'autres phénomènes sont trop complexes pour être intégrés directement dans les équations dynamiques. On utilise des modèles simplifiés pour chaque phénomène physique, afin de calculer les effets moyens sur la maille en tenant compte de façon plus ou moins réaliste des boucles de rétroaction (rayonnement vs couverture nuageuse...). Voici une liste non exhaustive des principaux phénomènes physiques paramétrisés :
  • apports de chaleur par le rayonnement du soleil, et à l'inverse la perte nocturne par rayonnement infrarouge.
  • la condensation de la vapeur d'eau, ou d'une manière plus générale la formation des nuages. Cela génère des précipitation, et libère une grande quantité de chaleur latente qu'il faut prendre en compte dans les calculs. Il existe également une boucle de rétroaction avec le rayonnement solaire.
  • l'influence de la surface et de son type (urbain, forêt, champs...), qui freine le vent, réchauffe l'atmosphère par la base, et génère des turbulences (couche limite atmosphérique, couche limite de surface)
  • le sol et les océans sont également pris en compte. On calcule les échanges de chaleur, en prenant en compte plusieurs couches de sol et d'océan, et on modélise également de manière simplifiée l'hydrologie pour prendre en compte la saturation des sols et l'écoulement des eaux de pluie, car cela modifie l'interaction sol/atmosphère (évapotranspiration, etc...).
  • l'influence du relief qui bloque ou dévie les vents, provoque d'importants forçages verticaux
Dans WRF, il est possible de choisir entre plusieurs modèles pour différents paramètres physiques. Nous n'avons fait jusqu'ici qu'utiliser les modèles standard, qui fonctionnent certes très bien, mais qui ne donnent pas forcément le résultat escompté dans certaines situations. C'est le cas par exemple du modèle de cumulus. Notre fichier namelist.input comportait la ligne suivante :

cu_physics = 2,

Cela signifie que l'on utilise le modèle Betts Miller Janjic, un modèle qui ne traite pas explicitement les courants ascendants/descendants. En simulant dans des situations orageuses, on comprend assez vite que ce modèle rend compte assez mal des phénomènes convectifs. La part convective des précipitations produites est peu importante, et en-dessous de la réalité. Il faut donc trouver un modèle plus élaboré.

Il existe une dizaine de schémas de cumulus dans WRF, mais ils ne fonctionnent pas forcément avec tous les coeurs. J'en ai testé trois qui fonctionnent avec NMM :
  • 2 : Bretts Milles Janjic, qui vient par défaut. Il ne gère qu'un seul type de particule de nuages.
  • 14 : New SAS (Simplified Arakawa Schubert). Il gère la glace de nuage ainsi que les courants descendants.
  • 1 : Kain-Fritsch. Il gère la glace, la pluie et la neige, les courants descendants et ascendants.

Le type 2 donne un résultat assez peu réaliste en présence de convection profonde. Le type 14 donne un bien meilleur résultat, mais il apparaît beaucoup trop lissé, les cellules semblent isolées les unes des autres. Ce schéma semble ne pas gérer forcément bien la convection peu profonde. Le meilleur résultat selon moi est obtenu avec le type 1, d'une part parce qu'il est plus complet, et d'autre part parce qu'il simule de façon satisfaisante les différents niveaux de convection.

On voit donc l'influence que peut avoir le choix d'une paramétrisation physique sur la qualité de simulation. Sans faire de test chiffré, je n'ai pas vu de différence de temps de calcul significative avec ces différents modèles. Le choix peut donc se porter sur le réalisme de la simulation.

Je ne vais pas faire de discussion longue sur tous les types de schéma des autres phénomènes physiques, je laisse le lecteur expérimenter sur le sujet en utilisant les documentations du modèle.


Augmenter la résolution ou réduire les temps de calcul



Il est tentant pour augmenter la qualité de simulation d'augmenter la résolution du domaine. Avec WRF, il est possible de descendre sous les 5 km, mais il faut savoir que tout doublement de la taille horizontale de grille se traduit par une augmentation du temps de calcul par 4. De même, un doublement du nombre de niveaux verticaux entrainera un doublement du temps de calcul. En réalité, c'est même pire que cela : le pas de temps devra être réduit de moitié pour avoir un modèle numériquement stable, ce qui augmente encore d'un facteur 2 le nombre de calculs à effectuer.

Si l'on prend une grille de 5km, sur 30 niveaux, sur la France, soit une résolution de 256x256x30, il faut à un Intel core i7 3770k environ 1h30 pour calculer 48h de simulation. Un doublement de la résolution horizontale nécessiterais de traiter 4 fois plus de points de maille, et de le faire pour 2 fois plus de pas de temps en raison de la condition de stabilité. Le temps de calcul serait donc multiplié par 8 !

On comprend alors très vite qu'armé d'un simple PC, même puissant, il va devenir compliqué de faire ce genre de simulations dans un temps raisonnable. Pour augmenter la puissance de calcul, il n'y a pas 36 solutions : il faut plus de processeurs. Quand vous vous appelez Météo France, vous disposez d'un monstre équipé de milliers de processeurs capable de vous faire cela avant même que vous n'ayez eu le temps de vous servir un café. Mais tout le monde n'a pas le budget d'un centre de calcul... Une autre solution est de mettre en cluster plusieurs PC et de les faire calculer en parallèle.

Cela nécessite bien entendu un investissement matériel qui peut être conséquent, mais avec des processeurs plus grand public on peut arriver à de très bon résultats sans trop exploser le budget. Dans ce tutorial, on est quasiment prêt pour la mise en cluster. En effet, si vous avez compilé WRF avec l'option MPI, vous avez quasiment tout ce qu'il vous faut. Pour mettre en cluster, l'idée est donc la suivante :
  • Acheter et assembler plusieurs machines identiques, les mettre en réseau gigabit et prendre un bon abonnement chez EDF
  • Installer linux sur chacune des machines en respectant les pré-requis pour WRF,
  • installer NFS sur chaque machine. Partager votre dossier météo sur la machine maître pour toutes les autres.
  • Mettre en place la communication inter-machines (SSH, fichier host...)
  • Compiler WRF sur la machine maître (si vous avez suivi ce tutorial, c'est déjà fait), dans le dossier partagé sur NFS
  • Installer MPICH et configurer le fichier machine
  • Lancer WRF depuis la machine maître avec la commande MPI qui va bien...
Tout ceci ne représente que les grandes lignes, et n'a pas été testé par l'auteur de cet article (qui n'a pas de datacenter chez lui). Une procédure plus détaillée pour Ubuntu est accessible dans la documentation officielle de la distribution, pour les lecteurs qui seraient tentés d'essayer.


Se concentrer sur les régions d'intérêt


Pour des raisons de budget, il n'est pas toujours envisageable de monter un cluster de machines chez soi. Pour limiter la puissance de calcul nécessaire, on peut limiter la simulation à des domaines d'intérêt plus petits. Mais un domaine trop petit ne donne pas forcément de bonnes simulations, en raison de la vitesse de propagation de certains phénomènes. Un vent de 100km/h traverserait en 2h un domaine de 200km, ce qui peut poser des problèmes de précision aux limites car par exemple avec GFS on ne dispose que d'un fichier toutes les 3 heures.

On partirais plutôt dans l'optique de baisser la résolution sur notre domaine (la france entière) : une diminution par 3 divise en gros notre temps de calcul par 9. Dans notre cas, on passe de 0.05 degrés à 0.15 degrés, ce qui est grosso-modo trois fois mieux que GFS 0.5 degrés. Notre grille passe de 256x256 à 84x84, et  au lieu de 1H30 de calcul on passe à 10-15 minutes. On utilise cette solution rapide pour initialiser ensuite un domaine plus petit mais avec une résolution 3 fois plus fine, centré sur notre région d'intérêt. On prend une grille de 0.05 dgrés, faisant seulement 128x128 ce qui divise malgré tout le temps de calcul initial par 4 soit environ 25 minutes. Votre temps de calcul pour l'ensemble serait ramené à environ 40-45 minutes.

WRF permet de faire cela très facilement, il suffit pour cela de le recompiler avec l'option nesting. Reprenez pour cela la procédure de la partie 7 de ce tutorial, mais déclarez une nouvelle variable d'environnement WRF_NMM_NEST=1 avant de compiler. Ensuite, il faut déclarer les différents domaines dans les fichiers namelist.wps et namelist.input. Inspirez-vous du fichier d'origine fourni avec les sources, qui utilise cette fonctionnalité. A la sortie, WRF produira autant de fichiers que de domaines existants. Il faudra en tenir compte dans le script Unipost et dans les scripts de génération de cartes.

Exemple de nesting 0.15° et 0.05°C sur le champ de précipitations.
Exemple de nesting 0.15° (à gauche) et 0.05°C (à droite) sur le champ de précipitations.
Voici les portions importantes du fichier namelist.wps permettant de faire cela :

&share
 wrf_core = 'NMM',
 max_dom = 2,

&geogrid
 parent_id         =   1,    1
 parent_grid_ratio =   1,    3
 i_parent_start    =   1,    10   
 j_parent_start    =   1,    26
 e_we              =  42,     64
 e_sn              =  84,    128
 geog_data_res     = '5m', '2m'
 dx = 0.15,
 dy = 0.15,

Les variables à modifier sont :
  • max_dom : indique le nombre de domaines
  • parent_id : permet lier hiérarchiquement les domaines entre eux.
  •  i_parent_start, j_parent_start : indices de lignes et de colonnes permettant de positionner la grille "détaillée" dans la grille "large".
  • e_we, e_sn, geog_data_res : des paramètres connus qu'il faut définir pour chaque grille
  • dx, dy : la résolution de la grille la plus large.

Et pareil pour le fichier namelist.input :

&time_control
 start_year                          = 2014,   2014
 start_month                         = 06,     06,
 start_day                           = 30,     30,
 start_hour                          = 12,     12,
 start_minute                        = 00,      00,
 start_second                        = 00,       00,
 end_year                            = 2014,     2014,
 end_month                           = 07,        07,
 end_day                             = 02,        02,
 end_hour                            = 12,        12,
 end_minute                          = 00,        00,
 end_second                          = 00,        00,
 

...
 history_interval                    = 60,        60,
 frames_per_outfile                  = 1000,      1000,
...
 &domains
 e_we                                = 42,        64,
 e_sn                                = 84,        128,
 e_vert                              = 30,        30,
...

 dx                                  = 0.15,    0.05,
 dy                                  = 0.15,    0.05,
 grid_id                             = 1,        2,
 parent_id                           = 0,        1,
 i_parent_start                      = 1,        10,
 j_parent_start                      = 1,       26,
 parent_grid_ratio                   = 1,        3,
 parent_time_step_ratio              = 1,        3,
...

 &physics
 mp_physics                          = 5,        5,
 ra_lw_physics                       = 99,        99,
 ra_sw_physics                       = 99,        99,
 radt                      = 15,  5,
 nrads                               = 105,     105,
 nradl                               = 105,     105,
 co2tf                               = 1,
 sf_sfclay_physics                   = 2,        2,
 sf_surface_physics                  = 2,        2,
 bl_pbl_physics                      = 2,        2,
 nphs                                = 6,
 cu_physics                          = 1,        1,
 ncnvc                               = 6,
 tprec                               = 3,         3,
 theat                               = 6,         6,
 tclod                               = 6,         6,
 trdsw                               = 6,         6,
 trdlw                               = 6,         6,
 tsrfc                               = 6,         6,
 pcpflg                              = .false., .false.,
...
 &bdy_control
 spec_bdy_width                      = 1,
 spec_zone                           = 1,
 relax_zone                          = 4,
 specified                           = .true.,     .false,
 nested                              = .false., .true.,
 /
...



Sans rentrer dans tous les détails, on définit pour chaque grille les différents paramètres de taille de grille et de paramétrisations physiques. En principe, celles-ci sont identiques pour chaque domaine. Mais si vous descendez sous la barre des 3km sur un domaine, il peut être judicieux de se passer des modèles de cumulus car le modèle NMM simule alors explicitement la convection. Le domaine correspondant pourra avoir le cu_physics à 0.

L'économie de temps de calcul est substantielle, puisque je met deux fois moins de temps pour avoir une simulation sur ma région d'intérêt, sans perdre du réalisme de la simulation. Dans certains cas, j'ai même constaté une meilleure précision qu'avec un modèle qui part directement de GFS sur toute la France. En effet, le "nest" se trouve initialisé avec des données aux limites affinées par rapport au GFS d'origine, les forçages par exemple se trouvent mieux rendus.


Conclusion


Nous arrivons au terme de cette serie de tutoriaux sur le modèle WRF. Il y aurait certainement encore beaucoup de choses à dire, ce qui fera peut-être l'objet d'articles complémentaires futurs au gré de mes expérimentations. En espérant que cette série vous aura intéressé, que vous aurez trouvé les informations techniques utiles. Et surtout que vous prendrez plaisir autant que moi à expérimenter vous aussi vos propres simulations avec WRF.