program packgrib_example
	parameter (idim=145,jdim=37)
	integer grib_lbl(43)
	character pds_ext*256
	real gridpoints(idim,jdim)
	parameter (odim=50000)
	integer output_buffer(odim)
	parameter (iwordsz=4)
	parameter (xmissing=-999.)

c GRIB Edition number
	grib_lbl(1)=1
c parameter table version
	grib_lbl(3)=3
c center id
	grib_lbl(4)=60
c generating process ID number
	grib_lbl(5)=255
c include the GDS
	grib_lbl(7)=1
c parameter code
	grib_lbl(8)=11
c level type
	grib_lbl(9)=100
c level
	grib_lbl(10)=500
	grib_lbl(11)=0
c date
	grib_lbl(12)=2004
	grib_lbl(13)=6
	grib_lbl(14)=21
	grib_lbl(15)=1625
c time unit
	grib_lbl(16)=1
c periods of time
	grib_lbl(17)=0
	grib_lbl(18)=0
c time range indicator
	grib_lbl(19)=10
c for averaged grids
	grib_lbl(20)=0
	grib_lbl(21)=0
c sub-center id
	grib_lbl(22)=0
c decimal scale factor
	grib_lbl(23)=0
c binary data section flag
	grib_lbl(24)=0
c packing width of data points
	grib_lbl(25)=16
c length of GDS
	grib_lbl(26)=32
c grid (data representation) type
	grib_lbl(29)=0
c latitude/longitude grid
	grib_lbl(30)=idim
	grib_lbl(31)=jdim
	grib_lbl(32)=0
	grib_lbl(33)=0
	grib_lbl(34)=0
	grib_lbl(35)=90000
	grib_lbl(36)=360000
	grib_lbl(37)=2500
	grib_lbl(38)=2500
	grib_lbl(39)=64

c add a comment to the PDS
	pds_ext='This is a PDS comment'
	iplen=21;

c fill the gridpoints array with randomly generated temperatures
	ii=1
	do i=1,145
	  do j=1,37
	    ii=mod((ii*i*j),2999)
	    if (ii .eq. 0) ii=1
	    xx=ii*0.01
c make two of the data points missing
	    if (i .eq. 3 .and. j .eq. 17 .or. i .eq. 27 .and. j .eq. 4) then
		gridpoints(i,j)=xmissing
	    else
		gridpoints(i,j)=253.15+xx
	    endif
	  enddo
	enddo

c pack the grid
	call packgrib(grib_lbl,pds_ext,iplen,gridpoints,idim,xmissing,
     +output_buffer,odim,length)
c write the grid to an output file
	length=(length+iwordsz-1)/iwordsz
	write(12) (output_buffer(i),i=1,length)
	end