RPN standard files

The RPN standard file format is a file format developed by "Le Centre de recherche en prévision numérique" (RPN) to store weather and climate related data. Usually, there is one 2-D (latitude and longitude) record for each level, variable and timestep.

For the end-user, the appearance of the RPN is that of a flat-file database structure. A single file contains numerous records which hold fields on 1, 2 or 3D data grids and self-describing information about the fields (header).

A thorough set of FST-related routines are contained in the RPN Library (armnlib). These routines provide easy access to the file contents. They can be called from FORTRAN or C routines.

A complete description of the FST file format and the utilities which can be used to manipulate and view their contents is available online at www.cmc.ec.gc.ca/rpn/utilities.

Sometimes the files are archived in an cmc-archive file (*.ca). If this is the case they need to be unarchived first.

The maximum size of RPN file is a little more than 8 GB!!! If you are trying to write a file larger than the maximum size the file will not get closed. Meaning no application will be able to open it and it can be considered as "corrupted". The usual error message one get when trying to open such a file is "file (unit=...) currently used by another application in write mode".

Looking at RPN files

vl3, vl4, vl5 -> list variable with description

There is also a set of 'vl' (variable list) commands (one for each model version) which will give you a list of all the variables in an RPN file. The first description is always correct, but rarely contains the unit. For physics fields take the second description with a grain of salt. If it matches the first one, it's unit is correct as well. If it doesn't match, ignore it. However, to all units of accumulators a 'per second' needs to get added, for example:

precipitation is not in [m] but in [m/s]
radiation fields are not in [J/m²] but in [J/s/m²] so in [W/m²]

Usage:

vl5 filename

voir -> list the content of an RPN file

'voir' is the most simple way to list the content of an RPN file (also xvoir)
It prints the headers for all of the records of the given RPN standard file:

voir [-style] -iment rpn-file

Without the -style option the date displayed is for raw GEM model output the start date of the simulation, for averages/accumulators/minimum/maximum fields it is the start of the avg/accum/min/max period and levels are coded.

With the -style option the date displayed is for raw GEM model output the validity (actual) date of the fields, for averages/accumulators/minimum/maximum fields it is the end of the avg/accum/min/max period and levels are decoded.

-style key	levels	date
-style key	levels	instantaneous fields	avg/accum/min/max
without -style	coded	start date of simulation	start date of avg/accum/min/max period
with -style	decoded	actual date	end date of avg/accum/min/max period

You can find more information about the output of 'voir' in this PDF: voir_output_explanation.pdf

You can also use the following aliases:

alias vo='voir -iment'
alias vos='voir -style -iment'

Since the output of 'voir' can be rather long I usually pipe it into less, for example:

vo filename | less
vos filename | less

xrec -> graphical visualization tool

xrec is a X-window visualization program used to display 2D (horizontal or vertical) meteorological fields stored in FST files. Usage:

xrec

or

xrec -imflds list_of_RPN_files

You can also use the following alias:

alias xr='xrec -imflds'

The above command will open 3 windows (so make sure you connected with 'ssh -YC ...'!!!):

xrecControlPanel : Here you make general selections about what you want to display
xrecRecordSelector : Here all the records you can see with 'voir' (see above) are listed. You can double click on them to be displayed.
xrecDisplay : Results will be displayed in this window. You might want to change its size to match your domain

Do not trust the wind direction for any other vectors but 'UU' and 'VV'!!!

Usage

Among other things, xrec can get used to :

display 2-D fields
vertical cross sections and profiles
time series for selected points.
time animation

Check the full documentation (xrec-documentation) about how to do this or just click on the buttons and look at the menus.

Dictionary

xrec is using a dictionary to define the title, unit (including a multiplication factor) and contour intervals for each variable.
The name of the dictionary is:
~/.recrc
It needs to have exactly this name and it needs to be in your home. To have a dictionary to start from you are welcome to take a copy of mine ( cp ~winger/.recrc ~/. ) or ask a college if you can copy hers/his.
You can find more information about this dictionary by clicking on the following link: xrec-dictionary

Remove default contour lines

The first time you open any field with xrec you will see there are "Colors" and "Contours". If you do not want to have "Contours" by default (you can always add or remove them later) do the following:

Open xrec with any file
Click on "Display" and unselect "Contours"
Click on "File" and on "Quit and save configuration"
Open xrec with any file again
Select any file. If the contours are gone, you did the above correctly. If the background is now gray and you would like to have a white one again...
Click on "Options" and "Contours" - a new window will open. In that window...
Click on "Background" and select "Color:" "white"
Click on "File" and on "Quit and save configuration"

The next time you open xrec the default contour lines should be back and the background should be white again.

High resolution geography

If you want to use high resolution geography on our UQAM servers you need to export:

export GDB_PATH=~data/Xrec_high_res

