external/io_grib2/g2lib/getfield.F

   1       subroutine getfield(cgrib,lcgrib,ifldnum,igds,igdstmpl,igdslen,
   2      &                    ideflist,idefnum,ipdsnum,ipdstmpl,ipdslen,
   3      &                    coordlist,numcoord,ndpts,idrsnum,idrstmpl,
   4      &                    idrslen,ibmap,bmap,fld,ierr)
   5 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   6 !                .      .    .                                       .
   7 ! SUBPROGRAM:    getfield
   8 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
   9 !
  10 ! ABSTRACT: This subroutine returns the Grid Definition, Product Definition,
  11 !   Bit-map ( if applicable ), and the unpacked data for a given data
  12 !   field.  Since there can be multiple data fields packed into a GRIB2
  13 !   message, the calling routine indicates which field is being requested
  14 !   with the ifldnum argument.
  15 !
  16 ! PROGRAM HISTORY LOG:
  17 ! 2000-05-26  Gilbert
  18 !
  19 ! USAGE:    CALL getfield(cgrib,lcgrib,ifldnum,igds,igdstmpl,igdslen,
  20 !    &                    ideflist,idefnum,ipdsnum,ipdstmpl,ipdslen,
  21 !    &                    coordlist,numcoord,ndpts,idrsnum,idrstmpl,
  22 !    &                    idrslen,ibmap,bmap,fld,ierr)
  23 !   INPUT ARGUMENT LIST:
  24 !     cgrib    - Character array that contains the GRIB2 message
  25 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
  26 !     ifldnum  - Specifies which field in the GRIB2 message to return.
  27 !
  28 !   OUTPUT ARGUMENT LIST:
  29 !     igds     - Contains information read from the appropriate GRIB Grid
  30 !                Definition Section 3 for the field being returned.
  31 !                Must be dimensioned >= 5.
  32 !                igds(1)=Source of grid definition (see Code Table 3.0)
  33 !                igds(2)=Number of grid points in the defined grid.
  34 !                igds(3)=Number of octets needed for each
  35 !                            additional grid points definition.
  36 !                            Used to define number of
  37 !                            points in each row ( or column ) for
  38 !                            non-regular grids.
  39 !                            = 0, if using regular grid.
  40 !                igds(4)=Interpretation of list for optional points
  41 !                            definition.  (Code Table 3.11)
  42 !                igds(5)=Grid Definition Template Number (Code Table 3.1)
  43 !     igdstmpl - Contains the data values for the specified Grid Definition
  44 !                Template ( NN=igds(5) ).  Each element of this integer
  45 !                array contains an entry (in the order specified) of Grid
  46 !                Defintion Template 3.NN
  47 !                A safe dimension for this array can be obtained in advance
  48 !                from maxvals(2), which is returned from subroutine gribinfo.
  49 !     igdslen  - Number of elements in igdstmpl().  i.e. number of entries
  50 !                in Grid Defintion Template 3.NN  ( NN=igds(5) ).
  51 !     ideflist - (Used if igds(3) .ne. 0)  This array contains the
  52 !                number of grid points contained in each row ( or column ).
  53 !                (part of Section 3)
  54 !                A safe dimension for this array can be obtained in advance
  55 !                from maxvals(3), which is returned from subroutine gribinfo.
  56 !     idefnum  - (Used if igds(3) .ne. 0)  The number of entries
  57 !                in array ideflist.  i.e. number of rows ( or columns )
  58 !                for which optional grid points are defined.
  59 !     ipdsnum  - Product Definition Template Number ( see Code Table 4.0)
  60 !     ipdstmpl - Contains the data values for the specified Product Definition
  61 !                Template ( N=ipdsnum ).  Each element of this integer
  62 !                array contains an entry (in the order specified) of Product
  63 !                Defintion Template 4.N
  64 !                A safe dimension for this array can be obtained in advance
  65 !                from maxvals(4), which is returned from subroutine gribinfo.
  66 !     ipdslen  - Number of elements in ipdstmpl().  i.e. number of entries
  67 !                in Product Defintion Template 4.N  ( N=ipdsnum ).
  68 !     coordlist- Array containg floating point values intended to document
  69 !                the vertical discretisation associated to model data
  70 !                on hybrid coordinate vertical levels.  (part of Section 4)
  71 !                The dimension of this array can be obtained in advance
  72 !                from maxvals(5), which is returned from subroutine gribinfo.
  73 !     numcoord - number of values in array coordlist.
  74 !     ndpts    - Number of data points unpacked and returned.
  75 !     idrsnum  - Data Representation Template Number ( see Code Table 5.0)
  76 !     idrstmpl - Contains the data values for the specified Data Representation
  77 !                Template ( N=idrsnum ).  Each element of this integer
  78 !                array contains an entry (in the order specified) of Product
  79 !                Defintion Template 5.N
  80 !                A safe dimension for this array can be obtained in advance
  81 !                from maxvals(6), which is returned from subroutine gribinfo.
  82 !     idrslen  - Number of elements in idrstmpl().  i.e. number of entries
  83 !                in Data Representation Template 5.N  ( N=idrsnum ).
  84 !     ibmap    - Bitmap indicator ( see Code Table 6.0 )
  85 !                0 = bitmap applies and is included in Section 6.
  86 !                1-253 = Predefined bitmap applies
  87 !                254 = Previously defined bitmap applies to this field
  88 !                255 = Bit map does not apply to this product.
  89 !     bmap()   - Logical*1 array containing decoded bitmap. ( if ibmap=0 )
  90 !                The dimension of this array can be obtained in advance
  91 !                from maxvals(7), which is returned from subroutine gribinfo.
  92 !     fld()    - Array of ndpts unpacked data points.
  93 !                A safe dimension for this array can be obtained in advance
  94 !                from maxvals(7), which is returned from subroutine gribinfo.
  95 !     ierr     - Error return code.
  96 !                0 = no error
  97 !                1 = Beginning characters "GRIB" not found.
  98 !                2 = GRIB message is not Edition 2.
  99 !                3 = The data field request number was not positive.
 100 !                4 = End string "7777" found, but not where expected.
 101 !                6 = GRIB message did not contain the requested number of
 102 !                    data fields.
 103 !                7 = End string "7777" not found at end of message.
 104 !                9 = Data Representation Template 5.NN not yet implemented.
 105 !               10 = Error unpacking Section 3.
 106 !               11 = Error unpacking Section 4.
 107 !               12 = Error unpacking Section 5.
 108 !               13 = Error unpacking Section 6.
 109 !               14 = Error unpacking Section 7.
 110 !
 111 ! REMARKS: Note that subroutine gribinfo can be used to first determine
 112 !          how many data fields exist in a given GRIB message.
 113 !
 114 ! ATTRIBUTES:
 115 !   LANGUAGE: Fortran 90
 116 !   MACHINE:  IBM SP
 117 !
 118 !$$$
 119
 120       character(len=1),intent(in) :: cgrib(lcgrib)
 121       integer,intent(in) :: lcgrib,ifldnum
 122       integer,intent(out) :: igds(*),igdstmpl(*),ideflist(*)
 123       integer,intent(out) :: ipdsnum,ipdstmpl(*)
 124       integer,intent(out) :: idrsnum,idrstmpl(*)
 125       integer,intent(out) :: ndpts,ibmap,idefnum,numcoord
 126       integer,intent(out) :: ierr
 127       logical*1,intent(out) :: bmap(*)
 128       real,intent(out) :: fld(*),coordlist(*)
 129
 130       character(len=4),parameter :: grib='GRIB',c7777='7777'
 131       character(len=4) :: ctemp
 132       integer:: listsec0(2)
 133       integer iofst,ibeg,istart
 134       integer(4) :: ieee
 135       logical have3,have4,have5,have6,have7
 136
 137       have3=.false.
 138       have4=.false.
 139       have5=.false.
 140       have6=.false.
 141       have7=.false.
 142       ierr=0
 143       numfld=0
 144 !
 145 !  Check for valid request number
 146 !
 147       if (ifldnum.le.0) then
 148         print *,'getfield: Request for field number must be positive.'
 149         ierr=3
 150         return
 151       endif
 152 !
 153 !  Check for beginning of GRIB message in the first 100 bytes
 154 !
 155       istart=0
 156       do j=1,100
 157         ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
 158         if (ctemp.eq.grib ) then
 159           istart=j
 160           exit
 161         endif
 162       enddo
 163       if (istart.eq.0) then
 164         print *,'getfield:  Beginning characters GRIB not found.'
 165         ierr=1
 166         return
 167       endif
 168 !
 169 !  Unpack Section 0 - Indicator Section
 170 !
 171       iofst=8*(istart+5)
 172       call g2lib_gbyte(cgrib,listsec0(1),iofst,8)     ! Discipline
 173       iofst=iofst+8
 174       call g2lib_gbyte(cgrib,listsec0(2),iofst,8)     ! GRIB edition number
 175       iofst=iofst+8
 176       iofst=iofst+32
 177       call g2lib_gbyte(cgrib,lengrib,iofst,32)        ! Length of GRIB message
 178       iofst=iofst+32
 179       lensec0=16
 180       ipos=istart+lensec0
 181 !
 182 !  Currently handles only GRIB Edition 2.
 183 !
 184       if (listsec0(2).ne.2) then
 185         print *,'getfield: can only decode GRIB edition 2.'
 186         ierr=2
 187         return
 188       endif
 189 !
 190 !  Loop through the remaining sections keeping track of the
 191 !  length of each.  Also keep the latest Grid Definition Section info.
 192 !  Unpack the requested field number.
 193 !
 194       do
 195         !    Check to see if we are at end of GRIB message
 196         ctemp=cgrib(ipos)//cgrib(ipos+1)//cgrib(ipos+2)//cgrib(ipos+3)
 197         if (ctemp.eq.c7777 ) then
 198           ipos=ipos+4
 199           !    If end of GRIB message not where expected, issue error
 200           if (ipos.ne.(istart+lengrib)) then
 201             print *,'getfield: "7777" found, but not where expected.'
 202             ierr=4
 203             return
 204           endif
 205           exit
 206         endif
 207         !     Get length of Section and Section number
 208         iofst=(ipos-1)*8
 209         call g2lib_gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 210         iofst=iofst+32
 211         call g2lib_gbyte(cgrib,isecnum,iofst,8)         ! Get Section number
 212         iofst=iofst+8
 213         !print *,' lensec= ',lensec,'    secnum= ',isecnum
 214         !
 215         !   If found Section 3, unpack the GDS info using the
 216         !   appropriate template.  Save in case this is the latest
 217         !   grid before the requested field.
 218         !
 219         if (isecnum.eq.3) then
 220           iofst=iofst-40       ! reset offset to beginning of section
 221           call unpack3(cgrib,lcgrib,iofst,igds,igdstmpl,igdslen,
 222      &                 ideflist,idefnum,jerr)
 223           if (jerr.eq.0) then
 224             have3=.true.
 225           else
 226             ierr=10
 227             return
 228           endif
 229         endif
 230         !
 231         !   If found Section 4, check to see if this field is the
 232         !   one requested.
 233         !
 234         if (isecnum.eq.4) then
 235           numfld=numfld+1
 236           if (numfld.eq.ifldnum) then
 237             iofst=iofst-40       ! reset offset to beginning of section
 238             call unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,ipdslen,
 239      &                   coordlist,numcoord,jerr)
 240             if (jerr.eq.0) then
 241               have4=.true.
 242             else
 243               ierr=11
 244               return
 245             endif
 246           endif
 247         endif
 248         !
 249         !   If found Section 5, check to see if this field is the
 250         !   one requested.
 251         !
 252         if ((isecnum.eq.5).and.(numfld.eq.ifldnum)) then
 253           iofst=iofst-40       ! reset offset to beginning of section
 254           call unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
 255      &                 idrslen,jerr)
 256           if (jerr.eq.0) then
 257             have5=.true.
 258           else
 259             ierr=12
 260             return
 261           endif
 262         endif
 263         !
 264         !   If found Section 6, Unpack bitmap.
 265         !   Save in case this is the latest
 266         !   bitmap before the requested field.
 267         !
 268         if (isecnum.eq.6) then
 269           iofst=iofst-40       ! reset offset to beginning of section
 270           call unpack6(cgrib,lcgrib,iofst,igds(2),ibmap,bmap,jerr)
 271           if (jerr.eq.0) then
 272             have6=.true.
 273           else
 274             ierr=13
 275             return
 276           endif
 277         endif
 278         !
 279         !   If found Section 7, check to see if this field is the
 280         !   one requested.
 281         !
 282         if ((isecnum.eq.7).and.(numfld.eq.ifldnum)) then
 283           if (idrsnum.eq.0) then
 284             call simunpack(cgrib(ipos+5),lensec-6,idrstmpl,ndpts,fld)
 285             have7=.true.
 286           elseif (idrsnum.eq.2.or.idrsnum.eq.3) then
 287             call comunpack(cgrib(ipos+5),lensec-6,lensec,idrsnum,
 288      &                     idrstmpl,ndpts,fld,ier)
 289             if ( ier .ne. 0 ) then
 290                 ierr=14
 291                 return
 292             endif
 293             have7=.true.
 294           elseif (idrsnum.eq.50) then
 295             call simunpack(cgrib(ipos+5),lensec-6,idrstmpl,ndpts-1,
 296      &                     fld(2))
 297             ieee=idrstmpl(5)
 298             call rdieee(ieee,fld(1),1)
 299             have7=.true.
 300           else
 301             print *,'getfield: Data Representation Template ',idrsnum,
 302      &              ' not yet implemented.'
 303             ierr=9
 304             return
 305           endif
 306         endif
 307         !
 308         !   Check to see if we read pass the end of the GRIB
 309         !   message and missed the terminator string '7777'.
 310         !
 311         ipos=ipos+lensec                 ! Update beginning of section pointer
 312         if (ipos.gt.(istart+lengrib)) then
 313           print *,'getfield: "7777"  not found at end of GRIB message.'
 314           ierr=7
 315           return
 316         endif
 317
 318         if (have3.and.have4.and.have5.and.have6.and.have7) return
 319
 320       enddo
 321
 322 !
 323 !  If exited from above loop, the end of the GRIB message was reached
 324 !  before the requested field was found.
 325 !
 326       print *,'getfield: GRIB message contained ',numlocal,
 327      &        ' different fields.'
 328       print *,'getfield: The request was for the ',ifldnum,
 329      &        ' field.'
 330       ierr=6
 331
 332       return
 333       end
 334
 335
 336       subroutine unpack3(cgrib,lcgrib,iofst,igds,igdstmpl,
 337      &                   mapgridlen,ideflist,idefnum,ierr)
 338 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 339 !                .      .    .                                       .
 340 ! SUBPROGRAM:    unpack3
 341 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 342 !
 343 ! ABSTRACT: This subroutine unpacks Section 3 (Grid Definition Section)
 344 !   starting at octet 6 of that Section.
 345 !
 346 ! PROGRAM HISTORY LOG:
 347 ! 2000-05-26  Gilbert
 348 !
 349 ! USAGE:    CALL unpack3(cgrib,lcgrib,lensec,iofst,igds,igdstmpl,
 350 !    &                   mapgridlen,ideflist,idefnum,ierr)
 351 !   INPUT ARGUMENT LIST:
 352 !     cgrib    - Character array that contains the GRIB2 message
 353 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 354 !     iofst    - Bit offset of the beginning of Section 3.
 355 !
 356 !   OUTPUT ARGUMENT LIST:
 357 !     iofst    - Bit offset at the end of Section 3, returned.
 358 !     igds     - Contains information read from the appropriate GRIB Grid
 359 !                Definition Section 3 for the field being returned.
 360 !                Must be dimensioned >= 5.
 361 !                igds(1)=Source of grid definition (see Code Table 3.0)
 362 !                igds(2)=Number of grid points in the defined grid.
 363 !                igds(3)=Number of octets needed for each
 364 !                            additional grid points definition.
 365 !                            Used to define number of
 366 !                            points in each row ( or column ) for
 367 !                            non-regular grids.
 368 !                            = 0, if using regular grid.
 369 !                igds(4)=Interpretation of list for optional points
 370 !                            definition.  (Code Table 3.11)
 371 !                igds(5)=Grid Definition Template Number (Code Table 3.1)
 372 !     igdstmpl - Contains the data values for the specified Grid Definition
 373 !                Template ( NN=igds(5) ).  Each element of this integer
 374 !                array contains an entry (in the order specified) of Grid
 375 !                Defintion Template 3.NN
 376 !     mapgridlen- Number of elements in igdstmpl().  i.e. number of entries
 377 !                in Grid Defintion Template 3.NN  ( NN=igds(5) ).
 378 !     ideflist - (Used if igds(3) .ne. 0)  This array contains the
 379 !                number of grid points contained in each row ( or column ).
 380 !                (part of Section 3)
 381 !     idefnum  - (Used if igds(3) .ne. 0)  The number of entries
 382 !                in array ideflist.  i.e. number of rows ( or columns )
 383 !                for which optional grid points are defined.
 384 !     ierr     - Error return code.
 385 !                0 = no error
 386 !                5 = "GRIB" message contains an undefined Grid Definition
 387 !                    Template.
 388 !
 389 ! REMARKS: Uses Fortran 90 module gridtemplates.
 390 !
 391 ! ATTRIBUTES:
 392 !   LANGUAGE: Fortran 90
 393 !   MACHINE:  IBM SP
 394 !
 395 !$$$
 396
 397       use gridtemplates
 398
 399       character(len=1),intent(in) :: cgrib(lcgrib)
 400       integer,intent(in) :: lcgrib
 401       integer,intent(inout) :: iofst
 402       integer,intent(out) :: igds(*),igdstmpl(*),ideflist(*)
 403       integer,intent(out) :: ierr,idefnum
 404
 405       integer,allocatable :: mapgrid(:)
 406       integer :: mapgridlen,ibyttem
 407       logical needext
 408
 409       ierr=0
 410
 411       call g2lib_gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 412       iofst=iofst+32
 413       iofst=iofst+8     ! skip section number
 414
 415       call g2lib_gbyte(cgrib,igds(1),iofst,8)     ! Get source of Grid def.
 416       iofst=iofst+8
 417       call g2lib_gbyte(cgrib,igds(2),iofst,32)    ! Get number of grid pts.
 418       iofst=iofst+32
 419       call g2lib_gbyte(cgrib,igds(3),iofst,8)     ! Get num octets for opt. list
 420       iofst=iofst+8
 421       call g2lib_gbyte(cgrib,igds(4),iofst,8)     ! Get interpret. for opt. list
 422       iofst=iofst+8
 423       call g2lib_gbyte(cgrib,igds(5),iofst,16)    ! Get Grid Def Template num.
 424       iofst=iofst+16
 425       if (igds(1).eq.0) then
 426 !      if (igds(1).eq.0.OR.igds(1).eq.255) then  ! FOR ECMWF TEST ONLY
 427         allocate(mapgrid(lensec))
 428         !   Get Grid Definition Template
 429         call getgridtemplate(igds(5),mapgridlen,mapgrid,needext,
 430      &                       iret)
 431         if (iret.ne.0) then
 432           ierr=5
 433           return
 434         endif
 435       else
 436 !        igdstmpl=-1
 437         mapgridlen=0
 438         needext=.false.
 439       endif
 440       !
 441       !   Unpack each value into array igdstmpl from the
 442       !   the appropriate number of octets, which are specified in
 443       !   corresponding entries in array mapgrid.
 444       !
 445       ibyttem=0
 446       do i=1,mapgridlen
 447         nbits=iabs(mapgrid(i))*8
 448         if ( mapgrid(i).ge.0 ) then
 449           call g2lib_gbyte(cgrib,igdstmpl(i),iofst,nbits)
 450         else
 451           call g2lib_gbyte(cgrib,isign,iofst,1)
 452           call g2lib_gbyte(cgrib,igdstmpl(i),iofst+1,nbits-1)
 453           if (isign.eq.1) igdstmpl(i)=-igdstmpl(i)
 454         endif
 455         iofst=iofst+nbits
 456         ibyttem=ibyttem+iabs(mapgrid(i))
 457       enddo
 458       !
 459       !   Check to see if the Grid Definition Template needs to be
 460       !   extended.
 461       !   The number of values in a specific template may vary
 462       !   depending on data specified in the "static" part of the
 463       !   template.
 464       !
 465       if ( needext ) then
 466         call extgridtemplate(igds(5),igdstmpl,newmapgridlen,mapgrid)
 467         !   Unpack the rest of the Grid Definition Template
 468         do i=mapgridlen+1,newmapgridlen
 469           nbits=iabs(mapgrid(i))*8
 470           if ( mapgrid(i).ge.0 ) then
 471             call g2lib_gbyte(cgrib,igdstmpl(i),iofst,nbits)
 472           else
 473             call g2lib_gbyte(cgrib,isign,iofst,1)
 474             call g2lib_gbyte(cgrib,igdstmpl(i),iofst+1,nbits-1)
 475             if (isign.eq.1) igdstmpl(i)=-igdstmpl(i)
 476           endif
 477           iofst=iofst+nbits
 478           ibyttem=ibyttem+iabs(mapgrid(i))
 479         enddo
 480         mapgridlen=newmapgridlen
 481       endif
 482       !
 483       !   Unpack optional list of numbers defining number of points
 484       !   in each row or column, if included.  This is used for non regular
 485       !   grids.
 486       !
 487       if ( igds(3).ne.0 ) then
 488          nbits=igds(3)*8
 489          idefnum=(lensec-14-ibyttem)/igds(3)
 490          call g2lib_gbytes(cgrib,ideflist,iofst,nbits,0,idefnum)
 491          iofst=iofst+(nbits*idefnum)
 492       else
 493          idefnum=0
 494       endif
 495       if( allocated(mapgrid) ) deallocate(mapgrid)
 496       return    ! End of Section 3 processing
 497       end
 498
 499
 500       subroutine unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,mappdslen,
 501      &                   coordlist,numcoord,ierr)
 502 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 503 !                .      .    .                                       .
 504 ! SUBPROGRAM:    unpack4
 505 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 506 !
 507 ! ABSTRACT: This subroutine unpacks Section 4 (Product Definition Section)
 508 !   starting at octet 6 of that Section.
 509 !
 510 ! PROGRAM HISTORY LOG:
 511 ! 2000-05-26  Gilbert
 512 !
 513 ! USAGE:    CALL unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,mappdslen,
 514 !    &                   coordlist,numcoord,ierr)
 515 !   INPUT ARGUMENT LIST:
 516 !     cgrib    - Character array that contains the GRIB2 message
 517 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 518 !     iofst    - Bit offset of the beginning of Section 4.
 519 !
 520 !   OUTPUT ARGUMENT LIST:
 521 !     iofst    - Bit offset of the end of Section 4, returned.
 522 !     ipdsnum  - Product Definition Template Number ( see Code Table 4.0)
 523 !     ipdstmpl - Contains the data values for the specified Product Definition
 524 !                Template ( N=ipdsnum ).  Each element of this integer
 525 !                array contains an entry (in the order specified) of Product
 526 !                Defintion Template 4.N
 527 !     mappdslen- Number of elements in ipdstmpl().  i.e. number of entries
 528 !                in Product Defintion Template 4.N  ( N=ipdsnum ).
 529 !     coordlist- Array containg floating point values intended to document
 530 !                the vertical discretisation associated to model data
 531 !                on hybrid coordinate vertical levels.  (part of Section 4)
 532 !     numcoord - number of values in array coordlist.
 533 !     ierr     - Error return code.
 534 !                0 = no error
 535 !                5 = "GRIB" message contains an undefined Product Definition
 536 !                    Template.
 537 !
 538 ! REMARKS: Uses Fortran 90 module pdstemplates.
 539 !
 540 ! ATTRIBUTES:
 541 !   LANGUAGE: Fortran 90
 542 !   MACHINE:  IBM SP
 543 !
 544 !$$$
 545
 546       use pdstemplates
 547
 548       character(len=1),intent(in) :: cgrib(lcgrib)
 549       integer,intent(in) :: lcgrib
 550       integer,intent(inout) :: iofst
 551       real,intent(out) :: coordlist(*)
 552       integer,intent(out) :: ipdsnum,ipdstmpl(*)
 553       integer,intent(out) :: ierr,numcoord
 554
 555       real(4),allocatable :: coordieee(:)
 556       integer,allocatable :: mappds(:)
 557       integer :: mappdslen
 558       logical needext
 559
 560       ierr=0
 561
 562       call g2lib_gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 563       iofst=iofst+32
 564       iofst=iofst+8     ! skip section number
 565       allocate(mappds(lensec))
 566
 567       call g2lib_gbyte(cgrib,numcoord,iofst,16)    ! Get num of coordinate values
 568       iofst=iofst+16
 569       call g2lib_gbyte(cgrib,ipdsnum,iofst,16)    ! Get Prod. Def Template num.
 570       iofst=iofst+16
 571       !   Get Product Definition Template
 572       call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret)
 573       if (iret.ne.0) then
 574         ierr=5
 575         return
 576       endif
 577       !
 578       !   Unpack each value into array ipdstmpl from the
 579       !   the appropriate number of octets, which are specified in
 580       !   corresponding entries in array mappds.
 581       !
 582       do i=1,mappdslen
 583         nbits=iabs(mappds(i))*8
 584         if ( mappds(i).ge.0 ) then
 585           call g2lib_gbyte(cgrib,ipdstmpl(i),iofst,nbits)
 586         else
 587           call g2lib_gbyte(cgrib,isign,iofst,1)
 588           call g2lib_gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
 589           if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
 590         endif
 591         iofst=iofst+nbits
 592       enddo
 593       !
 594       !   Check to see if the Product Definition Template needs to be
 595       !   extended.
 596       !   The number of values in a specific template may vary
 597       !   depending on data specified in the "static" part of the
 598       !   template.
 599       !
 600       if ( needext ) then
 601         call extpdstemplate(ipdsnum,ipdstmpl,newmappdslen,mappds)
 602         !   Unpack the rest of the Product Definition Template
 603         do i=mappdslen+1,newmappdslen
 604           nbits=iabs(mappds(i))*8
 605           if ( mappds(i).ge.0 ) then
 606             call g2lib_gbyte(cgrib,ipdstmpl(i),iofst,nbits)
 607           else
 608             call g2lib_gbyte(cgrib,isign,iofst,1)
 609             call g2lib_gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
 610             if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
 611           endif
 612           iofst=iofst+nbits
 613         enddo
 614         mappdslen=newmappdslen
 615       endif
 616       !
 617       !   Get Optional list of vertical coordinate values
 618       !   after the Product Definition Template, if necessary.
 619       !
 620       if ( numcoord .ne. 0 ) then
 621         allocate (coordieee(numcoord))
 622         call g2lib_gbytes(cgrib,coordieee,iofst,32,0,numcoord)
 623         call rdieee(coordieee,coordlist,numcoord)
 624         deallocate (coordieee)
 625         iofst=iofst+(32*numcoord)
 626       endif
 627       if( allocated(mappds) ) deallocate(mappds)
 628       return    ! End of Section 4 processing
 629       end
 630
 631
 632       subroutine unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
 633      &                   mapdrslen,ierr)
 634 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 635 !                .      .    .                                       .
 636 ! SUBPROGRAM:    unpack5
 637 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 638 !
 639 ! ABSTRACT: This subroutine unpacks Section 5 (Data Representation Section)
 640 !   starting at octet 6 of that Section.
 641 !
 642 ! PROGRAM HISTORY LOG:
 643 ! 2000-05-26  Gilbert
 644 !
 645 ! USAGE:    CALL unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
 646 !                        mapdrslen,ierr)
 647 !   INPUT ARGUMENT LIST:
 648 !     cgrib    - Character array that contains the GRIB2 message
 649 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 650 !     iofst    - Bit offset of the beginning of Section 5.
 651 !
 652 !   OUTPUT ARGUMENT LIST:
 653 !     iofst    - Bit offset at the end of Section 5, returned.
 654 !     ndpts    - Number of data points unpacked and returned.
 655 !     idrsnum  - Data Representation Template Number ( see Code Table 5.0)
 656 !     idrstmpl - Contains the data values for the specified Data Representation
 657 !                Template ( N=idrsnum ).  Each element of this integer
 658 !                array contains an entry (in the order specified) of Data
 659 !                Representation Template 5.N
 660 !     mapdrslen- Number of elements in idrstmpl().  i.e. number of entries
 661 !                in Data Representation Template 5.N  ( N=idrsnum ).
 662 !     ierr     - Error return code.
 663 !                0 = no error
 664 !                7 = "GRIB" message contains an undefined Data
 665 !                    Representation Template.
 666 !
 667 ! REMARKS: None
 668 !
 669 ! ATTRIBUTES:
 670 !   LANGUAGE: Fortran 90
 671 !   MACHINE:  IBM SP
 672 !
 673 !$$$
 674
 675       use drstemplates
 676
 677       character(len=1),intent(in) :: cgrib(lcgrib)
 678       integer,intent(in) :: lcgrib
 679       integer,intent(inout) :: iofst
 680       integer,intent(out) :: ndpts,idrsnum,idrstmpl(*)
 681       integer,intent(out) :: ierr
 682
 683 C      integer,allocatable :: mapdrs(:)
 684       integer,allocatable :: mapdrs(:)
 685       integer :: mapdrslen
 686       logical needext
 687
 688       ierr=0
 689
 690       call g2lib_gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 691       iofst=iofst+32
 692       iofst=iofst+8     ! skip section number
 693       allocate(mapdrs(lensec))
 694
 695       call g2lib_gbyte(cgrib,ndpts,iofst,32)    ! Get num of data points
 696       iofst=iofst+32
 697       call g2lib_gbyte(cgrib,idrsnum,iofst,16)     ! Get Data Rep Template Num.
 698       iofst=iofst+16
 699       !   Gen Data Representation Template
 700       call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret)
 701       if (iret.ne.0) then
 702         ierr=7
 703         return
 704       endif
 705       !
 706       !   Unpack each value into array ipdstmpl from the
 707       !   the appropriate number of octets, which are specified in
 708       !   corresponding entries in array mappds.
 709       !
 710       do i=1,mapdrslen
 711         nbits=iabs(mapdrs(i))*8
 712         if ( mapdrs(i).ge.0 ) then
 713           call g2lib_gbyte(cgrib,idrstmpl(i),iofst,nbits)
 714         else
 715           call g2lib_gbyte(cgrib,isign,iofst,1)
 716           call g2lib_gbyte(cgrib,idrstmpl(i),iofst+1,nbits-1)
 717           if (isign.eq.1) idrstmpl(i)=-idrstmpl(i)
 718         endif
 719         iofst=iofst+nbits
 720       enddo
 721       !
 722       !   Check to see if the Data Representation Template needs to be
 723       !   extended.
 724       !   The number of values in a specific template may vary
 725       !   depending on data specified in the "static" part of the
 726       !   template.
 727       !
 728       if ( needext ) then
 729         call extdrstemplate(idrsnum,idrstmpl,newmapdrslen,mapdrs)
 730         !   Unpack the rest of the Data Representation Template
 731         do i=mapdrslen+1,newmapdrslen
 732           nbits=iabs(mapdrs(i))*8
 733           if ( mapdrs(i).ge.0 ) then
 734             call g2lib_gbyte(cgrib,idrstmpl(i),iofst,nbits)
 735           else
 736             call g2lib_gbyte(cgrib,isign,iofst,1)
 737             call g2lib_gbyte(cgrib,idrstmpl(i),iofst+1,nbits-1)
 738             if (isign.eq.1) idrstmpl(i)=-idrstmpl(i)
 739           endif
 740           iofst=iofst+nbits
 741         enddo
 742         mapdrslen=newmapdrslen
 743       endif
 744       if( allocated(mapdrs) ) deallocate(mapdrs)
 745       return    ! End of Section 5 processing
 746       end
 747
 748
 749       subroutine unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
 750 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 751 !                .      .    .                                       .
 752 ! SUBPROGRAM:    unpack6
 753 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 754 !
 755 ! ABSTRACT: This subroutine unpacks Section 6 (Bit-Map Section)
 756 !   starting at octet 6 of that Section.
 757 !
 758 ! PROGRAM HISTORY LOG:
 759 ! 2000-05-26  Gilbert
 760 !
 761 ! USAGE:    CALL unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
 762 !   INPUT ARGUMENT LIST:
 763 !     cgrib    - Character array that contains the GRIB2 message
 764 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 765 !     iofst    - Bit offset of the beginning of Section 6.
 766 !     ngpts    - Number of grid points specified in the bit-map
 767 !
 768 !   OUTPUT ARGUMENT LIST:
 769 !     iofst    - Bit offset at the end of Section 6, returned.
 770 !     ibmap    - Bitmap indicator ( see Code Table 6.0 )
 771 !                0 = bitmap applies and is included in Section 6.
 772 !                1-253 = Predefined bitmap applies
 773 !                254 = Previously defined bitmap applies to this field
 774 !                255 = Bit map does not apply to this product.
 775 !     bmap()   - Logical*1 array containing decoded bitmap. ( if ibmap=0 )
 776 !     ierr     - Error return code.
 777 !                0 = no error
 778 !                4 = Unrecognized pre-defined bit-map.
 779 !
 780 ! REMARKS: None
 781 !
 782 ! ATTRIBUTES:
 783 !   LANGUAGE: Fortran 90
 784 !   MACHINE:  IBM SP
 785 !
 786 !$$$
 787
 788       character(len=1),intent(in) :: cgrib(lcgrib)
 789       integer,intent(in) :: lcgrib,ngpts
 790       integer,intent(inout) :: iofst
 791       integer,intent(out) :: ibmap
 792       integer,intent(out) :: ierr
 793       logical*1,intent(out) :: bmap(ngpts)
 794
 795       integer :: intbmap(ngpts)
 796
 797       ierr=0
 798
 799       iofst=iofst+32    ! skip Length of Section
 800       iofst=iofst+8     ! skip section number
 801
 802       call g2lib_gbyte(cgrib,ibmap,iofst,8)    ! Get bit-map indicator
 803       iofst=iofst+8
 804
 805       if (ibmap.eq.0) then               ! Unpack bitmap
 806         call g2lib_gbytes(cgrib,intbmap,iofst,1,0,ngpts)
 807         iofst=iofst+ngpts
 808         do j=1,ngpts
 809           bmap(j)=.true.
 810           if (intbmap(j).eq.0) bmap(j)=.false.
 811         enddo
 812       elseif (ibmap.eq.254) then               ! Use previous bitmap
 813         return
 814       elseif (ibmap.eq.255) then               ! No bitmap in message
 815         bmap(1:ngpts)=.true.
 816       else
 817         print *,'unpack6: Predefined bitmap ',ibmap,' not recognized.'
 818         ierr=4
 819       endif
 820
 821       return    ! End of Section 6 processing
 822       end
 823