program pcGRIB
c
c PC MS_FORTRAN version 
c
c The maximum number of grids in the BDS section is limited by nnx,nny set
c in the subroutine genGRIB1_decode.  If nnx*nny exceeds the system CPM memory
c limit, this program will hang the system.  It is recommend to use the
c current nnx,nny settings (e.g. 144,73) for debugging purpose.  Although not
c all grids are decoded when number of grids exceeds (nnx*nny), the program
c will run through various grib records.
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
      parameter (mxrd=1024)                     ! number of bytes in each strmread
      parameter (nbit=32)                       ! number of bits in a word,
      parameter (nword=nbit/8)                  ! number of bytes in a word,
      parameter (mxcray=28000)                  ! max number of words for a GRIB on Cray
      parameter (mxlen=mxcray*8/nword)          ! max number of words for a GRIB on Sun
      parameter (mxone=mxrd/nword)              ! max length of ione
      dimension ingrib(mxlen), ione(mxone)
      character*80 infl
      integer*4 ii,iskp,ibloc(4),lengrib,lenpds,gribedi
      character one*1024,ach*4,bch*4,two*1024
      character asum*132
      equivalence (one,ione)
c
c if (isum.eq.1) print out one line summary of GRIB records
c reset isum=0 if you want to see more details
      isum=1
c
c open the file
c
      write(*,*) ' enter input GRIB file? <sample.pbi>'
      read(*,'(a)') infl
      lunin=30
      open(unit=lunin, file=infl, form='binary')
c
c get one character at a time
c
      ngrib=0
      ii=4
      read(lunin) ach(1:4)
   10 continue
      if(ach.eq.'GRIB') go to 30
      read(lunin,end=9998,err=9998) bch(4:4)
      bch(1:3)=ach(2:4)
      ach=bch
      ii=ii+1
      go to 10
   30 continue
      one=' '
c     write(*,*) ' found GRIB at ',ii-1
      one(1:4)=ach
c
c read 4 more bytes
c
      read(lunin) ach(1:4)
      one(5:8)=ach
      call swap4(ach,ach,4)
      call gbytes(ach,ibloc,0,8,0,4)
c     write(*,*) 'ibloc: ',(ibloc(k),k=1,4)
      if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then
        if(isum.eq.0) write(*,*) ' GRIB Edition 1 '
        call gbyte(ach,lengrib,0,24)
        call gbyte(ach,gribedi,24,8)
      elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and.
     &       ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then
        if(isum.eq.0) write(*,*) ' ECMWF GRIB Edition 0 '
        call gbyte(ach,lengrib,0,24)
        call gbyte(ach,gribedi,24,8)
      elseif(ibloc(1).eq.98) then
        if(isum.eq.0) write(*,*) ' ECMWF GRIB Edition X '
        lengrib=20
        gribedi=-1
      else
        write(*,*) ' UNKNOWN GRIB format, exit '
        go to 10
      endif
      if(isum.eq.0) write(*,*) 'lengrib,edition: ',lengrib,gribedi
c
c decode for lengths of gds,bms,bds
c
      if(gribedi.eq.0) then
        ioff=4
        lenpds=lengrib
        ibeg=4
        read(lunin) one(ibeg+5:ibeg+lenpds)
c       write(*,*) '....ECMWF edition 0'
        iskp=ioff*8
        call swap4(one,two,lenpds+ibeg)
        call gbyte(two,kgds,iskp+56,1)
        call gbyte(two,kbms,iskp+57,1)
c       write(*,*) ' kgds,kbms = ',kgds,kbms
        if(kgds.eq.0) then
          lengds=0
          go to 40
        endif
        read(lunin) ach
        one(lenpds+5:lenpds+8)=ach
        call swap4(ach,ach,4)
        call gbyte(ach,lengds,0,24)
c       write(*,*) ' GDS length: ',lengds
        ibeg=4+lenpds
        read(lunin) one(ibeg+5:ibeg+lengds)
   40   continue
        if(kbms.eq.0) then
          lenbms=0
          go to 50
        endif
        read(lunin) ach
        ibeg=4+lenpds+lengds
        one(ibeg+1:ibeg+4)=ach
        call swap4(ach,ach,4)
        call gbyte(ach,lenbms,0,24)
c       write(*,*) ' BMS length: ',lenbms
        read(lunin) one(ibeg+5:ibeg+lenbms)
   50   continue
        read(lunin) ach
        ibeg=4+lenpds+lengds+lenbms
        one(ibeg+1:ibeg+4)=ach
        call swap4(ach,ach,4)
        call gbyte(ach,lenbds,0,24)
c       write(*,*) '     bdslen: ',lenbds
        lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
        if(isum.eq.0) write(*,*) 'lenpds,gds,bms,bds: ',
     &       ioff,lenpds,lengds,lenbms,lenbds
c
      elseif(gribedi.eq.-1) then
c       write(*,*) '....ECMWF edition X,  no PDS length'
        ioff=4
        lenpds=20
        iskp=ioff*8
        read(lunin) one(9:lenpds+4)
        call swap4(one,two,lenpds+4)
        call gbyte(two,kgds,iskp+24+7,1)
        call gbyte(two,kbms,iskp+24+6,1)
