Skip Navigation Links www.nws.noaa.gov 
NOAA logo - Click to go to the NOAA home page National Weather Service   NWS logo - Click to go to the NWS home page
Climate Prediction Center
 
 

 
About Us
   Our Mission
   Who We Are

Contact Us
   CPC Information
   CPC Web Team

 
HOME > Monitoring_and_Data > Oceanic and Atmospheric Data > Reanalysis: Atmospheric Data > ftn_wgrib2api
 

HOME: ftn_wgrib2api

The Easy Way of Reading and Writing Grib2 Using Fortran

Introduction

Can a set of fortran routines make reading and writing grib2 easy? Is it possible? Can you read grib2 without referring to a the WMO grib standard? Can you write grib2 without becoming an grib2 expert? Is reading and writing grib now easier than other popular file formats?

!       read 6 hour forecast from z500 field from file 'FILE.grb2'

        use wgrib2api
	real, allocatable :: grid(:,:)

	iret = grb2_mk_inv('FILE.grb2','FILE.inv')	              ! make index/inventory file, FILE.inv
	if (iret.ne.0) stop 1

!       search index file for terms: ':HGT:500 mb:' and ':6 hour fcst:'
	iret = grb2_inq('FILE,grb2','FILE.inv',':HGT:500 mb:',':6 hour fcst:',data2=grid)
	if (iret.ne.1) stop 2                                         ! error if not 1 match
	write(*,*)'It was easy!  Z500 fhr=6 grid(10,12) =', grid(10,12)

	stop
	end

Reading doesn't get much easier than that. You are allowed upto 20 search terms that are based on the wgrib2 inventory. You don't have to even allocate the variable grid. Now wgrib2api uses optional parameters to make the grb2_inq(..) call both simple and powerful. We can expand the above program to get much more information about the field that was read. Do some testing, is it easier to read grib2, netcdf or hdf5?

!       read 6 hour forecast from z500 field from file 'FILE.grb2'
!       showing the use of some optional parameters

        use wgrib2api
	real, allocatable :: grid(:,:), lat(:,:), lon(:,:)
        character (len=300) :: metadata, grid_info
        integer :: nx, ny

	iret = grb2_mk_inv('FILE.grb2','FILE.inv')	              ! make index/inventory file, FILE.inv
	if (iret.ne.0) stop 1

!       search index file for terms: ':HGT:', ':500 mb:' and ':6 hour fcst:'
	iret = grb2_inq('FILE,grb2','FILE.inv',':HGT:',':500 mb:',':6 hour fcst:',data2=grid, &
          nx=nx, ny=ny, lat=lat, lon=lon, desc=metadata, grid_desc=grid_info)
        if (iret.ne.1) stop 2                                         ! error if not 1 match

        write(*,*)'It was easy!  Z500 fhr=6 grid(10,12) =', grid(10,12)
        write(*,*) 'at latitude=',lat(10,12),' longitude=',lon(10,12)
        write(*,*) 'grid size=(',nx,',',ny,')'
        write(*,*) 'metadata: ',trim(metadata)
        write(*,*) 'grid: ', trim(grid_info)

	stop
	end

Writing grib2 is also easy. First you need a template which is a sample grib2 file. The template has the correct grid and unchanging metadata such as the center and grid definition. The writing process consists of adding the gridded data, the new variable name, level, reference date and forecast information such as analysis, or 12 hour forecast and perhaps some ensemble information or some similar changing metadata.

Finding a template has gotten easier as grib2 has become more common. There are tools to change the grid information (wgrib2) and change specific metadata information (wgrib2). Suppose we already have a template. Writing a grib2 file is as easy as

!       template.grb2: grib2 file, message/record 1 is template
!       fort.11: 10 mb TMP binary data WE:SN order
!       OUT.grb2:  output file, 10 mb TMP in grib2 format

        use wgrib2api
        real, allocatable :: grid(:,:)
        integer, parameter :: nx=360, ny=181

        allocate (grid(nx,ny))
        read(11) grid
        i = grb2_wrt('OUT.grb2','template.grb2',1,data2=grid,meta='d=2001020100:TMP:10 mb:anl:packing=j')
!       writes out 10 mb temperature, reference date = 00Z01Feb2001, analysis, jpeg2000 packing
        if (i.ne.0) stop 8

        stop
        end

