addgrid.F90

Go to the documentation of this file.

 
 
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
 
  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_gbytec(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_gbytec(cgrib, ilen, iofst, 32)
     iofst = iofst + 32
     call g2_gbytec(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_sbytec(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_sbytec(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_sbytec(cgrib, one, iofst, 1)
        call g2_sbytec(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_sbytec(cgrib, lensec3, ibeg, 32)
 
 
  ! Update current byte total of message in Section 0.
  call g2_sbytec(cgrib, lencurr + lensec3, 96, 32)
end subroutine addgrid