getfield.F90

Go to the documentation of this file.

 
 
subroutine getfield(cgrib, lcgrib, ifldnum, igds, igdstmpl, &
     igdslen, ideflist, idefnum, ipdsnum, ipdstmpl, ipdslen, &
     coordlist, numcoord, ndpts, idrsnum, idrstmpl, idrslen,  &
     ibmap, bmap, fld, ierr)
 
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib, ifldnum
  integer, intent(out) :: igds(*), igdstmpl(*), ideflist(*)
  integer, intent(out) :: ipdsnum, ipdstmpl(*)
  integer, intent(out) :: idrsnum, idrstmpl(*)
  integer, intent(out) :: ndpts, ibmap, idefnum, numcoord
  integer, intent(out) :: ierr
  logical*1, intent(out) :: bmap(*)
  real, intent(out) :: fld(*), coordlist(*)
 
  character(len = 4), parameter :: grib = 'GRIB', c7777 = '7777'
  character(len = 4) :: ctemp
  integer :: listsec0(2)
  integer :: iofst, istart
  integer(4) :: ieee
  logical :: have3, have4, have5, have6, have7
 
  !implicit none additions
  integer, intent(out) :: igdslen, ipdslen, idrslen
  integer :: numfld, j, lengrib, lensec0, ipos
  integer :: lensec, isecnum, jerr, ier, numlocal
 
  have3 = .false.
  have4 = .false.
  have5 = .false.
  have6 = .false.
  have7 = .false.
  ierr = 0
  numfld = 0
  numlocal = 0
 
  !     Check for valid request number
  if (ifldnum .le. 0) then
     print *, 'getfield: Request for field number ' &
          ,'must be positive.'
     ierr = 3
     return
  endif
 
  !     Check for beginning of GRIB message in the first 100 bytes
  istart = 0
  do j = 1, 100
     ctemp = cgrib(j) // cgrib(j + 1) // cgrib(j + 2) // cgrib(j + 3)
     if (ctemp .eq. grib) then
        istart = j
        exit
     endif
  enddo
  if (istart .eq. 0) then
     print *, 'getfield:  Beginning characters GRIB not found.'
     ierr = 1
     return
  endif
 
  !     Unpack Section 0 - Indicator Section
  iofst = 8 * (istart + 5)
  call g2_gbytec(cgrib, listsec0(1), iofst, 8) ! Discipline
  iofst = iofst + 8
  call g2_gbytec(cgrib, listsec0(2), iofst, 8) ! GRIB edition number
  iofst = iofst + 8
  iofst = iofst + 32
  call g2_gbytec(cgrib, lengrib, iofst, 32) ! Length of GRIB message
  iofst = iofst + 32
  lensec0 = 16
  ipos = istart + lensec0
 
  !     Currently handles only GRIB Edition 2.
  if (listsec0(2) .ne. 2) then
     print *, 'getfield: can only decode GRIB edition 2.'
     ierr = 2
     return
  endif
 
  !     Loop through the remaining sections keeping track of the length of
  !     each. Also keep the latest Grid Definition Section info.  Unpack
  !     the requested field number.
  do
     !         Check to see if we are at end of GRIB message
     ctemp = cgrib(ipos) // cgrib(ipos + 1) // cgrib(ipos + 2) // &
          cgrib(ipos + 3)
     if (ctemp .eq. c7777) then
        ipos = ipos + 4
        !             If end of GRIB message not where expected, issue error
        if (ipos.ne.(istart + lengrib)) then
           print *, 'getfield: "7777" found, but not ' &
                ,'where expected.'
           ierr = 4
           return
        endif
        exit
     endif
     !         Get length of Section and Section number
     iofst = (ipos - 1) * 8
     call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section
     iofst = iofst + 32
     call g2_gbytec(cgrib, isecnum, iofst, 8) ! Get Section number
     iofst = iofst + 8
 
     !         If found Section 3, unpack the GDS info using the appropriate
     !         template. Save in case this is the latest grid before the
     !         requested field.
     if (isecnum .eq. 3) then
        iofst = iofst - 40    ! reset offset to beginning of section
        call unpack3(cgrib, lcgrib, iofst, igds, igdstmpl,  &
             igdslen, ideflist, idefnum, jerr)
        if (jerr .eq. 0) then
           have3 = .true.
        else
           ierr = 10
           return
        endif
     endif
 
     !         If found Section 4, check to see if this field is the one
     !         requested.
     if (isecnum .eq. 4) then
        numfld = numfld + 1
        if (numfld .eq. ifldnum) then
           iofst = iofst - 40 ! reset offset to beginning of section
           call unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
                ipdslen, coordlist, numcoord, jerr)
           if (jerr .eq. 0) then
              have4 = .true.
           else
              ierr = 11
              return
           endif
        endif
     endif
 
     !         If found Section 5, check to see if this field is the one
     !         requested.
     if ((isecnum .eq. 5) .and. (numfld .eq. ifldnum)) then
        iofst = iofst - 40    ! reset offset to beginning of section
        call unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum, &
             idrstmpl, idrslen, jerr)
        if (jerr .eq. 0) then
           have5 = .true.
        else
           ierr = 12
           return
        endif
     endif
 
     !         If found Section 6, Unpack bitmap. Save in case this is the
     !         latest bitmap before the requested field.
     if (isecnum .eq. 6) then
        iofst = iofst - 40    ! reset offset to beginning of section
        call unpack6(cgrib, lcgrib, iofst, igds(2), ibmap, bmap, &
             jerr)
        if (jerr .eq. 0) then
           have6 = .true.
        else
           ierr = 13
           return
        endif
     endif
 
     !         If found Section 7, check to see if this field is the one
     !         requested.
     if ((isecnum .eq. 7) .and. (numfld .eq. ifldnum)) then
        if (idrsnum .eq. 0) then
           call simunpack(cgrib(ipos + 5), lensec - 6, idrstmpl, &
                ndpts, fld)
           have7 = .true.
        elseif (idrsnum .eq. 2 .or. idrsnum .eq. 3) then
           call comunpack(cgrib(ipos + 5), lensec - 6, lensec, &
                idrsnum,idrstmpl, ndpts, fld, ier)
           if (ier .ne. 0) then
              ierr = 14
              return
           endif
           have7 = .true.
        elseif (idrsnum .eq. 50) then
           call simunpack(cgrib(ipos + 5), lensec - 6, idrstmpl, &
                ndpts - 1, fld(2))
           ieee = idrstmpl(5)
           call rdieee(ieee, fld(1), 1)
           have7 = .true.
        elseif (idrsnum .eq. 40 .or. idrsnum .eq. 40000) then
           call jpcunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
                ndpts, fld)
           have7 = .true.
        elseif (idrsnum .eq. 41 .or. idrsnum .eq. 40010) then
           call pngunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
                ndpts, fld)
           have7 = .true.
        else
           print *, 'getfield: Data Representation Template ', &
                idrsnum, ' not yet implemented.'
           ierr = 9
           return
        endif
     endif
 
     !         Check to see if we read pass the end of the GRIB message and
     !         missed the terminator string '7777'.
     ipos = ipos + lensec      ! Update beginning of section pointer
     if (ipos .gt. (istart + lengrib)) then
        print *, 'getfield: "7777"  not found at end' &
             ,' of GRIB message.'
        ierr = 7
        return
     endif
 
     if (have3 .and. have4 .and. have5 .and. have6 .and. have7) &
          return
 
  enddo
 
  !     If exited from above loop, the end of the GRIB message was reached
  !     before the requested field was found.
  print *, 'getfield: GRIB message contained ', numlocal,  &
       ' different fields.'
  print *, 'getfield: The request was for the ', ifldnum,  &
       ' field.'
  ierr = 6
 
