#! /usr/bin/csh
#
cd /ptmp/chifan
#
cat << "EOF" > ibmspGRIB.f
c
c This GRIB decoder is for IBM AIX workstation only.
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  If you want to see how to select a region out of 2.5x2.5 grids or 
c  gaussian (192x94) grids, look for 'select a region'.
c
c----------------------------------------------------------------------------
c
      program ibmGRIB
c
c This program works on SUN as well as IBM/RS6000.  
c
      parameter (mxrd=8)
      parameter (nbit=32)               ! number of bits in a word,
      parameter (nword=nbit/8)          ! number of bytes in a word,
      parameter (mxcray=60000)          ! max number of words for a GRIB on Cray
      parameter (mxlen=mxcray*8/nword)  ! max number of words for a GRIB on Suni
      character*72 inf, ach*4
      character one*8
      character*132 asum
      integer*4 ia,ibloc(4),ingrib(mxlen)
      integer*4 gribedi,lengrib
      logical endfile,newgrib
      equivalence (ach,ia)
c
      isum=1    ! if 0, no summary or set it to non-zero for printing out asum
c
c open the file
c
      infunit=1
      write(*,*) 'enter input GRIB filename? '
      read(*,'(a)') inf
      write(*,*) inf
      open(unit=1,file=inf,access='direct',
     &     recl=1,form='unformatted')
c
c     write(*,*) 'How many GRIB records to skip? <21>'
c     read(*,*) mgrib
      mgrib=0
c
      ioff=0
      endfile=.false.
      ngrib=0
c
    1 continue
      one=' '
      ioff=ioff+1
      read(1,rec=ioff,err=10)   one(1:1)
      read(1,rec=1+ioff,err=10) one(2:2)
      read(1,rec=2+ioff,err=10) one(3:3)
      read(1,rec=3+ioff,err=10) one(4:4)
      if(one(1:4).ne.'GRIB') go to 1
      ngrib=ngrib+1
      if (ngrib.le.mgrib) then
        ioff=3+ioff
        go to 1
      endif
      read(1,rec=4+ioff,err=10) one(5:5)
      read(1,rec=5+ioff,err=10) one(6:6)
      read(1,rec=6+ioff,err=10) one(7:7)
      read(1,rec=7+ioff,err=10) one(8:8)
cxxxxxxxxxxxxxxxxx
      ach(1:4)=one(5:8)
      write(*,*) 'found one GRIB at ',ioff
      call gbytes(ia,ibloc,0,8,0,4)
      write(*,*) ' ibloc 4: ',(ibloc(jk),jk=1,4)
      if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then
        write(*,*) ' GRIB Edition 1 '
        newgrib=.true.
        indxt=0
        gribedi=1
        call gbyte(ia,lengrib,0,24)
        write(*,*) ' lengrib = ',lengrib
      elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and. 
     &  ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then
        write(*,*) ' ECMWF GRIB Edition 0 '
        newgrib=.true.
        indxt=0
        gribedi=0
        call gbyte(ia,lenpds,0,24)
        write(*,*) ' lenpds = ',lenpds
      elseif(ibloc(1).eq.98) then
        write(*,*) ' ECMWF GRIB Edition X '
        newgrib=.true.
        indxt=0
        gribedi=-1
        lenpds=20
        write(*,*) ' lenpds = ',lenpds
      else
        newgrib=.false.
        write(*,*) ' FOUND a GRIB at wrong place'
        go to 999
      endif
c
      ingrib(2)=ia
      ach(1:4)=one(1:4)
      ingrib(1)=ia
      indxt=2
      if(gribedi.eq.1) then    ! read all GRIB message
        lenrd=4*(lengrib/4)
c       do i=ioff+8,ioff+lengrib-1,4
        do i=ioff+8,ioff+lenrd-1,4
          read(1,rec=i,err=10)   ach(1:1)
          read(1,rec=1+i,err=10) ach(2:2)
          read(1,rec=2+i,err=10) ach(3:3)
          read(1,rec=3+i,err=10) ach(4:4)
          indxt=indxt+1
          ingrib(indxt)=ia
        enddo
        lenleft=mod(lengrib,4)
        if(lenleft.ne.0) then
          ach=' '
          do i=1,lenleft
            read(1,rec=(i+ioff+lenrd-1),err=10) ach(i:i)
          enddo
          indxt=indxt+1
          ingrib(indxt)=ia
        endif
        ioff=ioff+lengrib-1
      else
        do i=ioff+8,ioff+lenpds+3,4     ! read PDS section first
          read(1,rec=i,err=10)   ach(1:1)
          read(1,rec=1+i,err=10) ach(2:2)
          read(1,rec=2+i,err=10) ach(3:3)
          read(1,rec=3+i,err=10) ach(4:4)
          indxt=indxt+1
          ingrib(indxt)=ia
        enddo
        ioff=ioff+lenpds+3
        if(gribedi.eq.0) then
          call gbyte(ingrib,kgds,32+56,1)
          call gbyte(ingrib,kbms,32+57,1)
        else
          call gbyte(ingrib,kgds,32+31,1)
          call gbyte(ingrib,kbms,32+30,1)
        endif
        write(*,*) ' kgds,kbms = ',kgds,kbms
        if(kgds.eq.0) then
          lengds=0
          go to 40
        endif
c
        read(1,rec=1+ioff,err=10) ach(1:1)
        read(1,rec=2+ioff,err=10) ach(2:2)
        read(1,rec=3+ioff,err=10) ach(3:3)
        read(1,rec=4+ioff,err=10) ach(4:4)
        call gbyte(ach,lengds,0,24)
        write(*,*) ' GDS length: ',lengds,indxt
        do i=ioff+1,ioff+lengds,4
          read(1,rec=i) ach(1:1)
          read(1,rec=1+i,err=10) ach(2:2)
          read(1,rec=2+i,err=10) ach(3:3)
          read(1,rec=3+i,err=10) ach(4:4)
          indxt=indxt+1
          ingrib(indxt)=ia
        enddo
        ioff=ioff+lengds
   40   continue
        if(kbms.eq.0) then
          lenbms=0
          go to 50
        endif
        read(1,rec=1+ioff,err=10) ach(1:1)
        read(1,rec=2+ioff,err=10) ach(2:2)
        read(1,rec=3+ioff,err=10) ach(3:3)
        read(1,rec=4+ioff,err=10) ach(4:4)
        call gbyte(ach,lenbms,0,24)
        write(*,*) ' BMS length: ',lenbms,indxt
        do i=ioff+1,ioff+lenbms,4
          read(1,rec=i) ach(1:1)
          read(1,rec=1+i,err=10) ach(2:2)
          read(1,rec=2+i,err=10) ach(3:3)
          read(1,rec=3+i,err=10) ach(4:4)
          indxt=indxt+1
          ingrib(indxt)=ia
        enddo
        ioff=ioff+lenbms
   50   continue
        read(1,rec=1+ioff,err=10) ach(1:1)
        read(1,rec=2+ioff,err=10) ach(2:2)
        read(1,rec=3+ioff,err=10) ach(3:3)
        read(1,rec=4+ioff,err=10) ach(4:4)
        call gbyte(ach,lenbds,0,24)
        write(*,*) ' BDS length: ',lenbds,indxt
        lenrd=4*(lenbds/4)
        if(mod(lenbds,4).ne.0) lenrd=4*((lenrd/4)+1)
c       do i=ioff+1,ioff+lenbds+1,4
        do i=ioff+1,ioff+lenrd+1,4
          read(1,rec=i) ach(1:1)
          read(1,rec=1+i,err=10) ach(2:2)
          read(1,rec=2+i,err=10) ach(3:3)
          read(1,rec=3+i,err=10) ach(4:4)
          indxt=indxt+1
          ingrib(indxt)=ia
          if(mod(i,2000).eq.0) write(*,*) ' read bds ',i
        enddo
        lengrib=4+lenpds+lengds+lenbms+lenbds+4      ! 'GRIB+...+7777'
        write(*,*) 'lenpds,gds,bms,bds: ',
     &       ioff,lengrib,lenpds,lengds,lenbms,lenbds
        ioff=ioff+lenbds
      endif
      write(*,*) ' total words: indxt = ',indxt
      ngrib=ngrib+1
      asum=' '
      if(isum.ne.0) write(asum(1:6),'(i6)') ngrib
      call genGRIB1_decode(mxlen,nbit,nword,indxt,asum,ingrib)
      if(isum.eq.0) then
        write(*,*) 'enter 1 to continue ...'
        read(*,*) imore
        if(imore.ne.1) go to 999
      else
        write(asum(1:6),'(i6)') ngrib
        write(*,'(a)') asum(1:125)
      endif
      if(.not. endfile) go to 1
   10 continue
      write(*,*) 'total length: ',ioff,i
  999 close(infunit)
      stop
      end
c
c
      subroutine genGRIB1_decode(mxlen,nbit,nword,lnin,asum,ingrib)
c
c asum: one line summary of current GRIB record
c       if(asum.eq.' ') then print out certain information
c
c maximum byte array expected: 60000*8
c
c basic logic:
c  (1) read one complete GRIB record       (machine dependent)
c  (2) decode PDS
c  (3) decode GDS
c  (4) decode BMS if there
c  (5) decode BDS
c  if more processing, repeat steps (1) to (5)
c
      parameter (nnx=720,nny=361,nmax=nnx*nny)  ! global, 0.5 degree
c     parameter (nbit=32)                       ! number of bits in a word,
c     parameter (nword=nbit/8)                  ! number of bytes in a word,
c     parameter (mxcray=60000)        ! max number of words for a GRIB on Cray
c     parameter (mxlen=mxcray*8/nword)! max number of words for a GRIB on Sun
      dimension ingrib(mxlen)
      integer ione
      character*4 one, short*5
      character*132 asum
      equivalence (ione,one)
