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

ja
cd $TMPDIR

#set sub1 = 'anyGRIB1.f'

cat << EOF > tsrc.f
c--------------------------------------------------------------------------
c This is to demonstrate how to find and read pgb.f0092110100
c Step 1. Get the corresponding YEAR/VSN list from our anonymous ftp server
c         on ncardata.ucar.edu under
c           ftp/datasets/ds090.0/inventories
c         with filenames like 'YYYY_A.list', for example, '1992_A.list'.
c
c Step 2. Locate the file_type, VSN from the YEAR/VSN list.
c         Find '9211.pgb.f00' in '1992_A.list', and the VSN is A03662.
c
c Step 3. Each VSN is a Unix tar file with a block size of 32768 bytes.
c         To extract the data files from this VSN, you need to know the
c         correct filenames.  You may find this from the VSN/tar_list in
c         YYYYREABALtar.lst, i.e. '1992REANALtar.lst'.  Or if you are
c         familiar with the filenaming convention, you can construct it
c         yourself.  For pgb.f00 data, the filenames follow 'pgb.f00YYMMDDHH'
c         convention.  In this example, we'll extract 'pgb.f0092110100' from
c         A03662 (or /DSS/A03662 on the MSS).
c
c Step 4. Modify the section near the end of this program to extract the
c         data file wanted.
c
c Step 5. If you want to select specific variables and/or levels, search
c         for block of codes bounded by 'select GRIB' string.  Please read
c         the 'Note' below as well.
c
c Note:
c  This program is designed to scan through the whole individual data
c  file.  If you want to select certain variables/levels, you need to
c  know both Table 2 and Table 3 in GRIB documentation to modify this
c  program.
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 To use this program, you need to make anyGRIB1.f (as defined above) available.
c The anyGRIB1.f is copied over to CRAY in this job.
c
c This program will read pure binary GRIB files on Cray.  It will not
c decode GRIB cos-blocked files properly.  To read cos-blocked,
c variable record length GRIB files, please use 'GBbin2cos.qsub' and
c 'crayGRIB.qsub' files.
c
c CRAY version  tested on various versions of GRIB 
c
c assign -a grib.data -s u -s unblocked fort.30
c
      program main
c
c maximum byte array expected: 60000*8
c
c     parameter (nbit=32)                       ! number of bits in a word,
      parameter (nbit=64)                       ! number of bits in a word,
      parameter (mxrd=480000)                   ! max number of bytes
      parameter (nnx=720,nny=361,nmax=nnx*nny)  ! global, 0.5 degree
      parameter (nword=nbit/8)                  ! number of bytes in a word,
      parameter (mxlen=mxrd/nword)              ! max number of words for a GRIB
      parameter (IUN=1,MSGUN=6,LMWD=nbit,IRW=1,LWK=mxlen)
      character name*72
      dimension ione(mxlen),iwk(lwk),ibloc(4)
      character one*(mxrd),two*(mxrd+3),tmpa*(mxrd)         !(nword)
c
      integer*4 ibeg,laspos,gribedi
      logical   shortr
      character asum*132    ! one line summary of a GRIB record
c
      equivalence (one,ione)
      data iuni/30/
c     data iuno/20/ ! this is moved to subroutine genGRIB1_decode
c
      save
c
c if the first time, open it and check for recl
c
c     write(*,*) 'How many GRIB records to skip? <21>'
c     read(*,*) mgrib
      mgrib=0
c
c set isum=1 to print out summary of each GRIB record
      isum=1
c
c     write(*,*)' Enter binary filename to read? '
c     read(*,'(a)') name
c
c if you want to save the GRIB records into a pure binary GRIB file,
c activate the next assignment.  The output unit is 'iuno'.
c     call assign("assign -a grib.out -s u -s unblocked fort.20")
c
      ngrib=0
      itot=0
      laspos=0
      tmpa=' '
      lentmpa=0
      ired=25
   10 continue
      buffer in (iuni,0) (ione(1),ione(ired))
      lblk=length(iuni)
      if(lblk.eq.0) then
        write(*,*) 'END_OF_FILE during read, exit'
        go to 999
      endif
      lenrd=lblk*nword             ! bytes
      itot=itot+lenrd
      if(lentmpa.gt.0) then
        one(lentmpa+1:lentmpa+lenrd)=one(1:lenrd)
        one(1:lentmpa)=tmpa(1:lentmpa)
        lenrd=lenrd+lentmpa
      endif
c
      do i=1,lenrd-3
        if(one(i:i+3).eq.'GRIB') then
          ngrib=ngrib+1
          call gbytes(one,ibloc,(i+3)*8,8,0,4)
          if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then
            call gbyte(one,lengrib,(i+3)*8,24)
            call gbyte(one,gribedi,(i+3)*8+24,8)
c           write(*,*) ' GRIB Edition 1 ',lengrib,gribedi
            ibeg=i
          elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and.
     &       ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then
            call gbyte(one,lenpds,(i+3)*8,24)
            call gbyte(one,gribedi,(i+3)*8+24,8)
c           write(*,*) ' ECMWF GRIB Edition 0 ',lenpds,gribedi
            call gbyte(one,kgds,(i+10)*8,1)
            call gbyte(one,kbms,(i+10)*8+1,1)
c           write(*,*) ' kgds,kbms = ',kgds,kbms
            if(kgds.eq.0) then
              lengds=0
              go to 40
            endif
            call gbyte(one,lengds,(i+3+lenpds)*8,24)
c           write(*,*) ' GDS length: ',lengds
   40       continue
            if(kbms.eq.0) then
              lenbms=0
              go to 50
            endif
            iskp=(lengds+lenpds+i+3)*8
            call gbyte(one,lenbms,iskp,24)
c           write(*,*) ' BMS length: ',lenbms
   50       continue
            iskp=(lenbms+lengds+lenpds+i+3)*8
            call gbyte(one,lenbds,iskp,24)
c           write(*,*) '     bdslen: ',lenbds
            lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
            ibeg=i
c           write(*,*) 'lenpds,gds,bms,bds: ',
c    &        lenpds,lengds,lenbms,lenbds
          elseif(ibloc(1).eq.98) then
c           write(*,*) '....ECMWF edition X,  no PDS length'
            gribedi=-1
            lenpds=20
            iskp=(i+3+3)*8
            call gbyte(one,kgds,iskp+7,1)
            call gbyte(one,kbms,iskp+6,1)
c           write(*,*) ' kgds,kbms = ',kgds,kbms
            if(kgds.eq.0) then
              lengds=0
              go to 43
            endif
            iskp=(lenpds+i+3)*8
            call gbyte(one,lengds,iskp,24)
c           write(*,*) ' GDS length: ',lengds
   43       continue
            if(kbms.eq.0) then
              lenbms=0
              go to 53
            endif
            iskp=(lengds+lenpds+i+3)*8
            call gbyte(one,lenbms,iskp,24)
c           write(*,*) ' BMS length: ',lenbms
   53       continue
            iskp=(lenbms+lengds+lenpds+i+3)*8
            call gbyte(one,lenbds,iskp,24)
c           write(*,*) '     bdslen: ',lenbds
            lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
            ibeg=i
            ired=mxlen
c           write(*,*) 'lenpds,gds,bms,bds: ',
c    &        lenpds,lengds,lenbms,lenbds
          else
            gribedi=-9
          endif
c
c fill one GRIB record in two
c
          if(gribedi.eq.0 .or. gribedi.eq.1 .or.
     &       gribedi.eq.-1) then
            leftrd=lengrib-(lenrd-ibeg+1)
c           write(*,*) ' one GRIB, ',ngrib,lengrib,leftrd,ibeg,lenrd
            laspos=lenrd-ibeg+1
            two(1:laspos)=one(ibeg:lenrd)
c
c short record
c
            shortr=.false.
            if(leftrd.lt.0) then
c             if(ired.ne.mxlen) then
c               write(*,*) '****SHORT record in one read'
c               call exit(0)
c               shortr=.true.
c             else
                two(1:lengrib)=one(ibeg:(ibeg+lengrib-1))
                lentmpa=lenrd-ibeg-lengrib+1
                tmpa(1:lentmpa)=one((ibeg+lengrib):lenrd)
c             endif
            else
              ndum=leftrd/nword
              if(mod(leftrd,nword).ne.0) ndum=ndum+1
