g2index.F90

Go to the documentation of this file.

 
 
subroutine g2_create_index(lugb, lugi, idxver, filename, iret)
  implicit none
  
  integer, intent(in) :: lugb, lugi, idxver
  character*(*) :: filename
  integer, intent(out) :: iret
  
  integer (kind = 8) :: msk1, msk2
  parameter(msk1 = 32000_8, msk2 = 4000_8)
  character(len=1), pointer, dimension(:) :: cbuf
  integer :: numtot, nnum, nlen, mnum, kw
  integer :: irgi, iw, nmess
  
  interface
     subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
          nlen, nnum, nmess, iret)
       integer, intent(in) :: lugb
       integer (kind = 8), intent(in) :: msk1, msk2
       integer, intent(in) :: mnum, idxver
       character(len = 1), pointer, dimension(:) :: cbuf
       integer, intent(out) :: nlen, nnum, nmess, iret
     end subroutine getg2i2r
  end interface
 
  ! Assume success.
  iret = 0
  numtot = 0
  nlen = 0
 
  ! Generate index records for all messages in file, or until memory
  ! runs out.
  mnum = 0
  call getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
       nlen, nnum, nmess, irgi)
  if (irgi .gt. 1 .or. nnum .eq. 0 .or. nlen .eq. 0) then
     iret = 92
     return
  endif
  numtot = numtot + nnum
  mnum = mnum + nmess
 
  ! Write headers.
  call g2_write_index_headers(lugi, nlen, numtot, idxver, filename)
  iw = 162
 
  ! Write the index data we have so far.
  call bawrite(lugi, iw, nlen, kw, cbuf)
  iw = iw + nlen
 
  ! Extend index file if index buffer length too large to hold in memory.
  if (irgi .eq. 1) then
     do while (irgi .eq. 1 .and. nnum .gt. 0)
        if (associated(cbuf)) then
           deallocate(cbuf)
           nullify(cbuf)
        endif
        call getg2i2r(11, msk1, msk2, mnum, idxver, cbuf, &
             nlen, nnum, nmess, irgi)
        if (irgi .le. 1 .and. nnum .gt. 0) then
           numtot = numtot + nnum
           mnum = mnum + nmess
           call bawrite(lugi, iw, nlen, kw, cbuf)
           iw = iw + nlen
        endif
     enddo
     ! Go back and overwrite headers with new info.
     call g2_write_index_headers(lugi, iw, numtot, idxver, filename)
  endif
  deallocate(cbuf)  
 
end subroutine g2_create_index
 
subroutine g2_write_index_headers(lugi, nlen, nnum, idxver, filename)
  implicit none
 
  integer, intent(in) :: lugi, nlen, nnum, idxver
  character, intent(in) :: filename*(*)
  
  character cd8*8, ct10*10, hostname*15
#ifdef __GFORTRAN__
  integer istat
#else
  character hostnam*15
  integer hostnm
#endif
  character chead(2)*81
  integer :: kw
 
  !  fill first 81-byte header
  call date_and_time(cd8, ct10)
  chead(1) = '!GFHDR!'
  chead(1)(9:10) = ' 1'
  chead(1)(12:14) = '  1'
  write(chead(1)(16:20),'(i5)') 162
  chead(1)(22:31) = cd8(1:4) // '-' // cd8(5:6) // '-' // cd8(7:8)
  chead(1)(33:40) = ct10(1:2) // ':' // ct10(3:4) // ':' // ct10(5:6)
  chead(1)(42:47) = 'GB2IX1'
  chead(1)(49:54) = '      '
#ifdef __GFORTRAN__
  istat = hostnm(hostname)
  if (istat .eq. 0) then
     chead(1)(56:70) = '0000'
  else
     chead(1)(56:70) = '0001'
  endif
#else
  chead(1)(56:70) = hostnam(hostname)
#endif
  chead(1)(72:80) = 'grb2index'
  chead(1)(81:81) = char(10)
 
  !  fill second 81-byte header
  if (idxver .eq. 1) then
     chead(2) = 'IX1FORM:'
  else
     chead(2) = 'IX2FORM:'
  endif
  write(chead(2)(9:38),'(3i10)') 162, nlen, nnum
  chead(2)(41:80) = filename
  chead(2)(81:81) = char(10)
 
  !  write headers at beginning of index file
  call bawrite(lugi, 0, 162, kw, chead)
 
  return
end subroutine g2_write_index_headers
 