c
      integer*4 idum,ioff,iskp,ibloc(4)
c
c PDS 
c
      integer*4 lengrib,gribedi,pds(22),kgds,kbms,lvl2
c
c GDS used
c
      integer*4 lengds
c
c BMS used
c
      integer*4 lenbms, mapuse
      logical*1 mapb(nmax)
c
c BDS used
c
      integer*4 lenbds,npoints,ix(nmax)
      dimension datan(nmax)
c     dimension grd1(nmax), grd2(144,73)        ! for select a region use
c
      save
c
c reset iprnt=1 if you want details printed out
      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 found in this grib record.'
      return
c
   15 continue
      ioff=i-1
c ioff is the number of bytes before 'GRIB'
      write(*,*) ' One GRIB starts at: ',ioff
c
c check to see if NMC or ECMWF
c
      iskp=(ioff+4)*8
      call gbytes(ingrib,ibloc,iskp,8,0,4)  
      if(ibloc(4).eq.1 .and. ibloc(1).ne.98) then
        call gbyte(ingrib,lengrib,iskp,24)
        call gbyte(ingrib,gribedi,iskp+24,8)
        write(*,*) '.. GRIB Edition 1 '
        ioff=ioff+8
        iskp=ioff*8
        call gbyte(ingrib,lenpds,iskp,24)
      elseif(ibloc(1).eq.0 .and. ibloc(2).eq.0 .and.
     &       ibloc(3).eq.24 .and. ibloc(4).eq.0 ) then
        call gbyte(ingrib,lengrib,iskp,24)
        call gbyte(ingrib,gribedi,iskp+24,8)
        write(*,*) '.. ECMWF GRIB Edition 0 '
        lenpds=lengrib
        ioff=ioff+4
        iskp=ioff*8
      elseif(ibloc(1).eq.98) then
        write(*,*) '.. ECMWF GRIB Edition X '
        gribedi=-1
        lenpds=20
        ioff=4
        iskp=ioff*8
      else
        write(*,*) '..unknown GRIB record, exit ',lengrib,idum
        call exit(0)
      endif
      write(*,*) ' GRIB, length,edi= ',lengrib,gribedi
      idum=lenpds+ioff
c
c at this time, ioff is the number of bytes before PDS
      write(*,*) 'lengrib,lenpds,ioff: ',lengrib,lenpds,ioff
c
c begin to process GRIB
c to check BLOK, loop through the whole ingrib array
c
c
c ioff is the number of bytes before PDS
c iskp is the number of bits before PDS
c
c process PDS
c
      call pds_GRIB1(iprnt,gribedi,iskp,idum,
     &   kgds,kbms,dscalf,ierr,lvl2,pds,ingrib)
c------------------------------------------------ select GRIB
c step 1: Look up nmc_grib.tbl2 for the variables wanted.
c         For example, U is 33, V is 34, Z is 7, T is 11, RH is 52.
c
c
c step 2: Check the value of pds(7).  If it is not one the variables
c         wanted, return and continue to the next GRIB record in the file.
c
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
c step 3: If a certain pressure level wanted, check pds(8), pds(9) and lvl2.
c         For a field at the surface, the codition is pds(8)=1.
c         From nmc_grib.tbl3 and GRIB doc, it is 850mb data when
c           pds(8)=100, pds(9)=850.
c         For a layer between 1000mb and 500mb, the conditions are
c           pds(8)=101, pds(9)=1000, lvl2=500.
c
c     if((pds(8).ne.1)) return                      !not at sfc, return
c     if((pds(8).ne.100).or.(pds(9).ne.850)) return !(not isobaric) or (not 850mb) return
c     if((pds(8).ne.101).or.(pds(9).ne.1000).or.    !(not pressure layer) or (not 1000mb) or
c    &   (lvl2.ne.500)) return                      !    (not 500mb), return
c
c Or, you may combine steps 2 and 3 together:
c For example, if you want surface pressure, pds(7)=1 and pds(8)=1
c     if((pds(7).ne.1).or.(pds(8).ne.1)) return  !(not pressure) or (not at sfc), return
c
c
c step 4: Check for the reference time.
c
c     if(pds(10).ne.92) return                      ! 1992, year
c     if(pds(11).ne.11) return                      ! 11, month
c     if(pds(12).ne.1) return                       ! 1, day
c     if(pds(13).ne.0) return                       ! 0, UTC hour
c------------------------------------------------ select GRIB
c
c
c To handle old GRIB Edition from ECMWF
c     if(gribedi.eq.-1) then
c       call gbyte(ingrib,kgds,7*8+7,1)
c       call gbyte(ingrib,kbms,7*8+6,1)
c       write(*,*) '....... revised kgds,kbms .....'
c       write(*,*) '.... kgds = ',kgds
c       write(*,*) '.... kbms = ',kbms
c       call gbytes(ingrib,ibloc,10*8,8,0,2)
c       if(pds(8).eq.100 .or. pds(8).eq.103 .or.
c    &     pds(8).eq.105 .or. pds(8).eq.107 .or.
c    &     pds(8).eq.109 .or. pds(8).eq.111 .or.
c    &     pds(8).eq.113 .or. pds(8).eq.125 .or.
c    &     pds(8).eq.160 .or. pds(8).eq.200 .or.
c    &     pds(8).eq.201 ) then
c         pds(9)=ibloc(1)*32+ibloc(2)
c         write(*,*) '....... revised level.....'
c         write(*,*) '....level = ',pds(9)
c       endif
c     endif
c
c
c process GDS
c
      if(kgds.eq.0) then
        lengds=0
        go to 40
      endif
      iskp=(pds(1)+ioff)*8
      call gbyte(ingrib,lengds,iskp,24)
c     write(*,*) ' GDS length: ',lengds
      call gds_GRIB1(iprnt,gribedi,lnin,iskp,
     &    kgds,ierr,nx,ny,la1,lo1,pds,ingrib)
c
c BMS section
c
   40 continue
      mapuse=0
      if(kbms.eq.0) then
        lenbms=0
c       go to 50
      else
        iskp=(lengds+pds(1)+ioff)*8
        call gbyte(ingrib,lenbms,iskp,24)
      endif
      call bms_GRIB1(iprnt,lnin,nmax,iskp,kbms,nx,
     &  ny,la1,lo1,ierr,numgd,mapuse,pds,ingrib,mapb)
c     write(*,*) ' BMS length: ',lenbms
c
c BDS section
c
   50 continue
      iskp=(lenbms+lengds+pds(1)+ioff)*8
      call gbyte(ingrib,lenbds,iskp,24)
c     write(*,*) '     bdslen: ',lenbds
c
c check the total length with GRIB header
c
      if(pds(3).eq.98 .and. gribedi.ne.1) then      ! ECMWF edition 0 or X '
        lenall=(4+lenpds+lengds+lenbms+lenbds+4)/1  ! 'GRIB+...+7777'
      else
        lenall=(8+pds(1)+lengds+lenbms+lenbds+4)/1
      endif
      if(iprnt.ne.0) write(*,*) 'lenall,pds(1),gds,bms,bds: ',
     &     ioff,pds(1),lengds,lenbms,lenbds,lenall
c     write(*,*) 'lnin = ',lnin
      short=' '
      if((lnin*nword).lt.lenall) then
        write(*,*) ' *** SHORT GRIB record ',lnin*nword,lenall
        if((lenall-(lnin*nword)).gt.4) then
          write(*,*) ' *** no bds_GRIB1, return '
          return
        endif
      endif
      if(asum.eq.' ') then
        write(*,*) 'enter 1 to continue ...'
        read(*,*) imore
        if(imore.ne.1) go to 9998
      endif
c     write(*,*) '..reading BDS, iskp = ',iskp
c
      call bds_GRIB1(iprnt,gribedi,lnin,iskp,dscalf,ierr,nmax,
     &  npoints,rc00,rmax,rmin,numgd,mapuse,mapb,pds,ingrib,
     &  ix,datan)
      if(ierr.eq.-3) then
        write(*,*) '****Too many grids than specified..'
        ierr=0
        return
      endif
      if(ierr.eq.-8) then
        write(*,*) '****Constant field, all grids = ',datan(1)
        write(*,*) '    No number of grids can be computed from bds'
        write(*,*) '    Use number of grids from gds section'
        ierr=0
        return
      endif
      if(ierr.eq.-11) then
        write(*,*) '****UNKNOWN predefined grid type, '
        write(*,*) '    No grids decoded'
        ierr=0
        return
      endif
      if(ierr.eq.-12) then
        write(*,*) '****COMPLEX PACKING of spherical harmonic '
        write(*,*) '    coefficients, to be developed!'
        write(*,*) '    No grids decoded'
        ierr=0
        return
      endif
      if(asum.eq.' ') write(*,83) (datan(k),k=1,15)
   83 format(4(E15.8,1X),E15.8)

 9998 continue
      if(asum.eq.' ')
     &  write(*,*) ' finish reading input grib record....'