c             write(*,*) ' read 2a: ',ndum,lenrd,leftrd,nword
              buffer in (iuni,0) (ione(1),ione(ndum))
              lblk=length(iuni)
              lenrd=lblk*nword
              itot=itot+lenrd
c             write(*,*) ' read 2b: ',ngrib,ndum,lblk,lenrd,leftrd
              if(lenrd.gt.leftrd) then     ! read one extra word
                two(laspos+1:laspos+leftrd)=one(1:leftrd)
                lentmpa=lenrd-leftrd
                tmpa(1:lentmpa)=one(leftrd+1:lenrd)
              elseif(lenrd.lt.leftrd) then  ! short
                write(*,*) '****SHORT record, exit ',lenrd,leftrd
c               call exit(0)
                buffer in (iuni,0) (ione(lblk+1),ione(ndum))
                lblk2=length(iuni)
                lenrd2=lblk2*nword
                itot=itot+lenrd2
c               write(*,*) ' read 3: ',ngrib,ndum-lblk,lblk2,lenrd2
                if(lblk2.eq.0) then
                  write(*,*) '****UNEXPECTED EOF ',lblk2
                  call exit(0)
                endif
              else
                two(laspos+1:laspos+leftrd)=one(1:leftrd)
                lentmpa=0
              endif
            endif
c           write(*,*) '.. a complete GRIB in two ',lengrib,laspos+leftrd
            one(1:lengrib)=two(1:lengrib)
            mdum=lengrib/nword
            if(mod(lengrib,nword).ne.0) then
              mdum=mdum+1
              one(lengrib+1:mdum*nword)=' '
            endif
c
            asum=' '
            if(isum.ne.0) write(asum(1:6),'(i6)') ngrib
            call genGRIB1_decode(mxlen,nbit,nword,mdum,asum,ione)
c           write(*,*) '.... save a complete GRIB ',ngrib,mdum,lengrib
            if(isum.ne.0 .and. asum(7:125).ne.' ')            ! 960429
     &         write(*,'(a)') asum(1:125)
            if(shortr) then
              write(*,*) '.... shortr ',ngrib,mdum,lengrib,
     &                   leftrd,mxlen,ired
              one(1:(laspos-lengrib))=two(lengrib+1:laspos)
              lenrd=laspos-lengrib
            endif
          endif
          go to 10
        endif
      enddo
      if(lenrd.lt.7) then
        write(*,*) '****This should not happen, but it did! ',ngrib,lenrd
        call exit(0)
      endif
      lentmpa=8
      tmpa(1:lentmpa)=one((lenrd-7):lenrd)
c
c     write(*,*) 'enter 1 for more? '
c     read(*,*) imore
c     if(imore.eq.1) go to 10
c     if(ngrib.ge.50) go to 999
      go to 10
c
c if all done, close the input file and exit
c
  999 continue
      write(*,*) '....total GRIB recs ',ngrib,itot
      stop
      end
c
c
      subroutine genGRIB1_decode(mxlen,nbit,nword,lnin,asum,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, asum*132, short*5
      equivalence (ione,one), (itwo,two)
      data iuno/20/
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 kgrib1=0, NMC; =1,ECMWF
c
      save
c
c reset iprnt=1 if you want details printed out
      iprnt=0
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                             !970712
          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'
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)
      if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then
        call gbyte(ingrib,lengrib,iskp,24)
        call gbyte(ingrib,gribedi,iskp+24,8)
c       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)
c       write(*,*) '.. ECMWF GRIB Edition 0 '
        lenpds=lengrib
        ioff=ioff+4
        iskp=ioff*8
      elseif(ibloc(1).eq.98) then
