Skip to content

WRF workflows

Specific tasks

Before running the model

Defining the vertical grid

Customizing model orography

Defining a new geographical database

Using ECMWF data as IC/BC

The long story made short is: you should link grib1 files and process them with ungrib.exe using Vtable.ECMWF_sigma.

More in detail, since a few years ECMWF has been distributing a mixture of grib2 and grib1 files. Namely:

  • grib1 files for surface and soil model levels.
  • grib2 files for atmospheric model levels.

The WPS has a predefined Vtable for grib1 files from ECMWF, so the easiest way to process ECMWF data is to:

  1. convert model-level grib2 files to grib1
  2. if necessary, for every time stamp, concatenate the model-level and surface grib1 files into a single file. This is only necessary if the grib1 and grib2 data were downloaded as separate sets of GRIB files.
  3. process the resulting files with ungrib after linking ungrib/Variable_Tables/Vtable.ECMWF_sigma as Vtable

In detail:

  1. Conversion to grib1 (needs the grib_set utility from eccodes):

    convert to grib1
    1
2
3
4
5
    for i in det.CROSSINN.mlv.20190913.0000.f*.grib2; 
do
  j=`basename $i .grib2`; 
  grib_set -s deletePV=1,edition=1 ${i} ${j}; 
done

  2. Concatenation of grib files (two sets of files, *mlv* and *sfc*, with names ending with "grib1" yield a new set of files with names ending with "grib"; everything is grib1):

    concatenate grib files
    1
2
3
4
5
6
    for i in det.CROSSINN.mlv.20190913.0000.f*.grib1; 
do 
  j=`echo $i|sed 's/.mlv./.sfc./'`; 
  k=`echo $i|sed 's/.mlv././'|sed 's/.grib1/.grib/'`; 
  cat $i $j > $k; 
done

  3. In the WPS main directory:

    link grib files and convert
    1
2
3
    link_grib.csh /data/GRIB_IC_for_LAM/ECMWF/20190913_CROSSINN_IOP8/det.CROSSINN.20190913.0000.f*.grib
ln -s ungrib/Variable_Tables/Vtable.ECMWF_sigma Vtable
./ungrib.exe

An alternative procedure would be to convert everything to grib2 instead of grib1. Then, one has to use a Vtable with grib2 information for the surface fields, for instance the one included here at the bottom. But: Data from the bottom soil level will not be read correctly with this Vtable, because the Level2 value for the bottom level is actually MISSING in grib2 files (at the moment of writing, 6 May 2022; this may be fixed in the future).

Text Only
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
GRIB1| Level| From |  To  | metgrid  | metgrid  | metgrid                                  |GRIB2|GRIB2|GRIB2|GRIB2|
Param| Type |Level1|Level2| Name     | Units    | Description                              |Discp|Catgy|Param|Level|
-----+------+------+------+----------+----------+------------------------------------------+-----------------------+
 130 | 109  |   *  |      | TT       | K        | Temperature                              |  0  |  0  |  0  | 105 |
 131 | 109  |   *  |      | UU       | m s-1    | U                                        |  0  |  2  |  2  | 105 |
 132 | 109  |   *  |      | VV       | m s-1    | V                                        |  0  |  2  |  3  | 105 |
 133 | 109  |   *  |      | SPECHUMD | kg kg-1  | Specific humidity                        |  0  |  1  |  0  | 105 |
 152 | 109  |   *  |      | LOGSFP   | Pa       | Log surface pressure                     |  0  |  3  |  25 | 105 |
 129 | 109  |   *  |      | SOILGEO  | m        | Surface geopotential                     |  0  |  3  |  4  |  1  |
     | 109  |   *  |      | SOILHGT  | m        | Terrain field of source analysis         |  0  |  3  |  5  |  1  |
 134 | 109  |   1  |      | PSFCH    | Pa       |                                          |  0  |  3  |  0  |  1  |
 157 | 109  |   *  |      | RH       | %        | Relative Humidity                        |  0  |  1  |  1  | 105 |
 165 |  1   |   0  |      | UU       | m s-1    | U                                        |  0  |  2  |  2  | 103 |
 166 |  1   |   0  |      | VV       | m s-1    | V                                        |  0  |  2  |  3  | 103 |
 167 |  1   |   0  |      | TT       | K        | Temperature                              |  0  |  0  |  0  | 103 |
 168 |  1   |   0  |      | DEWPT    | K        |                                          |  0  |  0  |  6  | 103 |
 172 |  1   |   0  |      | LANDSEA  | 0/1 Flag | Land/Sea flag                            |  2  |  0  |  0  |  1  |
 151 |  1   |   0  |      | PMSL     | Pa       | Sea-level Pressure                       |  0  |  3  |  0  | 101 |
 235 |  1   |   0  |      | SKINTEMP | K        | Sea-Surface Temperature                  |  0  |  0  | 17  |  1  |
  34 |  1   |   0  |      | SST      | K        | Sea-Surface Temperature                  |  10 |  3  |  0  |  1  |
 139 | 112  |     0|   700| ST000007 | K        | T of 0-7 cm ground layer                 | 192 | 128 | 139 | 106 |
 170 | 112  |   700|  2800| ST007028 | K        | T of 7-28 cm ground layer                | 192 | 128 | 170 | 106 |
 183 | 112  |  2800| 10000| ST028100 | K        | T of 28-100 cm ground layer              | 192 | 128 | 183 | 106 |
 236 | 112  | 10000|     0| ST100289 | K        | T of 100-289 cm ground layer             | 192 | 128 | 236 | 106 |
  39 | 112  |     0|   700| SM000007 | fraction | Soil moisture of 0-7 cm ground layer     | 192 | 128 |  39 | 106 |
  40 | 112  |   700|  2800| SM007028 | fraction | Soil moisture of 7-28 cm ground layer    | 192 | 128 |  40 | 106 |
  41 | 112  |  2800| 10000| SM028100 | fraction | Soil moisture of 28-100 cm ground layer  | 192 | 128 |  41 | 106 |
  42 | 112  | 10000|     0| SM100289 | fraction | Soil moisture of 100-289 cm ground layer | 192 | 128 |  42 | 106 |