subroutine getidx(lugb, lugi, cindex, nlen, nnum, iret)
  implicit none
 
  integer, intent(in) :: lugb, lugi
  character(len = 1), pointer, dimension(:) :: cindex
  integer, intent(out) :: nlen, nnum, iret
  integer :: idxver
 
  interface
     subroutine getidx2(lugb, lugi, idxver, cindex, nlen, nnum, iret)
       integer, intent(in) :: lugb, lugi
       integer, intent(inout) :: idxver
       character(len = 1), pointer, dimension(:) :: cindex
       integer, intent(out) :: nlen, nnum, iret
     end subroutine getidx2
  end interface
 
  ! When getidx() is called, always use index version 1. Call
  ! getidx2() for a chance to set the index version.
  idxver = 1
  call getidx2(lugb, lugi, idxver, cindex, nlen, nnum, iret)
 
  ! If index version 2 is being used, return error.
  if (iret .eq. 0 .and. idxver .eq. 2) iret = 97
 
end subroutine getidx
 
subroutine getidx2(lugb, lugi, idxver, cindex, nlen, nnum, iret)
  implicit none
 
  integer, intent(in) :: lugb, lugi
  integer, intent(inout) :: idxver
  character(len = 1), pointer, dimension(:) :: cindex
  integer, intent(out) :: nlen, nnum, iret
  
  integer, parameter :: maxidx = 10000
  integer (kind = 8), parameter :: msk1 = 32000_8, msk2 = 4000_8
  integer :: lux
  integer :: irgi, mskp, nmess, i
 
  type gindex
     integer :: nlen
     integer :: nnum
     integer :: idxver
     character(len = 1), pointer, dimension(:) :: cbuf
  end type gindex
 
  type(gindex), save :: idxlist(10000)
 
  data lux/0/
 
  !  declare interfaces (required for cbuf pointer)
  interface
     subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
       integer, intent(in) :: lugi
       character(len=1), pointer, dimension(:) :: cbuf
       integer, intent(out) :: idxver, nlen, nnum, iret
     end subroutine getg2i2
     subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
          nlen, nnum, nmess, iret)
       integer, intent(in) :: lugb
       integer (kind = 8), intent(in) :: msk1, msk2
       integer, intent(in) :: mnum, idxver
       character(len = 1), pointer, dimension(:) :: cbuf
       integer, intent(out) :: nlen, nnum, nmess, iret
     end subroutine getg2i2r
  end interface
 
  ! Free all associated memory and exit.
  if (lugb .eq. 0) then
     !print *, 'getidx: Freeing all memory'
     do i = 1, 10000
        if (associated(idxlist(i)%cbuf)) then
           !print *, 'deallocating ', loc(idxlist(i)%cbuf)
           deallocate(idxlist(i)%cbuf)
           nullify(idxlist(i)%cbuf)
        endif
     end do
     iret = 0
     return
  endif
 
  ! Determine whether index buffer needs to be initialized.
  lux = 0
  iret = 0
  if (lugb .le. 0 .or. lugb .gt. 9999) then
     print *, ' file unit number out of range'
     print *, ' use unit numbers in range: 0 - 9999 '
     iret = 90
     return
  endif
 
  ! Force regeneration of index from GRIB2 file.
  if (lugi .eq. lugb) then      
     if (associated(idxlist(lugb)%cbuf))  &
          deallocate(idxlist(lugb)%cbuf)
     !print *, 'Force regeneration'
     nullify(idxlist(lugb)%cbuf)
     idxlist(lugb)%nlen = 0
     idxlist(lugb)%nnum = 0
     lux = 0
  endif
 
  ! Force re-read of index from indexfile.
  if (lugi .lt. 0) then      
     ! associated with unit abs(lugi)
     if (associated(idxlist(lugb)%cbuf))  &
          deallocate(idxlist(lugb)%cbuf)
     !print *, 'Force re-read'
     nullify(idxlist(lugb)%cbuf)
     idxlist(lugb)%nlen = 0
     idxlist(lugb)%nnum = 0
     lux = abs(lugi)
  endif
 
  !  Check if index already exists in memory.
  if (associated(idxlist(lugb)%cbuf)) then
     !print *, 'Index exists in memory!'
     cindex => idxlist(lugb)%cbuf
     nlen = idxlist(lugb)%nlen
     nnum = idxlist(lugb)%nnum
     idxver = idxlist(lugb)%idxver
     return
  endif
 
  ! Either read index from index file, or generate it from the GRIB2
  ! file. This is an index for all messages in the file, each message
  ! gets an index record, all stuffed into idxlist(lugbb)%cbuf.
  irgi = 0
  if (lux .gt. 0) then
     call getg2i2(lux, idxlist(lugb)%cbuf, idxver, nlen, nnum, irgi)
  elseif (lux .le. 0) then
     mskp = 0
     call getg2i2r(lugb, msk1, msk2, mskp, idxver, idxlist(lugb)%cbuf, &
          nlen, nnum, nmess, irgi)
  endif
 
  ! Handle errors.
  if (irgi .ne. 0) then
     nlen = 0
     nnum = 0
     print *, ' error reading index file '
     iret = 96
     return
  endif
 
  ! Fill these values.
  cindex => idxlist(lugb)%cbuf
  idxlist(lugb)%nlen = nlen
  idxlist(lugb)%nnum = nnum
  idxlist(lugb)%idxver = idxver
 