c       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
c     write(*,*) ' GRIB, length,edi= ',lengrib,gribedi
      idum=lenpds+ioff

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
      if((pds(7).ne.33) .and. (pds(7).ne.34) .and.   !U,V
     &   (pds(7).ne.1)  .and. (pds(7).ne.7)  .and.   !Pres,Z
     &   (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 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
      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
      if(iprnt.ne.0) 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..', ngrib
c
c make one line summary
c
      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,6i4,i8,i6,2f9.1,1x,a5)
c    nn 32   2  7 80  2   7 100 1000    0 91070600   1   0   0  10  20   0   1315 2 10512+999999.9+999999.9 short
c123456iiijjjj1231231234iiiijjjj1234512345kkkiijjii12341234mmmmnnnnkkkkjjjj12345678123456fffffffffgggggggggxaaaaa
c12345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
c0        1         2         3         4         5         6         7         8         9
c
c to write the GRIB record to unit iuno
c
c           buffer out (iuno,0) (ingrib(1),ingrib(lnin))
      return
      end
c
c
c
c This version of anyGRIB1.f has changed the arguments passed to/from
c the subroutine calls.  It works with the GRIB decoders in this directory.
c
c
c Revision History:
c  Jan 03, 1997 -- minor, fix P1,P2 unpacking in pds_GRIB1
c  Nov 13, 1996 -- reverse the order of lines 1349,1350 (search for 961113).
c  Feb 12, 1996 -- initialize rmax,rmin to be -999999.0,999999.0
c  Jan  5, 1996 -- fix iskp in BMS.
c  Nov 30, 1995 -- fix BMS section when using a bitmap.
c  Aug 31, 1995 -- major change, to unpack THINNED grids.
c  Jul 10, 1995 -- minor change, to fix bit map section.
c  Feb 17, 1995 -- minor change, to correct GDS LOV computation.
c
c Systems Tested:
c  SUN,CRAY,PC/MS-DOS,DEC/ALPHA,IBM/RS6000
c==================================================================
c The subroutines below are to deocode GRIB format
c
      subroutine bds_GRIB1(iprnt,gribedi,nword,inoff,dscalf,ierr,
     &  nmax,npoints,rc00,rmax,rmin,numgd,mapuse,mapb,pds,ingrib,
     &  ix,datan)
c
c iprnt:   print out flag, 0=NO printing
c gribedi: GRIB edition number, 1, 0, -1
c nword:   length of ingrib
c inoff:   offset number of bits 
c dscalf:  decimal scale factor, from PDS
c ierr: error code, 
c       0=no error, 
c      -1=number of bits per data value exceeds word length or maximum allowed
c      -2=input array holds fewer coded values than expected,
c      -3=output array too small to hold decoded data values
c      -4=currently undefined GRIB Edition 1 option selected,
c      -8=constant field, npoints cannot be computed, use number of grids from 
c         gds to fill datan, the reference value is in datan(1) 
c      -9=blank record
c     -11=unknown predefined grid
c     -12=complex packing of spher harmonic coeff, to be developed
c       1=no '7777' at the end
c nmax:    maximum dimension of datan
c npoints: actual number of array unpacked
c rc00:
c rmax:    maximum value of array datan
c rmin:    minimum value of array datan
c numgd:   number of actual grids used, counted from BMP bitmap
c mapuse:  if 0, do not use mapb; otherwise, mapb(mapuse)
c mapb:    bit map from BMS or generated
c pds:     PDS array
c ingrib:  integer array containing GRIB record
c ix:      integer array of unpacked array, of dimension (nnx*nny)
c datan:   returned array of length (nmax=nnx*nny)
c
c bdslen: length of BDS
c bdsflg: BDS flag (byte 4 bits 1-4 and byte 14)
c bdsub:  number of unused bits in BDS
c bdse:   BDS scale factor
c bdsref: BDS reference value
c bdsn:   number of bits used for packing
c
      dimension ingrib(nword),datan(nmax)
      logical*1 mapb(nmax)
      real*4    refval, bscalf
      integer*4 bdslen,bdsflg(8),bdsub,bdse,bdsn,bdsn1
      integer*4 pds(22),ix(nmax),refa,refb,gribedi
      integer*4 N,P,JL,KL,ML
      integer*4 ifov(361),kbdsn2(361)  ! 361 is the maximum nny
      character*72 bdsnam(6)
      data ifirst/0/
c
c     save
c
      if(ifirst.eq.0) then
        bdsnam(1)=' BDS length'                       ! byte 1-3
        bdsnam(2)=' BDS flag'                         ! byte 4, bits 1-4 & byte 14
        bdsnam(3)=' BDS number of unused bits'        ! byte 4, bits 5-8
        bdsnam(4)=' BDS scale factor'                 ! byte 5-6
        bdsnam(5)=' BDS reference value'              ! byte 7-10
        bdsnam(6)=' BDS number of bits for packing'   ! byte 11
        ifirst=9
      endif
      scalten = 10**dscalf
c
c unpack BDS section
c
      iskp=inoff
      call gbyte(ingrib,bdslen,iskp,24)
      if(iprnt.ne.0) write(*,*) '     bdslen: ',bdslen
      do i=1,8
        bdsflg(i)=0
      enddo
      iskp=iskp+24
      if(gribedi.eq.-1) then
        call gbyte(ingrib,bdsflg(1),iskp+3,1)
        call gbyte(ingrib,bdsflg(2),iskp+2,1)
        call gbyte(ingrib,bdsflg(3),iskp+1,1)
        call gbyte(ingrib,bdsflg(4),iskp,1)
        call gbyte(ingrib,bdsub,iskp+4,4)
        if(bdsflg(4).eq.1) then
          call gbyte(ingrib,bdsflg(5),iskp+10*8,1)       ! should be zero
          call gbyte(ingrib,bdsflg(6),iskp+10*8+1,1)
          call gbyte(ingrib,bdsflg(7),iskp+10*8+2,1)
          call gbyte(ingrib,bdsflg(8),iskp+10*8+3,1)
        endif
      else
        call gbyte(ingrib,bdsflg(1),iskp,1)
        call gbyte(ingrib,bdsflg(2),iskp+1,1)
        call gbyte(ingrib,bdsflg(3),iskp+2,1)
        call gbyte(ingrib,bdsflg(4),iskp+3,1)
        call gbyte(ingrib,bdsub,iskp+4,4)
        if(bdsflg(4).eq.1) then
          call gbyte(ingrib,bdsflg(5),iskp+10*8,1)       ! should be zero
          call gbyte(ingrib,bdsflg(6),iskp+10*8+1,1)
          call gbyte(ingrib,bdsflg(7),iskp+10*8+2,1)
          call gbyte(ingrib,bdsflg(8),iskp+10*8+3,1)
        endif
      endif
      if(iprnt.ne.0) 
     &  write(*,*) '     bdsflg: ',(bdsflg(i),i=1,8)
      if(iprnt.ne.0) write(*,*) '     bdsub:  ',bdsub
c
c scale factor      ! byte 5-6
c
      iskp=iskp+8
      call gbyte(ingrib,bdse,iskp,16)
      if(bdse.gt.32768) bdse=32768-bdse
      bscalf=float(bdse)
      scaltwo=2**(bscalf)
      if(iprnt.ne.0) write(*,*) '     bdse:   ',bdse
c
c reference value   ! byte 7-10
c
      iskp=iskp+16
      call gbyte(ingrib,isign,iskp,1)
      call gbyte(ingrib,refa,iskp+1,7)
      call gbyte(ingrib,refb,iskp+8,24)
      rdum=refb
      refval=(-1)**isign * rdum * 16**(refa-64) / 2**24
c
c if bdsflg(3) == 1, refval needs to be round into an integer
c
      if(iprnt.ne.0) 
     &  write(*,*) '     bdsref: ',isign,refa,refb,refval
c
c number of bits used for packing   ! byte 11
c if bdsn == 0, constatn field
c
      iskp=iskp+32
      call gbyte(ingrib,bdsn,iskp,8)
      if(iprnt.ne.0) then
        write(*,*) '     bdsn:   ',bdsn
        write(*,*) ' refval,bscalf,dscalf: ',refval,bscalf,dscalf
      endif
c
c compute or find out how many points used
c
      if(bdsflg(4).eq.0) then        ! no flag in byte 14
c
c if (bdslen-11) is odd, the last byte may be padded
c
        idum=(bdslen-11)*8-bdsub
        if(bdsn.eq.0) then           ! constant field when bdsn=0
c
c Ierr: -8=constant field, 
c          npoints cannot be computed, use number of grids from 
c          gds to fill datan, the reference value is in datan(1) 
c
          npoints=1
          datan(1)=refval
          rmax=refval
          rmin=refval
          ijrmax=-9999
          ijrmin=-9999
          ierr=-8
          return
        else
          if(bdsflg(7).eq.0 .and. mapuse.gt.0) then                   !951130
            npoints=numgd                                             !951130
          else                                                        !951130
            npoints=idum/bdsn                                         !951130
          endif                                                       !951130
          kdum=bdsn                                                   !970712
          if(mod(idum,kdum).ne.0) then                                !970712
c                                                                     !950710
c fix NMC GRIB encoder error for 1x1 degree (360x181=65160) grids     !950710
c                                                                     !950710
            if(abs(npoints-65160).le.2) then                          !950710
              npoints=65160                                           !950710
            else                                                      !950710
              write(*,*) 'WARNING! length of packed bits can not ',   !950710
     &         'be divided by bdsn ',idum,bdsn,npoints                !950710
            endif                                                     !950710
          endif
        endif
      else                          ! extended flags in byte 14
c
c   Extended BDS flags:
c                        NAME        BYTE LOCATION
c                        bdsn1       11
c                        n1          12-13
c                                    14
c                        n2          15-16
c                        np1         17-18
c                        np2         19-20
c                                    21
c                        bdsn2(np1)  22-(xx-1)  one byte each  xx=np1+22
c                        bmap2(np2)  xx-(n1-1)
c                        mp1(np1)    n1-(n2-1)
c                        mp2(np2)    n2-...
c
        bdsn1=bdsn
        iskp=iskp+8
        call gbyte(ingrib,n1,iskp,16)
        iskp=iskp+16+8
        call gbyte(ingrib,n2,iskp,16)
        iskp=iskp+16
        call gbyte(ingrib,np1,iskp,16)
        iskp=iskp+16
        call gbyte(ingrib,np2,iskp,16)
        iskp=iskp+16+8
        if(iprnt.ne.0) then
          write(*,*) '    bdsn1:   ',bdsn1
          write(*,*) '       n1:   ',n1
          write(*,*) '       n2:   ',n2
          write(*,*) '      np1:   ',np1
          write(*,*) '      np2:   ',np2
        endif
        if(bdsflg(7).eq.0 .and. mapuse.gt.0) then
          npoints=numgd
        else
          npoints=-9999  ! to be determined counting secondary bit map
        endif
      endif
c
c All flags decoded, unpack the data
c
c bdsflg(1) = 0, grid point data
c             1, spherical harmonic coefficients
c bdsflg(2) = 0, simple packing
c             1, second order ("complex") packing
c bdsflg(3) = 0, original data were floating point values
c             1, original data were integer values
c bdsflg(4) = 0, no additional flags at oct 14
c             1, octet 14 contains flag bits 5-12
c bdsflg(5) = 0, reserved, should be 0
c bdsflg(6) = 0, single datum at each grid point
c             1, matrix of values at each grid point
c bdsflg(7) = 0, no secondary bit maps
c             1, secondary bit maps present
c bdsflg(8) = 0, second order values have constant width
c             1, second order values have different width
c
      if(bdsflg(2).eq.0) then        ! simple packing
        if(iprnt.ne.0) 
     &    write(*,*) ' total number of grids: ',npoints, nmax
        if(npoints.gt.nmax) then
          ierr=-3
          write(*,*) 'ERROR in bds_GRIB1, ierr = ',ierr
          return
        endif
        if(bdsflg(1).eq.0) then      ! simple packing, grid point data
          rc00=0.
          iskp=11*8+inoff
          call gbytes(ingrib,ix,iskp,bdsn,0,npoints)
          do i=1,npoints
            datan(i)=(refval+ix(i)*scaltwo)/scalten
          enddo
        elseif(bdsflg(1).eq.1) then  ! simple packing, spher harmonic coeff
          iskp=11*8+inoff
          call gbyte(ingrib,isign,iskp,1)
          call gbyte(ingrib,refa,iskp+1,7)
          call gbyte(ingrib,refb,iskp+8,24)
          rdum=refb
          rc00=(-1)**isign * rdum * 16**(refa-64) / 2**24
          if(iprnt.ne.0)
     &      write(*,*) ' real part of coeff (0,0): ',rc00
          iskp=iskp+32
          call gbytes(ingrib,ix,iskp,bdsn,0,npoints-1)
          datan(1)=rc00
          do i=2,npoints
            datan(i)=(refval+ix(i-1)*scaltwo)/scalten
          enddo
        else
        endif
      elseif(bdsflg(2).eq.1) then    ! complex packing
        if(bdsflg(1).eq.0) then      ! complex packing, grid point data
c
c octet
c 22 -- (xx-1): widths in bits of the 2nd order packed values
c xx -- (n1-1): secondary bit maps, padded to an even number octet
c n1 -- (n2-1): np1 local minima, each is bdsn1 wide, padded
c n2 -- bdslen: np2 second order packed values, padded at the end
c
c if bdsflg(7)==0, [22--(n1-1)]: widths in bits of 2nd order packed values
c                  this is "row by row (or column by column) packing"
c
c if bdsflg(8)==0, [22--22]: width of 2nd order packed values
c if bdsflg(8)==1, but width==0, this sub-section is constant
c
          if(bdsflg(7).ne.0) then
            write(*,*) '** Secondary Bit Map in use, no unpacking yet!'
            npoints=-9999
            return
          endif
          k1off=inoff+(n1-1)*8
          do k1=1,np1
            call gbyte(ingrib,ifov(k1),k1off,bdsn1)
c           if(k1.ge.(np1-10))write(*,*) ' ifov: ',k1,k1off,ifov(k1)
            k1off=k1off+bdsn1
          enddo
c
c unpack second order value
c
          k2off=inoff+21*8         ! width starts from byte 22
          do k2=1,np1
            call gbyte(ingrib,kbdsn2(k2),k2off,8)
            k2off=k2off+8
c           if(k2.ge.(np1-10))
c    &        write(*,*) ' kbdsn2: ',k2,k2off,kbdsn2(k2)
          enddo
c
c         write(*,*) ' unpacking 2nd order value: ',numgd,mapuse
          k3off=inoff+(n2-1)*8
          nt=0
          kone=mapuse/np1
          do k2=1,mapuse
            kk=((k2-1)/kone)+1
            if(mapb(k2)) then
              if(kbdsn2(kk).eq.0) then
                isov=0
              else
                call gbyte(ingrib,isov,k3off,kbdsn2(kk))
                k3off=k3off+kbdsn2(kk)
              endif
              nt=nt+1
              datan(nt)=(refval+float(ifov(kk)+isov)*scaltwo)/scalten
c             if(nt.ge.130 .and. nt.le.140) 
c    &          write(*,*) ' isov: ',k2,nt,datan(nt),kk,ifov(kk),isov
            endif
          enddo
c         write(*,*) ' unpacking 2nd done ',numgd,nt,k2off,k3off
          if(numgd.ne.nt) write(*,*) '**WARNING, number of grids',
     &      'mismatch in 2nd order unpacking! ',numgd,nt
          npoints=nt
        elseif(bdsflg(1).eq.1) then  ! complex packing, spher harmonic coeff
          ierr=-12
          if((pds(3).ne.98) .or. (gribedi.eq.1)) then
            write(*,*) '**Complex packing of non-ECMWF spectral'
            write(*,*) '**or GRIB Edition 1, NO unpacking yet!'
            npoints=-9999
            return
          else
c
c                  Second Order Packing Formats:
c
c               GRIB 1              ECMWF(version 0, x)
c              NAME  LOCATION      LOCATION   NAME
c                N1   12-13         12-13     N
c                     14            14-15     P
c                N2   15-16         16        JL
c                P1   17-18         17        KL
c                P2   19-20         18        ML
c                     21            19-(N-1)  binary data
c         BIT WIDTH   22-(P1-1)     N-        packed binary data
c Secondary BIT MAP   xx-(N1-1)
c          MP1(NP1)   N1-(N2-1)
c          MP2(NP2)   N2-
c
c ECMWF spherical harmonic complex packing
c
            iskp=11*8+inoff
            call gbyte(ingrib,N,iskp,16)
            iskp=iskp+16
            call gbyte(ingrib,P,iskp,16)
            iskp=iskp+16
            call gbyte(ingrib,JL,iskp,8)
            iskp=iskp+8
            call gbyte(ingrib,KL,iskp,8)
            iskp=iskp+8
            call gbyte(ingrib,ML,iskp,8)
            iskp=iskp+8
c
c compute how many values to unpack
c
            idum=bdslen-19
            num=idum/4
            if(mod(idum,4).ne.0) then
              write(*,*) 'UNEVEN ECMWF spherical harmonic unpacking!'
              write(*,*) '    ',bdslen,idum,num
              return
            endif
            if(iprnt.ne.0) write(*,*) 'ECMWF spherical harmonic: ',num
c
c unpack every 4 bytes into real
c
            do i=1,num
              iskp=iskp+(i-1)*32
              call gbyte(ingrib,isign,iskp,1)
              call gbyte(ingrib,refa,iskp+1,7)
              call gbyte(ingrib,refb,iskp+8,24)
              rdum=refb
              datan(i)=(-1)**isign * rdum * 16**(refa-64) / 2**24
            enddo
            write(*,*) '**PARTIALLY DONE of '
            write(*,*) '**ECMWF spherical harmonic complex unpacking!'
            write(*,*) '**DO NOT use yet!'
            npoints=-9999
            return
          endif
        else
        endif
      else
        write(*,*) '****ERROR, unknown packing, grid/spectral!'
        stop'UNKNOWN PACKING ERROR'
      endif
c
c after unpacking, locate the maximum and minimum among datan
c
      rmax=-999999.
      rmin=999999.
      do i=1,npoints
        if(datan(i).gt.rmax) then
          ijrmax=i
          rmax=datan(i)
        endif
        if(datan(i).lt.rmin) then
          rmin=datan(i)
          ijrmin=i
        endif
      enddo
      if(iprnt.ne.0) then
        write(*,*) ' rmax, at ij = ',rmax,ijrmax
        write(*,*) ' rmin, at ij = ',rmin,ijrmin
      endif
      ierr=0
 9998 return
      end
c
c
      subroutine bms_GRIB1(iprnt,nword,nmax,inoff,
     &  kbms,nx,ny,la1,lo1,ierr,numb,mapuse,pds,
     &  ingrib,mapb)
c
c This is a generic GRIB Edition 1 BMS decoder
c
c iprnt:  print out flag, 0=NO printing
c lenpds: input, length of PDS in bytes
c ipds:   input, integer array containing GRIB PDS
c nword:  input, number of words used in ipds
c kbms:   input, GDS exists if kbms=1
c numb:   actual number of grids counted from bitmap
c mapuse: number of grids expected from grid specification
c ierr:   returned error code, 
c         0=no error, -11=unknown predefined grid
c
      dimension ingrib(nword),igrd37(73)
      logical*1 mapb(nmax)
      character*72 bmsnam(3)
      integer*4 kbms,iskp,idum,pds(22)
      integer*4 bmslen,bmsub,bmstbl
      integer*4 nx,ny,la1,lo1,mapuse
c
      data ifirst/0/
c
c igrd37 is good for grid_id 37-44
      data igrd37/73, 73, 73, 
     &    73, 73, 73, 73, 73, 72, 72, 72, 71, 71,
     &    71, 70, 70, 69, 69, 68, 67, 67, 66, 65, 
     &    65, 64, 63, 62, 61, 60, 60, 59, 58, 57,
     &    56, 55, 54, 52, 51, 50, 49, 48, 47, 45, 
     &    44, 43, 42, 40, 39, 38, 36, 35, 33, 32, 
     &    30, 29, 28, 26, 25, 23, 22, 20, 19, 17,
     &    16, 14, 12, 11,  9,  8,  6,  5,  3,  2/                 
c
      ierr=0
      if(ifirst.eq.0) then
        bmsnam(1)=' 1 = BMS length'                               !byte 1-3
        bmsnam(2)=' 2 = number of unused bits in BMS'             !byte 4
        bmsnam(3)=' 3 = BMS Table Reference'                      !byte 5-6
        ifirst=9
      endif
c
c process BMS
c
      mapuse=0
      numb=0
      if(kbms.ne.0) then
        iskp=inoff
        call gbyte(ingrib,bmslen,iskp,24)
        iskp=iskp+24
        call gbyte(ingrib,bmsub,iskp,8)
        iskp=iskp+8
        call gbyte(ingrib,bmstbl,iskp,16)                         !950710
        iskp=iskp+16                                           !960105
        if(iprnt.ne.0) then
          write(*,*) ' BMS length: ',bmslen
          write(*,*) ' BMS number of unused bits: ',bmsub
          write(*,*) ' BMS Table Reference: ',bmstbl
          write(*,*) ' BMS LA1,LO1: ',la1,lo1
        endif
c                                                                 !950710
c if bmstbl =0, a bit map follows, otherwise see Center Specs     !950710
c                                                                 !950710
        if(bmstbl.eq.0) then                                      !950710
          mapuse=bmslen*8-6*8-bmsub                               !950710
c         call gbytes(ingrib,mapb,iskp,1,0,numb)                  !950710
          do i=1,mapuse                                           !951130
            mapb(i)=.false.
            call gbyte(ingrib,idum,iskp,1)
            if(idum.eq.1) then                                    !951130
              mapb(i)=.true.                                      !951130
              numb=numb+1                                         !951130
            endif                                                 !951130
            iskp=iskp+1                                           !951130
          enddo
c         write(*,*) ' BMS Numb bitmap:     ',numb,mapuse
          return
        endif
      endif
c
c pds(3) == 07, US NMC
c           08, US NWS Telecomms Gateway
c           09, US NWS Field Stations
c           34, Japan Meteorological Agency
c           52, US NHC, Miami
c           54, Canadian
c           57, US Air Forece, GWC
c           58, US FNOC
c           59, US NOAA FSL
c           74, UKMET
c           85, French Weather Service - Toulouse
c           97, European Space Agency (ESA)
c           98, ECMWF
c           99, DeBilt, Netherlands
c
c     if(pds(3).ne.7) return
c
c NMC's predefined grids: make an rectangular mapb
c
      if(pds(5).eq.21 .or. pds(5).eq.22) then
        mapuse=37*37
        numb=1333
        m=0
        k=0
        do j=1,ny+1
          do i=1,nx
            m=m+1
            mapb(m)=.true.
            if(j.eq.(ny+1) .and. i.ge.2) mapb(m)=.false.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-1 ',k,numb
        if(m.ne.mapuse) write(*,*) '**WARNING! in bms_GRIB1-2 ',m,mapuse
      elseif(pds(5).eq.23 .or. pds(5).eq.24) then
        mapuse=37*37
        numb=1333
        m=0
        k=0
        do j=1,ny+1
          do i=1,nx
            m=m+1
            mapb(m)=.true.
            if(j.eq.1 .and. i.ge.2) mapb(m)=.false.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-3 ',k,numb
        if(m.ne.mapuse) write(*,*) '**WARNING! in bms_GRIB1-4 ',m,mapuse
      elseif(pds(5).eq.25) then
        numb=1297
        mapuse=72*19
        m=0
        k=0
        do j=1,ny+1
          do i=1,nx
            m=m+1
            mapb(m)=.true.
            if(j.eq.(ny+1) .and. i.ge.2) mapb(m)=.false.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '****WARNING! in bms_GRIB1-5 ',k,numb
        if(m.ne.mapuse) write(*,*) '**WARNING! in bms_GRIB1-6 ',m,mapuse
      elseif(pds(5).eq.26) then
        numb=1297
        m=0
        k=0
        do j=1,ny+1
          do i=1,nx
            m=m+1
            mapb(m)=.true.
            if(j.eq.1 .and. i.ge.2) mapb(m)=.false.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-7 ',k,numb
        if(m.ne.mapuse) write(*,*) '**WARNING! in bms_GRIB1-8 ',m,mapuse
      elseif(pds(5).ge.37 .and. pds(5).le.44) then
        numb=3447
        mapuse=73*73
        m=0
        k=0
        do j=1,73
          jj=j
          if(la1.eq.-90000) jj=74-j
          do i=1,73
            m=m+1
            mapb(m)=.false.
            if(i.le.igrd37(jj)) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-37 ',k,numb
        if(m.ne.mapuse) 
     &    write(*,*) '**WARNING! in bms_GRIB1-38 ',m,mapuse
      elseif(pds(5).eq.50) then
        numb=964
        mapuse=36*33
        m=0
        k=0
        do j=1,4
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.8 .and. i.le.29) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=5,8
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.7 .and. i.le.30) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=9,12
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.6 .and. i.le.31) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=13,16
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.5 .and. i.le.32) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=17,20
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.4 .and. i.le.33) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=21,24
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.3 .and. i.le.34) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=25,28
          do i=1,36
            m=m+1
            mapb(m)=.false.
            if(i.ge.2 .and. i.le.35) mapb(m)=.true.
            if(mapb(m)) k=k+1
          enddo
        enddo
        do j=29,33
          do i=1,36
            m=m+1
            mapb(m)=.true.
            k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-50 ',k,numb
        if(m.ne.mapuse) write(*,*)'**WARNING! in bms_GRIB1-51 ',m,mapuse
      elseif(pds(5).eq.61 .or. pds(5).eq.62) then
        numb=4096
        mapuse=91*46
        m=0
        k=0
        do j=1,ny+1
          do i=1,nx
            m=m+1
            mapb(m)=.true.
            if(j.eq.(ny+1) .and. i.ge.2) mapb(m)=.false.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-9 ',k,numb
        if(m.ne.mapuse) write(*,*)'**WARNING! in bms_GRIB1-10 ',m,mapuse
      elseif(pds(5).eq.63 .or. pds(5).eq.64) then
        numb=4096
        mapuse=91*46
        m=0
        k=0
        do j=1,ny+1
          do i=1,nx
            m=m+1
            mapb(m)=.true.
            if(j.eq.1 .and. i.ge.2) mapb(m)=.false.
            if(mapb(m)) k=k+1
          enddo
        enddo
        if(k.ne.numb) write(*,*) '**WARNING! in bms_GRIB1-11 ',k,numb
        if(m.ne.mapuse) write(*,*)'**WARNING! in bms_GRIB1-12 ',m,mapuse
      else
        numb=0
        ierr=-11
      endif
      return
      end
