g2unpack.F90

Go to the documentation of this file.

 
 
subroutine gf_unpack1(cgrib, lcgrib, iofst, ids, idslen, ierr)
  implicit none
 
  character(len=1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib
  integer, intent(inout) :: iofst
  integer, pointer, dimension(:) :: ids
  integer, intent(out) :: ierr, idslen
  integer, dimension(:) :: mapid(13)
  integer :: i, istat, lensec, nbits
 
  data mapid /2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1/
 
  ierr = 0
  idslen = 13
  nullify(ids)
 
  call g2_gbytec(cgrib, lensec, iofst, 32)        ! Get Length of Section
  iofst = iofst + 32
  iofst = iofst + 8     ! skip section number
 
  ! Unpack each value into array ids from the the appropriate number
  ! of octets, which are specified in corresponding entries in array
  ! mapid.
  istat = 0
  allocate(ids(idslen), stat = istat)
  if (istat .ne. 0) then
     ierr = 6
     nullify(ids)
     return
  endif
 
  do i = 1, idslen
     nbits = mapid(i) * 8
     call g2_gbytec(cgrib, ids(i), iofst, nbits)
     iofst = iofst + nbits
  enddo
end subroutine gf_unpack1
 
subroutine gf_unpack2(cgrib, lcgrib, iofst, lencsec2, csec2, ierr)
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib
  integer, intent(inout) :: iofst
  integer, intent(out) :: lencsec2
  integer, intent(out) :: ierr
  character(len = 1), pointer, dimension(:) :: csec2
  integer :: lensec, istat, isecnum, ipos
 
  ierr = 0
  lencsec2 = 0
  nullify(csec2)
 
  call g2_gbytec(cgrib, lensec, iofst, 32)        ! Get Length of Section
  iofst = iofst + 32    
  lencsec2 = lensec-5
  call g2_gbytec(cgrib, isecnum, iofst, 8)         ! Get Section Number
  iofst = iofst + 8     
  ipos = (iofst / 8) + 1
 
  if (isecnum .ne. 2) then
     ierr = 6
     print *, 'gf_unpack2: Not Section 2 data. '
     return
  endif
 
  allocate(csec2(lencsec2), stat = istat)
  if (istat .ne. 0) then
     ierr = 6
     nullify(csec2)
     return
  endif
 
  csec2(1:lencsec2) = cgrib(ipos:ipos + lencsec2 - 1)
  iofst = iofst + (lencsec2 * 8)
 
end subroutine gf_unpack2
 
subroutine gf_unpack3(cgrib, lcgrib, iofst, igds, igdstmpl, &
     mapgridlen, ideflist, idefnum, ierr)
 
  use gridtemplates
  use re_alloc              !  needed for subroutine realloc
  implicit none
 
  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
 
  integer, allocatable :: mapgrid(:)
  integer, intent(out) :: mapgridlen
  integer :: ibyttem
  logical needext
  integer :: lensec, istat, i, nbits, isign, newmapgridlen, iret
 
  ierr = 0
  nullify(igdstmpl, ideflist)
 
  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 .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
        if (allocated(mapgrid)) deallocate(mapgrid)
        return
     endif
  else
     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.
  istat = 0
  if (mapgridlen .gt. 0) allocate(igdstmpl(mapgridlen), stat = istat)
  if (istat .ne. 0) then
     ierr = 6
     nullify(igdstmpl)
     if (allocated(mapgrid)) deallocate(mapgrid)
     return
  endif
  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.
     call realloc(igdstmpl, mapgridlen, newmapgridlen, istat)
     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
  if (allocated(mapgrid)) deallocate(mapgrid)
 
  ! 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)
     istat = 0
     if (idefnum .gt. 0) allocate(ideflist(idefnum), stat = istat)
     if (istat .ne. 0) then
        ierr = 6
        nullify(ideflist)
        return
     endif
     call g2_gbytesc(cgrib, ideflist, iofst, nbits, 0, idefnum)
     iofst = iofst + (nbits * idefnum)
  else
     idefnum = 0
     nullify(ideflist)
  endif
end subroutine gf_unpack3
 
subroutine gf_unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
     mappdslen, coordlist, numcoord, ierr)
  use pdstemplates
  use re_alloc              !  needed for subroutine realloc
  implicit none
 
  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
 
  real(4), allocatable :: coordieee(:)
  integer, allocatable :: mappds(:)
  integer :: mappdslen
  logical needext
  integer :: lensec, nbits, newmappdslen
  integer :: istat1, istat, isign, iret, i
 
  ierr = 0
  nullify(ipdstmpl, coordlist)
 
  ! Get Length of Section.
  call g2_gbytec(cgrib, lensec, iofst, 32)
  iofst = iofst + 32
  iofst = iofst + 8     ! skip section number
  allocate(mappds(lensec))
 
  ! Get num of coordinate values.
  call g2_gbytec(cgrib, numcoord, iofst, 16)
  iofst = iofst + 16
  ! Get Prod. Def Template num.
  call g2_gbytec(cgrib, ipdsnum, iofst, 16)
  iofst = iofst + 16
  ! Get Product Definition Template.
  call getpdstemplate(ipdsnum, mappdslen, mappds, needext, iret)
  if (iret.ne.0) then
     ierr = 5
     if (allocated(mappds)) deallocate(mappds)
     return
  endif
 
  ! Unpack each value into array ipdstmpl from the the appropriate
  ! number of octets, which are specified in corresponding entries in
  ! array mappds.
  istat = 0
  if (mappdslen.gt.0) allocate(ipdstmpl(mappdslen), stat = istat)
  if (istat.ne.0) then
     ierr = 6
     nullify(ipdstmpl)
     if (allocated(mappds)) deallocate(mappds)
     return
  endif
  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)
     call realloc(ipdstmpl, mappdslen, newmappdslen, istat)
     ! 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
  if (allocated(mappds)) deallocate(mappds)
 
  ! Get Optional list of vertical coordinate values
  ! after the Product Definition Template, if necessary.
  nullify(coordlist)
  if (numcoord .ne. 0) then
     allocate (coordieee(numcoord), stat = istat1)
     allocate(coordlist(numcoord), stat = istat)
     if ((istat1 + istat).ne.0) then
        ierr = 6
        nullify(coordlist)
        if (allocated(coordieee)) deallocate(coordieee)
        return
     endif
     call g2_gbytescr(cgrib, coordieee, iofst, 32, 0, numcoord)
     call rdieee(coordieee, coordlist, numcoord)
     deallocate (coordieee)
     iofst = iofst + (32 * numcoord)
  endif