end subroutine getidx2
 
subroutine getg2i(lugi, cbuf, nlen, nnum, iret)
  implicit none
  
  integer, intent(in) :: lugi
  character(len=1), pointer, dimension(:) :: cbuf
  integer, intent(out) :: nlen, nnum, iret
  integer :: idxver
 
  interface
     subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
       integer, intent(in) :: lugi
       character(len=1), pointer, dimension(:) :: cbuf
       integer, intent(out) :: idxver, nlen, nnum, iret
     end subroutine getg2i2
  end interface
 
  ! Call the version of this subroutine that handles version 1 and
  ! version 2. If getg2i() is called on a version 2 index, return
  ! error.
  call getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
  if (idxver .eq. 2) iret = 5
 
end subroutine getg2i
  
subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
  implicit none
  
  integer, intent(in) :: lugi
  character(len=1), pointer, dimension(:) :: cbuf
  integer, intent(out) :: idxver, nlen, nnum, iret
  
  character chead*162
  integer :: ios, istat, lbuf, lhead, nskp
 
  nullify(cbuf)
  nlen = 0
  nnum = 0
  iret = 4
  call baread(lugi, 0, 162, lhead, chead)
  if (lhead .eq. 162 .and. chead(42:47) .eq. 'GB2IX1') then
     read(chead(82:162), '(2x, i1, 5x, 3i10, 2x, a40)', iostat = ios) idxver, nskp, nlen, nnum
     if (ios .eq. 0) then
        allocate(cbuf(nlen), stat = istat)    ! Allocate space for cbuf.
        if (istat .ne. 0) then
           iret = 2
           return
        endif
        iret = 0
        call baread(lugi, nskp, nlen, lbuf, cbuf)
        if (lbuf .ne. nlen) iret = 3
     endif
  endif
end subroutine getg2i2
 
subroutine getg2ir(lugb, msk1, msk2, mnum, cbuf, nlen, nnum, nmess, iret)
  use re_alloc              ! needed for subroutine realloc
  implicit none
 
  integer, intent(in) :: lugb
  integer, intent(in) :: msk1, msk2
  integer, intent(in) :: mnum
  character(len = 1), pointer, dimension(:) :: cbuf
  integer, intent(out) :: nlen, nnum, nmess, iret
 
  integer (kind = 8) :: msk1_8, msk2_8
  
  interface
     subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
          nlen, nnum, nmess, iret)
       integer, intent(in) :: lugb
       integer (kind = 8), intent(in) :: msk1, msk2
       integer, intent(in) :: mnum, idxver
       character(len = 1), pointer, dimension(:) :: cbuf
       integer, intent(out) :: nlen, nnum, nmess, iret
     end subroutine getg2i2r
  end interface
 
  msk1_8 = msk1
  msk2_8 = msk2
  call getg2i2r(lugb, msk1_8, msk2_8, mnum, 1, cbuf, nlen, nnum, nmess, iret)
end subroutine getg2ir
     
subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, nlen, nnum, nmess, iret)
  use re_alloc              ! needed for subroutine realloc
  implicit none
 
  integer, intent(in) :: lugb
  integer (kind = 8), intent(in) :: msk1, msk2
  integer, intent(in) :: mnum, idxver
  character(len = 1), pointer, dimension(:) :: cbuf
  integer, intent(out) :: nlen, nnum, nmess, iret
  
  character(len = 1), pointer, dimension(:) :: cbuftmp
  integer :: nbytes, newsize, next, numfld, m, mbuf
  integer (kind = 8) :: iseek, lskip, lgrib
  integer :: istat, init, iret1, lgrib4
  parameter(init = 50000, next = 10000)
 
  interface      ! required for cbuf pointer
     subroutine ix2gb2(lugb, lskip8, idxver, lgrib, cbuf, numfld, mlen, iret)
       integer :: lugb
       integer (kind = 8) :: lskip8
       integer :: idxver, lgrib
       character(len = 1), pointer, dimension(:) :: cbuf
       integer :: numfld, mlen, iret
     end subroutine ix2gb2
  end interface
 
  ! Initialize.
  iret = 0
  nullify(cbuf)
  mbuf = init
  allocate(cbuf(mbuf), stat = istat)    ! allocate initial space for cbuf.
  if (istat .ne. 0) then
     iret = 2
     return
  endif
 
  ! Search for first grib message.
  iseek = 0_8
  call skgb8(lugb, iseek, msk1, lskip, lgrib)
  do m = 1, mnum
     if (lgrib .gt. 0) then
        iseek = lskip + lgrib
        call skgb8(lugb, iseek, msk2, lskip, lgrib)
     endif
  enddo
 
  ! Get index records for every grib message found.
  nlen = 0
  nnum = 0
  nmess = mnum
  do while (iret .eq. 0 .and. lgrib .gt. 0)
     lgrib4 = int(lgrib, kind(4))
     call ix2gb2(lugb, lskip, idxver, lgrib4, cbuftmp, numfld, nbytes, iret1)
     if (iret1 .ne. 0) print *, ' SAGT ', numfld, nbytes, iret1
     if((nbytes + nlen) .gt. mbuf) then             ! Allocate more space, if necessary.
        newsize = max(mbuf + next, mbuf + nbytes)
        call realloc(cbuf, nlen, newsize, istat)
        if (istat .ne. 0) then
           iret = 1
           return
        endif
        mbuf = newsize
     endif
 
     ! If index records were returned in cbuftmp from ixgb2, 
     ! copy cbuftmp into cbuf, then deallocate cbuftmp when done.
     if (associated(cbuftmp)) then
        cbuf(nlen + 1 : nlen + nbytes) = cbuftmp(1 : nbytes)
        deallocate(cbuftmp, stat = istat)
        if (istat .ne. 0) then
           print *, ' deallocating cbuftmp ... ', istat
           iret = 3
           return
        endif
        nullify(cbuftmp)
        nnum = nnum + numfld
        nlen = nlen + nbytes
        nmess = nmess + 1
     endif
 
     ! Look for next grib message.
     iseek = lskip + lgrib
     call skgb8(lugb, iseek, msk2, lskip, lgrib)
  enddo
end subroutine getg2i2r
 
subroutine getgb2s(cbuf, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
     jgdt, k, gfld, lpos, iret)
  use grib_mod
  implicit none
 
  character(len = 1), intent(in) :: cbuf(nlen)
  integer, intent(in) :: nlen, nnum, j, jdisc
  integer, dimension(:) :: jids(*)
  integer, intent(in) :: jpdtn
  integer, dimension(:) :: jpdt(*)
  integer, intent(in) :: jgdtn
  integer, dimension(:) :: jgdt(*)
  integer, intent(out) :: k
  type(gribfield), intent(out) :: gfld
  integer, intent(out) :: lpos, iret
 
  interface
     subroutine getgb2s2(cbuf, idxver, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
          jgdt, k, gfld, lpos, iret)
       import gribfield
       character(len = 1), intent(in) :: cbuf(nlen)
       integer, intent(in) :: idxver, nlen, nnum, j, jdisc
       integer, dimension(:) :: jids(*)
       integer, intent(in) :: jpdtn
       integer, dimension(:) :: jpdt(*)
       integer, intent(in) :: jgdtn
       integer, dimension(:) :: jgdt(*)
       integer, intent(out) :: k
       type(gribfield), intent(out) :: gfld
       integer, intent(out) :: lpos, iret
     end subroutine getgb2s2
  end interface
 
  ! When getgb2s() is called, always use index version 1. Call
  ! getgb2s2() to handle version 1 or 2.
  call getgb2s2(cbuf, 1, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
       jgdt, k, gfld, lpos, iret)
 
end subroutine getgb2s
 
