Post processing with gl

Introduction

gl ( as in griblist ) is a multi purpose tool for file manipulation and conversion. It uses ECMWF's ecCodes library, and can be compiled with and without support for HARMONIE FA/LFI or NETCDF files. The gl package also includes software for extraction for verification, fldextr, and field comparison, xtool.

 USAGE: gl file [-n namelist_file] [-o output_file] -[lfgmicp(nc)sdtq] [-lbc CONF]

 gl [-f] file, list the content of a file, -f for FA/lfi files  
 -c    : Convert a FA/lfi file to grib ( -f implicit )          
 -p    : Convert a FA file to grib output without extension zone
         (-c and -f implicit )                                  
 -nc   : Convert a FA/lfi file to NetCDF ( -f implicit )        
 -musc : Convert a MUSC FA file ASCII ( -c implicit )           
 -ufn  : Use the FA name as output name in NETCDF               
 -lbc ARG : Convert a CONF file to HARMONIE input               
            where CONF is ifs or hir as in ECMWF/HIRLAM data    
         climate_aladin assumed available                       
 -d    : Together with -lbc it gives a (bogus) NH boundary file   
         climate_aladin assumed available                       
 -s    : Work as silent as possible                             
 -g    : Prints ksec/cadre/lfi info                             
 -m    : Prints min,mean,max of the fields                      
 -i    : Prints the namelist options (useless)                  
 -tp   : Prints the GRIB parameter usage                        
 -t    : Prints the FA/lfi/GRIB table (useful)                  
 -wa   : Prints the atmosphere FA/NETCDF/GRIB table in wiki fmt 
 -ws   : Prints the surfex FA/NETCDF/GRIB table in wiki fmt     
 -q    : Cross check the FA/lfi/GRIB table (try)                
 -pl X : Give polster_projlat in degrees                        

 gl file -n namelist_file : interpolates file according to      
                            namelist_file                       
 gl -n namelist_file : creates an empty domain according to     
                       specifications in namelist_file          
 -igd  : Set lignore_duplicates=T                               
 -igs  : Set lignore_shortname=T. Use indicatorOfParameter      
             instead of shortName for selection                 

ecCodes definition tables

Since ecCodes has replaced grib_api as the ECMWF primary software package to handle GRIB, we will hereafter only refer to ecCodes but same similar settings applies for grib_api as well. With the change to ecCodes we heavily rely on the shortName key for identification. To get the correct connection between the shortnames and the GRIB1/GRIB2 identifiers we have defined specific tables for harmonie. These tables can be found in /util/gl/definitions. To use these tables you have to define the ECCODES_DEFINITION_PATH environment variable as

export ECCODES_DEFINITION_PATH=SOME_PATH/gl/definitions:PATH_TO_YOUR_ECCODES_INSTALLATION

If this is not set correctly the interpretation of the fields may be wrong.

GRIB/FA/LFI file listing

Listing of GRIB/ASIMOF/FA/LFI files.

 gl [-l] [-f] [-m] [-g] FILE

where FILE is in GRIB/ASIMOF/FA/LFI format

OptionDescription
-linput format is LFI
-finput format is FA
-l and -f are equivalent
-gprint GRIB/FA/LFI header
-mprint min/mean/max values

GRIB/FA/LFI file conversion

Output to GRIB1

gl [-c] [-p] FILE [ -o OUTPUT_FILE] [ -n NAMELIST_FILE]

where

-cconverts the full field (including extension zone) from FA to GRIB1
-pconverts field excluding the extension zone ("p" as in physical domain) from FA to GRIB1

The FA/LFI to GRIB mapping is done in a table defined by a util/gl/inc/trans_tab.h

To view the table:

gl -t
gl -tp

To check for duplicates in the table:

gl -q

The translation from FA/LFI to GRIB1 can be changed through a namelist like this one:

  &naminterp
    user_trans%full_name ='CLSTEMPERATURE',
    user_trans%t2v       = 253,
    user_trans%pid       = 123,
    user_trans%levtype   = 'heigthAboveGround',
    user_trans%level     = 002,
    user_trans%tri       = 000,
  /