c
c make one line summary
c
      if(asum.ne.' ')
     &  write(asum(7:132),3213) (pds(i),i=1,5),(pds(j),j=7,9),
     &    lvl2,(pds(i),i=10,13),(pds(k),k=15,18),pds(21),
     &    pds(22),lenbds,npoints,rmax,rmin,short
 3213 format(i3,i4,2i3,3i4,2i5,1x,4i2.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
c------------------------------------- select a region, say 40N-20S, 130W-70W
c 1. Check mapuse to detemine if a bitmap is used.
c    If a bitmap is used, npoints will be less than (nx*ny), so fill to grd1(nx*ny).
c 2. Remap grd1(1..npoints) to grd2(nx,ny).  If grd1 are gaussian grids, convert 
c    them to 2.5x2.5 lat,long grids.
c 3. Print the location of grd1(1), i.e. la1,lo1
c 4. Determine the starting,ending (i,j) on 2.5x2.5
c 5. Save the regional grids ((grd2(i,j),i=ibeg,iend),j=jbeg,jend)
c
c step 1:
c     If(mapuse.gt.0) then
c       if(mapuse.ne.(nx*ny)) then
c         write(*,*) '**** ERROR in select a region w/ bitmap ',mapuse,nx*ny,nx,ny
c         call exit(0)
c       endif
c       kcnt=0
c       do i=1,mapuse
c         if(mapb(i).eq.1) then
c            kcnt=kcnt+1
c            grd1(i)=datan(kcnt)
c         else
c            grd1(i)=-9999.
c         endif
c       enddo
c     else
c       if(npoints.ne.(nx*ny)) then
c         write(*,*) '**** ERROR in select a region, ',npoints,nx*ny,nx,ny
c         call exit(0)
c       endif
c       do i=1,npoints
c         grd1(i)=datan(i)
c       enddo
c     endif
c
c step 2:
c     if((nx*ny).ne.(144*73)) then
c       if((nx.ne.192) .or. (ny.ne.94)) then
c         write(*,*) '****ERROR, select a region only works for '
c         write(*,*) '           144x73 (2.5x2.5 degree) grids, or'
c         write(*,*) '           192x94 (gaussian) grids now......'
c         call exit(0)
c       endif
c       call interpolation(nx,ny,grd2,grd1)
c     else
c       do j=1,ny
c         do i=1,nx
c           grd2(i,j)=grd1((j-1)*nx+i)
c         enddo
c       enddo
c     endif
c
c step 3:
c     xlat1=float(la1)/1000.
c     xlon1=float(lo1)/1000.
c     write(*,*) '....grd2(1,1) is located at lat,long ',xlat1,xlon1
c
c step 4: 
c     if((la1 .ne. 90000) .or. (lo1 .ne. 0)) then
c       write(*,*) '****WARNING, adjust la1,lo1 to 90N,0E '
c       xlat1=90.
c       xlon1=0.
c     endif
c     slat=40.
c     elat=-20.
c     slong=360.-130.
c     elong=360.-70.
c     ibeg=(slong-xlon1)/2.5+1
c     iend=(elong-xlon1)/2.5+1
c     jbeg=(xlat1-slat)/2.5+1
c     jend=(xlat1-elat)/2.5+1
c
c step 5: save data
c     write(*,*) ((grd2(i,j),i=ibeg,iend),j=jbeg,jend)
c
c------------------------------------- select a region

      return
      end
c
c
c-- lines below are from ncep_int.f which is available in our anonymous ftp 
c   area on ncardata.ucar.edu in 'ftp/datasets/ds090.0/ncep_int.f' file.
c
      subroutine interpolation (nnx,nny,grdo2,grdi)
c
c This program demonstrates how to use the interpolation routines 
c (same as NCEP's) to get 2.5x2.5 grids from Gaussian (192x94) grids.
c
      parameter (mxx=192, mxy=94, imo=144,jmo=73)
      dimension grdi(mxx*mxy), grdo(imo*jmo),grdo2(imo,jmo)
      dimension gglat(mxy/2), bmapi(mxx*mxy),bmapo(mxx*mxy)
      real PK(mxy/2),PKM1(mxy/2)
      real SJN(mxy/2),WJN(mxy/2),CJN1(mxy/2)
      real AJN(mxy+1),AJO(jmo+1),XIO(imo+1),XJO(jmo+1)
      real WIO1(imo),WIO2(imo)
      integer JPO(jmo)
      real FPX(imo),FP0X(imo),WPX(imo),WTX(imo)
      character*72 inf, outf
c
c The file 'T62.gauss_lat' is available from anonymous ftp on ncardata.ucar.edu
c under 'ftp/datasets/ds090.0/' directory.  It can be generated here as well.
c
c     open(unit=21,file='T62.gauss_lat')
c     do i=1,14
c       read(21,*)
c     enddo
c     do i=1,47
c       read(21,*) idum,gglat(i)
c     enddo
c     close(unit=21)
c
      gglat(1)=0.
c
c     kin=1
c     write(*,*) 'enter input data file?'
c     read(*,'(a)') inf
c     open(kin,file=inf)
c     write(*,*) 'enter input data dimension? <192,94>'
c     read(*,*) nnx,nny
      if(nnx.gt.mxx .or. (nny.gt.mxy)) then
        write(*,*) '****input dimension too large, exit'
        write(*,*) nnx,nny,'   ',mxx,mxy
        call exit(0)
      endif
      rmax=-999999.0
      rmin= 999999.0
      k=0
      do j=1,nny
        do i=1,nnx
          k=k+1
c         read(kin,*) idum,grdi(k)
          if(grdi(k).ge.rmax) then
            rmax=grdi(k)
c           kmax=idum
          endif
          if(grdi(k).le.rmin) then
            rmin=grdi(k)
c           kmin=idum
          endif
        enddo
      enddo
c     close(kin)
c     write(*,*) 'last idum = ',idum
c     write(*,*) 'max: ',kmax,rmax,'  min: ',kmin,rmin
c
      mp=0 !0:scalar, 1:vector, 2:flag, 3:budget
      mb=0 !0:no bitmap
      imn=192
      ixn=1
      jmn=94
      jxn=192
      ixo=1
      jxo=144
      cio=0.    ! starting longitude in degrees East
      cjo=90.   ! starting latitude in degrees North
      dio=2.5   ! longitude increment in degrees East
      djo=-2.5  ! latitude increment in degrees North
      call gg2ll(mp,mb,nnx,ixn,nny,jxn,gglat,bmapi,grdi,
     &           imo,ixo,jmo,jxo,cio,cjo,dio,djo,bmapo,grdo,
     &           sjn,wjn,cjn1,ajn,ajo,xio,xjo,wio1,wio2,jpo,
     &           fpx,wpx,wtx,fp0x,PK,PKM1)
c
c     do j=1,jmn/2
c       write(*,*) ' Gaussian lat ',j,gglat(j)
c     enddo
      k=0
      do j=1,jmo
        do i=1,imo
          k=k+1
          grdo2(i,j)=grdo(k)
        enddo
      enddo
c
c save interpolated grids into a GrADS (or IEEE) data file
c
c     write(*,*) 'enter output data file?'
c     read(*,'(a)') outf
c     open(unit=23,file=outf,form='unformatted',access='direct',
c    &     recl=imo*jmo*4)
c     write(23,rec=1) ((grdo2(i,j),i=1,imo),j=jmo,1,-1)
c     close(unit=23)
c     stop
c
      return
      end
C-----------------------------------------------------------------------
CFPP$ NOCONCUR R
      SUBROUTINE GG2LL(MP,MB,IMN,IXN,JMN,JXN,CJN,LN,FN,
     &                 IMO,IXO,JMO,JXO,CIO,CJO,DIO,DJO,LO,FO,
     &                 sjn,wjn,cjn1,ajn,ajo,xio,xjo,wio1,wio2,jpo,
     &                 fpx,wpx,wtx,fp0x,PK,PKM1)
c
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    GG2LL       INTERPOLATE GAUSSIAN TO LAT-LON GRID.
C   PRGMMR: IREDELL          ORG: W/NMC23    DATE: 92-10-31
C
C ABSTRACT: A HORIZONTAL FIELD IS BILINEARLY INTERPOLATED
C   FROM A GLOBAL GAUSSIAN GRID TO A LATITUDE-LONGITUDE GRID.
C   THE GAUSSIAN LATITUDES MAY BE PASSED OR CALCULATED.
C   THE INPUT GLOBAL FIELD IS ASSUMED TO RUN EASTWARD AND SOUTHWARD,
C   STARTING AT GREENWICH AND THE NORTHERNMOST LATITUDE,
C   BUT THE OUTPUT FIELD ORIENTATION IS DEFINED BY PASSED ARGUMENTS.
C   THE FIELD MAY BE TAGGED A SCALAR, A VECTOR, A FLAG OR A BUDGET.
C   A FLAG FIELD SUCH AS THE LAND-SEA MASK IS NOT INTERPOLATED
C   BUT TAKEN FROM THE VALUES AT THE CLOSEST INPUT GRIDPOINTS.
C   A BUDGET FIELD SUCH AS PRECIPITATION IS INTERPOLATED
C   WHILE PRESERVING THE AREA INTEGRALS OF THE ORIGINAL FIELD.
C   POLAR SCALARS ARE THE AVERAGE OF THE CLOSEST LATITUDE CIRCLE VALUES.
C   POLAR VECTOR COMPONENTS ARE TAKEN FROM THE WAVENUMBER 1 COMPONENT
C   EXTRACTED FROM THE VALUES ON THE CLOSEST LATITUDE CIRCLE.
C   A LOGICAL BITMAP MAY BE PASSED TO MASK OUT PART OF THE INPUT FIELD.
C   IN THIS CASE, AN APPROPRIATE OUTPUT BITMAP IS CONSTRUCTED AND
C   THE OUTPUT FIELD IS COMPUTED ONLY WHERE THE OUTPUT BITMAP IS TRUE
C   AND ONLY USING DATA WHERE THE INPUT BITMAP IS TRUE.
C   WHERE THE OUTPUT BITMAP IS FALSE, THE OUTPUT IS SET TO ZERO.
C   THE BITMAP MAY ALTERNATIVELY BE CONSIDERED AS A TWO-WAY MASK,
C   IN WHICH CASE OUTPUT FIELD IS ALSO COMPUTED WHERE THE OUTPUT BITMAP
C   IS FALSE ONLY USING DATA WHERE THE INPUT BITMAP IS FALSE.
C   THE TWO-WAY MASK OPTION SHOULD NOT BE USED WITH VECTOR COMPONENTS.
C
C PROGRAM HISTORY LOG:
C   92-10-31  IREDELL
C   93-10-21  IREDELL   ALLOW VECTORS, FLAGS, BITMAP, NORTH/SOUTH FLIP
C   94-12-05  IREDELL   ALLOW BUDGETS, OPTIONALLY CALCULATE LATITUDES,
C                       PRECOMPUTE INDICES, GENERALIZE OUTPUT DOMAIN.
C   95-07-11  IREDELL   REWEIGHT VECTOR INTERPOLATION NEAR POLE.
C
C USAGE:    CALL GG2LL(MP,MB,IMN,IXN,JMN,JXN,CJN,LN,FN,
C    &                 IMO,IXO,JMO,JXO,CIO,CJO,DIO,DJO,LO,FO)
C   INPUT ARGUMENT LIST:
C     MP       - INTEGER FIELD PARAMETER IDENTIFIER
C                (0 FOR SCALAR, 1 FOR VECTOR, 2 FOR FLAG, 3 FOR BUDGET)
C     MB       - INTEGER BITMAP IDENTIFIER
C                (0 FOR NO BITMAP, 1 TO INTERPOLATE BITMAP,
C                 2 TO ALSO INTERPOLATE DATA TO BITMAP FALSES)
C     IMN      - INTEGER INPUT LONGITUDE DIMENSION
C     IXN      - INTEGER NUMBER TO SKIP BETWEEN INPUT LONGITUDES
C     JMN      - INTEGER INPUT LATITUDE DIMENSION (EVEN)
C     JXN      - INTEGER NUMBER TO SKIP BETWEEN INPUT LATITUDES
C     CJN      - REAL (JMN/2) GAUSSIAN LATITUDES IN DEGREES NORTH
C                (IF CJN(1)=0., THE LATITUDES ARE CALCULATED.)
C     LN       - LOGICAL ((IMN-1)*IXN+(JMN-1)*JXN+1) BITMAP IF MB=1
C     FN       - REAL ((IMN-1)*IXN+(JMN-1)*JXN+1) FIELD TO INTERPOLATE
C     IMO      - INTEGER OUTPUT LONGITUDE DIMENSION
C     IXO      - INTEGER NUMBER TO SKIP BETWEEN OUTPUT LONGITUDES
C     JMO      - INTEGER OUTPUT LATITUDE DIMENSION
C     JXO      - INTEGER NUMBER TO SKIP BETWEEN OUTPUT LATITUDES
C     CIO      - REAL START LONGITUDE IN DEGREES EAST
C     CJO      - REAL START LATITUDE IN DEGREES NORTH
C     DIO      - REAL LONGITUDE INCREMENT IN DEGREES EAST
C     DJO      - REAL LATITUDE INCREMENT IN DEGREES NORTH
C
C   OUTPUT ARGUMENT LIST:
C     LO       - LOGICAL ((IMO-1)*IXO+(JMO-1)*JXO+1) BITMAP IF MB=1
C     FO       - REAL ((IMO-1)*IXO+(JMO-1)*JXO+1) INTERPOLATED FIELD
C
C SUBPROGRAMS CALLED:
C   GLAT         COMPUTE GAUSSIAN LATITUDES
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN 77
C
C$$$
      REAL CJN(JMN/2)
      REAL FN((IMN-1)*IXN+(JMN-1)*JXN+1)
      REAL FO((IMO-1)*IXO+(JMO-1)*JXO+1)
      LOGICAL LN((IMN-1)*IXN+(JMN-1)*JXN+1)
      LOGICAL LO((IMO-1)*IXO+(JMO-1)*JXO+1)
      REAL SJN(JMN/2),WJN(JMN/2),CJN1(JMN/2)
      REAL PK(JMN/2),PKM1(JMN/2)
      REAL AJN(JMN+1),AJO(JMO+1),XIO(IMO+1),XJO(JMO+1)
      REAL WIO1(IMO),WIO2(IMO)
      INTEGER JPO(JMO)
      REAL FPX(IMO),FP0X(IMO),WPX(IMO),WTX(IMO)
      IJN(I,J)=(I-1)*IXN+(J-1)*JXN+1
      IJO(I,J)=(I-1)*IXO+(J-1)*JXO+1
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  PRECOMPUTE INDICES
      PI=ACOS(-1.)
      IF(CJN(1).EQ.0.) THEN
        write(*,*) '****No Gaussian latitudes, calculate them...'
        CALL GLAT(JMN/2,SJN,WJN,PK,PKM1)
        DO J1=1,JMN/2
          CJN1(J1)=(180/PI)*ASIN(SJN(J1))
          cjn(j1)=cjn1(j1)
        ENDDO
      ELSE
        DO J1=1,JMN/2
          CJN1(J1)=CJN(J1)
        ENDDO
      ENDIF
      DO J=1,JMO
        JPO(J)=0
      ENDDO
      IF(MP.NE.3) THEN
C  SCALAR OR VECTOR OR FLAG
        DO I=1,IMO
          XIO(I)=MOD((CIO+(I-1)*DIO)/360*IMN+2*IMN,FLOAT(IMN))+1
        ENDDO
        IF(MP.EQ.0) THEN
          DO I=1,IMO
            I1=XIO(I)
            WIO1(I)=I1+1-XIO(I)
            WIO2(I)=XIO(I)-I1
          ENDDO
        ELSEIF(MP.EQ.1) THEN
          DL=2*PI/IMN
          DO I=1,IMO
            I1=XIO(I)
            WIO1(I)=SIN((I1+1-XIO(I))*DL)/SIN(DL)
            WIO2(I)=SIN((XIO(I)-I1)*DL)/SIN(DL)
          ENDDO
        ENDIF
CDIR$ IVDEP
        DO J1=1,JMN/2
          AJN(J1)=CJN1(J1)
          AJN(JMN+1-J1)=-CJN1(J1)
        ENDDO
        DO J=1,JMO
          AJO(J)=CJO+(J-1)*DJO
        ENDDO
        IF(DJO.LT.0.) THEN
          JB=1
          JE=JMO
          JI=1
        ELSE
          JB=JMO
          JE=1
          JI=-1
        ENDIF
        J1=1
        DO J=JB,JE,JI
          IF(AJO(J).GE.AJN(1)) THEN
            XJO(J)=1
          ELSEIF(AJO(J).LE.AJN(JMN)) THEN
            XJO(J)=JMN
          ELSE
            DOWHILE(AJO(J).LT.AJN(J1+1))
              J1=J1+1
            ENDDO
            XJO(J)=J1+(AJN(J1)-AJO(J))/(AJN(J1)-AJN(J1+1))
          ENDIF
        ENDDO
        IF(XJO(1).EQ.1.) THEN
          JPO(1)=1
        ELSEIF(XJO(1).EQ.FLOAT(JMN)) THEN
          JPO(1)=JMN
        ENDIF
        IF(XJO(JMO).EQ.1.) THEN
          JPO(JMO)=1
        ELSEIF(XJO(JMO).EQ.FLOAT(JMN)) THEN
          JPO(JMO)=JMN
        ENDIF
      ELSE
C  BUDGET
        DO I=1,IMO+1
          XIO(I)=MOD((CIO+(I-1.5)*DIO)/360*IMN+2*IMN+0.5,FLOAT(IMN))+1
        ENDDO
        IBX=0
        DO I=1,IMO
          XIB=XIO(INT(I+0.5-SIGN(0.5,DIO)))
          XIE=XIO(INT(I+0.5+SIGN(0.5,DIO)))
          IB=XIB
          IE=XIE
          IF(IE.LT.IB) IE=IE+IMN
          IBX=MAX(IBX,IE-IB+1)
        ENDDO
        AJN(1)=1.
        AJN(JMN/2+1)=0.
        AJN(JMN+1)=-1.
CDIR$ IVDEP
        DO J1=2,JMN/2
          AJN(J1)=SIN((PI/180)*0.5*(CJN1(J1-1)+CJN1(J1)))
          AJN(JMN+2-J1)=-AJN(J1)
        ENDDO
        DO J=1,JMO+1
          PJO=CJO+(J-1.5)*DJO
          IF(PJO.GE.90.) THEN
            AJO(J)=1.
          ELSEIF(PJO.LE.-90.) THEN
            AJO(J)=-1.
          ELSE
            AJO(J)=SIN((PI/180)*PJO)
          ENDIF
        ENDDO
        IF(DJO.LT.0.) THEN
          JB=1
          JE=JMO+1
          JI=1
        ELSE
          JB=JMO+1
          JE=1
          JI=-1
        ENDIF
        J1=1
        DO J=JB,JE,JI
          IF(AJO(J).GE.AJN(1)) THEN
            XJO(J)=1
          ELSEIF(AJO(J).LE.AJN(JMN+1)) THEN
            XJO(J)=JMN+1
          ELSE
            DOWHILE(AJO(J).LT.AJN(J1+1))
              J1=J1+1
            ENDDO
            XJO(J)=J1+(AJN(J1)-AJO(J))/(AJN(J1)-AJN(J1+1))
          ENDIF
        ENDDO
        IF(XJO(1).EQ.1.) THEN
          JPO(1)=1
        ELSEIF(XJO(1).EQ.FLOAT(JMN)) THEN
          JPO(1)=JMN
        ENDIF
        IF(XJO(JMO+1).EQ.1.) THEN
          JPO(JMO)=1
        ELSEIF(XJO(JMO+1).EQ.FLOAT(JMN+1)) THEN
          JPO(JMO)=JMN
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  GET POLAR VALUES
      DO J=1,JMO,JMO-1
        IF(JPO(J).GT.0) THEN
          IF(MB.EQ.0) THEN
            IF(MP.EQ.0) THEN
C  FULL SCALAR
              J1=JPO(J)
              FP=0.
              DO I=1,IMN
                FP=FP+FN(IJN(I,J1))
              ENDDO
              FP=FP/IMN
              DO I=1,IMO
                FO(IJO(I,J))=FP
              ENDDO
            ELSEIF(MP.EQ.1) THEN
C  FULL VECTOR
              J1=JPO(J)
              FPC=0.
              FPS=0.
              DO I=1,IMN
                CI=COS(2*PI*(I-1)/IMN)
                SI=SIN(2*PI*(I-1)/IMN)
                FPC=FPC+CI*FN(IJN(I,J1))
                FPS=FPS+SI*FN(IJN(I,J1))
              ENDDO
              FPC=2*FPC/IMN
              FPS=2*FPS/IMN
              DO I=1,IMO
                CI=COS(2*PI*(I-1)/IMO)
                SI=SIN(2*PI*(I-1)/IMO)
                FO(IJO(I,J))=FPC*CI+FPS*SI
              ENDDO
            ELSEIF(MP.EQ.2) THEN
C  FULL FLAG
              J1=JPO(J)
              DO I=1,IMO
                FO(IJO(I,J))=FN(IJN(1,J1))
              ENDDO
            ELSEIF(MP.EQ.3) THEN
C  FULL BUDGET
              XJB=XJO(INT(J+0.5+SIGN(0.5,DJO)))
              XJE=XJO(INT(J+0.5-SIGN(0.5,DJO)))
              JB=XJB
              JE=XJE
              WJB=JB+1-XJB
              WJE=XJE-JE
              FP=0.
              WP=0.
              DO J1=JB,JE
                W=AJN(J1)-AJN(J1+1)
                IF(J1.EQ.JB) W=W*WJB
                IF(J1.EQ.JE) W=W*WJE
                DO I=1,IMN
                  FP=FP+W*FN(IJN(I,J1))
                ENDDO
                WP=WP+W*IMN
              ENDDO
              DO I=1,IMO
                FO(IJO(I,J))=FP/WP
              ENDDO
            ENDIF
          ELSE
            IF(MP.EQ.0) THEN
C  BITMAP OR TWOWAY SCALAR
              J1=JPO(J)
              FP=0.
              WP=0.
              IF(MB.EQ.2) FP0=0.
              DO I=1,IMN
                IF(LN(IJN(I,J1))) THEN
                  FP=FP+FN(IJN(I,J1))
                  WP=WP+1.
                ELSEIF(MB.EQ.2) THEN
                  FP0=FP0+FN(IJN(I,J1))
                ENDIF
              ENDDO
              DO I=1,IMO
                LO(IJO(I,J))=WP.GE.0.5*IMN
              ENDDO
              IF(LO(IJO(1,J))) THEN
                DO I=1,IMO
                  FO(IJO(I,J))=FP/WP
                ENDDO
              ELSEIF(MB.EQ.2.AND.WP.LT.IMN) THEN
                DO I=1,IMO
                  FO(IJO(I,J))=FP0/(IMN-WP)
                ENDDO
              ELSE
                DO I=1,IMO
                  FO(IJO(I,J))=0.
                ENDDO
              ENDIF
            ELSEIF(MP.EQ.1) THEN
C  BITMAP VECTOR
              J1=JPO(J)
              IP=0
              DO I=1,IMN
                IF(LN(IJN(I,J1))) IP=IP+1
              ENDDO
              DO I=1,IMO
                LO(IJO(I,J))=IP.EQ.IMN
              ENDDO
              IF(LO(IJO(1,J))) THEN
                FPC=0.
                FPS=0.
                DO I=1,IMN
                  CI=COS(2*PI*(I-1)/IMN)
                  SI=SIN(2*PI*(I-1)/IMN)
                  FPC=FPC+CI*FN(IJN(I,J1))
                  FPS=FPS+SI*FN(IJN(I,J1))
                ENDDO
                FPC=2*FPC/IMN
                FPS=2*FPS/IMN
                DO I=1,IMO
                  CI=COS(2*PI*(I-1)/IMO)
                  SI=SIN(2*PI*(I-1)/IMO)
                  FO(IJO(I,J))=FPC*CI+FPS*SI
                ENDDO
              ELSE
                DO I=1,IMO
                  FO(IJO(I,J))=0.
                ENDDO
              ENDIF
            ELSEIF(MP.EQ.2) THEN
C  BITMAP OR TWOWAY FLAG
              J1=JPO(J)
              DO I=1,IMO
                LO(IJO(I,J))=LN(IJN(1,J1))
                IF(LN(IJN(1,J1)).OR.MB.EQ.2) THEN
                  FO(IJO(I,J))=FN(IJN(1,J1))
                ELSE
                  FO(IJO(I,J))=0.
                ENDIF
              ENDDO
            ELSEIF(MP.EQ.3) THEN
C  BITMAP OR TWOWAY BUDGET
              XJB=XJO(INT(J+0.5+SIGN(0.5,DJO)))
              XJE=XJO(INT(J+0.5-SIGN(0.5,DJO)))
              JB=XJB
              JE=XJE
              WJB=JB+1-XJB
              WJE=XJE-JE
              FP=0.
              WP=0.
              WT=0.
              IF(MB.EQ.2) THEN
                FP0=0.
              ENDIF
              DO J1=JB,JE
                W=AJN(J1)-AJN(J1+1)
                IF(J1.EQ.JB) W=W*WJB
                IF(J1.EQ.JE) W=W*WJE
                DO I=1,IMN
                  IF(LN(IJN(I,J1))) THEN
                    FP=FP+W*FN(IJN(I,J1))
                    WP=WP+W
                  ELSEIF(MB.EQ.2) THEN
                    FP0=FP0+W*FN(IJN(I,J1))
                  ENDIF
                ENDDO
                WT=WT+W*IMN
              ENDDO
              DO I=1,IMO
                LO(IJO(I,J))=WP.GE.0.5*WT
              ENDDO
              IF(LO(IJO(1,J))) THEN
                DO I=1,IMO
                  FO(IJO(I,J))=FP/WP
                ENDDO
              ELSEIF(MB.EQ.2.AND.WP.LT.WT) THEN
                DO I=1,IMO
                  FO(IJO(I,J))=FP0/(WT-WP)
                ENDDO
              ELSE
                DO I=1,IMO
                  FO(IJO(I,J))=0.
                ENDDO
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  INTERPOLATE THE REST OF THE FIELD
      DO J=1,JMO
        IF(JPO(J).EQ.0) THEN
          IF(MB.EQ.0) THEN
            IF(MP.EQ.0.OR.MP.EQ.1) THEN
C  FULL SCALAR OR VECTOR
              J1=XJO(J)
              J2=MIN(J1+1,JMN)
              WJ2=XJO(J)-J1
              WJ1=1.-WJ2
              DO I=1,IMO
                I1=XIO(I)
                I2=MOD(I1,IMN)+1
                WI1=WIO1(I)
                WI2=WIO2(I)
                W11=WI1*WJ1
                W21=WI2*WJ1
                W12=WI1*WJ2
                W22=WI2*WJ2
                FO(IJO(I,J))=W11*FN(IJN(I1,J1))+W21*FN(IJN(I2,J1))+
     &                       W12*FN(IJN(I1,J2))+W22*FN(IJN(I2,J2))
              ENDDO
            ELSEIF(MP.EQ.2) THEN
C  FULL FLAG
              J1=NINT(XJO(J))
              DO I=1,IMO
                I1=MOD(NINT(XIO(I))-1,IMN)+1
                FO(IJO(I,J))=FN(IJN(I1,J1))
              ENDDO
            ELSEIF(MP.EQ.3) THEN
C  FULL BUDGET
              XJB=XJO(INT(J+0.5+SIGN(0.5,DJO)))
              XJE=XJO(INT(J+0.5-SIGN(0.5,DJO)))
              JB=XJB
              JE=XJE
              WJB=JB+1-XJB
              WJE=XJE-JE
              DO I=1,IMO
                FPX(I)=0.
                WPX(I)=0.
              ENDDO
              DO J1=JB,JE
                WJ=AJN(J1)-AJN(J1+1)
                IF(J1.EQ.JB) WJ=WJ*WJB
                IF(J1.EQ.JE) WJ=WJ*WJE
                DO I1X=1,IBX
                  DO I=1,IMO
                    XIB=XIO(INT(I+0.5-SIGN(0.5,DIO)))
                    XIE=XIO(INT(I+0.5+SIGN(0.5,DIO)))
                    IB=XIB
                    IE=XIE
                    WIB=IB+1-XIB
                    WIE=XIE-IE
                    IF(IE.LT.IB) IE=IE+IMN
                    I1=IB+(I1X-1)
                    IF(I1.LE.IE) THEN
                      W=WJ
                      IF(I1.EQ.IB) W=W*WIB
                      IF(I1.EQ.IE) W=W*WIE
                      I1M=MOD(I1-1,IMN)+1
                      FPX(I)=FPX(I)+W*FN(IJN(I1M,J1))
                      WPX(I)=WPX(I)+W
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
              DO I=1,IMO
                FO(IJO(I,J))=FPX(I)/WPX(I)
              ENDDO
            ENDIF
          ELSE
            IF(MP.EQ.0.OR.MP.EQ.1) THEN
C  BITMAP OR TWOWAY SCALAR OR VECTOR
              J1=XJO(J)
              J2=MIN(J1+1,JMN)
              WJ2=XJO(J)-J1
              WJ1=1.-WJ2
              DO I=1,IMO
                I1=XIO(I)
                I2=MOD(I1,IMN)+1
                WI1=WIO1(I)
                WI2=WIO2(I)
                W11=WI1*WJ1
                W21=WI2*WJ1
                W12=WI1*WJ2
                W22=WI2*WJ2
                WP=0.
                IJ11=IJN(I1,J1)
                IJ21=IJN(I2,J1)
                IJ12=IJN(I1,J2)
                IJ22=IJN(I2,J2)
                IF(LN(IJ11)) WP=WP+W11
                IF(LN(IJ21)) WP=WP+W21
                IF(LN(IJ12)) WP=WP+W12
                IF(LN(IJ22)) WP=WP+W22
                WT=W11+W21+W12+W22
                IJ=IJO(I,J)
                LO(IJ)=WP.GE.0.5*WT
                IF(LO(IJ)) THEN
                  FP=0.
                  IF(LN(IJ11)) FP=FP+W11*FN(IJ11)
                  IF(LN(IJ21)) FP=FP+W21*FN(IJ21)
                  IF(LN(IJ12)) FP=FP+W12*FN(IJ12)
                  IF(LN(IJ22)) FP=FP+W22*FN(IJ22)
                  FO(IJ)=FP*WT/WP
                ELSEIF(MB.EQ.2.AND.WP.LT.WT) THEN
                  FP=0.
                  IF(.NOT.LN(IJ11)) FP=FP+W11*FN(IJ11)
                  IF(.NOT.LN(IJ21)) FP=FP+W21*FN(IJ21)
                  IF(.NOT.LN(IJ12)) FP=FP+W12*FN(IJ12)
                  IF(.NOT.LN(IJ22)) FP=FP+W22*FN(IJ22)
                  FO(IJ)=FP*WT/(WT-WP)
                ELSE
                  FO(IJ)=0.
                ENDIF
              ENDDO
            ELSEIF(MP.EQ.2) THEN
C  BITMAP OR TWOWAY FLAG
              J1=NINT(XJO(J))
              DO I=1,IMO
                I1=MOD(NINT(XIO(I))-1,IMN)+1
                LO(IJO(I,J))=LN(IJN(I1,J1))
                IF(LO(IJO(I,J)).OR.MB.EQ.2) THEN
                  FO(IJO(I,J))=FN(IJN(I1,J1))
                ELSE
                  FO(IJO(I,J))=0.
                ENDIF
              ENDDO
            ELSEIF(MP.EQ.3) THEN
C  BITMAP OR TWOWAY BUDGET
              XJB=XJO(INT(J+0.5+SIGN(0.5,DJO)))
              XJE=XJO(INT(J+0.5-SIGN(0.5,DJO)))
              JB=XJB
              JE=XJE
              WJB=JB+1-XJB
              WJE=XJE-JE
              DO I=1,IMO
                FPX(I)=0.
                IF(MB.EQ.2) FP0X(I)=0.
                WPX(I)=0.
                WTX(I)=0.
              ENDDO
              DO J1=JB,JE
                WJ=AJN(J1)-AJN(J1+1)
                IF(J1.EQ.JB) WJ=WJ*WJB
                IF(J1.EQ.JE) WJ=WJ*WJE
                DO I1X=1,IBX
                  DO I=1,IMO
                    XIB=XIO(INT(I+0.5-SIGN(0.5,DIO)))
                    XIE=XIO(INT(I+0.5+SIGN(0.5,DIO)))
                    IB=XIB
                    IE=XIE
                    WIB=IB+1-XIB
                    WIE=XIE-IE
                    IF(IE.LT.IB) IE=IE+IMN
                    I1=IB+(I1X-1)
                    IF(I1.LE.IE) THEN
                      W=WJ
                      IF(I1.EQ.IB) W=W*WIB
                      IF(I1.EQ.IE) W=W*WIE
                      I1M=MOD(I1-1,IMN)+1
                      IF(LN(IJN(I1M,J1))) THEN
                        FPX(I)=FPX(I)+W*FN(IJN(I1M,J1))
                        WPX(I)=WPX(I)+W
                      ELSEIF(MB.EQ.2) THEN
                        FP0X(I)=FP0X(I)+W*FN(IJN(I1M,J1))
                      ENDIF
                      WTX(I)=WTX(I)+W
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
              DO I=1,IMO
                LO(IJO(I,J))=WPX(I).GE.0.5*WTX(I)
                IF(LO(IJO(I,J))) THEN
                  FO(IJO(I,J))=FPX(I)/WPX(I)
                ELSEIF(MB.EQ.2.AND.WPX(I).LT.WTX(I)) THEN
                  FO(IJO(I,J))=FP0X(I)/(WTX(I)-WPX(I))
                ELSE
                  FO(IJO(I,J))=0.
                ENDIF
              ENDDO
            ENDIF
          ENDIF
        ENDIF
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END

C-----------------------------------------------------------------------
      SUBROUTINE GLAT(JH,SLAT,WLAT,PK,PKM1)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    GLAT        COMPUTE GAUSSIAN LATITUDE FUNCTIONS
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 92-10-31
C
C ABSTRACT: COMPUTES SINES OF GAUSSIAN LATITUDE BY ITERATION.
C           THE GAUSSIAN WEIGHTS ARE ALSO COMPUTED.
C
C PROGRAM HISTORY LOG:
C   91-10-31  MARK IREDELL
C
C USAGE:    CALL GLAT(JH,SLAT,WLAT)
C
C   INPUT ARGUMENT LIST:
C     JH       - INTEGER NUMBER OF GAUSSIAN LATITUDES IN A HEMISPHERE
C
C   OUTPUT ARGUMENT LIST:
C     SLAT     - REAL (JH) SINES OF (POSITIVE) GAUSSIAN LATITUDE
C     WLAT     - REAL (JH) GAUSSIAN WEIGHTS FOR THE NH
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN
C
C$$$
      DIMENSION SLAT(JH),WLAT(JH)
c     PARAMETER(PI=3.14159265358979,C=(1.-(2./PI)**2)*0.25,EPS=1.E-14)
      PARAMETER(PI=3.14159265358979,C=(1.-(2./PI)**2)*0.25,EPS=1.E-7)
      PARAMETER(JBZ=50)
      DIMENSION PK(JH),PKM1(JH),BZ(JBZ)
      DATA BZ        / 2.4048255577,  5.5200781103,
     $  8.6537279129, 11.7915344391, 14.9309177086, 18.0710639679,
     $ 21.2116366299, 24.3524715308, 27.4934791320, 30.6346064684,
     $ 33.7758202136, 36.9170983537, 40.0584257646, 43.1997917132,
     $ 46.3411883717, 49.4826098974, 52.6240518411, 55.7655107550,
     $ 58.9069839261, 62.0484691902, 65.1899648002, 68.3314693299,
     $ 71.4729816036, 74.6145006437, 77.7560256304, 80.8975558711,
     $ 84.0390907769, 87.1806298436, 90.3221726372, 93.4637187819,
     $ 96.6052679510, 99.7468198587, 102.888374254, 106.029930916,
     $ 109.171489649, 112.313050280, 115.454612653, 118.596176630,
     $ 121.737742088, 124.879308913, 128.020877005, 131.162446275,
     $ 134.304016638, 137.445588020, 140.587160352, 143.728733573,
     $ 146.870307625, 150.011882457, 153.153458019, 156.295034268 /
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  ESTIMATE LATITUDES USING BESSEL FUNCTION
      R=1./SQRT((2*JH+0.5)**2+C)
      DO J=1,MIN(JH,JBZ)
        SLAT(J)=COS(BZ(J)*R)
      ENDDO
      DO J=JBZ+1,JH
        SLAT(J)=COS((BZ(JBZ)+(J-JBZ)*PI)*R)
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CONVERGE UNTIL ALL SINES OF GAUSSIAN LATITUDE ARE WITHIN EPS
      SPMAX=1.
      DO WHILE(SPMAX.GT.EPS)
        SPMAX=0.
        DO J=1,JH
          PKM1(J)=1.
          PK(J)=SLAT(J)
        ENDDO
        DO N=2,2*JH
          DO J=1,JH
            PKM2=PKM1(J)
            PKM1(J)=PK(J)
            PK(J)=((2*N-1)*SLAT(J)*PKM1(J)-(N-1)*PKM2)/N
          ENDDO
        ENDDO
        DO J=1,JH
          SP=PK(J)*(1.-SLAT(J)**2)/(2*JH*(PKM1(J)-SLAT(J)*PK(J)))
          SLAT(J)=SLAT(J)-SP
          SPMAX=MAX(SPMAX,ABS(SP))
        ENDDO
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE COSINES AND GAUSSIAN WEIGHTS
      DO J=1,JH
        WLAT(J)=2.*(1.-SLAT(J)**2)/(2*JH*PKM1(J))**2
      ENDDO
      RETURN
      END

C-----------------------------------------------------------------------
      SUBROUTINE MPFDEF(IPTV,MPF)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: MPFDEF         SETS DEFAULT PARAMETER FLAGS
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 92-10-31
C
C ABSTRACT: SETS FIELD IDENTIFIER DEFAULTS FOR VARIOUS PARAMETERS.
C   A FLAG OF 0 IS SCALAR, 1 IS VECTOR, 2 IS FLAG AND 3 IS BUDGET.
C   THESE IDENTIFIERS ARE USED IN INTERPOLATION.
C
C PROGRAM HISTORY LOG:
C   93-10-21  IREDELL
C   93-12-05  IREDELL   ADDED BUDGET FLAG
C
C USAGE:    CALL MPFDEF(IPTV,MPF)
C   INPUT ARGUMENTS:
C     IPTV         PARAMTER TABLE VERSION (ONLY 1 OR 2 IS RECOGNIZED)
C   OUTPUT ARGUMENTS:
C     MPF          INTEGER (255) FIELD PARAMETER IDENTIFIERS
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN
C
C$$$
      DIMENSION MPF(255)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      do i=1,255
        MPF(i)=0
      enddo
      IF(IPTV.EQ.1.OR.IPTV.EQ.2) THEN
        MPF(033)=1
        MPF(034)=1
        MPF(049)=1
        MPF(050)=1
        MPF(095)=1
        MPF(096)=1
        MPF(124)=1
        MPF(125)=1
        MPF(147)=1
        MPF(148)=1
        MPF(181)=1
        MPF(182)=1
        MPF(183)=1
        MPF(184)=1
        MPF(247)=1
        MPF(248)=1
        MPF(081)=2
        MPF(091)=2
        MPF(140)=2
        MPF(141)=2
        MPF(142)=2
        MPF(143)=2
        MPF(173)=2
        MPF(174)=2
        MPF(175)=2
        MPF(209)=2
        MPF(054)=3      ! PRECIPITABLE WATER (KG/M2)
        MPF(057)=3      ! EVAPORATION (KG/M2)
        MPF(058)=3      ! CLOUD ICE (KG/M2)
        MPF(059)=3      ! PRECIPITATION RATE (KG/M2/S)
        MPF(061)=3      ! TOTAL PRECIPITATION (KG/M2)
        MPF(062)=3      ! LARGE-SCALE PRECIPITATION (KG/M2)
        MPF(063)=3      ! CONVECTIVE PRECIPITATION (KG/M2)
        MPF(064)=3      ! WATER EQUIVALENT SNOWFALL RATE (KG/M2/S)
        MPF(065)=3      ! WATER EQUIVALENT OF SNOW DEPTH (KG/M2)
C       MPF(071)=3      ! TOTAL CLOUD COVER (PERCENT)
C       MPF(072)=3      ! CONVECTIVE CLOUD COVER (PERCENT)
        MPF(073)=3      ! LOW CLOUD COVER (PERCENT)
        MPF(074)=3      ! MIDDLE CLOUD COVER (PERCENT)
        MPF(075)=3      ! HIGH CLOUD COVER (PERCENT)
        MPF(076)=3      ! CLOUD WATER (KG/M2)
        MPF(079)=3      ! LARGE SCALE SNOW (KG/M2)
        MPF(086)=3      ! SOIL WETNESS (KG/M2)
        MPF(090)=3      ! RUNOFF (KG/M2)
        MPF(111)=3      ! NET SOLAR RADIATIVE FLUX AT SURFACE (W/M2)
        MPF(112)=3      ! NET LONGWAVE RADIATIVE FLUX AT SURFACE (W/M2)
        MPF(113)=3      ! NET SOLAR RADIATIVE FLUX AT TOP (W/M2)
        MPF(114)=3      ! NET LONGWAVE RADIATIVE FLUX AT TOP (W/M2)
        MPF(115)=3      ! NET LONGWAVE RADIATIVE FLUX (W/M2)
        MPF(116)=3      ! NET SOLAR RADIATIVE FLUX (W/M2)
        MPF(117)=3      ! TOTAL RADIATIVE FLUX (W/M2)
        MPF(121)=3      ! LATENT HEAT FLUX (W/M2)
        MPF(122)=3      ! SENSIBLE HEAT FLUX (W/M2)
        MPF(123)=3      ! BOUNDARY LAYER DISSIPATION (W/M2)
        MPF(144)=3      ! VOLUMETRIC SOIL MOISTURE CONTENT (FRACTION)
        MPF(145)=3      ! POTENTIAL EVAPORATION RATE (W/M2)
        MPF(146)=3      ! CLOUD WORKFUNCTION (J/KG)
        MPF(155)=3      ! GROUND HEAT FLUX (W/M2)
        MPF(156)=3      ! CONVECTIVE INHIBITION (W/M2)
        MPF(160)=3      ! CLEAR SKY UPWARD SOLAR FLUX (W/M2)
        MPF(161)=3      ! CLEAR SKY DOWNWARD SOLAR FLUX (W/M2)
        MPF(162)=3      ! CLEAR SKY UPWARD LONGWAVE FLUX (W/M2)
        MPF(163)=3      ! CLEAR SKY DOWNWARD LONGWAVE FLUX (W/M2)
        MPF(164)=3      ! CLOUD FORCING NET SOLAR FLUX (W/M2)
        MPF(165)=3      ! CLOUD FORCING NET LONGWAVE FLUX (W/M2)
        MPF(166)=3      ! VISIBLE BEAM DOWNWARD SOLAR FLUX (W/M2)
        MPF(167)=3      ! VISIBLE DIFFUSE DOWNWARD SOLAR FLUX (W/M2)
        MPF(168)=3      ! NEAR IR BEAM DOWNWARD SOLAR FLUX (W/M2)
        MPF(169)=3      ! NEAR IR DIFFUSE DOWNWARD SOLAR FLUX (W/M2)
        MPF(172)=3      ! MOMENTUM FLUX (N/M2)
        MPF(204)=3      ! DOWNWARD SOLAR RADIATIVE FLUX (W/M2)
        MPF(205)=3      ! DOWNWARD LONGWAVE RADIATIVE FLUX (W/M2)
        MPF(211)=3      ! UPWARD SOLAR RADIATIVE FLUX (W/M2)
        MPF(212)=3      ! UPWARD LONGWAVE RADIATIVE FLUX (W/M2)
        MPF(213)=3      ! NON-CONVECTIVE CLOUD COVER (PERCENT)
        MPF(214)=3      ! CONVECTIVE PRECIPITATION RATE (KG/M2/S)
        MPF(223)=3      ! PLANT CANOPY SURFACE WATER (KG/M2)
        MPF(228)=3      ! POTENTIAL EVAPORATION (KG/M2)
        MPF(229)=3      ! SNOW PHASE-CHANGE HEAT FLUX (W/M2)
        MPF(232)=3      ! DOWNWARD TOTAL RADIATION FLUX (W/M2)
        MPF(233)=3      ! UPWARD TOTAL RADIATION FLUX (W/M2)
        MPF(224)=3      ! BASEFLOW-GROUNDWATER RUNOFF (KG/M2)
        MPF(225)=3      ! STORM SURFACE RUNOFF (KG/M2)
        MPF(238)=3      ! SNOW COVER (PERCENT)
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      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  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
          if(mod(idum,bdsn).ne.0) then
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)
              datan(i)=(-1)**isign * refb * 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=-99999.
      rmin=99999.
      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,kpl
      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)