end subroutine getfield
 
 
subroutine unpack3(cgrib, lcgrib, iofst, igds, igdstmpl,  &
     mapgridlen, ideflist, idefnum, ierr)
 
  use gridtemplates
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib
  integer, intent(inout) :: iofst
  integer, intent(out) :: igds(*), igdstmpl(*), ideflist(*)
  integer, intent(out) :: ierr, idefnum
 
  integer, allocatable :: mapgrid(:)
  integer :: mapgridlen, ibyttem
  logical needext
 
  !implicit none additions
  integer :: lensec, iret, i, nbits, isign, newmapgridlen
 
  ierr = 0
 
  call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section
  iofst = iofst + 32
  iofst = iofst + 8             ! skip section number
 
  call g2_gbytec(cgrib, igds(1), iofst, 8) ! Get source of Grid def.
  iofst = iofst + 8
  call g2_gbytec(cgrib, igds(2), iofst, 32) ! Get number of grid pts.
  iofst = iofst + 32
  call g2_gbytec(cgrib, igds(3), iofst, 8) ! Get num octets for opt. list
  iofst = iofst + 8
  call g2_gbytec(cgrib, igds(4), iofst, 8) ! Get interpret. for opt. list
  iofst = iofst + 8
  call g2_gbytec(cgrib, igds(5), iofst, 16) ! Get Grid Def Template num.
  iofst = iofst + 16
  if (igds(1) .eq. 0) then
     !      if (igds(1).eq.0.OR.igds(1).eq.255) then  ! FOR ECMWF TEST ONLY
     allocate(mapgrid(lensec))
     !         Get Grid Definition Template
     call getgridtemplate(igds(5), mapgridlen, mapgrid, needext,  &
          iret)
     if (iret .ne. 0) then
        ierr = 5
        return
     endif
  else
     !        igdstmpl = -1
     mapgridlen = 0
     needext = .false.
  endif
 
  !     Unpack each value into array igdstmpl from the the appropriate
  !     number of octets, which are specified in corresponding entries in
  !     array mapgrid.
  ibyttem = 0
  do i = 1, mapgridlen
     nbits = iabs(mapgrid(i)) * 8
     if (mapgrid(i) .ge. 0) then
        call g2_gbytec(cgrib, igdstmpl(i), iofst, nbits)
     else
        call g2_gbytec(cgrib, isign, iofst, 1)
        call g2_gbytec(cgrib, igdstmpl(i), iofst + 1, nbits-1)
        if (isign .eq. 1) igdstmpl(i) = -igdstmpl(i)
     endif
     iofst = iofst + nbits
     ibyttem = ibyttem + iabs(mapgrid(i))
  enddo
 
  !      Check to see if the Grid Definition Template needs to be
  !      extended. The number of values in a specific template may vary
  !      depending on data specified in the "static" part of the template.
  if (needext) then
     call extgridtemplate(igds(5), igdstmpl, newmapgridlen, &
          mapgrid)
     !         Unpack the rest of the Grid Definition Template
     do i = mapgridlen + 1, newmapgridlen
        nbits = iabs(mapgrid(i)) * 8
        if (mapgrid(i) .ge. 0) then
           call g2_gbytec(cgrib, igdstmpl(i), iofst, nbits)
        else
           call g2_gbytec(cgrib, isign, iofst, 1)
           call g2_gbytec(cgrib, igdstmpl(i), iofst + 1, nbits - &
                1)
           if (isign .eq. 1) igdstmpl(i) = -igdstmpl(i)
        endif
        iofst = iofst + nbits
        ibyttem = ibyttem + iabs(mapgrid(i))
     enddo
     mapgridlen = newmapgridlen
  endif
 
  !     Unpack optional list of numbers defining number of points in each
  !     row or column, if included. This is used for non regular grids.
  if (igds(3) .ne. 0) then
     nbits = igds(3) * 8
     idefnum = (lensec - 14 - ibyttem) / igds(3)
     call g2_gbytesc(cgrib, ideflist, iofst, nbits, 0, idefnum)
     iofst = iofst + (nbits * idefnum)
  else
     idefnum = 0
  endif
  if (allocated(mapgrid)) deallocate(mapgrid)