or for the case where the level number is included in the FA name

  &naminterp
    user_trans%full_name='SNNNEZDIAG01',
    user_trans%cpar='S'
    user_trans%ctyp='EZDIAG01',
    ...
  /

Conversion can be refined to convert a selection of fields. Below is and example that will write out

  • T (shortname='t',pid=011), u (shortname='u',pid=033) andv (shortname='v',pid=034) on all (level=-1) model levels (levtype='hybrid')
  • T (shortname='t',pid=011) at 2m (lll=2) above the ground (levtype='heightAboveGround') [T2m]
  • Total precipitation (shortname='tp',pid=061,levtype='heightAboveGround',level=000)
  &naminterp
   readkey%shortname=   't',     'u',     'v',                't',               'tp',               'fg',
   readkey%levtype='hybrid','hybrid','hybrid','heightAboveGround','heightAboveGround','heightAboveGround',
   readkey%level=        -1,      -1,      -1,                  2,                  0,                 10,
   readkey%tri =          0,       0,       0,                  0,                  4,                  2,
  /

where

  • shortname is the ecCodes shortname of the parameter
  • levtype is the ecCodes level type
  • level is the GRIB level
  • tri means timeRangeIndicator and is set to distinguish between instantaneous, accumulated and min/max values.

The first three ones are well known to most users. The time range indicator is used in HARMONIE to distinguish between instantaneous and accumulated fields. Read more about the options here Note that for levtype hybrid setting level=-1 means all.

We can also pick variables using their FA/lfi name:

  &naminterp
    readkey%faname = 'SPECSURFGEOP','SNNNTEMPERATURE',
  /

Where SNNNTEMPERATURE means that we picks all levels.

Fields can be excluded from the conversion by name

  &naminterp
    exclkey%faname = 'SNNNTEMPERATURE'
  /

Output to GRIB2

To get GRIB2 files the format has to be set in the namelist as

  &naminterp
    output_format = 'GRIB2'
  /

The conversion from FA to GRIB2 is done in gl via the ecCodes tables. All translations are defined in util/gl/scr/harmonie_grib1_2_grib2.pm where we find all settings required to specify a parameter in GRIB1 and GRIB2.


  tmax => {
   editionNumber=> '2',
   comment=> 'Maximum temperature',
   discipline=> '0',
   indicatorOfParameter=> '15',
   paramId=> '253015',
   parameterCategory=> '0',
   parameterNumber=> '0',
   shortName=> 'tmax',
   table2Version=> '253',
   typeOfStatisticalProcessing=> '2',
   units=> 'K',
  },

To create ecCodes tables from this file run

   cd gl/scr
   ./gen_tables.pl harmonie_grib1_2_grib2

and copy the grib1/grib2 directories to gl/definitions.

Note that there are no GRIB2 transations yet defined for the SURFEX fields!

Output to NetCDF

gl -nc [-p] FILE [ -o OUTPUT_FILE] [ -n NAMELIST_FILE] [-ufn]

where

-pconverts field excluding the extension zone ("p" as in physical domain) from FA to NetCDF
-ooutput file name
-nnamelist file to be used in conversion
-ufnuse full name: use FA field names instead of NetCDF names

The FA/LFI to NetCDF mapping is done using tables defined by util/gl/inc/trans_tab.h and util/gl/inc/nc_tab.h. With the -ufn flag, FA field names are used instead of the NetCDF names defined in nc_tab.h.

Namelist options for NetCDF

The translation from FA/LFI to NetCDF can be changed through a namelist like this one:

  &naminterp
    user_nctrans%full_name ='SFX.SIC',
    user_nctrans%s_name    = ""       
    user_nctrans%l_name    = "Sea-Ice Area Percentage (Atmospheric Grid)",
    user_nctrans%unit      = "%"                 ,
  /

The unit entry can be used to do a limited set of unit conversions, in the example above SIC will be converted from the original units (fraction) to a percentage.

Other specific naminterp options for converting to netcdf:

