g2create.F90

Go to the documentation of this file.

 
 
subroutine gribcreate(cgrib, lcgrib, listsec0, listsec1, ierr)
  implicit none
 
  character(len = 1), intent(inout) :: cgrib(lcgrib)
  integer, intent(in) :: listsec0(*), listsec1(*)
  integer, intent(in) :: lcgrib
  integer, intent(out) :: ierr
 
  integer :: i, lensec1, nbits
  character(len = 4), parameter :: grib = 'GRIB'
  integer, parameter :: ZERO = 0, one = 1
  integer, parameter :: mapsec1len = 13
  integer, parameter :: mapsec1(mapsec1len) = (/ 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1 /)
  integer lensec0, iofst, ibeg
    
  interface
     subroutine g2_sbytec(out, in, iskip, nbits)
       character*1, intent(inout) :: out(*)
       integer, intent(in) :: in(*)
       integer, intent(in) :: iskip, nbits
     end subroutine g2_sbytec
     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
 
  ierr = 0
 
#ifdef LOGGING
  print *, 'gribcreate lcgrib ', lcgrib
#endif
 
  ! Currently handles only GRIB Edition 2.
  if (listsec0(2) .ne. 2) then
     print *, 'gribcreate: can only code GRIB edition 2.'
     ierr = 1
     return
  endif
 
  !  Pack Section 0 - Indicator Section (except for total length of
  !  GRIB message).
  cgrib(1) = grib(1:1) ! Beginning of GRIB message
  cgrib(2) = grib(2:2)
  cgrib(3) = grib(3:3)
  cgrib(4) = grib(4:4)
  call g2_sbytec1(cgrib, zero, 32, 16)           ! reserved for future use
  call g2_sbytec(cgrib, listsec0(1), 48, 8)     ! Discipline
  call g2_sbytec(cgrib, listsec0(2), 56, 8)     ! GRIB edition number
  lensec0 = 16      ! bytes (octets)
 
  ! Pack Section 1 - Identification Section.
  ibeg = lensec0 * 8        !   Calculate offset for beginning of section 1
  iofst = ibeg + 32         !   leave space for length of section
  call g2_sbytec1(cgrib, one, iofst, 8)     ! Store section number ( 1 )
  iofst = iofst + 8
 
  ! Pack up each input value in array listsec1 into the the
  ! appropriate number of octets, which are specified in corresponding
  ! entries in array mapsec1.
  do i = 1, mapsec1len
     nbits = mapsec1(i) * 8
     call g2_sbytec(cgrib, listsec1(i), iofst, nbits)
     iofst = iofst + nbits
  enddo
 
  ! Calculate length of section 1 and store it in octets 1-4 of
  ! section 1.
  lensec1 = (iofst - ibeg) / 8
  call g2_sbytec1(cgrib, lensec1, ibeg, 32)
 
  ! Put current byte total of message into Section 0.
  call g2_sbytec1(cgrib, zero, 64, 32)
  call g2_sbytec1(cgrib, lensec0 + lensec1, 96, 32)
end subroutine gribcreate
 