c
c
      subroutine gds_GRIB1(iprnt,gribedi,nword,inoff,
     &  kgds,ierr,nx,ny,la1,lo1,pds,igds)
c
c This is a generic GRIB Edition 1 GDS decoder
c
c iprnt:  print out flag, 0=NO printing
c lenpds: input, length of PDS in bytes
c igds:   input, integer array containing GRIB GDS section
c nword:  input, number of words used in igds
c kgds:   input, GDS exists if kgds=1
c ierr:   returned error code, 0=no error,
c
      dimension igds(nword)
      character*72 gdsnam(4), gds6name
      integer*4 kgds,iskp,pds(22),gribedi
      integer*4 gdslen,gdsnv,gdspv,gds6,shcrp,csmode
      integer*4 nx,ny,la1,lo1,gds17(8)
      integer*4 ksign,refa,refb
c
      integer*4 la2,lo2,di,dj,ngauss,latin,xssp,yssp,orient,nr
      integer*4 gds27,gds28(3),lov,dx,dy,latin2,latin1,latsp,lonsp
      integer*4 latpl,lonpl
      data ifirst/0/
c
      if(ifirst.eq.0) then
        gdsnam(1)=' 1 = length'                                       !byte 1-3
        gdsnam(2)=' 2 = NV, number of vertical coordinate parameter'  !byte 4
        gdsnam(3)=' 3 = PV, or PL or 255'                             !byte 5
        gdsnam(4)=' 4 = data representation type (table 6)'           !byte 6
        ifirst=9
      endif