Set min/max

You can use the Min-Max tool to fix a minimum and maximum value for a specific field.

    1) Open the Min-Max menu window by clicking on "Options" in the Control Panel and select "Min-Max / Missing values...".
    2) In the "Code" box search for and select the field for which you want to fix min/max.
    3) In the "Max" box replace "Auto" with your maximum value. Do the same for "Min".
    4) Now click on "Fixed"

Full documentation

You can find the full xrec documentation by clicking on the following link: xrec-documentation

r.fstliste -> list the content of an RPN file

r.fstliste is available if you need a bit more information than what is provided by voir.
It can be used in batch mode to extract specific or all of the record parameters of a RPN standard file record(s) based on a criteria given through the selection keys:

r.fstliste -izfst RPN-file [-datev -vdatev -ip1 -typvar -nomvar -col ...]

Code/Decode levels

In RPN files levels are coded. The different types of levels are:

0 - height [m] (metres)
1 - sigma [sg] (0.0->1.0)
2 - pressure [mb] (millibars)
3 - arbitrary code
4 - height [M] (metres) with respect to ground level
5 - hybrid coordinates [hy] (0.0->1.0)
6 - theta [th]

r.ip1 (Shell)

r.ip1 is a shell command that can be used to code and decode RPN levels as coded in RPN standard files. For more information click on the following link (login/password are the usual 'science'):

'r.ip1' documentation

Or simply execute the following command in a terminal:

r.ip1 -h

convertIp (Python)

'convertIp' is a function of the RPN Python module to code and decode RPN levels. To get access to it you first have to import the module, for example with:

import rpnpy.librmn.all as rmn

Once the module is loaded you can get a documentation of 'convertIp' with:

help(rmn.convertIp)

Example to convert pressure level to coded IP1:

import rpnpy.librmn.all as rmn # Module to read RPN files

p_lev = 500
ip1 = rmn.convertIp(rmn.CONVIP_ENCODE, p_lev, rmn.KIND_PRESSURE) # Code level 'p_lev' hPa into IP1 code
print ('Coded value for', p_lev, 'hPa is:', ip1)

Answer:
Coded value for 500 hPa is: 41394464

Example to decode IP1:

import rpnpy.librmn.all as rmn             # Module to read RPN files

# Read one record
filename = '...'    # Name of RPN file to read
varname = '...'                            # Name of variable to read
fid = rmn.fstopenall(filename,rmn.FST_RO)   # Open the file
rec = rmn.fstlir(fid,nomvar=varname)        # Read the full record of variable 'varname'
ip1 = rec['ip1']    # Read the IP1 code of the level

(val, kind) = rmn.convertIp(rmn.CONVIP_DECODE, ip1) # Decode IP1 into level value and unit
print ('Decoded value of', ip1, 'is:', val, rmn.kindToString(kind))

rmn.fstcloseall(fid) # Close the RPN file

Answer:
Decoded value of 41394464 is: 500.0 mb

convip (Fortran)

'convip' is a Fortran function to code and decode RPN levels. 'convip' is part of the RPN Fortran library. To get access to it you need to load the library 'librmn.a' when creating the executable. You can do that by added the key '-lrmn' to the compile command. For example:

s.f90 source_code.f90 -o executable.Abs -lrmn

For more information click on the following link (login/password are the usual 'science'):

'convip' documentation

Code/Decode time stamps

Dates and times are coded in RPN files. This also means that it is not possible to have any times in RPN files:

Between 1900 01 01 and 1980 01 01 the time resolution is 1 hour.
Between 1980 01 01 and 2236 01 01 the time resolution is 5 seconds.
Outside these time periods the time resolution is 3 hours.

r.date (Shell)

'r.date' is a shell command that can be used to code and decode RPN dates/times as coded in RPN standard files. For more information click on the following link (login/password are the usual 'science'):

'r.date' documentation

Or simply execute the following command in a terminal:

r.date -h

newdate (Python)

'newdate' is a Python function to code and decode RPN dates/times. To get access to it you first have to import the module, for example with:

import rpnpy.librmn.all as rmn

Once the module is loaded you can get a documentation of 'newdate' with:

help(rmn.newdate)

Example to convert date & time to coded timestamp:

import rpnpy.librmn.all as rmn # Module to read RPN files

date = 20241106 # YYYYMMDD
time = 13300000 # hhmmsshh => 13h30

timestamp = rmn.newdate(rmn.NEWDATE_PRINT2STAMP, date, time) # Code date and time into timestamp

print ('Coded value for date', date, 'at', time, 'is:', timestamp)

Answer:
Coded value of date 20241106 at 13300000 is: 477041750

Example to convert coded timestamp to date & time:

import rpnpy.librmn.all as rmn              # Module to read RPN files