c
c new, added on 20000511 for ECMWF ----------- begin
c
        if(pds(3).eq.98) then
          iskp=iskp+5*8
c         write(*,*) ' ***** before pv, iskp:',iskp
          do nn=1,gdsnv
            call gbyte(igds,isign,iskp,1)
            call gbyte(igds,refa,iskp+1,7)
            call gbyte(igds,refb,iskp+8,24)
            rdum=refb
            rpv=(-1)**isign * rdum * 16**(refa-64) / 2**24
            iskp=iskp+32
            write(*,*) 'PV: ',nn,rpv
          enddo
c         write(*,*) ' ***** before pl, iskp:',iskp
          msum=0
          do nn=1,ny
            call gbyte(igds,kpl,iskp,16)
            iskp=iskp+16
            write(*,*) 'PL: ',nn,kpl
            msum=msum+kpl
          enddo
          write(*,*) 'Sum of PL: ',msum
        endif
c
c new, added on 20000511 for ECMWF ----------- end
c
      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(ingrib,isign,iskp,1)
          call gbyte(ingrib,refa,iskp+1,7)
          call gbyte(ingrib,refb,iskp+8,24)
          angrot=(-1)**isign * refb * 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(ingrib,isign,iskp,1)
          call gbyte(ingrib,refa,iskp+1,7)
          call gbyte(ingrib,refb,iskp+8,24)
          strfac=(-1)**isign * refb * 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(ingrib,isign,iskp,1)
          call gbyte(ingrib,refa,iskp+1,7)
          call gbyte(ingrib,refb,iskp+8,24)
          angrot=(-1)**isign * refb * 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(ingrib,isign,iskp,1)
          call gbyte(ingrib,refa,iskp+1,7)
          call gbyte(ingrib,refb,iskp+8,24)
          strfac=(-1)**isign * refb * 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.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
        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
      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
      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)
        dscalf=float(pds(22))
        if(isign.eq.1) pds(22)=-pds(22)
      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
      SUBROUTINE GBYTE (IN,IOUT,ISKIP,NBYTE)
      CALL GBYTES (IN,IOUT,ISKIP,NBYTE,0,1)
      RETURN
      END

      SUBROUTINE SBYTE (IOUT,IN,ISKIP,NBYTE)
      CALL SBYTES (IOUT,IN,ISKIP,NBYTE,0,1)
      RETURN
      END

      SUBROUTINE GBYTES (IN,IOUT,ISKIP,NBYTE,NSKIP,N)