c
c process GDS
c
      if(kgds.eq.0) then
        gdslen=0
        return
      endif
      iskp=inoff
      call gbyte(igds,gdslen,iskp,24)
      iskp=iskp+24
      call gbyte(igds,gdsnv,iskp,8)
      iskp=iskp+8
      call gbyte(igds,gdspv,iskp,8)
      iskp=iskp+8
      call gbyte(igds,gds6,iskp,8)
      iskp=iskp+8
      call GDS_TAB6(pds(3),gribedi,gds6,gds6name)
      if(iprnt.ne.0) then
        write(*,*) ' GDS length: ',gdslen
        write(*,*) ' GDS NV: ',gdsnv
        write(*,*) ' GDS PV: ',gdspv
        write(*,*) ' GDS Data Repres Type:(table 6) ',gds6
        write(*,*) '   ',gds6name
      endif
c
c bytes 7-17: either "spherical harmonic" or "the rest"
c                      gds6 50-80
c
      if(gds6.ge.50 .and. gds6.le.80) then
        call gbyte(igds,nj,iskp,16)
        iskp=iskp+16
        call gbyte(igds,nk,iskp,16)
        iskp=iskp+16
        call gbyte(igds,nm,iskp,16)
        iskp=iskp+16
        call gbyte(igds,shcrp,iskp,8)
        iskp=iskp+8
        call gbyte(igds,csmode,iskp,8)
        iskp=iskp+8
        if(iprnt.ne.0) then
          write(*,*) ' GDS NJ: ',nj
          write(*,*) ' GDS NK: ',nk
          write(*,*) ' GDS NM: ',nm
          write(*,*) ' GDS shcrp: (table 9)',shcrp
          write(*,*) ' GDS Coeff Storage Mode: (table 10) ',csmode
        endif
      else
        call gbyte(igds,nx,iskp,16)
        if(nx.ge.65535) nx=-99999
        iskp=iskp+16
        call gbyte(igds,ny,iskp,16)
        iskp=iskp+16
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,la1,iskp+1,23)
        if(ksign.eq.1) la1=-la1
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,lo1,iskp+1,23)
        if(ksign.eq.1) lo1=-lo1
        iskp=iskp+24
        call gbytes(igds,gds17,iskp,1,0,8)
        if(iprnt.ne.0) then
          write(*,*) ' GDS Ni: ',nx
          write(*,*) ' GDS Nj: ',ny
          write(*,*) ' GDS La1: ',la1
          write(*,*) ' GDS Lo1: ',lo1
          write(*,*) ' GDS 17: (table 7) ',(gds17(i),i=1,8)
        endif
      endif