Reading and writing is based on text strings rather numbers in a table, so you don't have to refer to the WMO documentation. All the text strings are based on the wgrib2 inventories, so there should be no surprises as long as you are familiar with the wgrib2 inventories.

Here is a practical example, you want to calculate the 1000-500 mb thickness. We are assuming that the input file only contains one Z1000 and Z500 field. We are also assuming the first record can be used as a template.

!	sample program:
!       read z500, z1000
!       write 500-1000 mb  thickness  (z500-z1000)
!
	use wgrib2api

	real, allocatable :: z500(:,:), z1000(:,:)
	character (len=200) :: file, inv, metadata

	file = 'gep19.t00z.pgrb2af180'
	inv = 'gep19.t00z.pgrb2af180.inv'

	iret = grb2_mk_inv(file, inv)		! make an inventory file

!       read z500 and z1000
	iret = grb2_inq(file,inv,':HGT:1000 mb:', data2=z1000)
	if (iret.ne.1) stop 2
	iret = grb2_inq(file,inv,':HGT:500 mb:', data2=z500, desc=metadata)
	if (iret.ne.1) stop 3

        z1000 = z500 - z1000

!	the write needs metadata, going to use metadata from the z500 read

!	new style write, var=.. and level=.. override the z500 metadata
        iret = grb2_wrt('thick.grb',FILE,1,data2=z1000,meta=metadata,var='THICK',level='500-1000 mb')
	if (iret.ne.0) stop 4

!	old style write, convert z500 metadata -> thickness metadata
        write(*,*) 'metadata 0 ', trim(metadata)
        iret = grb2_set_substring(metadata,'THICK',2)
	if (iret.ne.0) stop 5
        iret = grb2_set_substring(metadata,'500-1000 mb',3)
	if (iret.ne.0) stop 6
        write(*,*) 'metadata 1 ', trim(metadata)
        iret = grb2_wrt('thick.grb',FILE,1,data2=z1000,meta=metadata)
	if (iret.ne.0) stop 7

        stop
        end

Tested Systems

  • RH6 gfortran
  • RH6 ifort
  • SUSE ifort
  • Ubuntu 14.04 gfortran
  • Cygwin64 (Windows) gfortran
  • MacOS ifort (beta v2.0.7)

    Source Code

    The source code is included with wgrib2. Please use the newest version of wgrib2 because wgrib2api is a new addition to wgrib2. You need to compile wgrib2api using the same fortran compiler that you will be using to compile your fortran programs. By default, wgrib2api/ftp_api is compiled with OpenMP. So your fortran needs to be compiled with the same OpenMP setting.

    compiling (making the ftn_api)
    Sample Code

    The default wgrib2/wgrib2api build includes the ipolates and netcdf3 libraries. If your fortran code uses a different ipolates or netcdf libraries, you will have to configure the wgrib2's makefile and rebuild the wgrib2 library without the ipolates and netcdf options. Removing the two libraries will not have an effect on the grb2_*(..) routines. However, if the missing libraries would adversely affect calls to wgrib2a(..) and wgrib2c(..) that use the -new_grid or -netcdf options.

    Documentation

    1. Introduction pptx
    2. Making an index file, in order to read grib: grb2_mk_inv(..)
    3. Reading grib: grb2_inq(..)
    4. Scanning grib: grb2_inq(..) with sequential reads
    5. Writing grib: grb2_wrt(..)
    6. Freeing files: grb2_free_file(..) does a flush
    7. undefined grid points
    8. HPC: introduction
    9. HPC: reading/writing memory buffers/files
    10. Sample fortran code: merge 2000 small grids
    11. Example: calculating the average surface to 700 mb relative humidity.

    Other Languages

    1. C_wgrib2api is in development and is based on ftn_wgrib2api
    2. python_wgrib2api is in development. Based on ctypes for ease of installation.
    3. python_wgrib2api first modes: ftn wgrib2api dopplegager
    4. python_wgrib2api second modes: more database approach

  • NOAA/ National Weather Service
    National Centers for Environmental Prediction
    Climate Prediction Center
    5830 University Research Court
    College Park, Maryland 20740
    Climate Prediction Center Web Team
    Page last modified: January 25, 2018, Aug 16, 2018
    Disclaimer Privacy Policy