C          Get bytes - unpack bits:  Extract arbitrary size values from a
C          packed bit string, right justifying each value in the unpacked
C          array.
      DIMENSION IN(*), IOUT(*)
C            IN    = packed array input
C            IO    = unpacked array output
C            ISKIP = initial number of bits to skip
C            NBYTE = number of bits to take
C            NSKIP = additional number of bits to skip on each iteration
C            N     = number of iterations
C************************************** MACHINE SPECIFIC CHANGES START HERE
C          Machine dependent information required:
C            LMWD   = Number of bits in a word on this machine
C            MASKS  = Set of word masks where the first element has only the
C                     right most bit set to 1, the second has the two, ...
C            LEFTSH = Shift left bits in word M to the by N bits
C            RGHTSH = Shift right
C            OR     = Logical OR (add) on this machine.
C            AND    = Logical AND (multiply) on this machine
C          This is for Sun UNIX Fortran, DEC Alpha, and RS6000
      PARAMETER (LMWD=32)
      DIMENSION MASKS(LMWD)
      SAVE      MASKS
      DATA      MASKS /'1'X,'3'X,'7'X,'F'X, '1F'X,'3F'X,'7F'X,'FF'X,
     +'1FF'X,'3FF'X,'7FF'X,'FFF'X, '1FFF'X,'3FFF'X,'7FFF'X,'FFFF'X,
     +'1FFFF'X,       '3FFFF'X,       '7FFFF'X,       'FFFFF'X,
     +'1FFFFF'X,      '3FFFFF'X,      '7FFFFF'X,      'FFFFFF'X,
     +'1FFFFFF'X,     '3FFFFFF'X,     '7FFFFFF'X,     'FFFFFFF'X,
     +'1FFFFFFF'X,    '3FFFFFFF'X,    '7FFFFFFF'X,    'FFFFFFFF'X/
