Fortran 77 Programs
for
Performing ECMWF EMOSLIB Standard Transformations

Data Support Section
Scientific Computing Division
National Center for Atmospheric Research1
Boulder, Colorado

5 July 2004


1The National Center for Atmospheric Research is sponsored by the National Science Foundation.

Contents


Introduction

This document provides Fortran code examples which make use of interpolation routines in the ECMWF EMOSLIB field interpolation software package, providing a basis of data transformations in four categories: The use of the term "interpolation" by ECMWF when referring to EMOSLIB routines is a very broad description indeed which covers the transformations listed above, and in some cases (e.g. spectral coefficients ->) may not involve interpolation at all.

In the code examples, input GRIB (GRIdded Binary) records are assumed, and the gridded output data are also encoded in GRIB. A typical processing sequence begins with opening an input GRIB file with pboben which returns a unique integer identifier for the input file that is subsequently used for the rest of the program. Similarly, an output GRIB file is opened and associated with its own unique integer identifier. The GRIB records of the input file are read sequentially, one record at a time, with pbgrib, which packs the GRIB record, still encoded in GRIB, into an integer array. It is the integer array returned by "pbgrib" that contains the bulk of the data and metadata necessary to the proper functioning of the interpolation routines. The remaining metadata is generated by invoking intout with sets of parameters describing the desired output grid. Once "intout" has been called, all the necessary preparations have been made to perform the transformation using intf3. The subroutine "intf" evaluates all available metadata, and generates output data encoded in GRIB which may be written to the output file using pbwrite. The input and output GRIB files are closed with pbclose after all records have been processed.

For the following programs, we have assumed a series of various input GRIB files containing 500 hPa temperature at 6-hourly resolution spanning September 1957, hence 120 GRIB records. The output is also in GRIB. In addition, generic 1 megabyte4 integer and real arrays are declared to handle all GRIB input, processing, and output. The actual number of bytes that are read, processed, to be written, and so on are signaled by integer counts returned and/or utilized by the EMOSLIB routines. We have found that 1 megabyte arrays are more than sufficient to handle all ERA-40 GRIB records, and alleviates potential difficulties if the arrays are declared with insufficient lengths. Finally, the EMOSLIB library must be installed on your system.


2Reduced Gaussian grids refer to Gaussian grids in which the number of grid points in the east-west direction vary as a function of latitude. Namely, the number of east-west grid points decrease as one approaches the poles. These grids are also known as "thinned" Gaussian grids.

3If one starts with spectral vorticity and divergence, and the desired output is U and V (eastward and northward wind components), then intuvp is used instead of "intf".

41 megabyte, i.e., a length of 262,144 4-byte words

Contents


Spectral coefficients ->

Please note: The following three programs use EMOSLIB to generate U and V from spectral vorticity and divergence. However, we have found that EMOSLIB generates slightly erroneous U and V (eastward and northward wind components, necessarily a vector pair) starting from spectral vorticity and divergence. Thus, use the following three programs with sufficient discretion. Please see the section "Outstanding ECMWF software issues" below for further details. The dependency "subr_check_error_status.f" is required for checking the status of error codes returned by EMOSLIB routines. See below.

Important Notes:

As stated previously, we have found that EMOSLIB generates slightly erroneous U and V (eastward and northward wind components, necessarily a vector pair) starting from spectral vorticity and divergence.

In addition, according to the ERA-40 Archive Plan (2002, see "References" below),

"Basic T159 spherical-harmonic fields will be truncated to T63 before the 2.5° grid values are computed."
To ensure this is achieved, we use "prog_sphhar_trsphh.f" to truncate spectral coefficients from T159 to T63, then "prog_sphhar_2p5reg.f" to transform truncated spectral coefficients to 2.5° grid values for scalars, or "prog_sphvod_uv_2p5reg.f" to transform truncated spectral vorticity and divergence to 2.5° grid values for U and V.

Please see the section "Outstanding ECMWF software issues" below for further details.

Contents


Regular latitude-longitude grids ->

Note that the dependency "subr_check_error_status.f" is required for checking the status of error codes returned by EMOSLIB routines. See below.

Contents


Regular Gaussian grids ->

Note that the dependency "subr_check_error_status.f" is required for checking the status of error codes returned by EMOSLIB routines. See below.

