program test_readgrib2 ! !----------------------------------------------------------------- ! ABSTRACT: This routine is to write out a new grib2 file ! Mar 2013 J.Wang !----------------------------------------------------------------- ! implicit none ! integer, parameter :: max_bytes=20000000 integer, parameter :: nx=192,ny=94 ! integer listsec0(2) integer listsec1(13) integer igds(5) integer igdstmpllen integer ipdstmpllen integer idrstmpllen integer idrsnum,ibmap,numcoord,ipdsnum,idefnum integer,dimension(100) :: igdstmpl integer,dimension(100) :: ipdstmpl integer,dimension(100) :: idrstmpl ! integer ideflist(1) real(4) coordlist(1) ! character(255) :: cout character*1 cgrib(max_bytes) ! logical*1 bmap(nx,ny) ! real(4),dimension(nx*ny) :: fld integer ifilw,i,j,lengrib,lonstt,lonlst,latstt,latlst,ierr integer yy,mm,dd,hh,mn,sc real(4) :: dxval ! ! code start !----------------------------------------------------------------- ! !-- set file unit ifilw=52 ! !-- get file name call getarg(1,cout) ! !-- Open GRIB2 file call baopen(ifilw,trim(cout),ierr) print *,'cout=',trim(cout),'ierr=',ierr ! !-- section 0: listsec0(1)=0 ! Discipline: table 0.0 listsec0(2)=2 ! grib edition number ! !-- section 1: listsec1(1)=7 ! Identification of orginating center (Table 0) (7:ncep) listsec1(2)=4 ! Identification of orginating subcenter (ON388-Table C) (4:emc) listsec1(3)=11 ! GRIB master tables version number (Table 1.0) (11:May 2013 version) listsec1(4)=1 ! Version number of GRIB local tables used to augment Master Tables (Table 1.1) listsec1(5)=1 ! Significance of reference time (Table 1.2) (0:ana 1:fcst 2:vrfy) listsec1(6)=2013 ! Reference time - Year (4 digits) listsec1(7)=07 ! Reference time - Month listsec1(8)=15 ! Reference time - Day listsec1(9)=00 ! Reference time - Hour listsec1(10)=00 ! Reference time - Minute listsec1(11)=00 ! Reference time - Second listsec1(12)=0 ! Production status of data (Table 1.3) (0:opn products 1:opn test products) listsec1(13)=1 ! Type of processed data (Table 1.4) (0:ana products 1:fcst products 2:ana & fcst 3: cntl fcst) call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr) print*,'gribcreate status=',ierr ! !-- section 3: grid definition section igds(1)=0 ! Source of grid definition (Table 3.0) (0:specified in the code) igds(2)=nx*ny ! Number of grid points in the defined grid igds(3)=0 ! Number of octets for optional list of numbers defining number of points igds(4)=0 ! Interpretation of list of numbers defining number of points !-- example: Gaussian Lat/lon igds(5)=40 ! Grid definition template number (Table 3.1) (0:Lat/lon, 30:Lambert 40:Gaussian) if( igds(5)==40) then igdstmpllen=19 igdstmpl=0 ! !-- set up grid definition template 3.40 igdstmpl=0 igdstmpl(1)=6 ! Shape of the Earth (Table 3.2) (6:Shape of the Earth = 6,371,229.0 m) igdstmpl(8)=nx ! Ni . number of points along a paralell igdstmpl(9)=ny ! Nj . number of points along a meridian igdstmpl(10)=0 ! Basic angle of the initial production domain igdstmpl(11)=0 ! Subdivisions of basic angle used to define extreme longitudes and latitudes, and direction increments latstt=88541961 lonstt=0. latlst=-88541961 lonlst=358125000 dxval=1875000 igdstmpl(12)=latstt ! La1 - latitude of first grid point igdstmpl(13)=lonstt ! Lo1 - longitude of first grid point igdstmpl(14)=48 ! Resolution and component flags (Table 3.3) igdstmpl(15)=latlst ! La2 - latitude of last grid point igdstmpl(16)=lonlst ! Lo2 - longitude of last grid point igdstmpl(17)=dxval ! Di - i direction increment igdstmpl(18)=ny/2 ! N - number of paralells between a pole and the equator igdstmpl(19)=0 ! Scanning mode (Table 3.4) (0:Points in the first row or column scan in the +i (+x) direction) endif ! idefnum=1 ! Used if igds(3) .ne. 0. The number of entries in array ideflist ideflist=0 ! Used if igds(3) .ne. 0. number of grid points contained in each row ( or column ), Dummy array otherwise call addgrid(cgrib,max_bytes,igds,igdstmpl,igdstmpllen,ideflist,idefnum,ierr) print*,'addgrid status=',ierr ! !-- section 4: product definition section ipdstmpl=0 ipdsnum=0 ! Product Definition Template Number (Table 4.0) (0: Analysis or forecast at a horizontal level or in a horizontal layer at a point in time) ipdstmpllen=15 ! pdt template length ipdstmpl(1)=0 ! catogory ipdstmpl(2)=0 ! parameter ipdstmpl(3)=2 ! Type of generating process (Table 4.3) (0:ana, 1:ic, 2:fcst) ipdstmpl(4)=0 ! Background generating process identifier ipdstmpl(5)=81 ! Analysis or forecast generating process identified (ON388TableA) ipdstmpl(6)=0 ! Hours of observational data cutoff after reference time ipdstmpl(7)=0 ! Minutes of observational data cutoff after reference time ipdstmpl(8)=1 ! Indicator of unit of time range (Table 4.4) (0:minute, 1:hour 2:day) ipdstmpl(9)=0 ! Forecast time in units defined by ipdstmpl(8) ipdstmpl(10)=100 ! Type of first fixed surface (see Code table 4.5) (100:isobaric level) ipdstmpl(11)=0 ! Scale factor of first fixed surface ipdstmpl(12)=85000 ! Scaled value of first fixed surface ipdstmpl(13)=255 ! Type of first second surface (see Code table 4.5) (100:isobaric level) ipdstmpl(14)=0 ! Scale factor of second fixed surface ipdstmpl(15)=0 ! Scaled value of second fixed surface ! numcoord=0 ! Number of coordinate values after template coordlist=0. ! Optional list of coordinate values ! !-- section 5: Data Representation Section idrstmpl=0 idrsnum=3 ! Data representation section template number (Table 5.0) (3:Grid Point Data - Complex Packing and Spatial Differencing) idrstmpllen=18 ! Length of Data representation section idrstmpl(2)=0 ! Binary scale factor idrstmpl(3)=2 ! Decimal scale factor idrstmpl(4)=16 ! Decimal scale factor idrstmpl(7)=0 ! Missing value management used (see Code Table 5.5) idrstmpl(8)=0 ! Primary missing value substitute idrstmpl(9)=0 ! Secondary missing value substitute idrstmpl(17)=2 ! Order of spatial difference (see Code Table 5.6) ! !-- section 6: ibmap=255 ! Bit-map indicator (Table 6.0) (0:A bit map applies, 255:A bit map does not apply) ! !-- test data do j=1,ny do i=1,nx fld(i+(j-1)*nx)=270.+i*0.01+(j-1)*0.1 enddo enddo print *,'fld=',maxval(fld(1:nx*ny)),minval(fld(1:nx*ny)) ! call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl,ipdstmpllen, & coordlist,numcoord,idrsnum,idrstmpl, & idrstmpllen,fld,nx*ny,ibmap,bmap,ierr) print*,'addfield status=',ierr !-- finalize GRIB message after all section !-- adds the End Section ( "7777" ) call gribend(cgrib,max_bytes,lengrib,ierr) print*,'gribend status=',ierr print*,'length of the final GRIB2 message in octets =',lengrib ! call wryte(ifilw, lengrib, cgrib) ! !-- generate wind grib message: cgrib='' ! !-- call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr) print*,'gribcreate status=',ierr call addgrid(cgrib,max_bytes,igds,igdstmpl,igdstmpllen,ideflist,idefnum,ierr) print*,'addgrid status=',ierr !u wind do j=1,ny do i=1,nx fld(i+(j-1)*nx)=i*0.01+(j-1)*0.1 enddo enddo ipdstmpl(1)=2 ! catogory (Momentum:2) ipdstmpl(2)=2 ! parameter (ugrd:2) ! idrstmpl(3)=4 ! decimal scale idrstmpl(4)=14 ! Decimal scale factor call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl,ipdstmpllen, & coordlist,numcoord,idrsnum,idrstmpl, & idrstmpllen,fld,nx*ny,ibmap,bmap,ierr) !v wind do j=1,ny do i=1,nx fld(i+(j-1)*nx)=i*0.1+(j-1)*0.01 enddo enddo ipdstmpl(1)=2 ! catogory (Momentum:2) ipdstmpl(2)=3 ! parameter (vgrd:3) ! idrstmpl(3)=4 ! decimal scale idrstmpl(4)=14 ! Decimal scale factor call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl,ipdstmpllen, & coordlist,numcoord,idrsnum,idrstmpl, & idrstmpllen,fld,nx*ny,ibmap,bmap,ierr) call gribend(cgrib,max_bytes,lengrib,ierr) ! call wryte(ifilw, lengrib, cgrib) print*,'after wrt cgrib2, lengrib=',lengrib ! return end