end subroutine unpack3
 
subroutine unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
     mappdslen, coordlist, numcoord, ierr)
 
  use pdstemplates
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib
  integer, intent(inout) :: iofst
  real, intent(out) :: coordlist(*)
  integer, intent(out) :: ipdsnum, ipdstmpl(*)
  integer, intent(out) :: ierr, numcoord
 
  real(4), allocatable :: coordieee(:)
  integer, allocatable :: mappds(:)
  integer :: mappdslen
  logical needext
 
  !implicit none additions
  integer :: lensec, iret, i, nbits, isign, newmappdslen
 
  ierr = 0
 
  call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section
  iofst = iofst + 32
  iofst = iofst + 8             ! skip section number
  allocate(mappds(lensec))
 
  call g2_gbytec(cgrib, numcoord, iofst, 16) ! Get num of coordinate values
  iofst = iofst + 16
  call g2_gbytec(cgrib, ipdsnum, iofst, 16) ! Get Prod. Def Template num.
  iofst = iofst + 16
  !     Get Product Definition Template.
  call getpdstemplate(ipdsnum, mappdslen, mappds, needext, iret)
  if (iret.ne.0) then
     ierr = 5
     return
  endif
 
  !     Unpack each value into array ipdstmpl from the the appropriate
  !     number of octets, which are specified in corresponding entries in
  !     array mappds.
  do i = 1, mappdslen
     nbits = iabs(mappds(i))*8
     if (mappds(i).ge.0) then
        call g2_gbytec(cgrib, ipdstmpl(i), iofst, nbits)
     else
        call g2_gbytec(cgrib, isign, iofst, 1)
        call g2_gbytec(cgrib, ipdstmpl(i), iofst + 1, nbits-1)
        if (isign.eq.1) ipdstmpl(i) = -ipdstmpl(i)
     endif
     iofst = iofst + nbits
  enddo
 
  !     Check to see if the Product Definition Template needs to be
  !     extended. The number of values in a specific template may vary
  !     depending on data specified in the "static" part of the template.
  if (needext) then
     call extpdstemplate(ipdsnum, ipdstmpl, newmappdslen, mappds)
 
     !         Unpack the rest of the Product Definition Template
     do i = mappdslen + 1, newmappdslen
        nbits = iabs(mappds(i))*8
        if (mappds(i).ge.0) then
           call g2_gbytec(cgrib, ipdstmpl(i), iofst, nbits)
        else
           call g2_gbytec(cgrib, isign, iofst, 1)
           call g2_gbytec(cgrib, ipdstmpl(i), iofst + 1, nbits-1)
           if (isign.eq.1) ipdstmpl(i) = -ipdstmpl(i)
        endif
        iofst = iofst + nbits
     enddo
     mappdslen = newmappdslen
  endif
 
  !     Get Optional list of vertical coordinate values after the Product
  !     Definition Template, if necessary.
  if (numcoord .ne. 0) then
     allocate (coordieee(numcoord))
     call g2_gbytesc(cgrib, coordieee, iofst, 32, 0, numcoord)
     call rdieee(coordieee, coordlist, numcoord)
     deallocate (coordieee)
     iofst = iofst + (32*numcoord)
  endif
  if (allocated(mappds)) deallocate(mappds)
