#
# lcase
# '#' @
# QSUB -eo
# QSUB -s /bin/csh
# QSUB -r printc -q reg -lT 80 -lM 10Mw

ja
cd $TMPDIR

set sub1 = 'anyGRIB1.f'

cat << EOF > tsrc.f
c
c To use this program, you need to make $sub1 (as defined above) available.
c The new_anyGRIB1.f is copied over to CRAY in this job.
c
c To read one GRIB1 record from COS blocked file using "buffer in"
c CRAY only
c
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
      program crayGRIB
c
c maximum byte array expected: 60000*8
c
      parameter (nbit=64)              ! number of bits in a word,
      parameter (nword=nbit/8)         ! number of bytes in a word,
      parameter (mxcray=480000)        ! 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=mxcray/nword)
      dimension ingrib(mxlen),ione(mxone),itwo(mxone)
      integer*4 lengrib,gribedi,idum,kkin,iskp, ibloc(4)
      character one*480000,two*480000
      equivalence (one,ione),(itwo,two)
c
      call assign("assign -a grib.dat -b 4 fort.10")
c     write(*,*) 'How many GRIB records to skip? <21>'
c     read(*,*) mgrib
      mgrib=0
c
c set iprint 1 to print out information in each section
      iprint=0
      ngrib=0
   11 continue
      buffer in (10,0) (ione(1),ione(mxone))
      klen=length(10)
      if(klen.eq.0) go to 999
      lenrd=klen*8
      do i=1,lenrd
        if(one(i:i+3).eq.'GRIB') then
c
c check to see if this is another new GRIB record
c
          call gbytes(one,ibloc,(i+3)*8,8,0,4)
          write(*,*) ' ibloc: ',(ibloc(kk),kk=1,4)
          if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then
c           write(*,*) ' GRIB Edition 1 '
            gribedi=1
            go to 25
          elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and.
     &       ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then
c           write(*,*) ' ECMWF GRIB Edition 0 '
            gribedi=0
            go to 25
          elseif(ibloc(1).eq.98) then
c           write(*,*) ' ECMWF GRIB Edition X '
            gribedi=-1
            go to 25
          else
          endif
        endif
      enddo
      go to 11 
   25 continue
      ngrib=ngrib+1
      if(ngrib.le.mgrib) go to 11
      ibeg=i
      do j=ibeg,lenrd
        ij=j-ibeg+1
        two(ij:ij)=one(j:j)
      enddo
      lcount=(lenrd-ibeg+1)/nword
      if((lcount*nword).ne.(lenrd-ibeg+1)) lcount=lcount+1
      kkin=0
      do j=1,lcount
        ingrib(j)=itwo(j)
      enddo
      kkin=lcount
      ioff=4
      iskp=ioff*8
      call gbyte(ingrib,lengrib,iskp,24)
c     call gbyte(ingrib,gribedi,iskp+24,8)
c     call gbyte(ingrib,idum,iskp+32,8)
c     write(*,*) ' klen = ',klen,ibeg,lengrib,gribedi,idum
      if(gribedi.eq.0) then
        write(*,*) '....ECMWF edition 0'
        lenpds=lengrib
        call gbyte(ingrib,kgds,iskp+56,1)
        call gbyte(ingrib,kbms,iskp+57,1)
        write(*,*) ' kgds,kbms = ',kgds,kbms
        if(kgds.eq.0) then
          lengds=0
          go to 40
        endif
        iskp=(lenpds+ioff)*8
        call gbyte(ingrib,lengds,iskp,24)
c       write(*,*) ' GDS length: ',lengds
   40   continue
        if(kbms.eq.0) then
          lenbms=0
          go to 50
        endif
        iskp=(lengds+lenpds+ioff)*8
        call gbyte(ingrib,lenbms,iskp,24)
c       write(*,*) ' BMS length: ',lenbms
   50   continue
        iskp=(lenbms+lengds+lenpds+ioff)*8
        call gbyte(ingrib,lenbds,iskp,24)
c       write(*,*) '     bdslen: ',lenbds
        lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