Contents


Reduced Gaussian grids ->

Note that the dependency "subr_check_error_status.f" is required for checking the status of error codes returned by EMOSLIB routines. See below.

Contents


Generic ->

Since a GRIB record used as input contains all the essential metadata describing grid conventions, variables, and the like, it is only necessary to provide information about an output grid that may potentially differ from the input grid. (As we have seen, this is accomplished with "intout" which sets up the desired output grid generated by the principal interpolation routines "intf" and "intuvp".) Thus, a sequence of input GRIB records may be mixed, for example consisting of data represented by spherical harmonics and reduced Gaussian grids, but due to their individual metadata the mixed records may be interpolated to a common grid by simply activating the appropriate set of parameters using "intout".

We provide three programs for generating 2.5° regular latitude-longitude grids, regular N80 Gaussian grids, or reduced N80 Gaussian grids, from generic input GRIB records.

Note that the dependency "subr_check_error_status.f" is required for checking the status of error codes returned by EMOSLIB routines. See below.

Contents


Diagnostics of ECMWF EMOSLIB error codes

Performing diagnostics of error codes returned by subprograms is an essential element of good programming practice. Error diagnostics in the programs listed above are performed by a subroutine which we provide, "check_error_status", which takes as its first argument the name of the EMOSLIB routine (in the form of a character string) just called, and as its second argument the integer error code returned by the EMOSLIB routine. Given these two arguments, the subroutine checks for a nonzero error code, and, if an error is encountered, prints a diagnostic message and stops execution of the program.

Thus far, "check_error_status" handles any errors returned by the EMOSLIB routines "pbopen", "pbwrite", "pbclose", "pbgrib", "intf", "intuvp", "intin", and "intout".

Contents


Decoding GRIB records to get real field data

Frequently it is desirable to decode GRIB records to get real field data that can be used for further processing and manipulation, or to write output in a format other than GRIB. We provide an example program which reads GRIB records, and invokes a subroutine (provided here) to decode the GRIB into real field data. The program also illustrates a simple technique to assign the field data to a two-dimensional latitude-longitude array. In addition, the EMOSLIB function "intf" returns not only an encoded GRIB record representing the transformed grid, but also the corresponding real field data as well. In this case, it is straightforward to assign the field data to a two-dimensional latitude-longitude array:

      integer length_array
      integer nlon
      integer nlat
c
      parameter(length_array=262144)
      parameter(nlon=73)
      parameter(nlon=144)
c
      integer iarray_in(length_array)
      integer iarray_out(length_array)
      integer rarray_in(length_array)
      integer rarray_out(length_array)
c
      real latlon_array(nlat,nlon)
c
      integer intout, intf
      external intout, intf
c
      integer ierror
      integer words_read_into_iarray_out
      integer ilat
      integer jlon
      integer k
      integer ival(80)
      real rval(80)
c
c     ...
c

      ierror = intout( 'form', ival, rval, 'unpacked' )
      call check_error_status('intout',ierror)
      ierror = intf( iarray_in,  length_array,               rarray_in
     +               iarray_out, words_read_into_iarray_out, rarray_out  )
      call check_error_status('intf', ierror)

      k = 0

      do 200 ilat = 1, nlat
        do 300 jlon = 1, nlon
          k = k + 1
          latlon_array(ilat,jlon) = rarray_out(k)
300     continue
200   continue
  
Similarly, the same technique can be applied to the real arrays returned by "intuvp".

Note also that the dependency "subr_check_error_status.f" is required for checking the status of error codes returned by EMOSLIB routines. See above.

Contents


Makefile for Fortran 77 programs for performing ECMWF EMOSLIB standard transformations

We provide a simple makefile that can be used to compile all of the Fortran 77 programs and dependencies listed in this document for performing ECMWF EMOSLIB standard transformations (assuming that the EMOSLIB routines are installed as a static library on your platform). For example, to compile the program for going from spherical harmonics to a 2.5° regular latitude-longitude grid, make sure the dependency "subr_check_error_status.f" (for checking the error status returned by EMOSLIB routines) is in the same directory as "prog_sphhar_2p5reg.f", and do "make sphhar_2p5reg". The executable "sphhar_2p5reg" is then generated in your working directory.

