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