c
c pds(3)  gds6    bytes 18-  : 
c         0,4    (1) lat/long grids (including Gaussian lat/long)          7-32
c        201     (2) Arakawa semi-staggered grid on rotated lat/long grids 7-32
c        202     (3) Arakawa filled E-grid on rotated lat/long grids       7-32
c          5     (4) Polar Stereographic grids                             7-32
c ????     2     (5) Gnomonic Projection or Stereographic grids            7-32
c          3     (6) Lambert Conformal secant or tangent cone grids        7-42
c          1     (7) Mercator grids                                        7-42
c         90     (8) Space View Perspective or Orthographic grids          7-42
c         50     (9) Spherical Harmonic Coefficients                       7-32
c  98     10    (10) ECMWF Rotated lat/long grids                          7-42
c  98     14    (11) ECMWF Rotated Gaussian grids                          7-42
c  98     60    (12) ECMWF Rotated Spherical Harmonics                     7-42
c  98     20    (13) ECMWF Stretched lat/long grids                        7-42
c  98     24    (14) ECMWF Stretched Gaussian grids                        7-42
c  98     70    (15) ECMWF Stretched Spherical Harmonics                   7-42
c  98     30    (16) ECMWF Stretched and Rotated lat/long grids            7-52
c  98     34    (17) ECMWF Stretched and Rotated Gaussian grids            7-52
c  98     80    (18) ECMWF Stretched and Rotated Spherical Harmonics       7-52
c
c bytes 18-
c
      iskp=iskp+8
      if(gds6.eq.0 .or. gds6.eq.4 .or.
     &   gds6.eq.10 .or. gds6.eq.14) then  !lat/long grids (including Gaussian)
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,la2,iskp+1,23)
        if(ksign.eq.1) la2=-la2
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,lo2,iskp+1,23)
        if(ksign.eq.1) lo2=-lo2
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,di,iskp+1,15)
        if(di.eq.32767) then
          di=-99999
        else
          if(ksign.eq.1) di=-di
        endif
        if(iprnt.ne.0) then
          write(*,*) ' GDS La2: ',la2
          write(*,*) ' GDS Lo2: ',lo2
          write(*,*) ' GDS Di: ',di
        endif
        iskp=iskp+16
        if(gds6.eq.4 .or. gds6.eq.14) then
          call gbyte(igds,ngauss,iskp,16)
          if(iprnt.ne.0) write(*,*) ' GDS N gaussian: ',ngauss
        else
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,dj,iskp+1,15)
          if(dj.eq.32767) then
            dj=-99999
          else
            if(ksign.eq.1) dj=-dj
          endif
          if(iprnt.ne.0) write(*,*) ' GDS Dj: ',dj
        endif
        iskp=iskp+16
        call gbytes(igds,gds28,iskp,1,0,3)
        if(iprnt.ne.0) 
     &    write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
      elseif(gds6.eq.201) then    !Arakawa semi-staggered rotated lat/long
        call gbyte(igds,la2,iskp,24)
        iskp=iskp+24
        call gbyte(igds,lo2,iskp,24)
        if(iprnt.ne.0) then
          write(*,*) ' GDS La2: ',la2
          write(*,*) ' GDS Lo2: ',lo2
        endif
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,di,iskp+1,15)
        if(di.eq.32767) then
          di=-99999
        else
          if(ksign.eq.1) di=-di
        endif
        if(iprnt.ne.0) write(*,*) ' GDS Di: ',di
        iskp=iskp+16
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,dj,iskp,15)
        if(dj.eq.32767) then
          dj=-99999
        else
          if(ksign.eq.1) dj=-dj
        endif
        if(iprnt.ne.0) write(*,*) ' GDS Dj: ',dj
        iskp=iskp+16
        call gbytes(igds,gds28,iskp,1,0,3)
        if(iprnt.ne.0) 
     &    write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
      elseif(gds6.eq.202) then    !Arakawa filled E-grid on rotated lat/long
        call gbyte(igds,la2,iskp,24)
        iskp=iskp+24
        call gbyte(igds,lo2,iskp,24)
        if(iprnt.ne.0) then
          write(*,*) ' GDS La2: ',la2
          write(*,*) ' GDS Lo2: ',lo2
        endif
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,di,iskp+1,15)
        if(di.eq.32767) then
          di=-99999
        else
          if(ksign.eq.1) di=-di
        endif
        if(iprnt.ne.0) write(*,*) ' GDS Di: ',di
        iskp=iskp+16
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,dj,iskp+1,15)
        if(dj.eq.32767) then
          dj=-99999
        else
          if(ksign.eq.1) dj=-dj
        endif
        if(iprnt.ne.0) write(*,*) ' GDS Dj: ',dj
        iskp=iskp+16
        call gbytes(igds,gds28,iskp,1,0,3)
        if(iprnt.ne.0) 
     &    write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
      elseif(gds6.eq.5) then      !Polar Stereographic grids
c
        call gbyte(igds,ksign,iskp,1)      ! added 02/17/95
        call gbyte(igds,lov,iskp+1,23)     ! added 02/17/95
        if(ksign.eq.1) lov=-lov            ! added 02/17/95
c       call gbyte(igds,lov,iskp,24)
c
        iskp=iskp+24
        call gbyte(igds,dx,iskp,24)
        iskp=iskp+24
        call gbyte(igds,dy,iskp,24)
        iskp=iskp+24
        call gbyte(igds,gds27,iskp,1)
        iskp=iskp+8
        call gbytes(igds,gds28,iskp,1,0,3)
c       call gbyte(igds,gds28(1),iskp,1)
c       call gbyte(igds,gds28(2),iskp+1,1)
c       call gbyte(igds,gds28(3),iskp+2,1)
        if(iprnt.ne.0) then
          write(*,*) ' GDS Lov: ',lov
          write(*,*) ' GDS Dx: ',dx
          write(*,*) ' GDS Dy: ',dy
          write(*,*) ' GDS 27: (note 5) ',gds27
          write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
        endif
      elseif(gds6.eq.2) then      !Gnomonic Projection or Stereographic 
        write(*,*) '****Gnomonic Projection,  not ready yet.'
      elseif(gds6.eq.3) then      !Lambert Conformal secant or tangent cone
c
        call gbyte(igds,ksign,iskp,1)      ! added 02/17/95
        call gbyte(igds,lov,iskp+1,23)     ! added 02/17/95
        if(ksign.eq.1) lov=-lov            ! added 02/17/95