subroutine addfield(cgrib, lcgrib, ipdsnum, ipdstmpl, ipdstmplen, &
     coordlist, numcoord, idrsnum, idrstmpl, &
     idrstmplen, fld, ngrdpts, ibmap, bmap, ierr)
  use pdstemplates
  use drstemplates
  implicit none
 
  logical :: match
  character(len=1), intent(inout) :: cgrib(lcgrib)
  integer, intent(in) :: ipdsnum, ipdstmpl(*)
  integer, intent(in) :: idrsnum, numcoord, ipdstmplen, idrstmplen
  integer, intent(in) :: lcgrib, ngrdpts, ibmap
  real, intent(in) :: coordlist(numcoord)
  real(kind = 4) :: coordlist_4(numcoord)
  real, target, intent(in) :: fld(ngrdpts)
  integer, intent(out) :: ierr
  integer, intent(inout) :: idrstmpl(*)
  logical*1, intent(in) :: bmap(ngrdpts)
 
  character(len=4), parameter :: grib='GRIB', c7777='7777'
  character(len=4):: ctemp
  character(len=1), allocatable :: cpack(:)
  real, pointer, dimension(:) :: pfld
  real(4) :: coordieee(numcoord), re00, tmpre00(1)
  integer(4) :: ire00, allones
  integer :: mappds(ipdstmplen), intbmap(ngrdpts), mapdrs(idrstmplen)
  integer, parameter :: zero=0, one=1, four=4, five=5, six=6, seven=7
  integer, parameter :: minsize=50000
  integer iofst, ibeg, lencurr, len, mappdslen, mapdrslen, lpos3
  integer width, height, ndpts
  integer lensec3, lensec4, lensec5, lensec6, lensec7
  logical issec3, needext, isprevbmap
  integer :: nbits, newlen, nsize, lcpack, left
  integer :: ibmprev, ilen, ioctet, iscan, isecnum, itemp
  integer :: i, jj, kk, mm
  integer :: iret, istat
  real (kind = 4) :: tmpfld(1)
 
  interface
     subroutine g2_gbytec(in, iout, iskip, nbits)
       character*1, intent(in) :: in(*)
       integer, intent(inout) :: iout(*)
       integer, intent(in) :: iskip, nbits
     end subroutine g2_gbytec
     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_gbytec81(in, siout, iskip, nbits)
       character*1, intent(in) :: in(*)
       integer (kind = 8), intent(inout) :: siout
       integer, intent(in) :: iskip, nbits
       integer (kind = 8) :: iout(1)
     end subroutine g2_gbytec81
     subroutine g2_sbytec(out, in, iskip, nbits)
       character*1, intent(inout) :: out(*)
       integer, intent(in) :: in(*)
       integer, intent(in) :: iskip, nbits
     end subroutine g2_sbytec
     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
 
  allones = int(z'FFFFFFFF')
  ierr = 0
 
  ! Check to see if beginning of GRIB message exists.
  match = .true.
  do i = 1, 4
     if (cgrib(i) /= grib(i:i)) then
        match = .false.
     endif
  enddo
  if (.not. match) then
     print *, 'addfield: GRIB not found in given message.'
     print *, 'addfield: Call to routine gribcreate required to initialize GRIB messge.'
     ierr = 1
     return
  endif
 
  ! Get current length of GRIB message.
  call g2_gbytec1(cgrib, lencurr, 96, 32)
 
  ! Check to see if GRIB message is already complete.
  ctemp = cgrib(lencurr-3) // cgrib(lencurr - 2) // cgrib(lencurr - 1) // cgrib(lencurr)
  if (ctemp .eq. c7777) then
     print *, 'addfield: GRIB message already complete.  Cannot add new section.'
     ierr = 2
     return
  endif
 
  ! Loop through all current sections of the GRIB message to find the
  ! last section number.
  issec3 = .false.
  isprevbmap = .false.
  len = 16 ! length of Section 0
  do
     ! Get number and length of next section
     iofst = len * 8
     call g2_gbytec1(cgrib, ilen, iofst, 32)
     iofst = iofst + 32
     call g2_gbytec1(cgrib, isecnum, iofst, 8)
     iofst = iofst + 8
     
     ! Check if previous Section 3 exists and save location of
     ! the section 3 in case needed later.
     if (isecnum .eq. 3) then
        issec3 = .true.
        lpos3 = len + 1
        lensec3 = ilen
     endif
     
     ! Check if a previous defined bitmap exists
     if (isecnum .eq. 6) then
        call g2_gbytec1(cgrib, ibmprev, iofst, 8)
        iofst = iofst + 8
        if ((ibmprev .ge. 0) .and. (ibmprev .le. 253)) isprevbmap = .true.
     endif
     len = len + ilen
     
     ! Exit loop if last section reached
     if (len .eq. lencurr) exit
     
     ! If byte count for each section does not match current
     ! total length, then there is a problem.
     if (len .gt. lencurr) then
        print *, 'addfield: Section byte counts don''t add to total.'
        print *, 'addfield: Sum of section byte counts = ', len
        print *, 'addfield: Total byte count in Section 0 = ', lencurr
        ierr = 3
        return
     endif
  enddo
 
  ! Sections 4 through 7 can only be added after section 3 or 7.
  if ((isecnum .ne. 3) .and. (isecnum .ne. 7)) then
     print *, 'addfield: Sections 4-7 can only be added after  Section 3 or 7.'
     print *, 'addfield: Section ', isecnum, ' was the last found in', &
          ' given GRIB message.'
     ierr = 4
     return
 
     ! Sections 4 through 7 can only be added if section 3 was previously defined.
  elseif (.not. issec3) then
     print *, 'addfield: Sections 4-7 can only be added if Section', &
          ' 3 was previously included.'
     print *, 'addfield: Section 3 was not found in given GRIB message.'
     print *, 'addfield: Call to routine addgrid required', &
          ' to specify Grid definition.'
     ierr = 6
     return
  endif
 
  ! Add Section 4 - Product Definition Section.
  ibeg = lencurr * 8 ! Calculate offset for beginning of section 4
  iofst = ibeg + 32 ! leave space for length of section
  call g2_sbytec1(cgrib, four, iofst, 8) ! Store section number (4)
  iofst = iofst + 8
  call g2_sbytec1(cgrib, numcoord, iofst, 16) ! Store num of coordinate values
  iofst = iofst + 16
  call g2_sbytec1(cgrib, ipdsnum, iofst, 16) ! Store 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
 
  ! Extend the Product Definition Template, if necessary.
  ! 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, mappdslen, mappds)
  endif
 
  ! Pack up each input value in array ipdstmpl into 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).or.(ipdstmpl(i) .ge. 0)) then
        call g2_sbytec(cgrib, ipdstmpl(i), iofst, nbits)
     else
        call g2_sbytec1(cgrib, one, iofst, 1)
        call g2_sbytec1(cgrib, iabs(ipdstmpl(i)), iofst + 1, nbits - 1)
     endif
     iofst = iofst+nbits
  enddo
 
  ! Add Optional list of vertical coordinate values
  ! after the Product Definition Template, if necessary.
  if (numcoord .ne. 0) then
     do i = 1, numcoord
        coordlist_4(i) = real(coordlist(i), 4)
     end do
     call mkieee(coordlist_4, coordieee, numcoord)
     call g2_sbytescr(cgrib, coordieee, iofst, 32, 0, numcoord)
     iofst = iofst + (32 * numcoord)
  endif
 
  ! Calculate length of section 4 and store it in octets
  ! 1-4 of section 4.
  lensec4 = (iofst-ibeg)/8
  call g2_sbytec1(cgrib, lensec4, ibeg, 32)
 
  ! Pack Data using appropriate algorithm
 
  ! Get Data Representation Template
  call getdrstemplate(idrsnum, mapdrslen, mapdrs, needext, iret)
  if (iret .ne. 0) then
     ierr = 5
     return
  endif
 
  ! contract data field, removing data at invalid grid points,
  ! if bit-map is provided with field.
  if (ibmap .eq. 0 .OR. ibmap .eq. 254) then
     allocate(pfld(max(2, ngrdpts)))
     ndpts = 0;
     do jj = 1, ngrdpts
        intbmap(jj) = 0
        if (bmap(jj)) then
           intbmap(jj) = 1
           ndpts = ndpts + 1
           pfld(ndpts) = fld(jj);
        endif
     enddo
     if (ndpts == 0 .and. ngrdpts > 0) then
        pfld(1) = 0
     endif
  else
     ndpts = ngrdpts;
     pfld=>fld;
  endif
  lcpack = 0
  nsize = ndpts*4
  if (nsize .lt. minsize) nsize = minsize
  allocate(cpack(nsize), stat = istat)
  if (idrsnum.eq.0) then ! Simple Packing
     call simpack(pfld, ndpts, idrstmpl, cpack, lcpack)
  elseif (idrsnum.eq.2.or.idrsnum.eq.3) then ! Complex Packing
     call cmplxpack(pfld, ndpts, idrsnum, idrstmpl, cpack, lcpack)
  elseif (idrsnum.eq.50) then ! Sperical Harmonic Simple Packing
     call simpack(pfld(2), ndpts-1, idrstmpl, cpack, lcpack)
     tmpfld(1) = real(pfld(1), 4)
     call mkieee(tmpfld, tmpre00, 1) ! ensure RE(0,0) value is IEEE format
     re00 = tmpre00(1)
     ire00 = transfer(re00, ire00)
     idrstmpl(5) = ire00
  elseif (idrsnum.eq.51) then ! Sperical Harmonic Complex Packing
     call getpoly(cgrib(lpos3), lensec3, jj, kk, mm)
     if (jj.ne.0 .AND. kk.ne.0 .AND. mm.ne.0) then
        call specpack(pfld, ndpts, jj, kk, mm, idrstmpl, cpack, lcpack)
     else
        print *, 'addfield: Cannot pack DRT 5.51.'
        ierr=9
        return
     endif
 
  elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then ! JPEG2000 encoding
     if (ibmap.eq.255) then
        call getdim(cgrib(lpos3), lensec3, width, height, iscan)
        if (width.eq.0 .OR. height.eq.0) then
           width=ndpts
           height=1
        elseif (width.eq.allones .OR. height.eq.allones) then
           width=ndpts
           height=1
        elseif (ibits(iscan, 5, 1) .eq. 1) then ! Scanning mode: bit 3
           itemp=width
           width=height
           height=itemp
        endif
     else
        width=ndpts
        height=1
     endif
     if (width<1 .or. height<1) then
        ! Special case: bitmask off everywhere.
        write(0, *) 'Warning: bitmask off everywhere.'
        write(0, *) '   Pretend one point in jpcpack to avoid crash.'
        width=1
        height=1
     endif
     lcpack=nsize
     ! print *, 'w, h=', width, height
     call jpcpack(pfld, width, height, idrstmpl, cpack, lcpack)
 
  elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then ! PNG encoding
     if (ibmap.eq.255) then
        call getdim(cgrib(lpos3), lensec3, width, height, iscan)
        if (width.eq.0 .OR. height.eq.0) then
           width=ndpts
           height=1
        elseif (width.eq.allones .OR. height.eq.allones) then
           width=ndpts
           height=1
        elseif (ibits(iscan, 5, 1) .eq. 1) then ! Scanning mode: bit 3
           itemp=width
           width=height
           height=itemp
        endif
     else
        width=ndpts
        height=1
     endif
     !print *, 'png size ', width, height
     call pngpack(pfld, width, height, idrstmpl, cpack, lcpack)
     !print *, 'png packed'