variabledescriptiondefault
lacc2fluxT: Convert accumulated fields (tri=4) to fluxes by dividing by the length of the interval FF-F0. Names, units and tri are adapted..FALSE.
lmergelevsT: write all levels of a variable to the same file; F: each level in a separate file.FALSE.
lclimate_fieldsT: don't add a time dimension and variable, useful for climate fields like land mask.FALSE.
lverticesT: add vertices (corner points) to netcdf file, only possible for newly created files.FALSE.
lhistoryT: add history global attribute to netcdf file.FALSE.
ikindncNetCDF version (34), 3: larger files, but faster; 4: compressed, but slow
ref_dateReference date, used to generate relative time axis19500101
ref_hourReference hour, used to generate relative time axis0
ctimeistime refers to "start", "middle", or "end" of interval for non-instantaneous fields. If writing several variables to 1 file, that don't have the same timing (e.g. accumulated vs. instantaneous), then "end" is probably the only safe option!end
csepseparator in derived netcdf file name_
cdatefnameformat for date in derived netcdf file name, if not recognized as format, use whatever is passedYYYYMMDDHH
cfidenused in derived netcdf file name to indicate origin (e.g. his, sfx, fp)
cfreqused in derived netcdf file name and as "frequency" global attribute (e.g. 1hr, 3hr, day, mon)
chm_revHARMONIE version, used as "model_id" global attribute
cdomaindomain name, used in derived netcdf file name and as "domain" global attribute
cexperimentexperiment id, used in derived netcdf file name and as "experiment_id" global attribute
cinstituteused as "institute_id" global attribute
chostmodused as "drivingmodelid" global attribute
Setting fstart for min-max fields

Min and max fields, with tri=2 (time range indicator) are valid for a certain period. By default the period is 3h, but this can be changed via variables FREQ_RESET_TEMP and FREQ_RESET_GUST in ecf/config_exp.h, for example to 1 to store min/max temperature over an hour. By default gl doesn’t have info on this frequency and it is assumed they are valid since the start of the run. Use the namelist option fstart to assign the appropriate starting value, e.g.:

fstart(15) = 3

for max t2m (parameter code 15). The fstart value is then used in the time_bnds. This value needs to be updated with FREQ_RESET. In Makegrib_gribex and convertFA there are examples of how to do this. Note that this works in the same way for NetCDF and GRIB.

Derived file name

If no output file name is supplied (-o flag) an output file name is derived from available info:

nc_fld_name[_levinfo][_cfiden][_cdomain][_cexperiment][_cfreq][_timeinfo].nc

with: | | | | –- | –- | | ncfldname | name of netcdf variable as defined in nctab.h | | levinfo | indication of level info, with *lev if all levels are written to the same file (lsplitlev), or a level number/height otherwise | | cfiden | identifier of input file (e.g. his, sfx), set via namelist | | cdomain | domain name, set via namelist | | cexperiment | experiment name, set via namelist | | cfreq | frequency, set via namelist | | timeinfo | indicator of file date/time, format controlled via cdatefname namelist variable, not used if lclimatefields |

If elements are not set, either via the namelist or by a default value, they are excluded from the name. The separator is a _ by default but can be changed via the csep namelist variable.

Time axis

A relative time axis (days since …) is created using the refdate and refhour as reference. In existing files the current time step is looked up, and if it already exists it is overwritten. If it doesn’t exist yet, it is appended to the end of the file. Note that if time steps are not converted in order the time axis will not be consecutive. For non-instantaneous fields a time bounds variable is added. The start of the interval is taken from fstart (from tri=2) or outkey%F0, or just from the input file (usually 0). The end of the interval is the current time step. Whether the time variables refers to the beginning, middle or end of the interval can be controlled with the ctimeis namelist variable. If instantaneous and non-instantaneous files are written to the same file, it is best to use ctimeis=end. Start may also work, but this should be tested first. With the namelist variable cdatefname you can write output from multiple cycles to the same file. For example by setting it to YYYYMM the derived file name will contain year and month info, but not day and hour, so all cycles from a month are written to the same file. Be careful with the first time step of a cycle when using cdatefname as gl will overwrite the last time step of the previous cycle with those of the first step of the new cycle. You can decide to skip the first time step, or multiple steps, if cycles overlap more than 1 step.