c       call gbyte(igds,lov,iskp,24)
c
        iskp=iskp+24
        call gbyte(igds,dx,iskp,24)
        iskp=iskp+24
        call gbyte(igds,dy,iskp,24)
        iskp=iskp+24
        call gbyte(igds,gds27,iskp,1)
        iskp=iskp+8
        call gbytes(igds,gds28,iskp,1,0,3)
        iskp=iskp+8
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,latin1,iskp+1,23)
        if(ksign.eq.1) latin1=-latin1
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,latin2,iskp+1,23)
        if(ksign.eq.1) latin2=-latin2
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,latsp,iskp+1,23)
        if(ksign.eq.1) latsp=-latsp
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,lonsp,iskp+1,23)
        if(ksign.eq.1) lonsp=-lonsp
        iskp=iskp+24
        if(iprnt.ne.0) then
          write(*,*) ' GDS Lov: ',lov
          write(*,*) ' GDS Dx: ',dx
          write(*,*) ' GDS Dy: ',dy
          write(*,*) ' GDS 27: (note 5) ',gds27
          write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
          write(*,*) ' GDS latin 1: ',latin1
          write(*,*) ' GDS latin 2: ',latin2
          write(*,*) ' GDS lat SP:  ',latsp
          write(*,*) ' GDS lon SP:  ',lonsp
        endif
      elseif(gds6.eq.1) then      !Mercator grids
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,la2,iskp+1,23)
        if(ksign.eq.1) la2=-la2
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,lo2,iskp+1,23)
        if(ksign.eq.1) lo2=-lo2
        iskp=iskp+24
        call gbyte(igds,ksign,iskp,1)
        call gbyte(igds,latin,iskp+1,23)
        if(ksign.eq.1) latin=-latin
        iskp=iskp+24
        call gbyte(igds,gds27,iskp,1)
        iskp=iskp+8
        call gbytes(igds,gds28,iskp,1,0,3)
        iskp=iskp+8
        call gbyte(igds,di,iskp,24)
        iskp=iskp+24
        call gbyte(igds,dj,iskp,24)
        if(iprnt.ne.0) then
          write(*,*) ' GDS La2:   ',la2
          write(*,*) ' GDS Lo2:   ',lo2
          write(*,*) ' GDS Latin: ',latin
          write(*,*) ' GDS 27: (note 5) ',gds27
          write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
          write(*,*) ' GDS Di: ',di
          write(*,*) ' GDS Dj: ',dj
        endif
      elseif(gds6.eq.90) then     !Space View Perspective or Orthographic grids
        call gbyte(igds,dx,iskp,24)
        iskp=iskp+24
        call gbyte(igds,dy,iskp,24)
        iskp=iskp+24
        call gbyte(igds,xssp,iskp,24)
        iskp=iskp+24
        call gbyte(igds,yssp,iskp,24)
        iskp=iskp+24
        call gbytes(igds,gds28,iskp,1,0,3)
        iskp=iskp+8
        call gbyte(igds,orient,iskp,24)
        iskp=iskp+24
        call gbyte(igds,nr,iskp,24)
        iskp=iskp+24
        if(iprnt.ne.0) then
          write(*,*) ' GDS Dx:     ',dx
          write(*,*) ' GDS Dy:     ',dy
          write(*,*) ' GDS Xssp:   ',xssp
          write(*,*) ' GDS Yssp:   ',yssp
          write(*,*) ' GDS 28: (table 8) ',(gds28(i),i=1,3)
          write(*,*) ' GDS Orient: ',orient
          write(*,*) ' GDS Nr:     ',nr
        endif
      elseif(gds6.eq.50) then     !Spherical Harmonic Coefficients

      else
      endif
c
      if(gds6.eq.10 .or. gds6.eq.14 .or. gds6.eq.60) then   !ECMWF Rotated only
        if(pds(3).eq.98) then
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,latsp,iskp+1,23)
          if(ksign.eq.1) latsp=-latsp
          iskp=iskp+24
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,lonsp,iskp+1,23)
          if(ksign.eq.1) lonsp=-lonsp
          iskp=iskp+24
          call gbyte(igds,isign,iskp,1)
          call gbyte(igds,refa,iskp+1,7)
          call gbyte(igds,refb,iskp+8,24)
          rdum=refb
          angrot=(-1)**isign * rdum * 16**(refa-64) / 2**24
          if(iprnt.ne.0) then
            write(*,*) ' GDS Lat SP:            ',latsp
            write(*,*) ' GDS Lon SP:            ',lonsp
            write(*,*) ' GDS Angle of Rotation: ',angrot
          endif
        endif
      elseif(gds6.eq.20 .or. gds6.eq.24 .or. gds6.eq.70) then
c                                                        !ECMWF Stretched only
        if(pds(3).eq.98) then
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,latpl,iskp+1,23)
          if(ksign.eq.1) latpl=-latpl
          iskp=iskp+24
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,lonpl,iskp+1,23)
          if(ksign.eq.1) lonpl=-lonpl
          iskp=iskp+24
          call gbyte(igds,isign,iskp,1)
          call gbyte(igds,refa,iskp+1,7)
          call gbyte(igds,refb,iskp+8,24)
          rdum=refb
          strfac=(-1)**isign * rdum * 16**(refa-64) / 2**24
          if(iprnt.ne.0) then
            write(*,*) ' GDS Lat Pole:          ',latpl
            write(*,*) ' GDS Lon Pole:          ',lonpl
            write(*,*) ' GDS Stretching Factor: ',strfac
          endif
        endif
      elseif(gds6.eq.30 .or. gds6.eq.34 .or. gds6.eq.80) then
c                                         !ECMWF Stretched & Rotated
        if(pds(3).eq.98) then
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,latsp,iskp+1,23)
          if(ksign.eq.1) latsp=-latsp
          iskp=iskp+24
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,lonsp,iskp+1,23)
          if(ksign.eq.1) lonsp=-lonsp
          iskp=iskp+24
          call gbyte(igds,isign,iskp,1)
          call gbyte(igds,refa,iskp+1,7)
          call gbyte(igds,refb,iskp+8,24)
          rdum=refb
          angrot=(-1)**isign * rdum * 16**(refa-64) / 2**24
          if(iprnt.ne.0) then
            write(*,*) ' GDS Lat SP:            ',latsp
            write(*,*) ' GDS Lon SP:            ',lonsp
            write(*,*) ' GDS Angle of Rotation: ',angrot
          endif
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,latpl,iskp+1,23)
          if(ksign.eq.1) latpl=-latpl
          iskp=iskp+24
          call gbyte(igds,ksign,iskp,1)
          call gbyte(igds,lonpl,iskp+1,23)
          if(ksign.eq.1) lonpl=-lonpl
          iskp=iskp+24
          call gbyte(igds,isign,iskp,1)
          call gbyte(igds,refa,iskp+1,7)
          call gbyte(igds,refb,iskp+8,24)
          rdum=refb
          strfac=(-1)**isign * rdum * 16**(refa-64) / 2**24
          if(iprnt.ne.0) then
            write(*,*) ' GDS Lat Pole:          ',latpl
            write(*,*) ' GDS Lon Pole:          ',lonpl
            write(*,*) ' GDS Stretching Factor: ',strfac
          endif
        endif
      else
      endif
      ierr=0
      return
      end
c
c
      subroutine GDS_TAB6(idcenter,gribedi,indx,name)
      parameter    (maxx=255)
      character*72 code(maxx),name
      integer*4    idcenter,gribedi
      data ifirst/0/
c
      if(ifirst.ne.0) go to 10
      do i=1,maxx
        code(i)=' '
      enddo
c
c               123456789012345678901234567890123456789012345678901234567890
      code(1)='latitude/longitude grid'
      code(2)='Mercator projection'
      code(3)='Gnomonic Porjection grid'
      code(4)='Lambert Conformal projection'
      code(5)='Gaussian latitude/longitude grid'
      code(6)='Polar stereographic'
      code(14)='Oblique Lambert conformal projection'
      code(51)='spherical harmonic coefficients'
      code(91)='space view perspective or orthographic grid'
      code(202)='Arakawa semi-stagg E-grid on rotated lat/long grids'
      code(203)='Arakawa filled E-grid on rotated lat/long grids'
c
   10 continue
      if(idcenter.eq.98 .and. (gribedi.eq.0 .or. gribedi.eq.-1)) then
        code(3)='Stereographic projection'
        code(11)='rotated latitude/longitude grids'
        code(14)=' '
        code(15)='rotated Gaussian latitude/longitude grids'
        code(21)='stretched latitude/longitude grids'
        code(25)='stretched Gaussian latitude/longitude grids'
        code(31)='stretched and rotated latitude/longitude grids'
        code(35)='stretched and rotated Gaussian lat/long grids'
        code(61)='rotated spherical harmonic coefficients'
        code(71)='streched spherical harmonics'
        code(81)='stretched and rotated spher harmonic coeff'
        code(91)=' '
        code(202)=' '
        code(203)=' '
      endif
      idum=indx+1
      name=code(idum)
      return
      end
c
c
      subroutine pds_GRIB1(iprnt,gribedi,iskp,nword,
     &  kgds,kbms,dscalf,ierr,lvl2,pds,ipds)
