ungrib/src/ngl/g2/addgrid.f

   1       subroutine addgrid(cgrib,lcgrib,igds,igdstmpl,igdstmplen,
   2      &                   ideflist,idefnum,ierr)
   3 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   4 !                .      .    .                                       .
   5 ! SUBPROGRAM:    addgrid
   6 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-01
   7 !
   8 ! ABSTRACT: This subroutine packs up a Grid Definition Section (Section 3)
   9 !   and adds it to a GRIB2 message.
  10 !   This routine is used with routines "gribcreate", "addlocal", "addfield",
  11 !   and "gribend" to create a complete GRIB2 message.  Subroutine
  12 !   gribcreate must be called first to initialize a new GRIB2 message.
  13 !
  14 ! PROGRAM HISTORY LOG:
  15 ! 2000-05-01  Gilbert
  16 !
  17 ! USAGE:    CALL addgrid(cgrib,lcgrib,igds,igdstmpl,igdstmplen,
  18 !                        ideflist,idefnum,ierr)
  19 !   INPUT ARGUMENT LIST:
  20 !     cgrib    - Character array to contain the GRIB2 message
  21 !     lcgrib   - Maximum length (bytes) of array cgrib.
  22 !     igds     - Contains information needed for GRIB Grid Definition Section 3.
  23 !                Must be dimensioned >= 5.
  24 !                igds(1)=Source of grid definition (see Code Table 3.0)
  25 !                igds(2)=Number of grid points in the defined grid.
  26 !                igds(3)=Number of octets needed for each
  27 !                            additional grid points definition.
  28 !                            Used to define number of
  29 !                            points in each row ( or column ) for
  30 !                            non-regular grids.
  31 !                            = 0, if using regular grid.
  32 !                igds(4)=Interpretation of list for optional points
  33 !                            definition.  (Code Table 3.11)
  34 !                igds(5)=Grid Definition Template Number (Code Table 3.1)
  35 !     igdstmpl - Contains the data values for the specified Grid Definition
  36 !                Template ( NN=igds(5) ).  Each element of this integer
  37 !                array contains an entry (in the order specified) of Grid
  38 !                Defintion Template 3.NN
  39 !   igdstmplen - Max dimension of igdstmpl()
  40 !     ideflist - (Used if igds(3) .ne. 0)  This array contains the
  41 !                number of grid points contained in each row ( or column )
  42 !      idefnum - (Used if igds(3) .ne. 0)  The number of entries
  43 !                in array ideflist.  i.e. number of rows ( or columns )
  44 !                for which optional grid points are defined.
  45 !
  46 !   OUTPUT ARGUMENT LIST:
  47 !     cgrib    - Character array to contain the GRIB2 message
  48 !     ierr     - Error return code.
  49 !                0 = no error
  50 !                1 = GRIB message was not initialized.  Need to call
  51 !                    routine gribcreate first.
  52 !                2 = GRIB message already complete.  Cannot add new section.
  53 !                3 = Sum of Section byte counts doesn't add to total byte count.
  54 !                4 = Previous Section was not 1, 2 or 7.
  55 !                5 = Could not find requested Grid Definition Template.
  56 !
  57 ! REMARKS: Note that the Local Use Section ( Section 2 ) can only follow
  58 !          Section 1 or Section 7 in a GRIB2 message.
  59 !
  60 ! ATTRIBUTES:
  61 !   LANGUAGE: Fortran 90
  62 !   MACHINE:  IBM SP
  63 !
  64 !$$$
  65
  66       use gridtemplates
  67
  68       character(len=1),intent(inout) :: cgrib(lcgrib)
  69       integer,intent(in) :: igds(*),igdstmpl(*),ideflist(idefnum)
  70       integer,intent(in) :: lcgrib,idefnum,igdstmplen
  71       integer,intent(out) :: ierr
  72
  73       character(len=4),parameter :: grib='GRIB',c7777='7777'
  74       character(len=4):: ctemp
  75       integer:: mapgrid(igdstmplen)
  76       integer,parameter :: one=1,three=3
  77       integer lensec3,iofst,ibeg,lencurr,len,mapgridlen
  78       logical needext
  79
  80       ierr=0
  81 !
  82 !  Check to see if beginning of GRIB message exists
  83 !
  84       ctemp=cgrib(1)//cgrib(2)//cgrib(3)//cgrib(4)
  85       if ( ctemp.ne.grib ) then
  86         print *,'addgrid: GRIB not found in given message.'
  87         print *,'addgrid: Call to routine gribcreate required',
  88      &          ' to initialize GRIB messge.'
  89         ierr=1
  90         return
  91       endif
  92 !
  93 !  Get current length of GRIB message
  94 !
  95       call gbyte(cgrib,lencurr,96,32)
  96 !
  97 !  Check to see if GRIB message is already complete
  98 !
  99       ctemp=cgrib(lencurr-3)//cgrib(lencurr-2)//cgrib(lencurr-1)
 100      &      //cgrib(lencurr)
 101       if ( ctemp.eq.c7777 ) then
 102         print *,'addgrid: GRIB message already complete.  Cannot',
 103      &          ' add new section.'
 104         ierr=2
 105         return
 106       endif
 107 !
 108 !  Loop through all current sections of the GRIB message to
 109 !  find the last section number.
 110 !
 111       len=16    ! length of Section 0
 112       do
 113       !    Get section number and length of next section
 114         iofst=len*8
 115         call gbyte(cgrib,ilen,iofst,32)
 116         iofst=iofst+32
 117         call gbyte(cgrib,isecnum,iofst,8)
 118         len=len+ilen
 119       !    Exit loop if last section reached
 120         if ( len.eq.lencurr ) exit
 121       !    If byte count for each section doesn't match current
 122       !    total length, then there is a problem.
 123         if ( len.gt.lencurr ) then
 124           print *,'addgrid: Section byte counts don''t add to total.'
 125           print *,'addgrid: Sum of section byte counts = ',len
 126           print *,'addgrid: Total byte count in Section 0 = ',lencurr
 127           ierr=3
 128           return
 129         endif
 130       enddo
 131 !
 132 !  Section 3 can only be added after sections 1, 2 and 7.
 133 !
 134       if ( (isecnum.ne.1) .and. (isecnum.ne.2) .and.
 135      &     (isecnum.ne.7) ) then
 136         print *,'addgrid: Section 3 can only be added after Section',
 137      &          ' 1, 2 or 7.'
 138         print *,'addgrid: Section ',isecnum,' was the last found in',
 139      &          ' given GRIB message.'
 140         ierr=4
 141         return
 142       endif
 143 !
 144 !  Add Section 3  - Grid Definition Section
 145 !
 146       ibeg=lencurr*8        !   Calculate offset for beginning of section 3
 147       iofst=ibeg+32         !   leave space for length of section
 148       call sbyte(cgrib,three,iofst,8)     ! Store section number ( 3 )
 149       iofst=iofst+8
 150       call sbyte(cgrib,igds(1),iofst,8)     ! Store source of Grid def.
 151       iofst=iofst+8
 152       call sbyte(cgrib,igds(2),iofst,32)    ! Store number of data pts.
 153       iofst=iofst+32
 154       call sbyte(cgrib,igds(3),iofst,8)     ! Store number of extra octets.
 155       iofst=iofst+8
 156       call sbyte(cgrib,igds(4),iofst,8)     ! Store interp. of extra octets.
 157       iofst=iofst+8
 158       !   if Octet 6 is not equal to zero, Grid Definition Template may
 159       !   not be supplied.
 160       if ( igds(1).eq.0 ) then
 161         call sbyte(cgrib,igds(5),iofst,16)  ! Store Grid Def Template num.
 162       else
 163         call sbyte(cgrib,65535,iofst,16)   ! Store missing value as Grid Def Template num.
 164       endif
 165       iofst=iofst+16
 166       !
 167       !   Get Grid Definition Template
 168       !
 169       if (igds(1).eq.0) then
 170         call getgridtemplate(igds(5),mapgridlen,mapgrid,needext,
 171      &                       iret)
 172         if (iret.ne.0) then
 173           ierr=5
 174           return
 175         endif
 176         !
 177         !   Extend the Grid Definition Template, if necessary.
 178         !   The number of values in a specific template may vary
 179         !   depending on data specified in the "static" part of the
 180         !   template.
 181         !
 182         if ( needext ) then
 183           call extgridtemplate(igds(5),igdstmpl,mapgridlen,mapgrid)
 184         endif
 185       else
 186         mapgridlen=0
 187       endif
 188       !
 189       !   Pack up each input value in array igdstmpl into the
 190       !   the appropriate number of octets, which are specified in
 191       !   corresponding entries in array mapgrid.
 192       !
 193       do i=1,mapgridlen
 194         nbits=iabs(mapgrid(i))*8
 195         if ( (mapgrid(i).ge.0).or.(igdstmpl(i).ge.0) ) then
 196           call sbyte(cgrib,igdstmpl(i),iofst,nbits)
 197         else
 198           call sbyte(cgrib,one,iofst,1)
 199           call sbyte(cgrib,iabs(igdstmpl(i)),iofst+1,nbits-1)
 200         endif
 201         iofst=iofst+nbits
 202       enddo
 203       !
 204       !   If requested,
 205       !   Insert optional list of numbers defining number of points
 206       !   in each row or column.  This is used for non regular
 207       !   grids.
 208       !
 209       if ( igds(3).ne.0 ) then
 210          nbits=igds(3)*8
 211          call sbytes(cgrib,ideflist,iofst,nbits,0,idefnum)
 212          iofst=iofst+(nbits*idefnum)
 213       endif
 214       !
 215       !   Calculate length of section 3 and store it in octets
 216       !   1-4 of section 3.
 217       !
 218       lensec3=(iofst-ibeg)/8
 219       call sbyte(cgrib,lensec3,ibeg,32)
 220
 221 !
 222 !  Update current byte total of message in Section 0
 223 !
 224       call sbyte(cgrib,lencurr+lensec3,96,32)
 225
 226       return
 227       end
 228