Multilevel fields and fields on heights or pressure levels

All levels of multilevel fields can be written to one file (lmergelevs=.TRUE.) or to separate files (default). This is possible for model levels, pressure levels and height levels.

All levels in one file

If all levels are written to the same file and additional dimension is needed: lev for model levels, height for height levels. For model level fields, named SNNN, the number of levels is derived from the input file (glistnlev). For height level fields, HNNNNN, currently the heights must be set via the hlevlist namelist variable. The heights in this list are used to expand the HNNNNN (to H00010, H00250 etc) and are also used as coordinate variable. For pressure levels PNNNNN is used in the same way.

Single level fields on height

For single level fields on a specific height such as 2m temperature and 10m wind a height variable is added with that height. This is done if level type = 105 and level ≠ 0. This may not be appropriate in all cases. Note that for some fields level is abused (e.g. level 760 for the sea tile), which gives useless height. The same approach is used when outputting multilevel fields with lsplitlev=.TRUE. (default).

Don’t mix fields

There is a check that all fields on height levels in one file have the same height specification, because only 1 height variable can be specified at the moment. So t2m and w10m cannot be in the same file. Not sure if the check is foolproof. It may be possible to define multiple heights in the code, e.g. height, height2 etc., but this has not been implemented yet.

NetCDF 3 or NetCDF 4

With the ikindnc namelist options the netcdf format can be set. NetCDF 4 files are compressed in gl, however, this makes the conversion much slower. At the moment it seems better to let gl use the NetCDF 3 format and then convert them to NetCDF 4 after creation of the file has finished. This can be done with the following command:

nccopy -k 4 -d 1 -s $nc3_file $nc4_file

where:

-ddeflate level
-sshuffling (can improve compression, speed and size)

In scr/convertFA this can be done by setting nc3to4=yes (default). At the end of the script the files are then converted from netcdf3 to netcdf4-classic with compression.

Direction of fluxes

A new element positive was added to the nctrans derived type in moduletypez, and in nctab.h positive, can be empty, 'd' or 'u'. If it is not empty a positive attribute is added to the variable in the NetCDF file. If it has value 'u', the values of the variable are multiplied by −1 to change the direction from towards the surface to away from the surface.

Fill value

By default no missing value is added for atmospheric fields. For SURFEX fields either 999 (version ≤ 4) or 1+e20 (version ≥ 5) is used. It is possible to add a missing value via the namelist. To do so, in the namelist set variable lcheck_misval to .TRUE. and set rmisval to the correct value.

Adding new netcdf variables

If you get messages like:

No NETCDF conversion for ....

then you need to add the field to util/gl/inc/nc_tab.h, which contains the translation from FA to netcdf names. The file util/gl/inc/trans_tab.h contains the conversion to FA names to GRIB codes. If the field you would like to is absent there, it is probably best to add it in that file as well, as for example GRIB level types are used for functionality in the netcdf conversion as well. Remember to recompile.

postprocessing

gl can be used to produce postprocessed parameters possibly not available directly from the model.

  • Postprocessed parameters are defined in util/gl/grb/postprocess.f90 and util/gl/grb/postp_pressure_level.f90. Some more popular parameters are listed:
    • Pseudo satellite pictures
    • Total precipitation and snow
    • Wind (gust) speed and direction
    • Cloud base, cloud top, cloud mask and significant cloud top

For a comprehensive list please check the output information for each cycle. NOTE that all parameters may not be implemented in gl

  • To produce "postprocessed" MSLP and accumulated total precipitation and visibility use the following namelist, nam_FApp:

    &naminterp
     pppkey(1:3)%shortname='pres','tp','vis',
     pppkey(1:3)%levtype='heightAboveSea','heightAboveGround','heightAboveGround'
     pppkey(1:3)%level=  0, 0, 0,
     pppkey(1:3)%tri=  0, 4, 0,
     lwrite_pponly= .TRUE.,
    /
    gl -p ICMSHHARM+0003 -o output_pp.grib -n nam_FApp
    • Note:
      • Set lwrite_pponly as true to only write the postprocessed fields to file
      • Set lwrite_pponly as false write all fields will be written to the file, input fields as well as the postprocessed fields.