c
c This is a generic GRIB Edition 1 PDS decoder 
c
c SUN FORTRAN version 
c
c lenpds: input, length of PDS in bytes
c iskp:   input, number of bits to skip from the beginning
c ipds:   input, integer array containing GRIB PDS
c nword:  input, number of words used in ipds
c pds:    returned, PDS array
c kgds:   returned, GDS exists if kgds=1
c kbms:   returned, BMS exists if kbms=1
c ierr:   returned error code, 0=no error, 
c
      dimension ipds(nword)
      character pdsnam(22)*72
      integer*4 pds(22),kgds,kbms,iskp,gribedi,ibloc(2)
      data ifirst/0/
c
      save
c
      if(ifirst.eq.0) then                                              ! pds bytes
        pdsnam(1)=' 1=length'                                           !  1-3
        pdsnam(2)=' 2=parameter table number'                           !  4
        pdsnam(3)=' 3=center ID'                                        !  5
        pdsnam(4)=' 4=process ID'                                       !  6
        pdsnam(5)=' 5=grid ID'                                          !  7
        pdsnam(6)=' 6=flag'                                             !  8
        pdsnam(7)=' 7=parameter ID (table 2)'                           !  9 
        pdsnam(8)=' 8=type of level or layer (table 3 & 3a)'            ! 10
        pdsnam(9)=' 9=pressure or height of level or layer'             ! 11-12
        pdsnam(10)='10=year of centry'                                  ! 13
        pdsnam(11)='11=month of year'                                   ! 14
        pdsnam(12)='12=day of month'                                    ! 15
        pdsnam(13)='13=hour of day'                                     ! 16
        pdsnam(14)='14=minute of hour'                                  ! 17
        pdsnam(15)='15=forecast time units (table 4)'                   ! 18
        pdsnam(16)='16=P1'                                              ! 19
        pdsnam(17)='17=P2'                                              ! 20
        pdsnam(18)='18=time range indicator (table 5)'                  ! 21
        pdsnam(19)='19=number included in averge when pds(18)<>0'       ! 22-23
        pdsnam(20)='20=number missing from averages or accumulations'   ! 24
        pdsnam(21)='21=century of initial (reference) time'             ! 25
c                                                                       ! 26
        pdsnam(22)='22=decimal scale factor D'                          ! 27-28
c                                                                       ! 29-40
c                                                                       ! 41-
        ifirst=9
      endif
c
      if(gribedi.eq.-1) then
        pds(1)=20
        pds(2)=128
      else
        call gbyte(ipds,pds(1),iskp,24)
        iskp=iskp+24
        call gbyte(ipds,pds(2),iskp,8)
        iskp=iskp+8
      endif
      call gbyte(ipds,pds(3),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(4),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(5),iskp,8)
      iskp=iskp+8
      if(gribedi.eq.-1) then
        call gbyte(ipds,kgds,7*8+7,1)                   ! ECMWF Edition X
        call gbyte(ipds,kbms,7*8+6,1)                   ! ECMWF Edition X
      else
        call gbyte(ipds,kgds,iskp,1)
        call gbyte(ipds,kbms,iskp+1,1)
        pds(6)=(kgds*10+kbms)*1000000                   !950711
      endif
      iskp=iskp+8
      call gbyte(ipds,pds(7),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(8),iskp,8)
      iskp=iskp+8
      if(pds(8).eq.100 .or. pds(8).eq.103 .or.
     &   pds(8).eq.105 .or. pds(8).eq.107 .or.
     &   pds(8).eq.109 .or. pds(8).eq.111 .or.
     &   pds(8).eq.113 .or. pds(8).eq.125 .or. 
     &   pds(8).eq.115 .or. pds(8).eq.117 .or.      !1998MAY14
     &   pds(8).eq.119 .or.                         !1998MAY14
     &   pds(8).eq.160 .or. pds(8).eq.200 .or. 
     &   pds(8).eq.201 ) then
        if(gribedi.eq.-1) then                   ! ECMWF Edition X
          call gbytes(ipds,ibloc,10*8,8,0,2)
          pds(9)=ibloc(1)*32+ibloc(2)
        else
          call gbyte(ipds,pds(9),iskp,16)
        endif
        lvl2=-9999
        iskp=iskp+16
      else
        call gbyte(ipds,pds(9),iskp,8)
        iskp=iskp+8
        call gbyte(ipds,lvl2,iskp,8)
        iskp=iskp+8
      endif
      call gbyte(ipds,pds(10),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(11),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(12),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(13),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(14),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(15),iskp,8)
      iskp=iskp+8
      idum=iskp                                 !970103
      call gbyte(ipds,pds(16),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(17),iskp,8)
      iskp=iskp+8
      call gbyte(ipds,pds(18),iskp,8)
      iskp=iskp+8
      if(pds(18).eq.10) then                    !970103
        call gbyte(ipds,pds(16),idum,16)        !970103
        pds(17)=-9999                           !970103
      endif                                     !970103
      call gbyte(ipds,pds(19),iskp,16)
      iskp=iskp+16
      call gbyte(ipds,pds(20),iskp,8)
      iskp=iskp+8
c
c fix if pdslen < 28
c
      dscalf=0.
      pds(22)=-9999
      if(pds(1).le.24) pds(21)=-9999
      if(pds(1).ge.28) then
        call gbyte(ipds,pds(21),iskp,8)
        iskp=iskp+8
c
c skip byte 26
c
        iskp=iskp+8
        call gbyte(ipds,isign,iskp,1)
        call gbyte(ipds,pds(22),iskp+1,15)
        if(isign.eq.1) pds(22)=-pds(22)           ! 961113
        dscalf=float(pds(22))                     ! 961113
      endif
c
c print out PDS
c
      if(iprnt.ne.0) then
        do i=1,22
          if(i.eq.6) then
            write(*,101) pds(i),kgds,kbms,pdsnam(i)
  101       format(i10,i2,i2,a70)
          elseif(i.eq.9) then
            write(*,115) pds(i),lvl2,pdsnam(i)
  115       format(i7,i6,1x,a70)
          else
            write(*,102) pds(i),pdsnam(i)
  102       format(i10,4x,a70)
          endif
        enddo
      endif
c
c fill ib1par
c
c     ib1par(1)=pds(3)
c     ib1par(2)=pds(4)
c     ib1par(3)=pds(5)
c     ib1par(4)=pds(6)
c     ib1par(5)=pds(7)
c     ib1par(6)=pds(8)
c     ib1par(7)=pds(9)
c     ib1par(8)=lvl2
c     ib1par(9)=pds(10)
c     ib1par(10)=pds(11)
c     ib1par(11)=pds(12)
c     ib1par(12)=pds(13)
c     ib1par(13)=pds(14)
c     ib1par(14)=pds(15)
c     ib1par(15)=pds(16)
c     ib1par(16)=pds(17)
c     ib1par(17)=pds(18)
c     ib1par(18)=pds(19)
c     ib1par(19)=pds(20)
c     ib1par(20)=pds(21)
c
      ierr=0
      return
      end
c
c
      subroutine fill_ingrib(mxone,mxlen,nword,lenrd,lentot,
     &                       ione,ingrib)
      dimension ingrib(mxlen),ione(mxone)
c
      mword=lenrd/nword
      if(mod(lenrd,nword).ne.0) mword=mword+1 
      do i=1,mword
        lentot=lentot+1
        ingrib(lentot)=ione(i)
      enddo
      return
      end
c
c
EOF

#rcp huron.scd.ucar.edu:/home/crestone/ftp/libraries/grib/$sub1 $sub1
#
#cat $sub1 >> tsrc.f

f90 -c tsrc.f

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

if ( -x a.out ) then
  set inmss = 'A03662'
  set flwant = 'pgb.f0092110100'
  msread -f BI $inmss /DSS/$inmss
#
# just one file
# tar xvf $inmss ./$flwant
#
# or, concatenate any GRIB files into one
  tar xvf $inmss ./$flwant ./pgb.f0092110106 ./pgb.f0092110112 ./pgb.f0092110118
  cat pgb.f0092110106 pgb.f0092110112 pgb.f0092110118 >> $flwant
#
  assign -a $flwant -s u -s unblocked -b 4 fort.30
  a.out
  if ( -r grib.out ) then
    rcp grib.out huron.scd.ucar.edu:/huron/u2/tmp/chifan/REANL/test.grib
  endif
endif

ls -l
ja -s
exit