end subroutine gf_unpack4
 
subroutine gf_unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum, idrstmpl, &
     mapdrslen, ierr)
  use drstemplates
  use re_alloc              !  needed for subroutine realloc
  implicit none
 
  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
 
  integer, allocatable :: mapdrs(:)
  integer :: mapdrslen
  logical needext
  integer :: newmapdrslen, nbits, istat, isign, lensec, iret, i
 
  ierr = 0
  nullify(idrstmpl)
 
  call g2_gbytec(cgrib, lensec, iofst, 32)        ! Get Length of Section
  iofst = iofst + 32
  iofst = iofst + 8     ! skip section number
  allocate(mapdrs(lensec))
 
  ! Get num of data points.
  call g2_gbytec(cgrib, ndpts, iofst, 32)
  iofst = iofst + 32
  ! Get Data Rep Template Num.
  call g2_gbytec(cgrib, idrsnum, iofst, 16)
  iofst = iofst + 16
  ! Gen Data Representation Template.
  call getdrstemplate(idrsnum, mapdrslen, mapdrs, needext, iret)
  if (iret .ne. 0) then
     ierr = 7
     if (allocated(mapdrs)) deallocate(mapdrs)
     return
  endif
 
  ! Unpack each value into array ipdstmpl from the the appropriate
  ! number of octets, which are specified in corresponding entries in
  ! array mappds.
  istat = 0
  if (mapdrslen .gt. 0) allocate(idrstmpl(mapdrslen), stat = istat)
  if (istat .ne. 0) then
     ierr = 6
     nullify(idrstmpl)
     if (allocated(mapdrs)) deallocate(mapdrs)
     return
  endif
  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)
     call realloc(idrstmpl, mapdrslen, newmapdrslen, istat)
 
     ! 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 gf_unpack5
 
subroutine gf_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, pointer, dimension(:) :: bmap
  integer :: intbmap(ngpts)
  integer :: istat, j
 
  ierr = 0
  nullify(bmap)
 
  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
     istat = 0
     if (ngpts .gt. 0) allocate(bmap(ngpts), stat = istat)
     if (istat .ne. 0) then
        ierr = 6
        nullify(bmap)
        return
     endif
     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
  endif
 
end subroutine gf_unpack6
 
subroutine gf_unpack7(cgrib, lcgrib, iofst, igdsnum, igdstmpl,  &
     idrsnum, idrstmpl, ndpts, fld, ierr)
  implicit none
 
  character(len = 1), intent(in) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib, ndpts, igdsnum, idrsnum
  integer, intent(inout) :: iofst
  integer, pointer, dimension(:) :: igdstmpl, idrstmpl
  integer, intent(out) :: ierr
  real, pointer, dimension(:) :: fld
  integer :: ier,  ipos,  istat,  lensec
  real (kind = 4) :: ieee(1)
  real :: tmpfld(1)
 
  ierr = 0
  nullify(fld)
 
  call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section
  iofst = iofst + 32
  iofst = iofst + 8 ! skip section number
 
  ipos = (iofst/8) + 1
  istat = 0
  allocate(fld(ndpts), stat = istat)
  if (istat.ne.0) then
     ierr = 6
     return
  endif
 
  if (idrsnum .eq. 0) then
     call simunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld)
  elseif (idrsnum.eq.2.or.idrsnum.eq.3) then
     call comunpack(cgrib(ipos), lensec-5, lensec, idrsnum, idrstmpl, ndpts, fld, ier)
     if (ier .ne. 0) then
        ierr = 7
        return
     endif
  elseif (idrsnum .eq. 50) then ! Spectral simple
     call simunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts-1, fld(2))
     ieee = transfer(idrstmpl(5), ieee, 1)
     call rdieee(ieee, tmpfld, 1)
     fld(1) = tmpfld(1)
  elseif (idrsnum .eq. 51) then ! Spectral complex
     if (igdsnum.ge.50.AND.igdsnum.le.53) then
        call specunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts,  &
             igdstmpl(1), igdstmpl(2), igdstmpl(3), fld)
     else
        print *, 'gf_unpack7: Cannot use GDT 3.', igdsnum, ' to unpack Data Section 5.51.'
        ierr = 5
        nullify(fld)
        return
     endif
  elseif (idrsnum .eq. 40 .OR. idrsnum .eq. 40000) then
     call jpcunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld)
  elseif (idrsnum .eq. 41 .OR. idrsnum .eq. 40010) then
     call pngunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld)
  else
     print *, 'gf_unpack7: Data Representation Template ', idrsnum, ' not yet implemented.'
     ierr = 4
     nullify(fld)
     return
  endif
 
  iofst = iofst + (8 * lensec)
 
end subroutine gf_unpack7