mercredi 14 août 2013

Mise en oeuvre des modèles météo WRF-ARW et WRF-NMM - Partie 5 - Notre premier run WRF

Dans l'article Partie 4 - Compilation de WRF-ARW et du WPS vous avez compilé le logiciel WRF et son pré-processeur. Nous allons maintenant apprendre comment réaliser une simulation du temps sur 24h pour la France avec WRF-ARW, en se basant sur un run du modèle GFS pour l'initialisation. La simulation aura une résolution horizontale de 10km.

Obtenir des données


Vous pourrez obtenir un jeu de données de test sur la page suivante (GFS GRIB 2 Data) :

http://www.mmm.ucar.edu/wrf/OnLineTutorial/DATA/GFS/

Vous pourrez également utiliser les données opérationnelles de GFS.

Créez un dossier DATA dans lequel vous téléchargerez les fichiers GRIB. Pour les données opérationnelle GFS, prenez bien les fichiers "pgrb2f", à ne pas confondre avec les "pgrb2bf" qui sont similaires mais ne contiennent pas toutes les données. Prenez bien les fichiers pour toutes les échéances requises (jusqu'à +24 pour 24h de simulation et ainsi de suite suivant vos besoins), et surtout n'oubliez pas le +00 qui correspond à l'analyse, indispensable. Nous ferons pour notre part une simulation sur 24h, ce qui correspond aux fichiers suivants :

gfs.t00z.pgrb2f00 gfs.t00z.pgrb2f03 gfs.t00z.pgrb2f06
gfs.t00z.pgrb2f09 gfs.t00z.pgrb2f12 gfs.t00z.pgrb2f15
gfs.t00z.pgrb2f18 gfs.t00z.pgrb2f21 gfs.t00z.pgrb2f24


Une fois les fichiers téléchargés, allez dans le dossier WPS et tapez la commande :

# ./link_grib.csh ../DATA

Ceci va créer des liens symboliques vers tous vos fichiers de données (attention à ne pas mettre autre chose que des fichiers GRIB dans le dossier DATA, sinon ils seront liés aussi...) :

GRIBFILE.AAA
GRIBFILE.AAB
etc...

Déterminer le domaine


Pour simplifier, nous utiliserons un seul domaine qui couvrira grossièrement la France, avec une grille de 100x100 pour une résolution horizontale de 10km. Ces paramétrages se font à l'aide du fichier namelist.wps, à éditer avec votre éditeur préféré. Je vous recommande vivement de faire une copie du fichier par défaut, pour le cas où voudriez suivre les exemples du tutoriel officiel ou pour chercher ce que vous auriez pu changer qui fait planter le programme (car on ne paramètre pas n'importe quoi n'importe comment, vous vous en rendrez vite compte !).

Je vous fournis à la fin de cette section mon fichier namelist pour mes fichiers d'exemple. Je décortique ici les différentes variables qui nous intéressent :
  • le type de coeur dynamique utilisé (variable wrf_core) : ARW ou NMM. Ne cherchons pas les ennuis, utilisons pour le moment ARW.
  • les dates de simulation (variables start_date et end_date). Il est impératif que vos fichiers de données couvrent tout cet intervalle.
  • l'intervalle entre deux fichiers de données. Pour GFS, c'est 3h soit 10800 secondes (variable interval_seconds)
  • debug_level : mettez cette valeur à 1000 afin d'avoir le plus de messages possibles, qui vous seront utiles pour résoudre tout problème éventuel.
  • une projection de la grille de type Lambert conforme. (variable map_proj). D'autres valeurs sont possibles, suivant le domaine que vous voulez couvrir et le type de carte que vous voudrez obtenir (polaire, etc...).
  • les paramètres de la projection (variables truelat1, truelat2, stand_lon). Les paramètres indiqués permettent d'avoir une bonne projection pour la France.
  • les variables ref_lat et ref_lon définissent habituellement la latitude et la longitude du centre de notre domaine. Les coordonnées indiquées dans le fichier correspondent au centre géographique officiel de la France.
  • les variables dx et dy définissent la taille de grille. Pour ARW, spécifiez une valeur en mètres (10000), pour NMM en degrés (0,1). On reviendra aux paramètres spécifiques NMM de manière plus détaillée ultérieurement.
  • les variables e_we et e_sn définissent l'extension de votre domaine dans la direction est-ouest et nord-sud. Cela définit concrètement le nombre de points de grille, en l'occurrence 100 dans les deux directions. On couvre alors un "rectangle" de 1000 km de côté centré sur notre ref_lat et ref_lon.
  • geog_data_path définit le chemin où vous avez extrait vos données géographiques.
  • geog_data_res définit la source de données géographique utilisée. Ces données sont disponibles en 3 résolution : 10m, 2m et 30s (pour 10 million, etc...). Nous prenons la résolution la plus fine...
Voici une illustration des concepts pour la définition du domaine.Vous remarquerez le positionnement du domaine par rapport à la carte dû à notre projection de type Lambert :

Définition du domaine de travail de WRF
Définition du domaine de travail de WRF

Voici finalement le fichier namelist.wps :

&share
wrf_core = 'ARW',
max_dom = 1,
start_date = '2013-07-28_00:00:00',
end_date = '2013-07-29_00:00:00',
interval_seconds = 10800
io_form_geogrid = 2,
debug_level = 1000,
/

&geogrid
parent_id = 1,
parent_grid_ratio = 1,
i_parent_start = 1,
j_parent_start = 1,
e_we = 100,
e_sn = 100,
geog_data_res = '30s',
dx = 10000,
dy = 10000,
map_proj = 'lambert',
ref_lat = 46.80,
ref_lon = 2.3372,
truelat1 = 45.8989,
truelat2 = 47.6960,
stand_lon = 2.3372,
geog_data_path = '/home/nicolas/Meteo/geog'
/

&ungrib
out_format = 'WPS',
prefix = 'FILE',
/

&metgrid
fg_name = 'FILE'
io_form_metgrid = 2,
/
Vous aurez sans doute remarqué que dans le fichier namelist.wps d'origine il y a plusieurs colonnes : ceci permet de créer des domaines imbriqués (nesting). Nous ne nous en servons pas ici, je les ai donc enlevées.

Les étapes du run


Nous allons maintenant pouvoir procéder aux différentes étapes du run. Nous devrons tout d'abord préparer des données à la résolution spatio-temporelle adaptée à notre simulation, et avec les variables requises par WRF. Puis l'étape suivante concernera le calcul proprement dit. De manière concrète et résumée, dans l'ordre :
  • On utilise l'utilitaire geogrid pour extraire les données géographiques à la bonne résolution pour notre domaine.
  • L'utilitaire ungrib sert à extraire les variables requises des fichiers GRIB (vent, pression, humidité...) et les mettre dans un format intermédiaire utilisé par WRF.
  • L'utilitaire metgrid sert à créer la grille du domaine proprement dite. Les variables meteo issues des deux étapes précédentes sont interpolées sur chaque point de grille horizontale du domaine. Les niveaux verticaux sont laissés inchangés.
  • L'utilitaire real sert quand à lui à créer les niveaux verticaux de la simulation. Il va donc interpoler verticalement les données créées par metgrid.
  • L'utilitaire wrf est le programme qui effectue la simulation proprement dite.

Geogrid


Depuis le dossier WPS, il suffit de lancer la commande :

# ./geogrid.exe 
 
Parsed 20 entries in GEOGRID.TBL
Processing domain 1 of 1
INFORM: Using 30s interpolator sequence for HGT_M.
INFORM: Using 30s data source for HGT_M.
INFORM: Using 30s interpolator sequence for LANDUSEF.
INFORM: Using 30s data source for LANDUSEF.
...
etc
...
Processing OL2
Processing OL3
Processing OL4
Processing VAR_SSO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Successful completion of geogrid. !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


A la fin, cela produit un fichier geo_em.d01.nc.
Ce n'est pas l'étape la plus difficile, mais en cas d'erreur, voici quelques causes possibles :
  • geog_data_path mal défini
  • des données manquent : avez-vous bien téléchargé toutes les données géographiques?
  • Paramètres de définition du domaine incorrects dans votre namelist. Vérifiez et re-vérifiez, comparez avec le fichier d'origine...

Ungrib


Assurez-vous d'avoir bien créés les liens symboliques vers vos GRIB avec link_grib, et de ne pas avoir de fichier "parasite". 

Il nous reste à indiquer à Ungrib comment sont définies les variables dans les fichiers GRIB GFS. En effet, chaque centre météo code différemment les différentes informations, il nous faut donc indiquer une table de correspondance. Pour cela, nous devons créer un lien symbolique vers la table de variables GFS :

# ln -s ungrib/Variable_Tables/Vtable.GFS Vtable

Tapez ensuite la commande :
 
# ./ungrib.exe
 
*** Starting program ungrib.exe ***
Start_date = 2013-07-28_00:00:00 , End_date = 2013-07-29_00:00:00
output format is WPS
INFORM: Interval value: 10800 seconds or 3.000000 hours
Path to intermediate files is ./
INFORM: Linked in png and jpeg libraries for Grib Edition 2
INFORM: Reading Grib Edition 2
...
Beaucoup de blahblah sur le contenu des données
...
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Successful completion of ungrib. !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Cela produit autant de fichier du style FILE:2013-07-28_* que de fichiers d'entrée.

Note : vous aurez probablement des messages à propos des variables du sol (SOIL) qu'il ne trouve par alors qu'elles sont bien dans le fichier. Ce message m'est arrivé sur ma VM en 32 bits. Je n'ai pas eu ce problème lors de mes essais en ubuntu 64 bits. Toutefois ce message n'est pas bloquant, vous n'aurez juste pas les températures et humidité du sol pour votre simulation. C'est dommage, mais pas trop gênant pour des simulations assez courtes, nous verrons comment nous en passer le moment venu.

Si vous avez des erreurs :
  • si ça vous dit que vos données ne sont pas en format GRIB, alors vous avez raté quelque chose à la compilation. Vérifiez que vous avez bien la bonne librairie Jasper, vos variables d'environnement, etc... Vous êtes bon pour recompiler. Mais avant de vous arracher les cheveux, vérifiez quand même que vous n'avez pas laissé trainer un fichier autre que GRIB dans le dossier DATA et qui aurait été linké par link_grib.
  • J'ai eu quelques erreurs, notamment concernant des dates ou d'autres erreurs, mais la plupart du temps c'était des problèmes dans mon namelist. Vérifiez, double-vérifiez !
  • Autre source d'erreur : le fichier Vtable. Avez-vous correctement lié ce fichier ? Est-ce bien le fichier GFS ?

Metgrid


Dernière étape du préprocesseur. La commande est :

# ./metgrid.exe

Processing domain 1 of 1
Processing 2013-07-28_00
FILE
...
Processing 2013-07-28_18
FILE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Successful completion of metgrid. !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Là encore, peut-être beaucoup de messages sur le manque de données sol (SOIL). Si vous avez eu des messages lors de l'étap Ungrib, alors c'est normal et non bloquant. Ces messages surviendront tout de même pour certains niveaux car GFS ne fournit pas tous les niveaux de sol que recherche metgrid.
Cette commande produit des fichiers du style met_em.d01.2013-07-28_00:00:00.nc (un par fichier d'entrée).

Vous ne devriez pas avoir trop d'erreurs bloquantes avec cette commande, si tout s'est bien passé auparavant. La seule que je me souvient avoir eu concerne le grid staggering. Si vous avez ce genre d'erreurs, c'est que vous avez modifié votre fichier METGRID.TBL et donc que vous avez déjà joué avec NMM ! On en reparlera donc plus tard.


Configuration de WRF


Précédemment, nous avons configuré le pré-processeur. Nous allons maintenant configurer WRF. Allez dans le dossier WRFV3/test/em_real et éditez fichier namelist.input. Comme d'habitude, c'est une bonne idée d'en faire une copie de sauvegarde avant de recopier celui que je vous fournit gracieusement. Et voici une description des variables d'intérêt :
  • run_days, run_hours, run_minutes, run_seconds : la durée de la simulation. 24H dans notre cas.
  • start_year, start_month, start_day, start_hour, start_minute, start_second : date de début, correspondant à vos fichiers de données
  • end_year, end_month, end_day, end_hour, end_minute, end_second : date de fin de la simulation.
  • interval_seconds : intervalle entre les fichiers de données (10800, comme dans le namelist.wps).
  • history_interval : intervalle auquel produire une écriture des données simulées dans le fichier de sortie. On veut avoir l'évolution heure par heure : on met donc 60 minutes.
  • time_step : pas de temps en secondes de la simulation. C'est la vitesse d'avancée du modèle : plus c'est court plus il y aura d'itérations et donc de calculs, plus c'est long moins il y en aura et le modèle terminera plus vite mais... si c'est trop long votre modèle devient instable et plantera. Préférez un pas de temps inférieur ou égal à 6*dx (dx en km).
  • time_step_fract_num, time_step_fract_den : permet de définir un pas de temps plus précis, en ajoutant une fraction pour faire par exemple 60s 1/4. On ne s'en sert pas ici, on a pile 60s.
  • e_we et e_sn représentent l'extention horizontale de votre domaine, à faire correspondre à votre namelist.wps
  • e_vert donne le nombre de niveaux verticaux à utiliser. L'atmosphère est alors découpée en autant de couches, espacées régulièrement par défaut. Il détermine donc la résolution verticale de votre modèle, afin de définir une grille tridimensionnelle de dimensions e_we x e_se x e_vert.
  • dx, dy donne la taille de grille, à faire correspondre également à votre namelist.wps.
Voici le fameux fichier :
&time_control
run_days = 0,
run_hours = 24,
run_minutes = 0,
run_seconds = 0,
start_year = 2013,
start_month = 07,
start_day = 28,
start_hour = 00,
start_minute = 00,
start_second = 00,
end_year = 2013,
end_month = 07,
end_day = 29,
end_hour = 00,
end_minute = 00,
end_second = 00,
interval_seconds = 10800
input_from_file = .true.,
history_interval = 60,
frames_per_outfile = 1000,
restart = .false.,
restart_interval = 5000,
io_form_history = 2
io_form_restart = 2
io_form_input = 2
io_form_boundary = 2
debug_level = 0
/



&domains
time_step = 60,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 100,
e_sn = 100,
e_vert = 30,
p_top_requested = 5000,
num_metgrid_levels = 27,
num_metgrid_soil_levels = 4,
dx = 10000,
dy = 10000,
grid_id = 1,
parent_id = 0,
i_parent_start = 1,
j_parent_start = 1,
parent_grid_ratio = 1,
parent_time_step_ratio = 1,
feedback = 1,
smooth_option = 0
/

&physics
mp_physics = 3,
ra_lw_physics = 1,
ra_sw_physics = 1,
radt = 30,
sf_sfclay_physics = 1,
sf_surface_physics = 2,
bl_pbl_physics = 1,
bldt = 0,
cu_physics = 1,
cudt = 5,
isfflx = 1,
ifsnow = 1,
icloud = 1,
surface_input_source = 1,
num_soil_layers = 4,
sf_urban_physics = 0,
/

&fdda
/

&dynamics
w_damping = 0,
diff_opt = 1,
km_opt = 4,
diff_6th_opt = 0,
diff_6th_factor = 0.12,
base_temp = 290.
damp_opt = 0,
zdamp = 5000.,
dampcoef = 0.2,
khdif = 0,
kvdif = 0,
non_hydrostatic = .true.,
moist_adv_opt = 1,
scalar_adv_opt = 1,
/

&bdy_control
spec_bdy_width = 5,
spec_zone = 1,
relax_zone = 4,
specified = .true.,
nested = .false.,
/

&grib2
/

&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/

Real


Dernière étape pour créer les fichiers d'entrée de WRF. Il nous reste à créer des liens symboliques vers les fichiers produits par metgrid :

# ln -sf ../../../WPS/met_em.d01.2013-07* .

Si vous avez compilé avec l'option serial ou smpar, alors tapez simplement :

# ./real.exe

Si vous avez compilé avec l'option dmpar tapez alors (-np 4 correspond au nombre de processeurs à exploiter) :

# mpirun -np 4 ./real.exe 
 
La sortie ressemble à cela (dans le cas de MPI elle n'est pas visible, vous la trouverez dans les fichiers rsl.out.000x pour chaque processus) :

Namelist logging not found in namelist.input. Using registry defaults for variables in logging.
Namelist dfi_control not found in namelist.input. Using registry defaults for variables in dfi_control
Namelist tc not found in namelist.input. Using registry defaults for variables in tc
Namelist noah_mp not found in namelist.input. Using registry defaults for variables in noah_mp
Namelist scm not found in namelist.input. Using registry defaults for variables in scm
Namelist fire not found in namelist.input. Using registry defaults for variables in fire
Namelist diags not found in namelist.input. Using registry defaults for variables in diags
Ntasks in X 2 , ntasks in Y 2
--- WARNING: traj_opt is zero, but num_traj is not zero; setting num_traj to zero.
--- NOTE: sst_update is 0, setting io_form_auxinput4 = 0 and auxinput4_interval = 0 for all domains
--- NOTE: grid_fdda is 0 for domain 1, setting gfdda interval and ending time to 0 for that domain.
--- NOTE: both grid_sfdda and pxlsm_soil_nudge are 0 for domain 1, setting sgfdda interval and ending time to 0 fo
r that domain.
--- NOTE: obs_nudge_opt is 0 for domain 1, setting obs nudging interval and ending time to 0 for that domain.
bl_pbl_physics /= 4, implies mfshconv must be 0, resetting
--- NOTE: num_soil_layers has been set to 5
REAL_EM V3.5 PREPROCESSOR
*************************************
Parent domain
ids,ide,jds,jde 1 100 1 100
ims,ime,jms,jme -4 57 -4 57
ips,ipe,jps,jpe 1 50 1 50
*************************************
DYNAMICS OPTION: Eulerian Mass Coordinate
alloc_space_field: domain 1 , 89528792 bytes allocated
Time period # 1 to process = 2013-07-28_00:00:00.
Time period # 2 to process = 2013-07-28_03:00:00.
Time period # 3 to process = 2013-07-28_06:00:00.
Time period # 4 to process = 2013-07-28_09:00:00.
Time period # 5 to process = 2013-07-28_12:00:00.
Time period # 6 to process = 2013-07-28_15:00:00.
Time period # 7 to process = 2013-07-28_18:00:00.
Time period # 8 to process = 2013-07-28_21:00:00.
Time period # 9 to process = 2013-07-29_00:00:00.
Total analysis times to input = 9.
....


d01 2013-07-29_00:00:00 real_em: SUCCESS COMPLETE REAL_EM INIT

Cette dernière ligne est plutôt bon signe ! Les fichiers wrfbdy_d01 et wrfinput_d01 sont produits.
Il peut encore y avoir quelques erreurs :
  • mismatch de dx, dy ou autre variable : très facile à résoudre car le message est explicite. Assurez-vous que ces valeurs correspondent bien au fichier namelist.wps et le cas échéant refaites le cycle complet du pré-processeur (geogrid/ungrib/metgrib) pour régénérer les données.
  • niveaux manquants pour la physique du sol (SOIL). Ceci arrive si vous avez eu des problèmes avec ces données dans les étapes précédentes. Dans ce cas, vous devrez utiliser une autre physique pour le sol, sans niveau. Editez votre namelist.input et modifiez les variables :

    num_metgrid_soil_levels = 0
    num_soil_layers=0
    sf_surface_physics = 1

WRF


Nous y voilà enfin ! Allez sans plus attendre voici la commande à utiliser si vous avez utilisé serial ou smpar :
 
# ./wrf.exe

Pour dmpar la commande devient (-np 4 pour le nombre de processeurs à exploiter) :

# mpirun -np 4 ./wrf.exe

La sortie ressemble à cela (dans les fichiers rsl* si vous utilisez MPI) :

taskid: 0 hostname: ubuntu
Namelist logging not found in namelist.input. Using registry defaults for variables in logging.
Namelist dfi_control not found in namelist.input. Using registry defaults for variables in dfi_control
Namelist tc not found in namelist.input. Using registry defaults for variables in tc
Namelist noah_mp not found in namelist.input. Using registry defaults for variables in noah_mp
Namelist scm not found in namelist.input. Using registry defaults for variables in scm
Namelist fire not found in namelist.input. Using registry defaults for variables in fire
Namelist diags not found in namelist.input. Using registry defaults for variables in diags
Ntasks in X 2 , ntasks in Y 2
--- WARNING: traj_opt is zero, but num_traj is not zero; setting num_traj to zero.
--- NOTE: sst_update is 0, setting io_form_auxinput4 = 0 and auxinput4_interval = 0 for all domains
--- NOTE: grid_fdda is 0 for domain 1, setting gfdda interval and ending time to 0 for that domain.
--- NOTE: both grid_sfdda and pxlsm_soil_nudge are 0 for domain 1, setting sgfdda interval and ending time to 0 for that domain.
--- NOTE: obs_nudge_opt is 0 for domain 1, setting obs nudging interval and ending time to 0 for that domain.
bl_pbl_physics /= 4, implies mfshconv must be 0, resetting
--- NOTE: num_soil_layers has been set to 5
WRF V3.5 MODEL
*************************************
Parent domain
ids,ide,jds,jde 1 100 1 100
ims,ime,jms,jme -4 57 -4 57
ips,ipe,jps,jpe 1 50 1 50
*************************************
DYNAMICS OPTION: Eulerian Mass Coordinate
alloc_space_field: domain 1 , 80026424 bytes allocated
med_initialdata_input: calling input_input
Timing for processing wrfinput file (stream 0) for domain 1: 0.26179 elapsed seconds
INPUT LandUse = "USGS"
LANDUSE TYPE = "USGS" FOUND 33 CATEGORIES 2 SEASONS WATER CATEGORY = 16 SNOW CATEGORY = 24
Timing for Writing wrfout_d01_2013-07-28_00:00:00 for domain 1: 0.59684 elapsed seconds
Timing for processing lateral boundary for domain 1: 0.02816 elapsed seconds
Tile Strategy is not specified. Assuming 1D-Y
WRF TILE 1 IS 1 IE 50 JS 1 JE 50
WRF NUMBER OF TILES = 1
Timing for main: time 2013-07-28_00:01:00 on domain 1: 2.12459 elapsed seconds
...
Les différentes itérations
...
Timing for Writing wrfout_d01_2013-07-29_00:00:00 for domain 1: 0.30211 elapsed seconds
d01 2013-07-29_00:00:00 wrf: SUCCESS COMPLETE WRF

Si vous en arrivez là, alors :
!!!!! FELICITATIONS !!!!!

Vous avez alors votre précieux fichier de sortie au format NetCDF. Il portera un nom du style wrfout_d01_2013-07-28_00:00:00. C'est ce fichier qui contient les données que vous utiliserez pour générer vos cartes, avec NCL par exemple. 

Pour les moins chanceux, il peut encore y avoir des problèmes. Mais comme d'habitude, il n'y a pas de problèmes, que des solutions :
  • le programme commence à calculer un peu (timing for main time défile) puis il commence à bégayer des choses bizarres et plante avec un segmentation fault (ou core dump). Cela signifie que vous avez mis un pas de temps trop large et le modèle devient numériquement instable ! Essayez de baisser la variable time_step de votre namelist.input
  • un segmentation fault peut également provenir d'une limite de pile. Tapez ulimit -s unlimited avant de lancer WRF et voyez si ça résoud le problème.
  • Toute autre erreur résulte d'un problème de paramétrage de namelist ou de données d'entrées erronnées. Vérifiez bien les étapes précédentes et refaites-les si nécessaire
Dans la prochaine partie, nous ferons une pause plus ludique avant d'aborder le coeur NMM. Nous avons en effet bien mérité de visualiser le résultat de tous nos efforts : nous sortirons donc nos premières cartes de prévision. Prochain article : Partie 6 - Visualisation des données.