#ifdef USE_AEC
  elseif (idrsnum.eq.42) then ! AEC compression
    if (ibmap.eq.255) then
      call getdim(cgrib(lpos3), lensec3, width, height, iscan)
      if (width.eq.0 .OR. height.eq.0) then
          width=ndpts
          height=1
      elseif (width.eq.allones .OR. height.eq.allones) then
          width=ndpts
          height=1
      elseif (ibits(iscan, 5, 1) .eq. 1) then ! Scanning mode: bit 3
          itemp=width
          width=height
          height=itemp
      endif
    else
      width=ndpts
      height=1
    endif
    call aecpack(pfld, width, height, idrstmpl, cpack, lcpack);
#endif /* USE_AEC */
  else
     print *, 'addfield: Data Representation Template 5.', idrsnum, &
          ' not yet implemented.'
     ierr=7
     return
  endif
  if (ibmap.eq.0 .OR. ibmap.eq.254) then
     deallocate(pfld)
  endif
  if (lcpack .lt. 0) then
     if (allocated(cpack))deallocate(cpack)
     ierr=10
     return
  endif
 
  !     Add Section 5 - Data Representation Section
  ibeg=iofst ! Calculate offset for beginning of section 5
  iofst=ibeg + 32 ! leave space for length of section
  call g2_sbytec1(cgrib, five, iofst, 8) ! Store section number (5)
  iofst=iofst + 8
  call g2_sbytec1(cgrib, ndpts, iofst, 32) ! Store num of actual data points
  iofst=iofst + 32
  call g2_sbytec1(cgrib, idrsnum, iofst, 16) ! Store Data Repr. Template num.
  iofst = iofst + 16
 
  ! Pack up each input value in array idrstmpl into the
  ! the appropriate number of octets, which are specified in
  ! corresponding entries in array mapdrs.
  do i=1, mapdrslen
     nbits =iabs(mapdrs(i)) * 8
     if ((mapdrs(i) .ge. 0) .or. (idrstmpl(i) .ge. 0)) then
        call g2_sbytec(cgrib, idrstmpl(i), iofst, nbits)
     else
        call g2_sbytec1(cgrib, one, iofst, 1)
        call g2_sbytec1(cgrib, iabs(idrstmpl(i)), iofst+1, nbits-1)
     endif
     iofst = iofst + nbits
  enddo
 
  ! Calculate length of section 5 and store it in octets
  ! 1-4 of section 5.
  lensec5 = (iofst - ibeg) / 8
  call g2_sbytec1(cgrib, lensec5, ibeg, 32)
 
  ! Add Section 6 - Bit-Map Section.
  ibeg = iofst ! Calculate offset for beginning of section 6
  iofst = ibeg + 32 ! leave space for length of section
  call g2_sbytec1(cgrib, six, iofst, 8) ! Store section number (6)
  iofst=iofst + 8
  call g2_sbytec1(cgrib, ibmap, iofst, 8) ! Store Bit Map indicator
  iofst = iofst + 8
 
  ! Store bitmap, if supplied
  if (ibmap.eq.0) then
     call g2_sbytesc(cgrib, intbmap, iofst, 1, 0, ngrdpts) ! Store BitMap
     iofst=iofst+ngrdpts
  endif
 
  ! If specifying a previously defined bit-map, make sure
  ! one already exists in the current GRIB message.
  if ((ibmap.eq.254).and.(.not.isprevbmap)) then
     print *, 'addfield: Requested previously defined bitmap, ', &
          ' but one does not exist in the current GRIB message.'
     ierr=8
     return
  endif
 
  ! Calculate length of section 6 and store it in octets
  ! 1-4 of section 6. Pad to end of octect, if necessary.
  left = 8 - mod(iofst, 8)
  if (left .ne. 8) then
     call g2_sbytec1(cgrib, zero, iofst, left) ! Pad with zeros to fill Octet
     iofst = iofst + left
  endif
  lensec6 = (iofst - ibeg) / 8
  call g2_sbytec1(cgrib, lensec6, ibeg, 32)
 
  ! Add Section 7 - Data Section.
  ibeg = iofst ! Calculate offset for beginning of section 7
  iofst = ibeg + 32 ! leave space for length of section
  call g2_sbytec1(cgrib, seven, iofst, 8) ! Store section number (7)
  iofst = iofst + 8
 
  ! Store Packed Binary Data values, if non-constant field
  if (lcpack .ne. 0) then
     ioctet = iofst / 8
     cgrib(ioctet + 1:ioctet + lcpack) = cpack(1:lcpack)
     iofst = iofst + (8 * lcpack)
  endif
 
  ! Calculate length of section 7 and store it in octets
  ! 1-4 of section 7.
  lensec7 = (iofst - ibeg) / 8
  call g2_sbytec1(cgrib, lensec7, ibeg, 32)
 
  if (allocated(cpack)) deallocate(cpack)
 
  ! Update current byte total of message in Section 0.
  newlen = lencurr + lensec4 + lensec5 + lensec6 + lensec7
  call g2_sbytec1(cgrib, newlen, 96, 32)
 
  return