c       write(*,*) 'lenpds,gds,bms,bds: ',
c    &       ioff,lenpds,lengds,lenbms,lenbds
c
      elseif(gribedi.eq.-1) then
        write(*,*) '....ECMWF edition X,  no PDS length'
        lenpds=20
        call gbyte(ingrib,kgds,iskp+24+7,1)
        call gbyte(ingrib,kbms,iskp+24+6,1)
        write(*,*) ' kgds,kbms = ',kgds,kbms
        if(kgds.eq.0) then
          lengds=0
          go to 43
        endif
        iskp=(lenpds+ioff)*8
        call gbyte(ingrib,lengds,iskp,24)
        write(*,*) ' GDS length: ',lengds
   43   continue
        if(kbms.eq.0) then
          lenbms=0
          go to 53
        endif
        iskp=(lengds+lenpds+ioff)*8
        call gbyte(ingrib,lenbms,iskp,24)
        write(*,*) ' BMS length: ',lenbms
   53   continue
        iskp=(lenbms+lengds+lenpds+ioff)*8
        call gbyte(ingrib,lenbds,iskp,24)
        write(*,*) '     bdslen: ',lenbds
        lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
        write(*,*) 'lenpds,gds,bms,bds: ',
     &       ioff,lenpds,lengds,lenbms,lenbds
      elseif(gribedi.eq.1) then
        write(*,*) '....Normal GRIB 1 format '
      else
        write(*,*) '....UNKNOWN GRIB format, exit'
        call exit(0)
      endif
      leftrd=lengrib-(lenrd-ibeg+1)
      mloop=leftrd/(klen*8)
      write(*,*) ' leftrd...: ',leftrd,lengrib,lenrd,klen,kkin,mloop
      ii=0
      if(mloop.gt.0) then
        do i=1,mloop
          ii=kkin+(i-1)*klen
          buffer in (10,0) (ingrib(ii+1),ingrib(ii+klen))
        enddo
      endif
      ii=ii+klen
      if(mod(leftrd,(klen*8)).gt.0) then
        ilast=(leftrd-klen*8*mloop)/8
        if((ilast*8).ne.(leftrd-klen*8*mloop)) ilast=ilast+1
        buffer in (10,0) (ingrib(ii+1),ingrib(ii+ilast))
      endif
c
c double check for short GRIB
c
      write(*,*) 'one GRIB in ingrib, length=',lengrib,ii+ilast
      lenword=ii+ilast
      call genGRIB1_decode(mxlen,nbit,nword,lenword,ingrib)
c     write(*,*) 'enter 1 for more? '
c     read(*,*) imore
c     if(imore.eq.1) go to 11
      go to 11
  999 stop
      end
c
c
      subroutine genGRIB1_decode(mxlen,nbit,nword,lnin,ingrib)
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
      dimension ingrib(mxlen)
      integer*4  ione
      integer    itwo
      character one*4, two*8
      equivalence (ione,one), (itwo,two)
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
      data ifirst/0/
c
c kgrib1=0, NMC; =1,ECMWF
c
      save
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
        if(nbit.eq.64) then
          itwo=ingrib(i)
          if(two(1:4).eq.'GRIB') go to 15
          if(two(5:8).eq.'GRIB') go to 15
        else
          ione=ingrib(i)
          if(one.eq.'GRIB') go to 15
        endif
      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 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
        go to 50
      endif
      iskp=(lengds+pds(1)+ioff)*8
      call gbyte(ingrib,lenbms,iskp,24)
c     write(*,*) ' BMS length: ',lenbms
      call bms_GRIB1(iprnt,lnin,nmax,iskp,
     &  kbms,nx,ny,la1,lo1,ierr,numgd,mapuse,pds,
     &  ingrib,mapb)

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
      write(*,*) 'pds(1),gds,bms,bds: ',
     &     ioff,pds(1),lengds,lenbms,lenbds
      write(*,*) 'lnin = ',lnin
      if((lnin*nword).lt.lenall) then
        write(*,*) ' *** SHORT GRIB record ',lnin*nword,lenall
        write(*,*) ' *** no bds_GRIB1, return '
        return
      endif
c     write(*,*) 'enter 1 to continue ...'
c     read(*,*) imore
c     if(imore.ne.1) go to 9998
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(rmax.lt.1) then
        write(*,81) (datan(k),k=1,300)
   81   format(15f8.3)
      else
        write(*,83) (datan(k),k=1,300)
   83   format(15f8.0)
      endif
 9998 continue
      write(*,*) ' finish reading input grib record....'
      return
c
 9999 write(*,*) '**** ERROR during read, stop'
      call exit(0)
      end
c
c
EOF

rcp huron.scd.ucar.edu:/huron/ftp/libraries/grib/$sub1 $sub1

cat $sub1 >> tsrc.f

#cft77 -e h tsrc.f
f90 -c tsrc.f

segldr tsrc.o -L /usr/local/lib -l ncarm,ncaro,mss

if ( -x a.out ) then
  rcp huron.scd:/huron/u2/tmp/chifan/GRIB/sample.5cos grib.dat
# assign -a grib.dat -b 4 fort.10
  a.out
#
  rcp huron.scd:/huron/u2/tmp/chifan/GRIB/y04746.pcos grib.dat
  a.out
endif
ja -s
exit