subroutine getgb2s2(cbuf, idxver, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
     jgdt, k, gfld, lpos, iret)
  use grib_mod
  implicit none
 
  character(len = 1), intent(in) :: cbuf(nlen)
  integer, intent(in) :: idxver, nlen, nnum, j, jdisc
  integer, dimension(:) :: jids(*)
  integer, intent(in) :: jpdtn
  integer, dimension(:) :: jpdt(*)
  integer, intent(in) :: jgdtn
  integer, dimension(:) :: jgdt(*)
  integer, intent(out) :: k
  type(gribfield), intent(out) :: gfld
  integer, intent(out) :: lpos, iret
 
  integer :: kgds(5)
  logical :: match1, match3, match4
  integer :: i, icnd, inlen, iof, ipos, jpos, lsec1, lsec3, lsec4, lsec5, numgdt, numpdt, inc
 
  interface
     subroutine gf_unpack1(cgrib, lcgrib, iofst, ids, idslen, ierr)
       character(len = 1), intent(in) :: cgrib(lcgrib)
       integer, intent(in) :: lcgrib
       integer, intent(inout) :: iofst
       integer, pointer, dimension(:) :: ids
       integer, intent(out) :: ierr, idslen
     end subroutine gf_unpack1
     subroutine gf_unpack3(cgrib, lcgrib, iofst, igds, igdstmpl, &
          mapgridlen, ideflist, idefnum, ierr)
       character(len = 1), intent(in) :: cgrib(lcgrib)
       integer, intent(in) :: lcgrib
       integer, intent(inout) :: iofst
       integer, pointer, dimension(:) :: igdstmpl, ideflist
       integer, intent(out) :: igds(5)
       integer, intent(out) :: ierr, idefnum
     end subroutine gf_unpack3
     subroutine gf_unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
          mappdslen, coordlist, numcoord, ierr)
       character(len = 1), intent(in) :: cgrib(lcgrib)
       integer, intent(in) :: lcgrib
       integer, intent(inout) :: iofst
       real, pointer, dimension(:) :: coordlist
       integer, pointer, dimension(:) :: ipdstmpl
       integer, intent(out) :: ipdsnum
       integer, intent(out) :: ierr, numcoord
     end subroutine gf_unpack4
     subroutine gf_unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum, &
          idrstmpl, mapdrslen, ierr)
       character(len = 1), intent(in) :: cgrib(lcgrib)
       integer, intent(in) :: lcgrib
       integer, intent(inout) :: iofst
       integer, intent(out) :: ndpts, idrsnum
       integer, pointer, dimension(:) :: idrstmpl
       integer, intent(out) :: ierr
     end subroutine gf_unpack5
  end interface
 
  ! Initialize.
  k = 0
  lpos = 0
  iret = 1
  ipos = 0
  nullify(gfld%idsect, gfld%local)
  nullify(gfld%list_opt, gfld%igdtmpl, gfld%ipdtmpl)
  nullify(gfld%coord_list, gfld%idrtmpl, gfld%bmap, gfld%fld)
  if (idxver .eq. 1) then
     inc = 0
  else
     ! Add the extra 4 bytes in the version 2 index record, starting
     ! at byte 9.
     inc = 4
  endif
 
  ! Search for request.
  do while(iret .ne. 0 .and. k .lt. nnum)
     k = k + 1
     ! Get length of current index record.
     call g2_gbytec(cbuf, inlen, ipos * 8, 4 * 8)    
     if (k .le. j) then           ! skip this index
        ipos = ipos + inlen
        cycle
     endif
 
     ! Check if grib2 discipline is a match.
     call g2_gbytec(cbuf, gfld%discipline, (ipos + inc + 41) * 8, 1 * 8)
     if (jdisc .ne. -1 .and. jdisc .ne. gfld%discipline) then
        ipos = ipos + inlen
        cycle
     endif
 
     ! Check if identification section is a match.
     match1 = .false.
     ! Get length of ids.
     call g2_gbytec(cbuf, lsec1, (ipos + inc + 44) * 8, 4 * 8)  
     iof = 0
     call gf_unpack1(cbuf(ipos + inc + 45), lsec1, iof, gfld%idsect, gfld%idsectlen, icnd)
     if (icnd .eq. 0) then
        match1 = .true.
        do i = 1, gfld%idsectlen
           if (jids(i) .ne. -9999 .and. jids(i) .ne. gfld%idsect(i)) then
              match1 = .false.
              exit
           endif
        enddo
     endif
     if (.not. match1) then
        deallocate(gfld%idsect)
        ipos = ipos + inlen
        cycle
     endif
 
     ! Check if grid definition template is a match.
     jpos = ipos + 44 + inc + lsec1
     match3 = .false.
     call g2_gbytec(cbuf, lsec3, jpos * 8, 4 * 8)  ! get length of gds 
     if (jgdtn .eq. -1) then
        match3 = .true.
     else
        call g2_gbytec(cbuf, numgdt, (jpos + 12) * 8, 2 * 8)  ! get gdt template no.
        if (jgdtn .eq. numgdt) then
           iof = 0
           call gf_unpack3(cbuf(jpos + 1), lsec3, iof, kgds, gfld%igdtmpl, &
                gfld%igdtlen, gfld%list_opt, gfld%num_opt, icnd)
           if (icnd .eq. 0) then
              match3 = .true.
              do i = 1, gfld%igdtlen
                 if (jgdt(i) .ne. -9999 .and. jgdt(i).ne.gfld%igdtmpl(i)) then
                    match3 = .false.
                    exit
                 endif
              enddo
           endif
        endif
     endif
     if (.not. match3) then
        if (associated(gfld%idsect)) deallocate(gfld%idsect)
        if (associated(gfld%igdtmpl)) deallocate(gfld%igdtmpl)
        if (associated(gfld%list_opt)) deallocate(gfld%list_opt)
        ipos = ipos + inlen
        cycle
     else
        gfld%griddef = kgds(1)
        gfld%ngrdpts = kgds(2)
        gfld%numoct_opt = kgds(3)
        gfld%interp_opt = kgds(4)
        gfld%igdtnum = kgds(5)
     endif
 
     ! Check if product definition template is a match.
     jpos = jpos + lsec3
     match4 = .false.
 
     ! Get length of pds.      
     call g2_gbytec(cbuf, lsec4, jpos * 8, 4 * 8)  
     if (jpdtn .eq. -1) then
        match4 = .true.
     else
        ! Get pdt template no.        
        call g2_gbytec(cbuf, numpdt, (jpos + 7) * 8, 2 * 8)  
        if (jpdtn .eq. numpdt) then
           iof = 0
           call gf_unpack4(cbuf(jpos + 1), lsec4, iof, gfld%ipdtnum, &
                gfld%ipdtmpl, gfld%ipdtlen, gfld%coord_list, gfld%num_coord, icnd)
           if (icnd .eq. 0) then
              match4 = .true.
              do i = 1, gfld%ipdtlen
                 if (jpdt(i) .ne. -9999 .and. jpdt(i) .ne. gfld%ipdtmpl(i)) then
                    match4 = .false.
                    exit
                 endif
              enddo
           endif
        endif
     endif
     if (.not. match4) then
        if (associated(gfld%idsect)) deallocate(gfld%idsect)
        if (associated(gfld%ipdtmpl)) deallocate(gfld%ipdtmpl)
        if (associated(gfld%coord_list)) deallocate(gfld%coord_list)
     endif
 
     ! If request is found set values for derived type gfld and return.
     if (match1 .and. match3 .and. match4) then
        lpos = ipos + 1
        call g2_gbytec(cbuf, gfld%version, (ipos + inc + 40) * 8, 1 * 8)
        call g2_gbytec(cbuf, gfld%ifldnum, (ipos + inc + 42) * 8, 2 * 8)
        gfld%unpacked = .false.
        jpos = ipos + 44 + inc + lsec1
        if (jgdtn .eq. -1) then     ! unpack gds, if not done before
           iof = 0
           call gf_unpack3(cbuf(jpos + 1), lsec3, iof, kgds, gfld%igdtmpl, &
                gfld%igdtlen, gfld%list_opt, gfld%num_opt, icnd)
           gfld%griddef = kgds(1)
           gfld%ngrdpts = kgds(2)
           gfld%numoct_opt = kgds(3)
           gfld%interp_opt = kgds(4)
           gfld%igdtnum = kgds(5)
        endif
        jpos = jpos + lsec3
        if (jpdtn .eq. -1) then     ! unpack pds, if not done before
           iof = 0
           call gf_unpack4(cbuf(jpos + 1), lsec4, iof, gfld%ipdtnum, &
                gfld%ipdtmpl, gfld%ipdtlen, gfld%coord_list, gfld%num_coord, icnd)
        endif
        jpos = jpos + lsec4
        call g2_gbytec(cbuf, lsec5, jpos * 8, 4 * 8)  ! get length of drs 
        iof = 0
        call gf_unpack5(cbuf(jpos + 1), lsec5, iof, gfld%ndpts, &
             gfld%idrtnum, gfld%idrtmpl, gfld%idrtlen, icnd)
        jpos = jpos + lsec5
        call g2_gbytec(cbuf, gfld%ibmap, (jpos + 5) * 8, 1 * 8)  ! get ibmap
        iret = 0
     else      ! pdt did not match
        ipos = ipos + inlen
     endif
  enddo