C    +'1FFFFFFFF'X,   '3FFFFFFFF'X,   '7FFFFFFFF'X,   'FFFFFFFFF'X,
C    +'1FFFFFFFFF'X,  '3FFFFFFFFF'X,  '7FFFFFFFFF'X,  'FFFFFFFFFF'X,
C    +'1FFFFFFFFFF'X, '3FFFFFFFFFF'X, '7FFFFFFFFFF'X, 'FFFFFFFFFFF'X,
C    +'1FFFFFFFFFFF'X,'3FFFFFFFFFFF'X,'7FFFFFFFFFFF'X,'FFFFFFFFFFFF'X,
C    +'1FFFFFFFFFFFF'X,   '3FFFFFFFFFFFF'X,   '7FFFFFFFFFFFF'X,
C    +                                        'FFFFFFFFFFFFF'X,
C    +'1FFFFFFFFFFFFF'X,  '3FFFFFFFFFFFFF'X,  '7FFFFFFFFFFFFF'X,
C                                             'FFFFFFFFFFFFFF'X,
C    +'1FFFFFFFFFFFFFF'X, '3FFFFFFFFFFFFFF'X, '7FFFFFFFFFFFFFF'X,
C                                             'FFFFFFFFFFFFFFF'X,
C    +'1FFFFFFFFFFFFFFF'X,'3FFFFFFFFFFFFFFF'X,'7FFFFFFFFFFFFFFF'X,
C                                             'FFFFFFFFFFFFFFFF'X/
C          IBM PC using Microsoft Fortran uses different syntax:
C     DATA MASKS/16#1,16#3,16#7,16#F,16#1F,16#3F,16#7F,16#FF,
C    + 16#1FF,16#3FF,16#7FF,16#FFF,16#1FFF,16#3FFF,16#7FFF,16#FFFF,
C    + 16#1FFFF,16#3FFFF,16#7FFFF,16#FFFFF,16#1FFFFF,16#3FFFFF,
C    + 16#7FFFFF,16#FFFFFF,16#1FFFFFF,16#3FFFFFF,16#7FFFFFF,16#FFFFFFF,
C    + 16#1FFFFFFF,16#3FFFFFFF,16#7FFFFFFF,16#FFFFFFFF/
      INTEGER RGHTSH, OR, AND
      LEFTSH(M,N) = ISHFT(M,N)
      RGHTSH(M,N) = ISHFT(M,-N)