# Read one record
filename = '...'                            # Name of RPN file to read
varname = '...'                            # Name of variable to read

fid = rmn.fstopenall(filename,rmn.FST_RO)   # Open the file
rec = rmn.fstlir(fid,nomvar=varname)        # Read the full record of variable 'varname'
timestamp = rec['datev']                    # Read the code of the validity date

(date, time) = rmn.newdate(rmn.NEWDATE_STAMP2PRINT, timestamp) # Decodes timestamp to date (YYYYMMDD) and time (hhmmsshh)
print ('Decoded value of timestamp', timestamp, 'is:', date, 'at', time)

rmn.fstcloseall(fid)                        # Close the RPN file

Answer:
Decoded value of timestamp 477041750 is: 20241106 at 13300000

newdate (Fortran)

'newdate' is a Fortran function to code and decode RPN dates/times. 'newdate' is part of the RPN Fortran library. To get access to it you need to load the library 'librmn.a' when creating the executable. You can do that by added the key '-lrmn' to the compile command. For example:

s.f90 source_code.f90 -o executable.Abs -lrmn

For more information click on the following link (login/password are the usual 'science'):

'newdate' documentation

Manipulating RPN files

The two main tools to modify FST files are called editfst and r.diag. And of course, you can also write your own C or Fortran programs.

editfst

editfst (also called editfst2000) is a utility used for editing and copying records from RPN standard files into a new or existing RPN standard file. It can do a straight (complete) copy of the input file or it can copy records selectively as indicated from the standard input or from a file of directives named in the "-i" key. One can select (via desire), exclude (via exclure) or modify (zap) records by their:

TYPVAR: type of variable; 'A'(analysis),'P'(forecast),'X'
NOMVAR: 4 caracter variable name
ETIKET: 12 caracter label recorded from the original model run
DATE: CMC date timestamp
IP1: the level of the field (in pressure,eta,meters,other)
IP2: the hour of the forecast (rounded off if not exact)
IP3: value of 0 unless otherwise modified

Example:
You could use it to select i.e. temperature (TT) and winds (UU,VV) at 500 hPa:

echo "desire (-1,['UU','VV'], -1 ,-1 ,[500.,MBAR])" | editfst -s input_file -d output_file

Note:
editfst appends to its output file. This means that if editfst is used more than once writing into the same output file, the output of the second set of commands will be appended to the output of the first set.

r.diag

r.diag (also called r.diag2000) is a toolkit used to manipulate FST or CCCma files. As editfst, it can be used to select or exclude specific records but in addition to this it also contains a large set of functions which can be used to manipulate the files.

'r.diag' is always used in the way:

r.diag function input_file(s) [ output_file(s) key(s) ]

To get information about a specific function simply type:

r.diag function

Notes:

r.diag generally overwrites its output files. There are a few functions, eg. select, which have an option to append to the output file instead of overwriting it.
By default, output files will be written under ${TMPDIR}!!!
If there is no path written in front of the output file name and if the output file does not yet exist in the local directory, it will be written in one of the ${BIG_TMPDIR} or ${TMPDIR} directories. Also, if the input file is not in the local directory, r.diag will look for it in the ${BIG_TMPDIR} and ${TMPDIR} directories.
So if you want to have the output file in a certain directory you have to specify said directory name before the file name. Fir the current directory you can just preceed the file name by './'.

Examples:

'select' - Select fields according to name and level

For example, temperature (TT) and winds (UU,VV) at 500 hPa can get selected with 'select' with the following command:

r.diag select input_file output_file_tt output_file_uu output_file_vv -name TT UU VV -lv1 500

'gsapl' - Interpolate model to pressure levels

Fields on model levels can get interpolated to a set of pressure levels with 'gsapl'.

Input files needed

Files can have any names, but to facilitate the explanation the following names are chosen for the example below:

mod_lev : input file on model levels
pres : input file containing surface pressure for the same time steps which are in 'mod_lev'
niveau.dat : file containing list if pressure levels to interpolate to. See more below!!!

Also just to facilitate the example below, the following parameters should get set:

plv : parameter set to number of pressure levels in 'niveau.dat'
kind : parameter set to type of model levels in the input file. Depending on the model version with which the output got produced set:
          kind=GEM2 : sigma model level ('hyb' in gem_settings.nml) starting NOT with 0.000
          kind=GEM3 : hybrid model level ('hyb' in gem_settings.nml) starting with 0.000
          kind=GEM4 : GEM4 staggered fields with '!!' record
          kind=GEM5 : GEM5 staggered fields with '!!' record

File 'niveau.dat'

You need to create this file which needs to contain the list of all the pressure levels you want to interpolate to. This file needs to have the following format:

max 16 levels per row
5 digits per level
right bound
all integer numbers!!!
Here are some example to put pressure levels smaller than 1:
     '-1100' stands for '0.1'
     '-1200' stands for '0.2'
        :
   ' -100' is the same as '1'