end subroutine getgb2s2
 
subroutine ixgb2(lugb, lskip, lgrib, cbuf, numfld, mlen, iret)
  use re_alloc              ! needed for subroutine realloc
  implicit none
 
  integer :: lugb, lskip, lgrib
  character(len = 1), pointer, dimension(:) :: cbuf
  integer :: numfld, mlen, iret
  integer (kind = 8) :: lskip8
 
  interface
     subroutine ix2gb2(lugb, lskip8, idxver, lgrib, cbuf, numfld, mlen, iret)
       integer :: lugb
       integer (kind = 8) :: lskip8
       integer :: idxver, lgrib
       character(len = 1), pointer, dimension(:) :: cbuf
       integer :: numfld, mlen, iret
     end subroutine ix2gb2
  end interface
 
  ! Always use index version 1 from this subroutine.
  lskip8 = lskip
  call ix2gb2(lugb, lskip8, 1, lgrib, cbuf, numfld, mlen, iret)
end subroutine ixgb2
 
subroutine ix2gb2(lugb, lskip8, idxver, lgrib, cbuf, numfld, mlen, iret)
  use re_alloc              ! needed for subroutine realloc
  implicit none
 
  integer :: lugb
  integer (kind = 8) :: lskip8
  integer :: idxver, lgrib
  character(len = 1), pointer, dimension(:) :: cbuf
  integer :: numfld, mlen, iret
  
  character cver, cdisc
  character(len = 4) :: ctemp
  integer loclus, locgds, locbms
  integer :: indbmp, numsec, next, newsize, g2_mova2i, mbuf, lindex
  integer :: linmax, ixskp
  integer :: mxspd, mxskp, mxsgd, mxsdr, mxsbm, mxlus
  integer :: mxlen, mxds, mxfld, mxbms
  integer :: init, ixlus, lskip
  integer :: ixsgd, ilndrs, ilnpds, istat, ixds
  integer (kind = 8) :: ibread8, lbread8, ibskip8, lengds8
  integer (kind = 8) :: ilnpds8, ilndrs8
  integer :: ixspd, ixfld, ixids, ixlen, ixsbm, ixsdr
  integer :: lensec, lensec1
  parameter(linmax = 5000, init = 50000, next = 10000)
  parameter(ixskp = 4, ixlus = 8, ixsgd = 12, ixspd = 16, ixsdr = 20, ixsbm = 24, &
       ixds = 28, ixlen = 36, ixfld = 42, ixids = 44)
  parameter(mxskp = 4, mxlus = 4, mxsgd = 4, mxspd = 4, mxsdr = 4, mxsbm = 4, &
       mxds = 4, mxlen = 4, mxfld = 2, mxbms = 6)
  character cbread(linmax), cindex(linmax)
  character cids(linmax), cgds(linmax)
  integer :: INT1_BITS, INT2_BITS, INT4_BITS, INT8_BITS
  parameter(int1_bits = 8, int2_bits = 16, int4_bits = 32, int8_bits = 64)
  integer :: mypos, inc
 
  if (idxver .eq. 1) then
     inc = 0
  else
     ! Add the extra 4 bytes in the version 2 index record, starting
     ! at byte 9.
     inc = 4
  endif
 
  loclus = 0
  iret = 0
  mlen = 0
  numfld = 0
  nullify(cbuf)
  mbuf = init
  allocate(cbuf(mbuf), stat = istat)    ! allocate initial space for cbuf
  if (istat .ne. 0) then
     iret = 1
     return
  endif
 
  ! Read sections 0 and 1 for GRIB version number and discipline.
  ibread8 = min(lgrib, linmax)
  call bareadl(lugb, lskip8, ibread8, lbread8, cbread)
  if (lbread8 .ne. ibread8) then
     iret = 2
     return
  endif
  if(cbread(8) .ne. char(2)) then          !  not grib edition 2
     iret = 3
     return
  endif
  cver = cbread(8)
  cdisc = cbread(7)
  call g2_gbytec(cbread, lensec1, 16 * 8, int4_bits)
  lensec1 = min(lensec1, int(ibread8, kind(lensec1)))
  cids(1:lensec1) = cbread(17:16 + lensec1)
  ibskip8 = lskip8 + 16_8 + int(lensec1, kind(8))
 
  ! Loop through remaining sections creating an index for each field.
  ibread8 = max(5, mxbms)
  do
     call bareadl(lugb, ibskip8, ibread8, lbread8, cbread)
     ctemp = cbread(1)//cbread(2)//cbread(3)//cbread(4)
     if (ctemp .eq. '7777') return        ! end of message found
     if (lbread8 .ne. ibread8) then
        iret = 2
        return
     endif
     call g2_gbytec(cbread, lensec, 0, int4_bits)
     call g2_gbytec(cbread, numsec, int4_bits, int1_bits)
 
     if (numsec .eq. 2) then                 ! save local use location
        loclus = int(ibskip8 - lskip8, kind(4))
     elseif (numsec .eq. 3) then                 ! save gds info
        lengds8 = lensec
        cgds = char(0)
        call bareadl(lugb, ibskip8, lengds8, lbread8, cgds)
        if (lbread8 .ne. lengds8) then
           iret = 2
           return
        endif
        locgds = int(ibskip8 - lskip8, kind(4))
     elseif (numsec .eq. 4) then                 ! found pds
        cindex = char(0)
        mypos = int4_bits
        if (idxver .eq. 1) then
           lskip = int(lskip8, kind(4))
           call g2_sbytec(cindex, lskip, mypos, int4_bits)    ! bytes to skip
           mypos = mypos + int4_bits
        else
           call g2_sbytec8(cindex, lskip8, mypos, int8_bits)    ! bytes to skip
           mypos = mypos + int8_bits
        endif
        call g2_sbytec(cindex, loclus, mypos, int4_bits)   ! location of local use
        mypos = mypos + int4_bits
        call g2_sbytec(cindex, locgds, mypos, int4_bits)   ! location of gds
        mypos = mypos + int4_bits
        call g2_sbytec(cindex, int(ibskip8 - lskip8, kind(4)), mypos, int4_bits)  ! location of pds
        mypos = mypos + int4_bits * 4 ! skip ahead in cbuf
        call g2_sbytec(cindex, lgrib, mypos, int8_bits)    ! len of grib2
        mypos = mypos + int8_bits
        cindex((mypos / 8) + 1) = cver
        mypos = mypos + int1_bits
        cindex((mypos / 8) + 1) = cdisc
        mypos = mypos + int1_bits
        call g2_sbytec(cindex, numfld + 1, mypos, int2_bits)   ! field num
        mypos = mypos + int2_bits
        cindex(ixids + 1 + inc:ixids + lensec1 + inc) = cids(1:lensec1)
        lindex = ixids + lensec1 + inc
        cindex(lindex + 1:lindex + lengds8) = cgds(1:lengds8)
        lindex = lindex + int(lengds8, kind(lindex))
        ilnpds = lensec
        ilnpds8 = ilnpds        
        call bareadl(lugb, ibskip8, ilnpds8, lbread8, cindex(lindex + 1))
        if (lbread8 .ne. ilnpds8) then
           iret = 2
           return
        endif
        lindex = lindex + ilnpds
     elseif (numsec .eq. 5) then                 ! found drs
        mypos = (ixsdr + inc) * int1_bits
        call g2_sbytec(cindex, int(ibskip8 - lskip8, kind(4)), mypos, int4_bits)  ! location of drs
        ilndrs = lensec
        ilndrs8 = ilndrs
        call bareadl(lugb, ibskip8, ilndrs8, lbread8, cindex(lindex + 1))
        if (lbread8 .ne. ilndrs8) then
           iret = 2
           return
        endif
        lindex = lindex + ilndrs
     elseif (numsec .eq. 6) then                 ! found bms
        indbmp = g2_mova2i(cbread(6))
        mypos = (ixsbm + inc) * int1_bits           
        if (indbmp.lt.254) then
           locbms = int(ibskip8 - lskip8, kind(4))
           call g2_sbytec(cindex, locbms, mypos, int4_bits)  ! loc. of bms
        elseif (indbmp.eq.254) then
           call g2_sbytec(cindex, locbms, mypos, int4_bits)  ! loc. of bms
        elseif (indbmp.eq.255) then
           call g2_sbytec(cindex, int(ibskip8 - lskip8, kind(4)), mypos, int4_bits)  ! loc. of bms
        endif
        cindex(lindex + 1:lindex + mxbms) = cbread(1:mxbms)
        lindex = lindex + mxbms
        call g2_sbytec(cindex, lindex, 0, int4_bits)    ! num bytes in index record
     elseif (numsec .eq. 7) then                 ! found data section
        mypos = (ixds + inc) * int1_bits           
        call g2_sbytec(cindex, int(ibskip8 - lskip8, kind(4)), mypos, int4_bits)   ! loc. of data sec.
        numfld = numfld + 1
        if ((lindex + mlen) .gt. mbuf) then ! allocate more space if necessary
           newsize = max(mbuf + next, mbuf + lindex)
           call realloc(cbuf, mlen, newsize, istat)
           if (istat .ne. 0) then
              numfld = numfld-1
              iret = 4
              return
           endif
           mbuf = newsize
        endif
        cbuf(mlen + 1:mlen + lindex) = cindex(1:lindex)
        mlen = mlen + lindex
     else                           ! unrecognized section
        iret = 5
        return
     endif
     ibskip8 = ibskip8 + lensec
  enddo
end subroutine ix2gb2
 
subroutine gf_finalize(iret)
  implicit none
 
  integer, intent(out) :: iret  
  character(len = 1), pointer, dimension(:) :: cindex
  integer :: nlen, nnum
 
  ! Declare interfaces (required for cbuf pointer).
  interface
     subroutine getidx(lugb, lugi, cbuf, nlen, nnum, irgi)
       character(len = 1), pointer, dimension(:) :: cbuf
       integer, intent(in) :: lugb, lugi
       integer, intent(out) :: nlen, nnum, irgi
     end subroutine getidx
  end interface
 
  ! Call getidx with 0 for the first parameter, ensuring that the
  ! internal memory is freed.
  call getidx(0, 0, cindex, nlen, nnum, iret)
 
end subroutine gf_finalize