program gribread c c New version as of June 1996 to handle the extended product definition c section which ECMWF began using in March 1996. The major c change is that many arrays are longer. c c This program must be modified for your own purposes. You must c do the following things: c c 1) Set the variable maxcount appropriately -- it limits the number c of output records -- for testing purposes. You will probably want c to eliminate the maxcount test once you have determined that the c code works and gives you what you want. c 2) Change the test in the if-block for parameter testing to c meet your needs. For example, you may want to test for c parm = 139 (surface temperature) or lev1 = 850, or test for c a range of dates. Parameter codes can be found in table 2 c of the ECMWF document "FM 92 GRIB." Contact Data Support c (datahelp@ncar.ucar.edu) for a copy of this document. c 3) Change or remove plotting code (calls to gks and subroutine c gridplot). As it stands, a quick-and-dirty contour plot is c made and saved to a file called 'gmeta', if you are using c NCAR graphics. This is commented out in this version. c 4) You may have to change the PARAMETER jpword to the number of c bits/word on your system. c 5) I/O is done by direct access, which reads records without c Fortran control words. Most Sun-like systems require opening c the file with the record length in bytes, but some, SGI and c DEC Ultrix among them, require the record length in words. c If you are having problems, you may need to try specifying c the length in words. If you have asynchronous I/O available c (check the "libraries/io" directory in the ftp area on ncardata) c you can use those routines. The Fortran unit for data input in c this program is set by a parameter statement, currently iunit=9. c 6) The record length is set by a parameter statement. The c various ECMWF/TOGA datasets have differing record lengths; c make sure that you have the record length (variable NRECL) c set correctly for the dataset you are reading. c 7) The routines GBYTES and GBYTE are needed for this program. c If you accessed this program by ftp, you can find these programs c in the subdirectory /ftp/libraries/gbytes. Several versions c are available, and you will need to determine which c version is appropriate for your local computer system. c On a Sun or Sun-like system, f77.f is the Fortran version c you should use (also see c.c for a C version). c 8) There are several files containing necessary subroutines c in the same directory as this file. You will need the file c gribsubs.f and the appropriate subs_.f for your system. c Most UNIX systems can use the Sun version. c c If you can't read the data, re-check steps 1-8 above to make sure c that you have made all the necessary modifications. If you are c still having problems, phone me at (303) 497-1214 or (preferably) c e-mail me at ilana@ncar.ucar.edu with a complete description of c your problem. c c WRITTEN BY ILANA STERN, DSS/SCD/NCAR c Subroutines supplied by ECMWF c c universal read code for ds111.x data sets (x=0,1,2,3,5) c set up for a Sun computer c PARAMETER (JPWORD=32) parameter (iunit=9) c PARAMETER (JPFLEN=60000, JPGLEN=60000, NRECL=?????) c c Values for NRECL: c c ECMWF/TOGA datasets (operational analyses) c c ds111.0: 8400 bytes through May 1985 c 23280 bytes Jun 1985 - 16 Sep 1991 c 11640 bytes 17 Sep 1991 - Dec 1995 c 24152 bytes Jan 1995 - present c ds111.1: 18720 bytes through May 1985 c 25680 bytes Jun 1985 - present c ds111.2: 10560 bytes through Feb 1996 c 10680 bytes Mar 1996 - Jun 1998 c 21152 bytes Jul 1998 - present c ds111.3: 12840 bytes through Dec 1989 c 25680 bytes Jan 1990 - present c ds111.5: 10560 bytes through Feb 1996 c 10680 bytes Mar 1996 - present c c ECMWF Reanalyses (2.5 deg) c c ds115.0: 21136 bytes c ds115.1: 21136 bytes c ds115.2: 21144 bytes c c ECMWF Reanalyses (high resolution) c c ds115.3: 17968 bytes c ds115.4: 23752 bytes c ds115.5: 17984 bytes c ds115.6: 18032 bytes c c Other ECMWF c ds115.7 21240 bytes c character*8 vsns(20), path*20, fmt*2, ofmt*1, hoper*1 character*24 filename dimension ksec1(60), ksec2(384), igrib(jpglen), f(jpflen), & finv(jpflen), cf(jpflen) dimension ksec0(2), ksec3(2), ksec4(42) dimension psec2(96), psec3(2) c ksec1 = iblk1, ksec2 = iblk2, psec4 = f c but note there is an extra variable at the front of ksec1! integer vers2, center, model, griddef, iflag1, parm, levtyp, & lev1, lev2, year, mon, day, hour, min, timeunit, & time1, time2, timeind, numavg common /head1/ vers2, center, model, griddef, iflag1, parm, & levtyp, lev1, lev2, year, mon, day, hour, min, & timeunit, time1, time2, timeind, numavg equivalence (ksec1(1), vers2) c note, was equivalence (iblk1(1), center) c c Header block 2 is set up for gaussian or lat/lon grid. c integer drtype, nlatpts, nlonpts, orglat, orglon, iflagr, & extlat, extlon, inclat, inclon, iflag8 common /head2/ drtype, nlatpts, nlonpts, orglat, orglon, iflagr, & extlat, extlon, inclat, inclon, iflag8 equivalence (ksec2(1), drtype) integer hdsave(7) c data fmt /'TR'/, ofmt /'U'/ c c ** ncar graphics c call opngks c ** ncar graphics write (*, 1000) 1000 format ('Name of file to read? ',$) read (*, 1010) filename 1010 format (A) irecg = 0 irec = 1 len = nrecl / (jpword / 8) icount = 0 maxcount = 6 open (unit=iunit, form='unformatted', file=filename, err=800, & recl=nrecl, access='direct', status='old', iostat=ierr) c c read in one complete grib, which may be several blocks long c 100 continue iblock = 1 istart = 1 110 continue read (iunit, rec=irec, err=700) & (igrib(i), i=istart,istart+len-1) c write (*,*) (igrib(i), i=1, 24) irec = irec + 1 c c for the first tape block of the grib, find the grib length c iend = istart + len - 1 if (iblock .eq. 1) then hoper = 'L' call gribex (ksec0, ksec1, ksec2, psec2, ksec3, psec3, & ksec4, f, jpflen, igrib, jpglen, idum, hoper, & ierr) igwds = ksec0(1) / (jpword / 8) end if c c If the length is unrealistic or lengrib returns -1 (error) c may be a dropped record c c if (igwds .gt. 20000) go to 100 if (igwds .eq. -1) go to 100 c c read the rest of the grib if there's more c if (igwds .gt. iend) then iblock = iblock + 1 istart = istart + len go to 110 end if c c decode the grib c irecg = irecg + 1 hoper = 'D' call gribex (ksec0, ksec1, ksec2, psec2, ksec3, psec3, & ksec4, f, jpflen, igrib, jpglen, idum, hoper, ierr) ivals = ksec4(1) if (ierr. ne. 0) then write (*, 1900) ierr, irec, vsns(i) 1900 format (' Grib conversion error ', I5, ' in record ', & I5, ' of file ', A) stop end if c c return from decogb: c ksec1 = product definition block (see /head1/) c ksec2 = grid description block (see /head2/) c f(j), j=1, ivals = grid values c c This is the parameter test if-block -- you will want to change c this to test for parameter, or date, or level, or some combination. c if (1 .eq. 1) then c c write out header blocks and sample data c call GRPRS0 (KSEC0) call GRPRS1 (KSEC0, ksec1) call GRPRS2 (KSEC0, ksec2, PSEC2) write (*, 9000) (f(j), j=1, 60) 9000 format (15F8.0) c c plot quantity c c ** ncar graphics c call gridplot (f, finv, nlatpts, nlonpts, ivals) c ** ncar graphics icount = icount + 1 end if c if (icount .le. maxcount) go to 100 c 700 continue close (iunit) go to 900 800 continue write (*, 8000) ierr, filename 8000 format ('Error ', i4, ' attempting to open ', A) c 900 continue c call clsgks stop end c ** ncar graphics subroutine gridplot (f, finv, nlatpts, nlonpts, ivals) dimension f(1), finv(1) do 300 j = 1, ivals finv(j) = f(ivals-j+1) 300 continue do 320 j = 1, nlonpts nlim = nlatpts/2 do 310 k = 1, nlim ftemp = finv((j-1)*nlatpts+k) finv((j-1)*nlatpts+k) = & finv((j-1)*nlatpts+(nlatpts-k+1)) finv((j-1)*nlatpts+(nlatpts-k+1)) = ftemp 310 continue 320 continue c c call cpezct (finv, nlatpts, nlonpts) c call cpcnrc (finv, nlatpts, nlatpts, nlonpts, c & 0., 0., 0., 0, 0, -682) c call frame return end c ** ncar graphics