Vertical interpolation

gl can be used to carry out vertical interpolation of parameters. Four types are available

  • HeightAboveSea, give height above sea in meters

  • HeightAboveGround, give height above ground in meters

  • HeightAboveGroundHighPrecision, give height above ground in centimeters

  • isobaricInHpa, give height above sea in hPa

  • To interpolation temperature to 1.40m (level 140 in cm) use the following namelist, nam_hl:

    &naminterp
     pppkey(1:1)%shortname='t',
     pppkey(1:1)%levtype='heightAboveGroundHighPrecision',
     pppkey(1:1)%level=  140,
     pppkey(1:1)%tri=  0,
     vint_z_order=1,
     lwrite_pponly= .TRUE.,
    /
    gl -p ICMSHHARM+0003 -o output_hl.grib -n nam_hl
    • Note:
      • Vertical interpolation to z levels is controlled by VINTZORDER: 0 is nearest level, 1 is linear interpolation
  • To height interpolation (Levls 500, 850 and 925 in hPa, type=100) use the following namelist, nam_pl:

    &naminterp
     pppkey(1:3)%shortname='t','t','t',
     pppkey(1:3)%levtype='isobaricInhPa','isobaricInhPa','isobaricInhPa',
     pppkey(1:3)%level=  500, 850, 925,
     pppkey(1:3)%tri=  0, 0, 0,
     vint_z_order=1,
     lwrite_pponly= .TRUE.,
    /
    gl -p ICMSHHARM+0003 -o output_pl.grib -n nam_pl

Horizontal interpolation

  • Interpolation/resampling between different geometries such as regular lat lon, Lambert conformal, Polar steregraphic, rotated lat lon and rotated Mercator is possible with gl

  • The interpolation methods available are:

    • nearest grid-point (order=-2)
    • most representative grid-point (order=-1)
    • nearest grid-point (order=0)
    • bi-linear (order=1)
    • bi-quadratic (order=2, mask not respected)
    • bi-cubic (order=3, mask not respected)
  • Example of (an Irish) rotated lat lon domain, nam_FArotll:

    &naminterp
     outgeo%nlon=50,
     outgeo%nlat=50,
     outgeo%nlev=-1,
     outgeo%gridtype='rotated_ll',
     outgeo%west=-2.5,
     outgeo%south=-2.5,
     outgeo%dlon=0.1,
     outgeo%dlat=0.1,
     outgeo%polon=-6.7,
     outgeo%polat=-36.2,
     order= 1,
    /

where DLON/DLAT are in degrees.The HIRLAM Domain Tool may be of use for viewing rotated lat lon domains.