end subroutine addfield
 
subroutine addgrid(cgrib, lcgrib, igds, igdstmpl, igdstmplen, &
     ideflist, idefnum, ierr)
  use gridtemplates
  implicit none
 
  character(len = 1), intent(inout) :: cgrib(lcgrib)
  integer, intent(in) :: igds(*), igdstmpl(*), ideflist(idefnum)
  integer, intent(in) :: lcgrib, idefnum, igdstmplen
  integer, intent(out) :: ierr
 
  character(len = 4), parameter :: grib = 'GRIB', c7777 = '7777'
  character(len = 4):: ctemp
  integer:: mapgrid(igdstmplen)
  integer, parameter :: ONE = 1, three = 3
  integer lensec3, iofst, ibeg, lencurr, len, mapgridlen
  logical needext
  integer :: i, ilen, iret, isecnum, nbits
 
  interface
     subroutine g2_gbytec(in, iout, iskip, nbits)
       character*1, intent(in) :: in(*)
       integer, intent(inout) :: iout(*)
       integer, intent(in) :: iskip, nbits
     end subroutine g2_gbytec
     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_gbytec81(in, siout, iskip, nbits)
       character*1, intent(in) :: in(*)
       integer (kind = 8), intent(inout) :: siout
       integer, intent(in) :: iskip, nbits
       integer (kind = 8) :: iout(1)
     end subroutine g2_gbytec81
     subroutine g2_sbytec(out, in, iskip, nbits)
       character*1, intent(inout) :: out(*)
       integer, intent(in) :: in(*)
       integer, intent(in) :: iskip, nbits
     end subroutine g2_sbytec
     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
 
  ierr = 0
 