c       write(*,*) ' kgds,kbms = ',kgds,kbms
        if(kgds.eq.0) then
          lengds=0
          go to 43
        endif
        read(lunin) ach
        ibeg=4+lenpds
        one(ibeg+1:ibeg+4)=ach
        call swap4(ach,ach,4)
        call gbyte(ach,lengds,0,24)
c       write(*,*) ' GDS length: ',lengds
        read(lunin) one(ibeg+5:ibeg+lengds)
   43   continue
        if(kbms.eq.0) then
          lenbms=0
          go to 53
        endif
        read(lunin) ach
        ibeg=4+lenpds+lengds
        one(ibeg+1:ibeg+4)=ach
        call swap4(ach,ach,4)
        call gbyte(ach,lenbms,0,24)
c       write(*,*) ' BMS length: ',lenbms
        read(lunin) one(ibeg+5:ibeg+lenbms)
   53   continue
        read(lunin) ach
        ibeg=4+lenpds+lengds+lenbms
        one(ibeg+1:ibeg+4)=ach
        call swap4(ach,ach,4)
        call gbyte(ach,lenbds,0,24)
c       write(*,*) '     bdslen: ',lenbds
        lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
        if(isum.eq.0) write(*,*) 'lenpds,gds,bms,bds: ',
     &       ioff,lenpds,lengds,lenbms,lenbds
c
      else
        ibeg=4
      endif
      if(isum.eq.0) then
        write(*,*) 'hit <CR> to continue....'
        read(*,*)
      endif
      kread=lengrib-(ibeg+4)
c     write(*,*) 'lengrib,kread: ',lengrib,kread
c finish reading one complete GRIB
c   length left to read is kread
c   starting from (ibeg+5)
c
c fill ingrib first
c
      ilast=ibeg+4
      ileft=mod((ibeg+4),nword)
      if(ileft.ne.0) then
        read(lunin) one(ibeg+5:ibeg+4+nword-ileft)
        ilast=ibeg+4+nword-ileft
      endif
      mloop=ilast/nword
      if(mloop.gt.mxone) then
        write(*,*) '****ERROR in fill ingrib first'
        go to 9998
      endif
      indxt=0
      if(mloop.ge.1) then
        do i=1,mloop
          indxt=indxt+1
          ingrib(indxt)=ione(i)
        enddo
      endif
c
c now read in the rest
c
      mloop=kread/mxrd
      if(mloop.ge.1) then
        do j=1,mloop
          one=' '
          read(lunin) one(1:mxrd)
          do i=1,mxone
            indxt=indxt+1
            ingrib(indxt)=ione(i)
          enddo
        enddo
      endif
      ileft=mod(kread,mxrd)
      if(ileft .gt.0) then
        one=' '
        read(lunin) one(1:ileft)
        mloop=ileft/nword
        if(mod(ileft,nword).ne.0) mloop=mloop+1
        do i=1,mloop
          indxt=indxt+1
          ingrib(indxt)=ione(i)
        enddo
      endif
      if(isum.eq.0) 
     &  write(*,*) 'TOTAL length of ingrib is ',indxt,indxt*nword
      ngrib=ngrib+1
c
c swap4 of the ingrib and decode it
c
      do i=2,indxt
        ione(1)=ingrib(i)
        call swap4(ione(1),ione(1),4)
        ingrib(i)=ione(1)
      enddo
c     call swap4(ingrib,ingrib,lengrib)
      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 9998
      else
        write(*,'(a)') asum(1:125)
      endif
      go to 10
c
   90 continue
      close(lunin)
 9998 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     parameter (nnx=720,nny=361,nmax=nnx*nny)  ! global, 0.5 degree
      parameter (nnx=144,nny=73,nmax=nnx*nny)
      dimension ingrib(mxlen)
      character*4 one
      character asum*132, short*5
      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
c     save
c
c to see more details, reset iprnt=1
      iprnt=0
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 in this grib record.',nbit
      return
c
   15 continue
      ioff=i-1
c ioff is the number of bytes before 'GRIB'
c     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)
c     write(*,*) ' ibloc: ',(ibloc(ka),ka=1,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)
        if(asum.eq.' ') 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)
        if(asum.eq.' ') write(*,*) '.. ECMWF GRIB Edition 0 '
        lenpds=lengrib
        ioff=ioff+4
        iskp=ioff*8
      elseif(ibloc(1).eq.98) then
        if(asum.eq.' ') write(*,*) '.. ECMWF GRIB Edition X '
        gribedi=-1
        lenpds=20
        ioff=4
        iskp=ioff*8
      else
        write(*,*) '..unknown GRIB record, exit ',lengrib,idum
        go to 9998  
      endif
c     write(*,*) ' GRIB, length,edi= ',lengrib,gribedi
      idum=lenpds+ioff
c
c at this time, ioff is the number of bytes before PDS
c     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
      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,numgd,mapuse
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      ! 'GRIB+...+7777'
      else
        lenall=8+pds(1)+lengds+lenbms+lenbds+4
      endif
      if(asum.eq.' ') write(*,*) 'pds(1),gds,bms,bds: ',
     &     ioff,pds(1),lengds,lenbms,lenbds
c     write(*,*) 'lnin = ',lnin
      short=' '
      if((lnin*nword).lt.lenall) then
        write(*,*) ' *** SHORT GRIB record ',lnin*nword,lenall
        short='SHORT'
        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