gl -p ICMSHHARM+0003 -n nam_FArotll  -o output.grib
  • Example of a lambert domain

    &naminterp
      outgeo%nlon       =  50 ,
      outgeo%nlat       =  50,
      outgeo%nlev       =  -1,
      outgeo%gridtype   =  'lambert',
      outgeo%west       =  15.0
      outgeo%south      =  50.0
      outgeo%dlon       = 10000.
      outgeo%dlat       = 10000.
      outgeo%projlat    =  60.
      outgeo%projlat2   =  60.
      outgeo%projlon    =  15.
    /

    where DLON/DLAT are in meters.The HIRLAM Domain Tool may be of use for viewing rotated lat lon domains.

  • Example polar stereographic projection

    &naminterp
      outgeo%nlon       =  50 ,
      outgeo%nlat       =  50,
      outgeo%nlev       =  -1,
      outgeo%gridtype   =  'polar_stereographic',
      outgeo%west       =  15.0
      outgeo%south      =  50.0
      outgeo%dlon       = 10000.
      outgeo%dlat       = 10000.
      outgeo%projlat    =  60.
      outgeo%projlon    =  15.
    /

    where DLON/DLAT are in meters.Note: the GRIB1 standard assumes that the projection plane is at 60 degrees north whereas HARMONIE assumes it is at 90 degrees north.

  • Example rotated Mercator

    &naminterp
      outgeo%nlon       =  50 ,
      outgeo%nlat       =  50,
      outgeo%nlev       =  -1,
      outgeo%projection =  11,
      outgeo%west       =  15.0
      outgeo%south      =  50.0
      outgeo%dlon       = 10000.
      outgeo%dlat       = 10000.
      outgeo%projlat    =  60.
      outgeo%projlon    =  15.
    /

    where DLON/DLAT are in metersNote: rotated Mercator is not supported in GRIB1.

  • Geographical points is a special case of projection 0 use namelist file, nam_FAgp:

    &naminterp
      outgeo%nlon=3 ,
      outgeo%nlat=1,
      outgeo%nlev=-1,
      outgeo%gridtype='regular_ll',
      outgeo%arakawa=  'a',
      order             =   0,
      readkey(1:3)%shortname='t','u','v',
      readkey(1:3)%levtype='heightAboveGround','heightAboveGround','heightAboveGround',
      readkey(1:3)%level=   2,  10, 10,
      readkey(1:3)%tri=  0, 0, 0,
      linterp_field     = f,
      gplat          = 57.375,57.35,57.60
      gplon          = 13.55,13.55,14.63
    /

    The result will be written to a ASCII file with the name gpYYYYMMDDHHLLL.

    gl -p ICMSHHARM+0003 -n nam_FAgp 
    cat gp20140702_1200+003

Extract (crop) a sub-domain

gl can be used to "cut out" a sub-domain from an input file using the namelist namCUT:

Crop using lower left and upper right coordinates

&naminterp
istart = 150
jstart = 150
istop = 350
jstop = 350
/

Use this command:

gl input.grib -n namCut -o cutout.grib

Another way of specifying your sub domain is to define how many points to exclude in the end

&naminterp
istart = 150
jstart = 150
istop = -10
jstop = -10
/

Crop using SW,NE corner and/or number of points

Here, you specify any of the SouthWest, NorthEast corners and/or the number of gridpoints

&naminterp
outgeo%gridtype = 'crop',
outgeo%nlat  = 200,
outgeo%nlon  = 300,
outgeo%south =  50.155,
outgeo%west  = -12.88,
/

Or

&naminterp
outgeo%gridtype = 'crop',
outgeo%north =  58.277,
outgeo%east  =  12.3,
outgeo%south =  50.155,
outgeo%west  = -12.88,
/

If you specify outgeo%gridtype as 'crop', the SouthWest corner will be translated to lower left grid coordinates and Nlat,Nlon will translate to upper right coordinates. You may specify any of SW, NE, nlat/nlon. Priority is given to SW, NE. The behaviour is as follows:

SW and NE have priority, these will anchor either corner. If a corner is not specified, Nlat/Nlon will extend from the other corner.

If only one coordinate is specified, the other corner becomes the corner of the input domain. So:

  1. Specify only SW and you get a crop from there to the NE corner of the input domain.
  1. Specify only NE and you get a crop from the SW corner of the input domain.

  2. Specify SW and NE and you get a crop between these corners

  1. Specify SW and Nlat/Nlon and you get Nlat x Nlon from SW corner.
  1. Specify NE and Nlat/Nlon and you get Nlat x Nlon south and west of NE corner.
  1. Specify SW, NE and Nlat/Nlon and you get a crop between SW/NE corners. Nlat/Nlon are ignored.
  1. Specify only Nlat/Nlon and you get Nlat x Nlon from SW corner of the input domain.

The crop must be within the original domain unless you set

ldemand_inside = .FALSE.

in the namelist. In this case, the crop will be adjusted to lie within the original domain and the output will be smaller than Nlat x Nlon. In the case where the requested crop lies entirely outside the original domain, the program will abort.

Rotating wind components