end subroutine unpack4
 
subroutine unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum,  &
     idrstmpl, mapdrslen, ierr)
 
  use drstemplates
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib
  integer, intent(inout) :: iofst
  integer, intent(out) :: ndpts, idrsnum, idrstmpl(*)
  integer, intent(out) :: ierr
 
  !      integer, allocatable :: mapdrs(:)
  integer, allocatable :: mapdrs(:)
  integer :: mapdrslen
  logical needext
 
  !implicit none additions
  integer :: lensec, i, nbits, isign, newmapdrslen, iret
 
  ierr = 0
 
  call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section
  iofst = iofst + 32
  iofst = iofst + 8             ! skip section number
  allocate(mapdrs(lensec))
 
  call g2_gbytec(cgrib, ndpts, iofst, 32) ! Get num of data points
  iofst = iofst + 32
  call g2_gbytec(cgrib, idrsnum, iofst, 16) ! Get Data Rep Template Num.
  iofst = iofst + 16
  !     Gen Data Representation Template
  call getdrstemplate(idrsnum, mapdrslen, mapdrs, needext, iret)
  if (iret.ne.0) then
     ierr = 7
     return
  endif
 
  !     Unpack each value into array ipdstmpl from the the appropriate
  !     number of octets, which are specified in corresponding entries in
  !     array mappds.
  do i = 1, mapdrslen
     nbits = iabs(mapdrs(i))*8
     if (mapdrs(i).ge.0) then
        call g2_gbytec(cgrib, idrstmpl(i), iofst, nbits)
     else
        call g2_gbytec(cgrib, isign, iofst, 1)
        call g2_gbytec(cgrib, idrstmpl(i), iofst + 1, nbits-1)
        if (isign.eq.1) idrstmpl(i) = -idrstmpl(i)
     endif
     iofst = iofst + nbits
  enddo
 
  !     Check to see if the Data Representation Template needs to be
  !     extended. The number of values in a specific template may vary
  !     depending on data specified in the "static" part of the template.
  if (needext) then
     call extdrstemplate(idrsnum, idrstmpl, newmapdrslen, mapdrs)
     !         Unpack the rest of the Data Representation Template
     do i = mapdrslen + 1, newmapdrslen
        nbits = iabs(mapdrs(i))*8
        if (mapdrs(i).ge.0) then
           call g2_gbytec(cgrib, idrstmpl(i), iofst, nbits)
        else
           call g2_gbytec(cgrib, isign, iofst, 1)
           call g2_gbytec(cgrib, idrstmpl(i), iofst + 1, nbits - 1)
           if (isign.eq.1) idrstmpl(i) = -idrstmpl(i)
        endif
        iofst = iofst + nbits
     enddo
     mapdrslen = newmapdrslen
  endif
  if (allocated(mapdrs)) deallocate(mapdrs)