-----+------+------+------+----------+----------+------------------------------------------+-----------------------+

Spinning up soil fields

After running the model

Converting model output to CF-compliant NetCDF

  1. To convert WRF output to CF-compliant NetCDF, use wrfout_to_cf.ncl (from https://sundowner.colorado.edu/wrfout_to_cf/overview.html):

    Text Only
    1
    ncl 'file_in="wrfinput_d01"' 'file_out="wrfpost.nc"' wrfout_to_cf.ncl

Interpolating model output to a new grid

  1. First convert to CF-compliant NetCDF (see above)

  2. Then use cdo to interpolate the CF-compliant WRF output:

    Text Only
    1
    cdo -remapnn,gridfile.lonlat.txt wrfpost.nc wrfpost_interpolated.nc

  3. In the code snippet above, -remapnn specifies the interpolation engine, in this case nearest-neighbour. See alternatives here: https://code.mpimet.mpg.de/projects/cdo/wiki/Tutorial#Horizontal-fields

  4. File gridfile.lonlat.txt contans the grid specifications, e.g.:

    Text Only
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
    gridtype  = lonlat
gridsize  = 721801
xsize     = 1201
ysize     = 601
xname     = lon
xlongname = "longitude" 
xunits    = "degrees_east" 
yname     = lat
ylongname = "latitude" 
yunits    = "degrees_north" 
xfirst    = 5.00
xinc      = 0.01
yfirst    = 43.00
yinc      = 0.01

Subsetting model output

Further compression of model output (data packing)

3D visualization

For 3D visualization of WRF output, it is recommended to use either Paraview or Mayavi.

  • Both softwares are based on the Visualization Toolkit (VTK) libraries, so visualizations are rather similar in the end.

  • Both sotwares can be used interactively from a graphical user interface or in batch mode (i.e., writing the visualization directives in a Python script).

  • While Paraview requires converting model data into one of a few supported formats, Mayavi supports direct rendering of Numpy objects, so it is easier to integrate it into Python code.

  • It is recommended to run 3D visualization software on GPUs. Running on a CPU (e.g., own laptop) is possible, but will be extremely slow. CPU is not the only bottleneck, because visualization software uses a lot of computer memory. Rendering 3D fields, in particular, is out of reach for normal laptops with 8GB or 16GB of RAM. Paraview is available on VSC5 and should be available soon on srvx8. Currently, Mayavi must be installed by individual users as a Python package.

Notes for readers/contributors: (1) Mayavi is untested yet. (2) It would be useful to add example batch scripts for both Paraview and Mayavi.

Paraview workflow

  1. Pre-requisite: download and install the Paraview application on your computer.

  2. Log in to VSC5 in a terminal window.

  3. On VSC5, convert the WRF output in a format that Paraview can ingest. One option is to use siso.

    Bash
    1
2
    siso -f vts ~/serafins/TEAMx_LES/output/100m/wrfout_d01_0001-01-01_00\:00\:00 > siso.out 2>&1 &
siso -f vts --planar ~/serafins/TEAMx_LES/output/100m/wrfout_d01_0001-01-01_00\:00\:00 > siso_sfc.out 2>&1 &
    The first and second statements handle respectively 3D and 2D WRF output. They process the native output from WRF in netcdf format and return collections of files in VTS format (the VTK format for structured grids). There will be two independent datasets (for 3D and 2D output).

  4. In the VSC5 terminal, request access to a GPU node. One of the private IMGW nodes has a GPU, and can be accessed with specific account/partition/quality of service directives.

    Bash
    1
2
3
4
5
6
7
    (zen3) [sserafin4@l50 ~]$ salloc -N 1 --gres=gpu:2 --account=p71386 -p zen3_0512_a100x2 -q p71386_a100dual
salloc: Pending job allocation 233600
salloc: job 233600 queued and waiting for resources
salloc: job 233600 has been allocated resources
salloc: Granted job allocation 233600
salloc: Waiting for resource configuration
salloc: Nodes n3072-006 are ready for job

  5. Once the GPU node becomes available, open up a new terminal session on your local machine, and set up an ssh tunnel to the GPU node through the login node.

    Bash
    1
    (mypy39) stefano@stefano-XPS-13-9370:~$ ssh -L 11111:n3072-006:11112 sserafin4@vsc5.vsc.ac.at
    This will redirect TCP/IP traffic from port 11111 of your local machine to port 11112 of the VSC5 GPU node, through the VSC5 login node. Port numbers are arbitary, but the remote port (11112) needs to match the Paraview server settings (see below).

  6. In the VSC5 terminal, log in to the GPU node:

    Bash
    1
2
3
4
5
    (zen3) [sserafin4@l50 ~]$ ssh n3072-006
Warning: Permanently added 'n3072-006,10.191.72.6' (ECDSA) to the list of known hosts.
sserafin4@n3072-006's password: 

(zen3) [sserafin4@n3072-006 ~]$

  7. In the VSC5 terminal on the GPU node, load the Paraview module and start the Paraview server:

    Bash
    1
2
3
4
5
    (zen3) [sserafin4@n3072-006 ~]$ module load paraview
(zen3) [sserafin4@n3072-006 ~]$ pvserver --force-offscreen-rendering --server-port=11112
Waiting for client...
Connection URL: cs://n3072-006:11112
Accepting connection(s): n3072-006:11112

  8. On your local machine, open the Paraview client (graphical user interface, GUI). Then select File > Connect and enter the url of the Paraview server (localhost:11111). Select the datasets you want to display and work on them in the GUI. Save the Paraview state to avoid repeating work at the next session. Paraview has extensive documentation, tutorials (one, two and three) and a wiki.

Mayavi workflow

Not tested yet.

Creating a video

Whether done with Paraview or with Mayavi, the visualization will result in a collection of png files, e.g., InnValley.%04d.png. There are several tools to convert invidual frames into movies. Among them, ffmpeg and apngasm. At the moment neither of them is available on IMGW servers (precompiled binaries are available through apt-get for Ubuntu).

The basic method to create an mp4 movie is:

Bash
1
ffmpeg -i InnValley.%04d.png -c:v libx264 -r 12 -pix_fmt yuv420p InnValley.mp4

The method above might return an error if frames have an odd number of pixels in one dimension:

Bash
1
[libx264 @ 0x5651e5f02980] height not divisible by 2 (1066x1083)

The fix is as follows:

Bash
1
ffmpeg -i InnValley.%04d.png -c:v libx264 -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" -r 12 -pix_fmt yuv420p InnValley.mp4

It is possible to add movie repetitions (similar to a loop). In this case, 3 additional loops are appended after the first one:

Bash
1
ffmpeg -stream_loop 3 -framerate 12 -i InnValley.%04d.png -c:v libx264 -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" -pix_fmt yuv420p InnValley.mp4

It also possible to generate movies in other formats, better suited for the web:

  • webp (most efficient compression for loops):

    Bash
    1
    ffmpeg -framerate 12 -i InnValley.%04d.png InnValley.webp

  • animated png (bigger in size):

    Bash
    1
    apngasm InnValley.png InnValley.0*png

  • gif (much bigger in size):

    Bash
    1
    ffmpeg -framerate 12 -i InnValley.%04d.png InnValley.gif

  • For the example dataset, the collection of raw png files takes 59 MB while the video file sizes range between 4.5 and 70 MB:

    Bash
    1
2
3
4
5
6
7
8
    (mypy39) stefano@stefano-XPS-13-9370:~/Desktop/Paraview_animation/anim$ du -hcs InnValley.0*png
59M total

(mypy39) stefano@stefano-XPS-13-9370:~/Desktop/Paraview_animation/anim$ du -hcs InnValley.[pgmw]*
70M InnValley.gif
14M InnValley.mp4
51M InnValley.png
4,5M    InnValley.webp

Useful tools

Last update: May 5, 2023
Created: May 5, 2023