#ifdef LOGGING
  print *, 'addgrid lcgrib ', lcgrib, ' igdstmplen ', igdstmplen, ' idefnum ', idefnum
#endif
 
  ! Check to see if beginning of GRIB message exists.
  do i = 1, 4
     if (cgrib(i) /= grib(i:i)) then
        print *, 'addgrid: GRIB not found in given message.'
        print *, 'addgrid: Call to routine gribcreate required', &
             ' to initialize GRIB messge.'
10      format('"', 4a1, '" /= "GRIB"')
        print 10, cgrib(1:4)
        ierr = 1
        return
     endif
  enddo
 
  ! Get current length of GRIB message.
  call g2_gbytec1(cgrib, lencurr, 96, 32)
 
  ! Check to see if GRIB message is already complete.
  ctemp = cgrib(lencurr - 3) // cgrib(lencurr - 2) // cgrib(lencurr &
       - 1) // cgrib(lencurr)
  if (ctemp .eq. c7777) then
     print *, 'addgrid: GRIB message already complete.  Cannot', &
          ' add new section.'
     ierr = 2
     return
  endif
 
  ! Loop through all current sections of the GRIB message to find the
  ! last section number.
  len = 16                  ! length of Section 0
  do
     ! Get length and section number of next section.
     iofst = len * 8
     call g2_gbytec1(cgrib, ilen, iofst, 32)
     iofst = iofst + 32
     call g2_gbytec1(cgrib, isecnum, iofst, 8)
     len = len + ilen
 
     ! Exit loop if last section reached.
     if (len .eq. lencurr) exit
 
     ! If byte count for each section doesn't match current
     ! total length, then there is a problem.
     if (len .gt. lencurr) then
        print *, 'addgrid: Section byte counts don''t add to total.'
        print *, 'addgrid: Sum of section byte counts = ', len
        print *, 'addgrid: Total byte count in Section 0 = ', lencurr
        ierr = 3
        return
     endif
  enddo
 
  ! Section 3 can only be added after sections 1, 2 and 7.
  if ((isecnum .ne. 1) .and. (isecnum .ne. 2) .and. &
       (isecnum .ne. 7)) then
     print *, 'addgrid: Section 3 can only be added after Section',  &
          ' 1, 2 or 7.'
     print *, 'addgrid: Section ', isecnum, &
          ' was the last found in given GRIB message.'
     ierr = 4
     return
  endif
 
  ! Add Section 3  - Grid Definition Section.
  ibeg = lencurr * 8        !   Calculate offset for beginning of section 3
  iofst = ibeg + 32         !   leave space for length of section
  call g2_sbytec1(cgrib, three, iofst, 8) ! Store section number (3)
  iofst = iofst + 8
  call g2_sbytec(cgrib, igds(1), iofst, 8) ! Store source of Grid def.
  iofst = iofst + 8
  call g2_sbytec(cgrib, igds(2), iofst, 32) ! Store number of data pts.
  iofst = iofst + 32
  call g2_sbytec(cgrib, igds(3), iofst, 8) ! Store number of extra octets.
  iofst = iofst + 8
  call g2_sbytec(cgrib, igds(4), iofst, 8) ! Store interp. of extra octets.
  iofst = iofst + 8
 
  ! If Octet 6 is not equal to zero, Grid Definition Template may not
  ! be supplied.
  if (igds(1) .eq. 0) then
     call g2_sbytec(cgrib, igds(5), iofst, 16) ! Store Grid Def Template num.
  else
     call g2_sbytec1(cgrib, 65535, iofst, 16) ! Store missing value as Grid Def Template num.
  endif
  iofst = iofst + 16
 
  ! Get Grid Definition Template.
  if (igds(1) .eq. 0) then
     call getgridtemplate(igds(5), mapgridlen, mapgrid, needext, &
          iret)
     if (iret .ne. 0) then
        ierr = 5
        return
     endif
 
     ! Extend the Grid Definition Template, if necessary. 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, mapgridlen, &
             mapgrid)
     endif
  else
     mapgridlen = 0
  endif
 
  ! Pack up each input value in array igdstmpl into the the
  ! appropriate number of octets, which are specified in corresponding
  ! entries in array mapgrid.
  do i = 1, mapgridlen
     nbits = iabs(mapgrid(i)) * 8
     if ((mapgrid(i) .ge. 0) .or. (igdstmpl(i) .ge. 0)) then
        call g2_sbytec(cgrib, igdstmpl(i), iofst, nbits)
     else
        call g2_sbytec1(cgrib, one, iofst, 1)
        call g2_sbytec1(cgrib, iabs(igdstmpl(i)), iofst + 1, nbits &
             - 1)
     endif
     iofst = iofst + nbits
  enddo
 
  ! If requested, insert optional list of numbers defining number of
  ! points in each row or column. This is used for non regular grids.
  if (igds(3) .ne. 0) then
     nbits = igds(3) * 8
     call g2_sbytesc(cgrib, ideflist, iofst, nbits, 0, idefnum)
     iofst = iofst + (nbits * idefnum)
  endif
 
  ! Calculate length of section 3 and store it in octets 1-4 of
  ! section 3.
  lensec3 = (iofst - ibeg) / 8
  call g2_sbytec1(cgrib, lensec3, ibeg, 32)
 
  ! Update current byte total of message in Section 0.
  call g2_sbytec1(cgrib, lencurr + lensec3, 96, 32)
