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)
  use g2logging
  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
 
#ifdef LOGGING
  ! Log results for debugging.
  write(g2_log_msg, '(a, i2, a, i2, a, i1)') 'getidx2: lugb ', lugb, ' lugi ', lugi, &
       ' idxver ', idxver
  call g2_log(1)
#endif
 
  ! 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)
  use g2logging
  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
 
#ifdef LOGGING
  write(g2_log_msg, *) 'getg2i2: lugi ', lugi
  call g2_log(1)
#endif
 
  nullify(cbuf)
  nlen = 0
  nnum = 0
  iret = 4
  call baread(lugi, 0, 162, lhead, chead)
#ifdef LOGGING
     write(g2_log_msg, *) 'lhead', lhead
     call g2_log(2)
#endif
  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
#ifdef LOGGING
     write(g2_log_msg, *) 'ios', ios, 'idxver', idxver, 'nskp', nskp, 'nlen', nlen, 'nnum', nnum
     call g2_log(2)
#endif 
    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 g2logging
  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, lgrib8, cbuf, numfld, mlen, iret)
       integer, intent(in) :: lugb
       integer (kind = 8), intent(in) :: lskip8
       integer, intent(in)  :: idxver
       integer (kind = 8), intent(in) :: lgrib8
       character(len = 1), pointer, intent(inout), dimension(:) :: cbuf
       integer, intent(out)  :: numfld, mlen, iret
     end subroutine ix2gb2
  end interface
 
#ifdef LOGGING
  ! Log results for debugging.
  write(g2_log_msg, *) 'getg2i2r: lugb ', lugb, ' msk1 ', msk1, ' msk2 ', msk2, 'idxver', idxver
  call g2_log(1)
#endif
 
  ! 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, lgrib, 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
  use g2logging
  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 g2_gbytec1(in, siout, iskip, nbits)
       character*1, intent(in) :: in(*)
       integer, intent(inout) :: siout
       integer, intent(in) :: iskip, nbits
     end subroutine g2_gbytec1
     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
 
#ifdef LOGGING
  ! Log results for debugging.
  write(g2_log_msg, *) 'getgb2s2: idxver ', idxver, ' nlen ', nlen, &
       ' nnum ', nnum, ' j ', j, ' jdisc ', jdisc
  call g2_log(1)
