c To compile and link this program, use command line c f77 ibmGRIB.f anyGRIB1.f f77.f c c The data array is in (datan(i),i=1,npoints) in subtoutine genGRIB1_decode. c Note: c This program is designed to scan through all GRIB records after skipping c the first (mgrib) records. If you want to select certain variables/levels, c you need to know both Table 2 and Table 3 in GRIB documentation to modify c this program. The selection is done in subroutine gebGRIB1_decode between c '------ select GRIB' lines. c c The decoded grid values are in array datan. The bitmap, if used, c is returned in array mapb. Both datan and mapb are returned to c subroutine genGRIB1_decode only. In genGRIB1_decode, c npoints = number of valid grid values in datan, (datan(i),i=1,npoints) c mapuse = number of valid bit map values, (mapb(j),j=1,mapuse) c Of course, if mapuse = 0, there is no bitmap used. c If npoints is less than the number of grid values expected, you need c to check the bitmap. c c If a bitmap is used, the relationship between bitmap and grid values c is shown below: c I, a counter 1 2 3 4 5 6 7 8 9 10 11 12 .... mapuse c mapb(I) value 1 0 0 1 1 1 0 0 1 1 0 1 c datan index 1 2 3 4 5 6 7 .... npoints c c---------------------------------------------------------------------------- c c c To use this ibmGRIB.f, you also need to ftp the following c files from our anonymous ftp on ncardata.ucar.edu c c libraries/grib/anyGRIB1.f c libraries/gbytes/f77.f c program ibmGRIB c c This program works on SUN as well as IBM/RS6000. c Due to the slow speed of using 'direct acess' with 'recl=1' on SUN, c it's better to use sunGRIB.f with stream I/O on SUN. c parameter (mxrd=8) parameter (nbit=32) ! number of bits in a word, parameter (nword=nbit/8) ! number of bytes in a word, parameter (mxcray=60000) ! max number of words for a GRIB on Cray parameter (mxlen=mxcray*8/nword) ! max number of words for a GRIB on Sun character*72 inf, ach*4 character one*8 character*132 asum integer*4 ia,ibloc(4),ingrib(mxlen) integer*4 gribedi,lengrib logical endfile,newgrib equivalence (ach,ia) c isum=0 ! if 0, no summary or set it to non-zero for printing out asum write(*,*) 'enter input file? ' read(*,'(a)') inf open(unit=1,file=inf,access='direct', & recl=1,form='unformatted') ioff=0 endfile=.false. infunit=1 1 continue one=' ' ioff=ioff+1 read(1,rec=ioff,err=10) one(1:1) read(1,rec=1+ioff,err=10) one(2:2) read(1,rec=2+ioff,err=10) one(3:3) read(1,rec=3+ioff,err=10) one(4:4) if(one(1:4).ne.'GRIB') go to 1 read(1,rec=4+ioff,err=10) one(5:5) read(1,rec=5+ioff,err=10) one(6:6) read(1,rec=6+ioff,err=10) one(7:7) read(1,rec=7+ioff,err=10) one(8:8) cxxxxxxxxxxxxxxxxx ach(1:4)=one(5:8) write(*,*) 'found one GRIB at ',ioff call gbytes(ia,ibloc,0,8,0,4) write(*,*) ' ibloc 4: ',(ibloc(jk),jk=1,4) if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then write(*,*) ' GRIB Edition 1 ' newgrib=.true. indxt=0 gribedi=1 call gbyte(ia,lengrib,0,24) write(*,*) ' lengrib = ',lengrib elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and. & ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then write(*,*) ' ECMWF GRIB Edition 0 ' newgrib=.true. indxt=0 gribedi=0 call gbyte(ia,lenpds,0,24) write(*,*) ' lenpds = ',lenpds elseif(ibloc(1).eq.98) then write(*,*) ' ECMWF GRIB Edition X ' newgrib=.true. indxt=0 gribedi=-1 lenpds=20 write(*,*) ' lenpds = ',lenpds else newgrib=.false. write(*,*) ' FOUND a GRIB at wrong place' go to 999 endif c ingrib(2)=ia ach(1:4)=one(1:4) ingrib(1)=ia indxt=2 if(gribedi.eq.1) then ! read all GRIB message lenrd=4*(lengrib/4) c do i=ioff+8,ioff+lengrib-1,4 do i=ioff+8,ioff+lenrd-1,4 read(1,rec=i,err=10) ach(1:1) read(1,rec=1+i,err=10) ach(2:2) read(1,rec=2+i,err=10) ach(3:3) read(1,rec=3+i,err=10) ach(4:4) indxt=indxt+1 ingrib(indxt)=ia enddo lenleft=mod(lengrib,4) if(lenleft.ne.0) then ach=' ' do i=1,lenleft read(1,rec=(i+ioff+lenrd-1),err=10) ach(i:i) enddo indxt=indxt+1 ingrib(indxt)=ia endif ioff=ioff+lengrib-1 else do i=ioff+8,ioff+lenpds+3,4 ! read PDS section first read(1,rec=i,err=10) ach(1:1) read(1,rec=1+i,err=10) ach(2:2) read(1,rec=2+i,err=10) ach(3:3) read(1,rec=3+i,err=10) ach(4:4) indxt=indxt+1 ingrib(indxt)=ia enddo ioff=ioff+lenpds+3 if(gribedi.eq.0) then call gbyte(ingrib,kgds,32+56,1) call gbyte(ingrib,kbms,32+57,1) else call gbyte(ingrib,kgds,32+31,1) call gbyte(ingrib,kbms,32+30,1) endif write(*,*) ' kgds,kbms = ',kgds,kbms if(kgds.eq.0) then lengds=0 go to 40 endif c read(1,rec=1+ioff,err=10) ach(1:1) read(1,rec=2+ioff,err=10) ach(2:2) read(1,rec=3+ioff,err=10) ach(3:3) read(1,rec=4+ioff,err=10) ach(4:4) call gbyte(ach,lengds,0,24) write(*,*) ' GDS length: ',lengds,indxt do i=ioff+1,ioff+lengds,4 read(1,rec=i) ach(1:1) read(1,rec=1+i,err=10) ach(2:2) read(1,rec=2+i,err=10) ach(3:3) read(1,rec=3+i,err=10) ach(4:4) indxt=indxt+1 ingrib(indxt)=ia enddo ioff=ioff+lengds 40 continue if(kbms.eq.0) then lenbms=0 go to 50 endif read(1,rec=1+ioff,err=10) ach(1:1) read(1,rec=2+ioff,err=10) ach(2:2) read(1,rec=3+ioff,err=10) ach(3:3) read(1,rec=4+ioff,err=10) ach(4:4) call gbyte(ach,lenbms,0,24) write(*,*) ' BMS length: ',lenbms,indxt do i=ioff+1,ioff+lenbms,4 read(1,rec=i) ach(1:1) read(1,rec=1+i,err=10) ach(2:2) read(1,rec=2+i,err=10) ach(3:3) read(1,rec=3+i,err=10) ach(4:4) indxt=indxt+1 ingrib(indxt)=ia enddo ioff=ioff+lenbms 50 continue read(1,rec=1+ioff,err=10) ach(1:1) read(1,rec=2+ioff,err=10) ach(2:2) read(1,rec=3+ioff,err=10) ach(3:3) read(1,rec=4+ioff,err=10) ach(4:4) call gbyte(ach,lenbds,0,24) write(*,*) ' BDS length: ',lenbds,indxt lenrd=4*(lenbds/4) if(mod(lenbds,4).ne.0) lenrd=4*((lenrd/4)+1) c do i=ioff+1,ioff+lenbds+1,4 do i=ioff+1,ioff+lenrd+1,4 read(1,rec=i) ach(1:1) read(1,rec=1+i,err=10) ach(2:2) read(1,rec=2+i,err=10) ach(3:3) read(1,rec=3+i,err=10) ach(4:4) indxt=indxt+1 ingrib(indxt)=ia if(mod(i,2000).eq.0) write(*,*) ' read bds ',i enddo lengrib=4+lenpds+lengds+lenbms+lenbds+4 ! 'GRIB+...+7777' write(*,*) 'lenpds,gds,bms,bds: ', & ioff,lengrib,lenpds,lengds,lenbms,lenbds ioff=ioff+lenbds endif write(*,*) ' total words: indxt = ',indxt asum=' ' if(isum.ne.0) write(asum(1:6),'(i6)') ngrib call genGRIB1_decode(mxlen,nbit,nword,indxt,asum,ingrib) if(isum.eq.0) then write(*,*) 'enter 1 to continue ...' read(*,*) imore if(imore.ne.1) go to 999 else c write(asum(1:6),'(i6)') ngrib write(*,'(a)') asum(1:125) endif if(.not. endfile) go to 1 10 continue write(*,*) 'total length: ',ioff,i 999 close(infunit) stop end c c subroutine genGRIB1_decode(mxlen,nbit,nword,lnin,asum,ingrib) c c asum: one line summary of current GRIB record c if(asum.eq.' ') then print out certain information c c maximum byte array expected: 60000*8 c c basic logic: c (1) read one complete GRIB record (machine dependent) c (2) decode PDS c (3) decode GDS c (4) decode BMS if there c (5) decode BDS c if more processing, repeat steps (1) to (5) c parameter (nnx=720,nny=361,nmax=nnx*nny) ! global, 0.5 degree c parameter (nbit=32) ! number of bits in a word, c parameter (nword=nbit/8) ! number of bytes in a word, c parameter (mxcray=60000) ! max number of words for a GRIB on Cray c parameter (mxlen=mxcray*8/nword)! max number of words for a GRIB on Sun dimension ingrib(mxlen) integer ione character*4 one, short*5 character*132 asum equivalence (ione,one) c integer*4 idum,ioff,iskp,ibloc(4) c c PDS c integer*4 lengrib,gribedi,pds(22),kgds,kbms,lvl2 c c GDS used c integer*4 lengds c c BMS used c integer*4 lenbms, mapuse logical*1 mapb(nmax) c c BDS used c integer*4 lenbds,npoints,ix(nmax) dimension datan(nmax) c save c c reset iprnt=1 if you want details printed out iprnt=1 c do i=1,nmax ix(i)=0 datan(i)=0. enddo c c read in one GRIB record until the next 'GRIB' c c To handle blank records of "NCAR NMC 2.5 degree global analysis grids", c use direct access and check for 'GRIB'. If no 'GRIB', ierr=-9. c ioff=0 do i=1,lnin ione=ingrib(i) if(one.eq.'GRIB') go to 15 enddo write(*,*) '****WARNING, no GRIB found in this grib record.' return c 15 continue ioff=i-1 c ioff is the number of bytes before 'GRIB' write(*,*) ' One GRIB starts at: ',ioff c c check to see if NMC or ECMWF c iskp=(ioff+4)*8 call gbytes(ingrib,ibloc,iskp,8,0,4) if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then call gbyte(ingrib,lengrib,iskp,24) call gbyte(ingrib,gribedi,iskp+24,8) write(*,*) '.. GRIB Edition 1 ' ioff=ioff+8 iskp=ioff*8 call gbyte(ingrib,lenpds,iskp,24) elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and. & ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then call gbyte(ingrib,lengrib,iskp,24) call gbyte(ingrib,gribedi,iskp+24,8) write(*,*) '.. ECMWF GRIB Edition 0 ' lenpds=lengrib ioff=ioff+4 iskp=ioff*8 elseif(ibloc(1).eq.98) then write(*,*) '.. ECMWF GRIB Edition X ' gribedi=-1 lenpds=20 ioff=4 iskp=ioff*8 else write(*,*) '..unknown GRIB record, exit ',lengrib,idum call exit(0) endif write(*,*) ' GRIB, length,edi= ',lengrib,gribedi idum=lenpds+ioff c c at this time, ioff is the number of bytes before PDS write(*,*) 'lengrib,lenpds,ioff: ',lengrib,lenpds,ioff c c begin to process GRIB c to check BLOK, loop through the whole ingrib array c c c ioff is the number of bytes before PDS c iskp is the number of bits before PDS c c process PDS c iprnt=1 call pds_GRIB1(iprnt,gribedi,iskp,idum, & kgds,kbms,dscalf,ierr,lvl2,pds,ingrib) c------------------------------------------------ select GRIB c if((pds(7).ne.33) .and. (pds(7).ne.34) .and. !U,V c & (pds(7).ne.1) .and. (pds(7).ne.7) .and. !Pres,Z c & (pds(7).ne.11) .and. (pds(7).ne.52)) return !T, RH c c if((pds(8).ne.100) .and. (pds(7).ne.1)) return !isobaric level or sfc Pres c if(pds(8).ne.100 .or. pds(9).ne.850) return ! isobaric level, 850mb c if(pds(9).ne.1000 .or. lvl2.ne.500) return ! layer between 1000,500 c if(pds(10).ne.92) return ! 1992, year c if(pds(11).ne.11) return ! 11, month c if(pds(12).ne.1) return ! 1, day c if(pds(13).ne.0) return ! 0, UTC hour c------------------------------------------------ select GRIB c c c To handle old GRIB Edition from ECMWF c if(gribedi.eq.-1) then c call gbyte(ingrib,kgds,7*8+7,1) c call gbyte(ingrib,kbms,7*8+6,1) c write(*,*) '....... revised kgds,kbms .....' c write(*,*) '.... kgds = ',kgds c write(*,*) '.... kbms = ',kbms c call gbytes(ingrib,ibloc,10*8,8,0,2) c if(pds(8).eq.100 .or. pds(8).eq.103 .or. c & pds(8).eq.105 .or. pds(8).eq.107 .or. c & pds(8).eq.109 .or. pds(8).eq.111 .or. c & pds(8).eq.113 .or. pds(8).eq.125 .or. c & pds(8).eq.160 .or. pds(8).eq.200 .or. c & pds(8).eq.201 ) then c pds(9)=ibloc(1)*32+ibloc(2) c write(*,*) '....... revised level.....' c write(*,*) '....level = ',pds(9) c endif c endif c c c process GDS c if(kgds.eq.0) then lengds=0 go to 40 endif iskp=(pds(1)+ioff)*8 call gbyte(ingrib,lengds,iskp,24) c write(*,*) ' GDS length: ',lengds call gds_GRIB1(iprnt,gribedi,lnin,iskp, & kgds,ierr,nx,ny,la1,lo1,pds,ingrib) c c BMS section c 40 continue mapuse=0 if(kbms.eq.0) then lenbms=0 c go to 50 else iskp=(lengds+pds(1)+ioff)*8 call gbyte(ingrib,lenbms,iskp,24) endif call bms_GRIB1(iprnt,lnin,nmax,iskp,kbms,nx, & ny,la1,lo1,ierr,numgd,mapuse,pds,ingrib,mapb) c write(*,*) ' BMS length: ',lenbms c c BDS section c 50 continue iskp=(lenbms+lengds+pds(1)+ioff)*8 call gbyte(ingrib,lenbds,iskp,24) c write(*,*) ' bdslen: ',lenbds c c check the total length with GRIB header c if(pds(3).eq.98 .and. gribedi.ne.1) then ! ECMWF edition 0 or X ' lenall=(4+lenpds+lengds+lenbms+lenbds+4)/1 ! 'GRIB+...+7777' else lenall=(8+pds(1)+lengds+lenbms+lenbds+4)/1 endif if(iprnt.ne.0) write(*,*) 'lenall,pds(1),gds,bms,bds: ', & ioff,pds(1),lengds,lenbms,lenbds,lenall c write(*,*) 'lnin = ',lnin short=' ' if((lnin*nword).lt.lenall) then write(*,*) ' *** SHORT GRIB record ',lnin*nword,lenall if((lenall-(lnin*nword)).gt.4) then write(*,*) ' *** no bds_GRIB1, return ' return endif endif if(asum.eq.' ') then write(*,*) 'enter 1 to continue ...' read(*,*) imore if(imore.ne.1) go to 9998 endif c write(*,*) '..reading BDS, iskp = ',iskp c call bds_GRIB1(iprnt,gribedi,lnin,iskp,dscalf,ierr,nmax, & npoints,rc00,rmax,rmin,numgd,mapuse,mapb,pds,ingrib, & ix,datan) if(ierr.eq.-3) then write(*,*) '****Too many grids than specified..' ierr=0 return endif if(ierr.eq.-8) then write(*,*) '****Constant field, all grids = ',datan(1) write(*,*) ' No number of grids can be computed from bds' write(*,*) ' Use number of grids from gds section' ierr=0 return endif if(ierr.eq.-11) then write(*,*) '****UNKNOWN predefined grid type, ' write(*,*) ' No grids decoded' ierr=0 return endif if(ierr.eq.-12) then write(*,*) '****COMPLEX PACKING of spherical harmonic ' write(*,*) ' coefficients, to be developed!' write(*,*) ' No grids decoded' ierr=0 return endif if(asum.eq.' ') write(*,83) (datan(k),k=1,15) 83 format(4(E15.8,1X),E15.8) 9998 continue if(asum.eq.' ') & write(*,*) ' finish reading input grib record....' c c make one line summary c if(asum.ne.' ') & write(asum(7:132),3213) (pds(i),i=1,5),(pds(j),j=7,9), & lvl2,(pds(i),i=10,13),(pds(k),k=15,18),pds(21), & pds(22),lenbds,npoints,rmax,rmin,short 3213 format(i3,i4,2i3,3i4,2i5,i3,3i2.2,2i4,i5,i4,2i5,i8,i6, & 2f9.1,1x,a5,18x) c3213 format(6x,i3,i4,2i3,3i4,2i5,i3,3i2.2,4i4,2i5,i8,i6,2f9.1,1x,a5) c nn 32 2 7 80 2 7 100 1000 0 91070600 1 0 0 10 20 0 13152 10512+999999.9+999999.9 short c123456iiijjjj1231231234iiiijjjj1234512345kkkiijjii12341234mmmmnnnnkkkkjjjj12345678123456fffffffffgggggggggxaaaaa c12345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890 c0 1 2 3 4 5 6 7 8 9 c return end c c