end subroutine addgrid
 
subroutine addlocal(cgrib, lcgrib, csec2, lcsec2, ierr)
  implicit none
  
  character(len = 1), intent(inout) :: cgrib(lcgrib)
  character(len = 1), intent(in) :: csec2(lcsec2)
  integer, intent(in) :: lcgrib, lcsec2
  integer, intent(out) :: ierr
 
  character(len = 4), parameter :: grib = 'GRIB', c7777 = '7777'
  character(len = 4):: ctemp
  integer, parameter :: two = 2
  integer :: lensec2, iofst, ibeg, lencurr, len
  integer :: ilen, isecnum, istart
 
  ierr = 0
 
  ! Check to see if beginning of GRIB message exists.
  ctemp = cgrib(1) // cgrib(2) // cgrib(3) // cgrib(4)
  if (ctemp .ne. grib) then
     print *, 'addlocal: GRIB not found in given message.'
     print *, 'addlocal: Call to routine gribcreate required', &
          ' to initialize GRIB messge.'
     ierr = 1
     return
  endif
 
  ! Get current length of GRIB message.
  call g2_gbytec1(cgrib, lencurr, 96, 32)
 
  ! Check to see if GRIB message is already complete
  ctemp = cgrib(lencurr - 3) // cgrib(lencurr - 2) // cgrib(lencurr - 1) // cgrib(lencurr)
  if (ctemp .eq. c7777) then
     print *, 'addlocal: GRIB message already complete.  Cannot add new section.'
     ierr = 2
     return
  endif
 
  ! Loop through all current sections of the GRIB message to find the
  ! last section number.
  len = 16    ! length of Section 0
  do 
     ! Get section number and length of next section.
     iofst = len * 8
     call g2_gbytec1(cgrib, ilen, iofst, 32)
     iofst = iofst + 32
     call g2_gbytec1(cgrib, isecnum, iofst, 8)
     len = len + ilen
     ! Exit loop if last section reached
     if (len .eq. lencurr) exit
     
     ! If byte count for each section doesn't match current
     ! total length, then there is a problem.
     if (len .gt. lencurr) then
        print *, 'addlocal: Section byte counts don''t add to total.'
        print *, 'addlocal: Sum of section byte counts = ', len
        print *, 'addlocal: Total byte count in Section 0 = ', lencurr
        ierr = 3
        return
     endif
  enddo
 
  ! Section 2 can only be added after sections 1 and 7.
  if ((isecnum .ne. 1) .and. (isecnum .ne. 7)) then
     print *, 'addlocal: Section 2 can only be added after Section 1 or Section 7.'
     print *, 'addlocal: Section ', isecnum, ' was the last found in given GRIB message.'
     ierr = 4
     return
  endif
 
  ! Add Section 2  - Local Use Section.
  ibeg = lencurr * 8  !    Calculate offset for beginning of section 2
  iofst = ibeg + 32   !    leave space for length of section
  call g2_sbytec1(cgrib, two, iofst, 8)     ! Store section number (2)
  istart = lencurr + 5
  cgrib(istart + 1:istart + lcsec2) = csec2(1:lcsec2)
 
  ! Calculate length of section 2 and store it in octets 1-4 of
  ! section 2.
  lensec2 = lcsec2 + 5 !  bytes
  call g2_sbytec1(cgrib, lensec2, ibeg, 32)
 
  ! Update current byte total of message in Section 0.
  call g2_sbytec1(cgrib, lencurr + lensec2, 96, 32)
 