end subroutine unpack5
 
subroutine unpack6(cgrib, lcgrib, iofst, ngpts, ibmap, bmap, ierr)
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib, ngpts
  integer, intent(inout) :: iofst
  integer, intent(out) :: ibmap
  integer, intent(out) :: ierr
  logical*1, intent(out) :: bmap(ngpts)
 
  integer :: intbmap(ngpts)
 
  !implicit none additions
  integer :: j
 
  ierr = 0
 
  iofst = iofst + 32            ! skip Length of Section
  iofst = iofst + 8             ! skip section number
 
  call g2_gbytec(cgrib, ibmap, iofst, 8) ! Get bit-map indicator
  iofst = iofst + 8
 
  if (ibmap.eq.0) then      ! Unpack bitmap
     call g2_gbytesc(cgrib, intbmap, iofst, 1, 0, ngpts)
     iofst = iofst + ngpts
     do j = 1, ngpts
        bmap(j) = .true.
        if (intbmap(j).eq.0) bmap(j) = .false.
     enddo
  elseif (ibmap.eq.254) then ! Use previous bitmap
     return
  elseif (ibmap.eq.255) then ! No bitmap in message
     bmap(1:ngpts) = .true.
  else
     print *, 'unpack6: Predefined bitmap ', ibmap, &
          ' not recognized.'
     ierr = 4
  endif
end subroutine unpack6