For example:

# Create file 'niveau.dat' (instead of with 'cat' you can also create it with an editor)
cat > niveau.dat <<+
-1100-1200-1400-1600 -100    2    4    6   10   20   30   50   70 100 150 200
250 300 400 500 600 700 850 925 1000
+

kind=GEM5   # When output is from GEM5
plv=25            # Number of pressure levels in niveau.dat

# As input, 'gsapl' needs the log of the surface pressure. So you need to calculate this first, for example:
# Create the log of the surface pressure
r.diag loge pres   log_pres
# Change variable name to 'LP' (required name by 'gsapl')
r.diag newnam log_pres lp -name LP
# Now the file 'lp' contains the log of the surface pressure and the variable has the name 'LP'.

# IF(!) there is more than 1 variable on model levels in your input file, you need to create a set of surface pressure for each variable!!!
# To do this, do the following. Otherwise skip the next four lines of code!
# Create a file 'mlp' containing as many records of the same ln(P0) as there are variables in the input file 'mod_lev'.
# And of course there must be one "set" of ln(P0) for each time step.
r.diag select cubet univeau -lv1 1.0
r.diag xlin   univeau un -a 0 -b 1 -name LP
r.diag gmlt   un lp mlp
mv mlp lp

# Interpolating all fields on model levels of the input file 'mod_lev' to pressure levels
# Make sure the input file 'mode_lev' contains the level record: 'HY' for GEM3 and '!!' for GEM4 and up.
r.diag gsapl mod_lev lp pres_lev -plv ${plv} -kind ${kind} < niveau.dat

The file 'pres_lev' should now contain the same fields and timesteps as the input file 'mod_lev' but on the pressure level set in 'niveau.dat'.

Writing your own C and Fortran programs

There is a number of functions which can be called in a C or Fortran program to read and write FST files, and they are documented in the RPN Standard File Functions. These functions are part of the library 'librmn'.

There are functions to:

open and close an FST file
search for a record
read a record
read the field descriptor of a record
write a recordcdf2rpn

Have a look here on how to compile.

Conversion: RPN ↔ NetCDF format

The tool 'cdf2rpn' which can get used to convert RPN to NetCDF format and vice versa. You can get the documentation of cdf2rpn by simply executing the command in a terminal:

cdf2rpn

The general usage is the following:

Convert from NetCDF to RPN format:

cdf2rpn -rpn RPN_standard_file -cdf NetCDF_file -attr attribut_list

Convert from RPN to NetCDF format - one needs to add the key '-dir' (for 'direction'):

cdf2rpn -rpn RPN_standard_file -cdf NetCDF_file -attr attribut_list -dir

Notes:

If you do not add a path to the output file name the file will end up in $TMPDIR!!!
All fields in a file must have the same levels.
You cannot have i.e. sigma, hybrid, arbitrary and pressure levels in the same file!!! Nor can there be, for example, a field on 20 pressure levels and another on only 1 pressure level. If your RPN input file contains fields with different types or numbers of levels you can select the fields your want to convert with 'r.diag' and/or 'editfst' - see above.
If converting RPN to NetCDF make sure your RPN file contains all necessary grid descriptors ('^^', '>>', '^>', 'HY','!!')
All variables you want to convert need to be in the attribute list.
You can find an example file on the servers under: ~winger/Scripts/NetCDF_converter/attribut_netcdf.dat
For each variable that is to be converted needs to be an entry starting with 'def_attribut' at the end of the file. Higher up in the file you find the documentation - see also below.
Make sure that the factors are correct!!! If one or more of your variables are not in this file, just take a copy and add the missing variables.

def_attribut[ name; ccname; mult; add;
               nom_attribut = valeur_attribut;
               nom_attribut = valeur_attribut; N;
               nom_attribut = valeur_attribut; valeur_attribut; N; ]

Where:

   name: name of netCDF variable - mandatory
   cccname: name of RPN variable - mandatory
   mult: multiplication factor for units (type real) DEFAUT 1.0
   add: additive factor for units (type real) DEFAUT 0.0
   nom_attribut: name of the attribute
   valeur_attribut: value associated to "nom_attribut"
   N: integer vale defining the type of "valeur_attribut" (2=character, 3=integer*2, 5=real*4)

Notes:

       - Unit conversion from RPN → netCDF:         unit(RPN) = mult*unit(netCDF)+add
       - Unit conversion from netCDF → RPN :    unit(netCDF)= (unit(RPN)-add)/mult
   - 'nom_attribut' and 'valeur_attribut' are only used for conversions to nsfCDF format.
       - Never put a space before the ";"!
       - "[", "]", "=" and ";" are reserved characters. Do not use them otherwise, for example in names or descriptions.