end subroutine addlocal
 
subroutine gribend(cgrib, lcgrib, lengrib, ierr)
  implicit none
 
  character(len = 1), intent(inout) :: cgrib(lcgrib)
  integer, intent(in) :: lcgrib
  integer, intent(out) :: lengrib, ierr
 
  integer ilen, isecnum
  character(len = 4), parameter :: grib = 'GRIB'
  character(len = 1), parameter :: c7777(4) = (/ '7', '7', '7', '7' /)
  character(len = 4):: ctemp
  integer iofst, lencurr, len
 
  ierr = 0
 
  ! Check to see if beginning of GRIB message exists.
  ctemp = cgrib(1) // cgrib(2) // cgrib(3) // cgrib(4)
  if (ctemp .ne. grib) then
     print *, 'gribend: GRIB not found in given message.'
     ierr = 1
     return
  endif
 
  ! Get current length of GRIB message.
  call g2_gbytec1(cgrib, lencurr, 96, 32)
 
  ! Loop through all current sections of the GRIB message to
  ! find the last section number.
  len = 16                    ! Length of Section 0
  do
     ! Get number and length of next section.
     iofst = len * 8
     call g2_gbytec1(cgrib, ilen, iofst, 32)
     iofst = iofst + 32
     call g2_gbytec1(cgrib, isecnum, iofst, 8)
     len = len + ilen
 
     ! Exit loop if last section reached.
     if (len .eq. lencurr) exit
 
     ! If byte count for each section doesn't match current total
     ! length, then there is a problem.
     if (len .gt. lencurr) then
        print *, 'gribend: Section byte counts don''t add ' &
             ,'to total.'
        print *, 'gribend: Sum of section byte counts = ', len
        print *, 'gribend: Total byte count in Section 0 = ', &
             lencurr
        ierr = 3
        return
     endif
  enddo
 
  ! Can only add End Section (Section 8) after Section 7.
  if (isecnum .ne. 7) then
     print *, 'gribend: Section 8 can only be added after Section 7.'
     print *, 'gribend: Section ', isecnum, &
          ' was the last found in',' given GRIB message.'
     ierr = 4
     return
  endif
 
  ! Add Section 8  - End Section.
  cgrib(lencurr + 1:lencurr + 4) = c7777
 
  ! Update current byte total of message in Section 0.
  lengrib = lencurr + 4
  call g2_sbytec1(cgrib, lengrib, 96, 32)
end subroutine gribend