C     OR(M,N)  = M.OR.N
C     AND(M,N) = M.AND.N
C     OR(M,N)  = IOR(M,N)
C     AND(M,N) = IAND(M,N)
C************************************** MACHINE SPECIFIC CHANGES END HERE
C          History:  written by Robert C. Gammill, jul 1972.


C          NBYTE must be less than or equal to LMWD
      ICON = LMWD-NBYTE
      IF (ICON.LT.0) RETURN
      MASK = MASKS (NBYTE)
C          INDEX  = number of words into IN before the next "byte" appears
C          II     = number of bits the "byte" is from the left side of the word
C          ISTEP  = number of bits from the start of one "byte" to the next
C          IWORDS = number of words to skip from one "byte" to the next
C          IBITS  = number of bits to skip after skipping IWORDS
C          MOVER  = number of bits to the right, a byte must be moved to be
C                   right adjusted
      INDEX = ISKIP/LMWD
      II    = MOD (ISKIP,LMWD)
      ISTEP = NBYTE+NSKIP
      IWORDS= ISTEP/LMWD
      IBITS = MOD (ISTEP,LMWD)

      DO 6 I=1,N                                                                
      MOVER = ICON-II
      IF (MOVER) 2,3,4

C          The "byte" is split across a word break.
    2 MOVEL = -MOVER
      MOVER = LMWD-MOVEL
      NP1 = LEFTSH (IN(INDEX+1),MOVEL)
      NP2 = RGHTSH (IN(INDEX+2),MOVER)
      IOUT(I) = AND (OR (NP1,NP2) , MASK)
      GO TO 5                                                                   