In the makefile notice that we have used the Portland Group Fortran 90 compiler (pgf90) and have linked the EMOS library available in /usr/local/emos-240/lib. Please modify these to identify the compiler and EMOSLIB library path that you are using. (The best practice in the case of EMOSLIB is to use the same compiler that was used to build the EMOSLIB library on your machine.)

Please note that if you cut and paste the makefile, then a tab must be inserted at the beginning of each command line, that is, prior to each instance of "${FORTRAN_COMPILER} ...". This is a cardinal convention of all makefiles.

Contents


Gaussian latitudes, weights, and points

If output Gaussian grids other than regular or reduced N80 Gaussian are required, we provide sufficient metadata in the form of Fortran data statements which initialize real arrays. These can then be substituted in your programs. Following ECMWF EMOSLIB conventions, the metadata is for the northern hemisphere only, arranged north to south, pole to equator. For example, if one requires T31 48x96 (latitude x longitude) Gaussian output, the following statements are in order:

      integer ngau
      integer stp
c
      parameter( ngau = 24 )
      parameter(  stp = 31 )
c
      integer intv(ngau)
      real realv(ngau)
      character*30 charv(ngau)
c
      real lats(ngau)
      real area(4)
c
      data lats /
     +   87.1591,  83.4789,  79.7770,  76.0702,  72.3616,  68.6520,
     +   64.9419,  61.2316,  57.5210,  53.8103,  50.0995,  46.3886,
     +   42.6776,  38.9666,  35.2556,  31.5445,  27.8334,  24.1223,
     +   20.4112,  16.7001,  12.9890,   9.2778,   5.5667,   1.8556 /
c
      data area / 4*0.0 /
c
c     ...
c
      ierror = intout('area',                  intv, area,  charv    )
      call check_error_status('intout', ierror)
c
      ierror = intout('user_regular_gaussian', ngau, realv, charv    )
      call check_error_status('intout', ierror)
c
      ierror = intout('g_lats',                intv, lats,  charv    )
      call check_error_status('intout', ierror)
c
      ierror = intout('truncation',            stp,  realv, charv    )
      call check_error_status('intout', ierror)
  

Contents


Fixed regular latitude-longitude grids other than 2.5°

A "fixed regular latitude-longitude grid" usually refers to a grid in which the first grid point represents 0°E and 90°N (other combinations such as 180°E and 90°S are possible, as long as the first grid point represents a longitude and latitude that are whole integers). Subsequent grid points are equally spaced in latitude and longitude. For example, in the programs listed above in which the output is a 73x144 2.5° regular latitude-longitude grid, the statement

      data grid / 2.5, 2.5 /
  
may be changed to

      data grid / 5.0, 5.0 /
  
to generate a 37x72 5.0° regular latitude-longitude output grid. Some common fixed regular latitude-longitude grids are (and the conventions adhered to in ECMWF EMOSLIB output):

  ------------------------------------------------------------------------------------------------
  latitude spacing   longitude spacing   latitudes   longitudes   first latitude   first longitude
  ------------------------------------------------------------------------------------------------
        1.0°               1.0°             181         360            90°N              0°E 
        1.125°             1.125°           161         320            90°N              0°E 
        2.5°               2.5°              73         144            90°N              0°E 
        5.0°               5.0°              37          72            90°N              0°E 
  ------------------------------------------------------------------------------------------------
  

Contents


Download tarfile

All code and makefiles described in this document are available for download in the form of a tarfile (133120 bytes).

Contents


Outstanding ECMWF software issues

In our work with ECMWF EMOSLIB software we have encountered three unresolved issues that may have a significant impact on users' data processing objectives:
5GNU Fortran 77 , part of the gcc compiler suite.

6See for example the Sun Fortran 77 Language Reference.

7Specifically, a Dell Precision 360 with Pentium IV processor running Red Hat Enterprise Linux WS Version 3.

Contents


References

European Centre for Medium-Range Weather Forecasts, 2002: The ERA-40 Archive. Reading, ECMWF, 40 pp. (On-line and A4 PostScript versions are available from the ECMWF ERA-40 Project Plan link.)

Contents


Web page written and maintained by David Stepaniak
davestep@ucar.edu
Last modified 23rd July 2004.

DSS, the Data Support Section in SCD at NCAR.

Contents