#endif
 
  ! 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 24 bytes in the version 2 index record, starting
     ! at byte 9.
     inc = 28
  endif
 
  ! Search for request.
  do while (iret .ne. 0 .and. k .lt. nnum)
     k = k + 1
     ! Get length of current index record.
     call g2_gbytec1(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_gbytec1(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_gbytec1(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_gbytec1(cbuf, lsec3, jpos * 8, 4 * 8)  ! get length of gds 
     if (jgdtn .eq. -1) then
        match3 = .true.
     else
        call g2_gbytec1(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_gbytec1(cbuf, lsec4, jpos * 8, 4 * 8)  
     if (jpdtn .eq. -1) then
        match4 = .true.
     else
        ! Get pdt template no.        
        call g2_gbytec1(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_gbytec1(cbuf, gfld%version, (ipos + inc + 40) * 8, 1 * 8)
        call g2_gbytec1(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_gbytec1(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_gbytec1(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, lgrib8
 
  interface
     subroutine ix2gb2(lugb, lskip8, idxver, lgrib8, cbuf, numfld, mlen, iret)
       integer, intent(in) :: lugb
       integer (kind = 8), intent(in) :: lskip8
       integer, intent(in)  :: idxver
       integer (kind = 8), intent(in) :: lgrib8
       character(len = 1), pointer, intent(inout), dimension(:) :: cbuf
       integer, intent(out)  :: numfld, mlen, iret
     end subroutine ix2gb2
  end interface
 
  ! Always use index version 1 from this subroutine.
  lskip8 = lskip
  lgrib8 = lgrib
  call ix2gb2(lugb, lskip8, 1, lgrib8, cbuf, numfld, mlen, iret)
  lgrib = int(lgrib8, 4)
end subroutine ixgb2
 
subroutine ix2gb2(lugb, lskip8, idxver, lgrib8, cbuf, numfld, mlen, iret)
  use re_alloc              ! needed for subroutine realloc
  use g2logging
  implicit none
 
  ! Subroutine parameters.
  integer, intent(in) :: lugb
  integer (kind = 8), intent(in) :: lskip8
  integer, intent(in)  :: idxver
  integer (kind = 8), intent(in) :: lgrib8
  character(len = 1), pointer, intent(inout), dimension(:) :: cbuf
  integer, intent(out)  :: numfld, mlen, iret
 
  character cver, cdisc
  character(len = 4) :: ctemp
  integer (kind = 8) :: loclus8, locgds8, locbms8
  integer locbms
  integer :: indbmp, numsec, newsize, g2_mova2i, mbuf, lindex
  integer :: lskip
  integer :: ilndrs, ilnpds, istat
  integer (kind = 8) :: ibread8, lbread8, ibskip8, lengds8
  integer (kind = 8) :: ilnpds8, ilndrs8
  integer :: lensec, lensec1
  integer :: mypos, inc
 
  ! Parameters.
  ! Size of the internal char buffers used in this subroutine.
  integer :: LINMAX
  parameter(linmax = 5000)
  ! Initial size of buffer to hold index record in cbuf array.
  integer :: INIT
  parameter(init = 50000)
  ! If more space is needed in cbuf, it will get this much more.
  integer :: NEXT
  parameter(next = 10000)
  ! Number of bytes of the BMS section put into index record.
  integer :: MXBMS
  parameter(mxbms = 6)
  integer :: IXDS1, IXDS2
  parameter(ixds1 = 28, ixds2 = 52)
  ! Bytes to skip in (version 1) index record to get to section 0.  
  integer :: IXIDS
  parameter(ixids = 44) 
  integer :: IXSDR
  parameter(ixsdr = 20)
  ! Bytes to skip in (version 1 and 2) index record to get to bms.    
  integer :: IXBMS1, IXBMS2
  parameter(ixbms1 = 24, ixbms2 = 44)
  ! Sizes of integers in bits.
  integer :: INT1_BITS, INT2_BITS, INT4_BITS, INT8_BITS
  parameter(int1_bits = 8, int2_bits = 16, int4_bits = 32, int8_bits = 64)
  ! Location of bytes to drs field in index version 1 and 2.
  integer :: IXDRS1, IXDRS2
  parameter(ixdrs1 = 20, ixdrs2 = 36)
 
  ! Buffers.
  character cbread(LINMAX), cindex(LINMAX)
  character cids(LINMAX), cgds(LINMAX)
 
  interface
     subroutine g2_gbytec1(in, siout, iskip, nbits)
       character*1, intent(in) :: in(*)
       integer, intent(inout) :: siout
       integer, intent(in) :: iskip, nbits
     end subroutine g2_gbytec1
     subroutine g2_sbytec81(out, sin, iskip, nbits)
       character*1, intent(inout) :: out(*)
       integer (kind = 8), intent(in) :: sin
       integer, intent(in) :: iskip, nbits
     end subroutine g2_sbytec81
     subroutine g2_sbytec1(out, in, iskip, nbits)
       character*1, intent(inout) :: out(*)
       integer, intent(in) :: in
       integer, intent(in) :: iskip, nbits
     end subroutine g2_sbytec1
  end interface
 
#ifdef LOGGING
  ! Log results for debugging.
  write(g2_log_msg, *) 'ix2gb2: lugb ', lugb, ' lskip8 ', lskip8, ' idxver ', idxver
  call g2_log(1)
#endif
 
  ! Are we using index version 1 (legacy), or version 2 (introduced to
  ! handle files > 2 GB).
  if (idxver .eq. 1) then
     inc = 0
  else
     ! Add the extra bytes in the version 2 index record, starting at
     ! byte 9. This is because some values early in the index record
     ! changed from 4-byte ints to 8-byte ints. This is the total
     ! extra bytes that were added to the beginning of the index
     ! record in version 2.
     inc = 28
  endif
 
  ! Initialize values and allocate buffer (at the user-provided cbuf
  ! pointer) where the index data will be written. When subroutine is
  ! complete, cbuf will hold either version 1 or 2 index data.
  loclus8 = 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 up to 5000 bytes from the file into buffer cbread. lskip8
  ! must be set to the beginning of a GRIB2 message in the file.
  ibread8 = min(lgrib8, linmax)
  call bareadl(lugb, lskip8, ibread8, lbread8, cbread)
  if (lbread8 .ne. ibread8) then
     iret = 2
     return
  endif
 
  ! Check GRIB version from section 0, must be 2.
  if (cbread(8) .ne. char(2)) then 
     iret = 3
     return
  endif
 
  ! Remember the GRIB version and discipline.
  cver = cbread(8)
  cdisc = cbread(7)
 
  ! Read the length of section 1 from the file data buffer.
  call g2_gbytec1(cbread, lensec1, 16 * int1_bits, int4_bits)
  lensec1 = min(lensec1, int(ibread8, kind(lensec1)))
 
  ! Copy section 1 values into cids array.
  cids(1:lensec1) = cbread(17:16 + lensec1)
  ! do i = 1, lensec1
  !    print *, i, ichar(cids(i))
  ! end do
 
  ! Skip past section 1 in the data buffer.
  ibskip8 = lskip8 + 16_8 + int(lensec1, kind(8))
 
  ! Loop through remaining sections creating an index for each
  ! field. This overwrites the cbread data buffer.
  ibread8 = max(5, mxbms)
  do
     ! Read 6 bytes from file into cbread buffer. (Why 6? Should this
     ! be 5?)
     call bareadl(lugb, ibskip8, ibread8, lbread8, cbread)
 
     ! Check if the first 4 bytes are '7777', indicating end of
     ! message. If we found end of message, we are done.
     ctemp = cbread(1)//cbread(2)//cbread(3)//cbread(4)
     if (ctemp .eq. '7777') return        ! end of message found
 
     ! If we did not find end of message, check that we read 6 bytes.
     if (lbread8 .ne. ibread8) then
        iret = 2
        return
     endif
 
     ! Read the 4-byte section length, and then the 1-byte section
     ! number.
     call g2_gbytec1(cbread, lensec, 0, int4_bits)
     call g2_gbytec1(cbread, numsec, int4_bits, int1_bits)
 
     ! Based on the section number, generate index data for each
     ! section.
     if (numsec .eq. 2) then
        ! Save the location of the local use section in the message.
        loclus8 = ibskip8 - lskip8
     elseif (numsec .eq. 3) then                 
        ! For the GDS section, read the whole section into the cgds
        ! buffer.
        lengds8 = lensec
        cgds = char(0)
        call bareadl(lugb, ibskip8, lengds8, lbread8, cgds)
        if (lbread8 .ne. lengds8) then
           iret = 2
           return
        endif
        ! Remember the GDS location in the message.        
        locgds8 = ibskip8 - lskip8
     elseif (numsec .eq. 4) then
        ! Having found the PDS, we write the beginning of the index
        ! record into the cindex buffer.
        cindex = char(0)
 
        ! The first 4 bytes are the length of the index record. We
        ! don't know that yet, so skip it for now,
        mypos = int4_bits
 
        ! Write the beginning of the index record. It contains bytes
        ! to skip in the file to reach the message, and offsets to
        ! each section within the message. Index version 1 uses 4-byte
        ! ints for these values, index version 2 uses 8-byte ints.
        if (idxver .eq. 1) then
           inc = 0
           lskip = int(lskip8, kind(4))
           call g2_sbytec1(cindex, lskip, mypos, int4_bits)    ! bytes to skip
           mypos = mypos + int4_bits
           call g2_sbytec1(cindex, int(loclus8, kind(4)), mypos, int4_bits)   ! location of local use
           mypos = mypos + int4_bits
           call g2_sbytec1(cindex, int(locgds8, kind(4)), mypos, int4_bits)   ! location of gds
           mypos = mypos + int4_bits
           call g2_sbytec1(cindex, int(ibskip8 - lskip8, kind(4)), mypos, int4_bits)  ! location of pds
#ifdef LOGGING
           write(g2_log_msg, *) ' writing pds location to index: mypos/8 ', mypos/8, &
                ' loc ', int(ibskip8 - lskip8, kind(4))
           call g2_log(4)
#endif
           mypos = mypos + int4_bits * 4 ! skip ahead in cbuf
        else
           inc = 28
           call g2_sbytec81(cindex, lskip8, mypos, int8_bits)    ! bytes to skip
           mypos = mypos + int8_bits
           call g2_sbytec81(cindex, loclus8, mypos, int8_bits)   ! location of local use
           mypos = mypos + int8_bits
           call g2_sbytec81(cindex, locgds8, mypos, int8_bits)   ! location of gds
           mypos = mypos + int8_bits
           call g2_sbytec81(cindex, ibskip8 - lskip8, mypos, int8_bits)  ! location of pds
           mypos = mypos + int8_bits * 4 ! skip ahead in cbuf
        endif
 
#ifdef LOGGING
        write(g2_log_msg, *) ' writing total len to index: mypos/8 ', mypos/8, lgrib8
        call g2_log(4)
#endif
        ! These ints are the same size in index version 1 and 2. The
        ! mypos variable contains the proper offset, which is
        call g2_sbytec81(cindex, lgrib8, mypos, int8_bits)    ! length 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_sbytec1(cindex, numfld + 1, mypos, int2_bits)   ! field num
        mypos = mypos + int2_bits
 
        ! Copy the section 1 values into the cindex buffer.
        cindex(ixids + 1 + inc:ixids + lensec1 + inc) = cids(1:lensec1)
        lindex = ixids + lensec1 + inc
 
        ! Copy the GDS values into the cindex buffer.
        cindex(lindex + 1:lindex + lengds8) = cgds(1:lengds8)
        lindex = lindex + int(lengds8, kind(lindex))
 
        ! Now read the PDS values from the file directly into cindex.
        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
        mypos = mypos + ilnpds
#ifdef LOGGING
        write(g2_log_msg, *) ' after writing pds location to index: mypos/8 ', mypos/8
        call g2_log(4)
#endif
     elseif (numsec .eq. 5) then
        ! Write the byte offset to the DRS section into the cindex buffer.
        !mypos = (IXSDR + inc) * INT1_BITS
#ifdef LOGGING
        write(g2_log_msg, *) ' before writing drs to index: ibskip8 - lskip8 ', ibskip8 - lskip8, ixdrs2
        call g2_log(4)
#endif
        ! Write the bytes to skip to the drs section into the index record.
        if (idxver .eq. 1) then
           call g2_sbytec1(cindex, int(ibskip8 - lskip8, kind(4)), ixdrs1 * int1_bits, int4_bits)
        else
           call g2_sbytec81(cindex, ibskip8 - lskip8, ixdrs2 * int1_bits, int8_bits)  ! location of drs
        endif
 
        ! Read the DRS section directly into the cindex buffer.
        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
        ! Write the location of the BMS section in the message into
        ! the cindex buffer.
        indbmp = g2_mova2i(cbread(6))
#ifdef LOGGING
        write(g2_log_msg, *) ' section 6: indbmp', indbmp
        call g2_log(4)
#endif
        
        if (indbmp .lt. 254) then
           if (idxver .eq. 1) then
              locbms = int(ibskip8 - lskip8, kind(4))
              call g2_sbytec1(cindex, locbms, ixbms1 * int1_bits, int4_bits)  ! loc. of bms
           else
              locbms8 = ibskip8 - lskip8
              call g2_sbytec81(cindex, locbms8, ixbms2 * int1_bits, int8_bits)  ! loc. of bms
           endif
#ifdef LOGGING
           write(g2_log_msg, *) ' section 6: locbms', locbms, 'locbms8', locbms8
           call g2_log(4)
#endif
        elseif (indbmp .eq. 254) then
           if (idxver .eq. 1) then
              call g2_sbytec1(cindex, locbms, ixbms1 * int1_bits, int4_bits)  ! loc. of bms
           else
              call g2_sbytec81(cindex, locbms8, ixbms2 * int1_bits, int8_bits)  ! loc. of bms
           endif
        elseif (indbmp .eq. 255) then
           if (idxver .eq. 1) then
              call g2_sbytec1(cindex, int(ibskip8 - lskip8, kind(4)), ixbms1 * int1_bits, int4_bits)  ! loc. of bms
           else
              call g2_sbytec81(cindex, ibskip8 - lskip8, ixbms2 * int1_bits, int8_bits)  ! loc. of bms
           endif
        endif
        
        ! Copy 6 bytes of the BMS from data buffer to the cindex buffer.
        cindex(lindex + 1:lindex + mxbms) = cbread(1:mxbms)
        lindex = lindex + mxbms
 
        ! The size of the index record is now known, so write it to
        ! the cindex buffer.
        call g2_sbytec1(cindex, lindex, 0, int4_bits)    ! num bytes in index record
     elseif (numsec .eq. 7) then                 ! found data section
 
#ifdef LOGGING
        write(g2_log_msg, *) ' writing  offset to the data in cindex: ibskip8 - lskip8 ', ibskip8 - lskip8, &
             'IXDS2', ixds2
        call g2_log(4)
#endif
        ! Write the offset to the data section in the cindex buffer.
        if (idxver .eq. 1) then
           call g2_sbytec1(cindex, int(ibskip8 - lskip8, kind(4)), ixds1 * int1_bits, int4_bits)
        else
           call g2_sbytec81(cindex, ibskip8 - lskip8, ixds2 * int1_bits, int8_bits)
        endif
           
        ! Increment the field count.
        numfld = numfld + 1
 
        ! Allocate more space in cbuf if necessary. The index record
        ! will be copied there.
        if (lindex + mlen .gt. mbuf) then 
           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
 
        ! Copy the index record into cbuf.
        cbuf(mlen + 1:mlen + lindex) = cindex(1:lindex)
        mlen = mlen + lindex
     else
        ! Unrecognized section.
        iret = 5
        return
     endif
 
     ! Skip past this section in the data buffer.
     ibskip8 = ibskip8 + lensec
  enddo ! next section
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