C          The "byte" is already right adjusted.
    3 IOUT(I) = AND (IN (INDEX+1) , MASK)
      GO TO 5                                                                   

C          Right adjust the "byte".
    4 IOUT(I) = AND (RGHTSH (IN (INDEX+1),MOVER) , MASK)

    5 II = II+IBITS
      INDEX = INDEX+IWORDS
      IF (II .LT. LMWD) GO TO 6
      II = II-LMWD
      INDEX = INDEX+1
    6 CONTINUE                                                                  

      RETURN                                                                    
      END                                                                       

      SUBROUTINE SBYTES (IOUT,IN,ISKIP,NBYTE,NSKIP,N)
C          Store bytes - pack bits:  Put arbitrary size values into a
C          packed bit string, taking the low order bits from each value
C          in the unpacked array.
      DIMENSION IN(*), IOUT(*)
C            IOUT  = packed array output
C            IN    = unpacked array input
C            ISKIP = initial number of bits to skip
C            NBYTE = number of bits to pack
C            NSKIP = additional number of bits to skip on each iteration
C            N     = number of iterations
C************************************** MACHINE SPECIFIC CHANGES START HERE
C          Machine dependent information required:
C            LMWD   = Number of bits in a word on this machine
C            MASKS  = Set of word masks where the first element has only the
C                     right most bit set to 1, the second has the two, ...
C            LEFTSH = Shift left bits in word M to the by N bits
C            RGHTSH = Shift right
C            OR     = Logical OR (add) on this machine
C            AND    = Logical AND (multiply) on this machine
C            NOT    = Logical NOT (negation) on this machine
C          This is for Sun UNIX Fortran
      PARAMETER (LMWD=32)
      DIMENSION MASKS(LMWD)
      SAVE      MASKS
      DATA      MASKS /'1'X,'3'X,'7'X,'F'X, '1F'X,'3F'X,'7F'X,'FF'X,
     +'1FF'X,'3FF'X,'7FF'X,'FFF'X, '1FFF'X,'3FFF'X,'7FFF'X,'FFFF'X,
     +'1FFFF'X,       '3FFFF'X,       '7FFFF'X,       'FFFFF'X,
     +'1FFFFF'X,      '3FFFFF'X,      '7FFFFF'X,      'FFFFFF'X,
     +'1FFFFFF'X,     '3FFFFFF'X,     '7FFFFFF'X,     'FFFFFFF'X,
     +'1FFFFFFF'X,    '3FFFFFFF'X,    '7FFFFFFF'X,    'FFFFFFFF'X/
      INTEGER RGHTSH, OR, AND
      LEFTSH(M,N) = ISHFT(M,N)
      RGHTSH(M,N) = ISHFT(M,-N)
C     OR(M,N)  = M.OR.N
C     AND(M,N) = M.AND.N
C     OR(M,N)  = IOR(M,N)
C     AND(M,N) = IAND(M,N)
C     NOT(M)   = .NOT.M
C***********************************************************************        

C          NBYTE must be less than or equal to LMWD
      ICON = LMWD-NBYTE
      IF (ICON .LT. 0) RETURN
      MASK = MASKS(NBYTE)
C          INDEX  = number of words into IOUT the next "byte" is to be stored
C          II     = number of bits in from the left side of the word to store it
C          ISTEP  = number of bits from the start of one "byte" to the next
C          IWORDS = number of words to skip from one "byte" to the next
C          IBITS  = number of bits to skip after skipping IWORDS
C          MOVER  = number of bits to the right, a byte must be moved to be
C                   right adjusted
      INDEX = ISKIP/LMWD
      II    = MOD(ISKIP,LMWD)
      ISTEP = NBYTE+NSKIP
      IWORDS = ISTEP/LMWD
      IBITS = MOD(ISTEP,LMWD)

      DO 6 I=1,N                                                                
      J = AND (MASK,IN(I))
      MOVEL = ICON-II
      IF (MOVEL) 2,3,4

C          The "byte" is to be split across a word break
    2 MSK = MASKS (NBYTE+MOVEL)
      IOUT(INDEX+1) = OR (AND(NOT(MSK),IOUT(INDEX+1)),RGHTSH(J,-MOVEL))
      ITEMP = AND (MASKS(LMWD+MOVEL),IOUT(INDEX+2))
      IOUT(INDEX+2) = OR(ITEMP,LEFTSH(J,LMWD+MOVEL))
      GO TO 5                                                                   

C          The "byte" is to be stored right-adjusted
    3 IOUT(INDEX+1) = OR ( AND (NOT(MASK),IOUT(INDEX+1)) , J)
      GO TO 5                                                                   

C          The "byte" is to be stored in middle of word, so shift left.
    4 MSK = LEFTSH(MASK,MOVEL)
      IOUT(INDEX+1) = OR(AND(NOT(MSK),IOUT(INDEX+1)),LEFTSH(J,MOVEL))

    5 II = II+IBITS
      INDEX = INDEX+IWORDS
      IF (II .LT. LMWD) GO TO 6
      II = II-LMWD
      INDEX = INDEX+1
    6 CONTINUE

      RETURN                                                                    
      END                                                                       
c
c
"EOF"

f77 ibmspGRIB.f

#use A21825 (p-A21825        6941 rcs, 227442688 bytes,  0209.pgb.f00)
msread A21825 /DSS/A21825
cosconvert -b A21825
tar xf A21825 pgb.f0002090100

echo 'pgb.f0002090100' >! dum.in
a.out < dum.in

exit