HARMONIE produces u and v wind components relative to the model grid. gl by default always outputs u and v relative to the output grid. So if no regridding is done and the output is still on the LCC grid, u and v will also still be relative to the LCC grid. But if the output is regridded to a regular lat-lon grid, then u and v will be rotated and will be relative to the regular lat-lon grid. Wind (from) direction (parameter 31), however, is always relative to a regular lat-lon grid. To rotate u and v to regular lat-lon, while retaining the data on the LCC grid set uvrelativetogrid=0 in the namelist. All u and v vectors that will be processed will be rotated to geographical E and N directions.

Output to several files

It is possible to let gl read data once and do processing loops with these data. Let us look at an example namelist

&naminterp
 OUTPUT_FORMAT='MEMORY'
/
&naminterp
 INPUT_FORMAT='MEMORY'
 OUTPUT_FORMAT='GRIB'
 OUTFILE='test1.grib'
/
&naminterp
 INPUT_FORMAT='MEMORY'
 OUTPUT_FORMAT='GRIB'
 OUTFILE='test2.grib'
 READKEY%FANAME='SNNNTEMPERATURE'
/
&naminterp
 INPUT_FORMAT='MEMORY'
 OUTPUT_FORMAT='GRIB'
 READKEY%FANAME='CLSTEMPERATURE'
 outgeo%nlon       =  50 ,
 outgeo%nlat       =  50,
 outgeo%nlev       =  -1,
 outgeo%gridtype   =  'polar_stereographic',
 outgeo%west       =  15.0
 outgeo%south      =  50.0
 outgeo%dlon       = 10000.
 outgeo%dlat       = 10000.
 outgeo%projlat    =  60.
 outgeo%projlon    =  15.
 OUTFILE='test3.grib'
/

In the first loop we read data and store it in memory. In the second look we read the data from memory and output to the file test1.grib. Then we make two more loops where we in the first one only output a subset and in the last one also do an interpolation to a new grid. The data in memory is however still untouched.

Input from several files

It's also possible to read several files and write them into one. This is used to gather the various FA fields written from the IO-server. A typical namelist would look like

&naminterp
 maxfl=28,
 output_format='MEMORY',
 output_type = 'APPEND',
 input_format='FA',
 infile='forecast/io_serv.000001.d/ICMSHHARM+0003.gridall',
/
&naminterp
 output_format='MEMORY',
 output_type = 'APPEND',
 input_format='FA',
 infile='forecast/io_serv.000002.d/ICMSHHARM+0003.gridall',
/
...
&naminterp
 input_format = 'MEMORY',
 output_format= 'GRIB'
 output_type  = 'NEW',
 outfile      = 'test.grib'
/

Where maxfl tells how many files that will be read.

domain_prop

domain_prop is used do extract various properties from a file.

Climate:  $MPPGL $BINDIR/domain_prop -DOMAIN_CHECK $LCLIMDIR/m$M1 -f || \

Check an existing domain with a namelist specification

domain_prop -DOMAIN_CHECK -f CLIMATE_FILE

The geometry is read from fort.10 and the program aborts if the new and old geometries differs. See scr/Climate for an example.

Check if Q is in gridpoint or spectral representation

domain_prop -f -QCHECK FAFILE

returns 1 if Q is spectral and 0 if Q is in gridpoint space.

Check if a specific field is present

domain_prop -f -CHECK_FIELD S001CLOUD_FRACTI

returns 1 if S001CLOUD_FRACTI is found, 0 otherwise

Check the number of levels in a file

domain_prop -f -NLEV FAFILE  

Check the geographical extension of the domain

domain_prop -f -MAX_EXT FAFILE  

This is used in several places to determine the domain to be extracted from MARS or limit the observations sample. Another way is to provide the projection parameters of your domain as input

domain_prop -MAX_EXTR \
-NLON $NLON -NLAT $NLAT \
-LATC $LATC -LONC $LONC \
-LAT0 $LAT0 -LON0 $LON0 \
-GSIZE $GSIZE

To get the geographical position of the lower left corner use

domain_prop -f -LOW_LEFT FAFILE  

To print out the important projection parameters in a file use:

domain_prop -f -4JB FAFILE

Get time information from a file

domain_prop -f -DATE FAFILE

fldextr and obsextr

Read about the verification extraction programs here