ungrib/src/ngl/g2/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 gbyte(cgrib,listsec0(1),iofst,8)     ! Discipline
 173       iofst=iofst+8
 174       call gbyte(cgrib,listsec0(2),iofst,8)     ! GRIB edition number
 175       iofst=iofst+8
 176       iofst=iofst+32
 177       call 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 gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 210         iofst=iofst+32
 211         call 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           elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then
 301             call jpcunpack(cgrib(ipos+5),lensec-5,idrstmpl,ndpts,fld)
 302             have7=.true.
 303           elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then
 304             call pngunpack(cgrib(ipos+5),lensec-5,idrstmpl,ndpts,fld)
 305             have7=.true.
 306           else
 307             print *,'getfield: Data Representation Template ',idrsnum,
 308      &              ' not yet implemented.'
 309             ierr=9
 310             return
 311           endif
 312         endif
 313         !
 314         !   Check to see if we read pass the end of the GRIB
 315         !   message and missed the terminator string '7777'.
 316         !
 317         ipos=ipos+lensec                 ! Update beginning of section pointer
 318         if (ipos.gt.(istart+lengrib)) then
 319           print *,'getfield: "7777"  not found at end of GRIB message.'
 320           ierr=7
 321           return
 322         endif
 323
 324         if (have3.and.have4.and.have5.and.have6.and.have7) return
 325
 326       enddo
 327
 328 !
 329 !  If exited from above loop, the end of the GRIB message was reached
 330 !  before the requested field was found.
 331 !
 332       print *,'getfield: GRIB message contained ',numlocal,
 333      &        ' different fields.'
 334       print *,'getfield: The request was for the ',ifldnum,
 335      &        ' field.'
 336       ierr=6
 337
 338       return
 339       end
 340
 341
 342       subroutine unpack3(cgrib,lcgrib,iofst,igds,igdstmpl,
 343      &                   mapgridlen,ideflist,idefnum,ierr)
 344 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 345 !                .      .    .                                       .
 346 ! SUBPROGRAM:    unpack3
 347 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 348 !
 349 ! ABSTRACT: This subroutine unpacks Section 3 (Grid Definition Section)
 350 !   starting at octet 6 of that Section.
 351 !
 352 ! PROGRAM HISTORY LOG:
 353 ! 2000-05-26  Gilbert
 354 !
 355 ! USAGE:    CALL unpack3(cgrib,lcgrib,lensec,iofst,igds,igdstmpl,
 356 !    &                   mapgridlen,ideflist,idefnum,ierr)
 357 !   INPUT ARGUMENT LIST:
 358 !     cgrib    - Character array that contains the GRIB2 message
 359 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 360 !     iofst    - Bit offset of the beginning of Section 3.
 361 !
 362 !   OUTPUT ARGUMENT LIST:
 363 !     iofst    - Bit offset at the end of Section 3, returned.
 364 !     igds     - Contains information read from the appropriate GRIB Grid
 365 !                Definition Section 3 for the field being returned.
 366 !                Must be dimensioned >= 5.
 367 !                igds(1)=Source of grid definition (see Code Table 3.0)
 368 !                igds(2)=Number of grid points in the defined grid.
 369 !                igds(3)=Number of octets needed for each
 370 !                            additional grid points definition.
 371 !                            Used to define number of
 372 !                            points in each row ( or column ) for
 373 !                            non-regular grids.
 374 !                            = 0, if using regular grid.
 375 !                igds(4)=Interpretation of list for optional points
 376 !                            definition.  (Code Table 3.11)
 377 !                igds(5)=Grid Definition Template Number (Code Table 3.1)
 378 !     igdstmpl - Contains the data values for the specified Grid Definition
 379 !                Template ( NN=igds(5) ).  Each element of this integer
 380 !                array contains an entry (in the order specified) of Grid
 381 !                Defintion Template 3.NN
 382 !     mapgridlen- Number of elements in igdstmpl().  i.e. number of entries
 383 !                in Grid Defintion Template 3.NN  ( NN=igds(5) ).
 384 !     ideflist - (Used if igds(3) .ne. 0)  This array contains the
 385 !                number of grid points contained in each row ( or column ).
 386 !                (part of Section 3)
 387 !     idefnum  - (Used if igds(3) .ne. 0)  The number of entries
 388 !                in array ideflist.  i.e. number of rows ( or columns )
 389 !                for which optional grid points are defined.
 390 !     ierr     - Error return code.
 391 !                0 = no error
 392 !                5 = "GRIB" message contains an undefined Grid Definition
 393 !                    Template.
 394 !
 395 ! REMARKS: Uses Fortran 90 module gridtemplates.
 396 !
 397 ! ATTRIBUTES:
 398 !   LANGUAGE: Fortran 90
 399 !   MACHINE:  IBM SP
 400 !
 401 !$$$
 402
 403       use gridtemplates
 404
 405       character(len=1),intent(in) :: cgrib(lcgrib)
 406       integer,intent(in) :: lcgrib
 407       integer,intent(inout) :: iofst
 408       integer,intent(out) :: igds(*),igdstmpl(*),ideflist(*)
 409       integer,intent(out) :: ierr,idefnum
 410
 411       integer,allocatable :: mapgrid(:)
 412       integer :: mapgridlen,ibyttem
 413       logical needext
 414
 415       ierr=0
 416
 417       call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 418       iofst=iofst+32
 419       iofst=iofst+8     ! skip section number
 420
 421       call gbyte(cgrib,igds(1),iofst,8)     ! Get source of Grid def.
 422       iofst=iofst+8
 423       call gbyte(cgrib,igds(2),iofst,32)    ! Get number of grid pts.
 424       iofst=iofst+32
 425       call gbyte(cgrib,igds(3),iofst,8)     ! Get num octets for opt. list
 426       iofst=iofst+8
 427       call gbyte(cgrib,igds(4),iofst,8)     ! Get interpret. for opt. list
 428       iofst=iofst+8
 429       call gbyte(cgrib,igds(5),iofst,16)    ! Get Grid Def Template num.
 430       iofst=iofst+16
 431       if (igds(1).eq.0) then
 432 !      if (igds(1).eq.0.OR.igds(1).eq.255) then  ! FOR ECMWF TEST ONLY
 433         allocate(mapgrid(lensec))
 434         !   Get Grid Definition Template
 435         call getgridtemplate(igds(5),mapgridlen,mapgrid,needext,
 436      &                       iret)
 437         if (iret.ne.0) then
 438           ierr=5
 439           return
 440         endif
 441       else
 442 !        igdstmpl=-1
 443         mapgridlen=0
 444         needext=.false.
 445       endif
 446       !
 447       !   Unpack each value into array igdstmpl from the
 448       !   the appropriate number of octets, which are specified in
 449       !   corresponding entries in array mapgrid.
 450       !
 451       ibyttem=0
 452       do i=1,mapgridlen
 453         nbits=iabs(mapgrid(i))*8
 454         if ( mapgrid(i).ge.0 ) then
 455           call gbyte(cgrib,igdstmpl(i),iofst,nbits)
 456         else
 457           call gbyte(cgrib,isign,iofst,1)
 458           call gbyte(cgrib,igdstmpl(i),iofst+1,nbits-1)
 459           if (isign.eq.1) igdstmpl(i)=-igdstmpl(i)
 460         endif
 461         iofst=iofst+nbits
 462         ibyttem=ibyttem+iabs(mapgrid(i))
 463       enddo
 464       !
 465       !   Check to see if the Grid Definition Template needs to be
 466       !   extended.
 467       !   The number of values in a specific template may vary
 468       !   depending on data specified in the "static" part of the
 469       !   template.
 470       !
 471       if ( needext ) then
 472         call extgridtemplate(igds(5),igdstmpl,newmapgridlen,mapgrid)
 473         !   Unpack the rest of the Grid Definition Template
 474         do i=mapgridlen+1,newmapgridlen
 475           nbits=iabs(mapgrid(i))*8
 476           if ( mapgrid(i).ge.0 ) then
 477             call gbyte(cgrib,igdstmpl(i),iofst,nbits)
 478           else
 479             call gbyte(cgrib,isign,iofst,1)
 480             call gbyte(cgrib,igdstmpl(i),iofst+1,nbits-1)
 481             if (isign.eq.1) igdstmpl(i)=-igdstmpl(i)
 482           endif
 483           iofst=iofst+nbits
 484           ibyttem=ibyttem+iabs(mapgrid(i))
 485         enddo
 486         mapgridlen=newmapgridlen
 487       endif
 488       !
 489       !   Unpack optional list of numbers defining number of points
 490       !   in each row or column, if included.  This is used for non regular
 491       !   grids.
 492       !
 493       if ( igds(3).ne.0 ) then
 494          nbits=igds(3)*8
 495          idefnum=(lensec-14-ibyttem)/igds(3)
 496          call gbytes(cgrib,ideflist,iofst,nbits,0,idefnum)
 497          iofst=iofst+(nbits*idefnum)
 498       else
 499          idefnum=0
 500       endif
 501       if( allocated(mapgrid) ) deallocate(mapgrid)
 502       return    ! End of Section 3 processing
 503       end
 504
 505
 506       subroutine unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,mappdslen,
 507      &                   coordlist,numcoord,ierr)
 508 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 509 !                .      .    .                                       .
 510 ! SUBPROGRAM:    unpack4
 511 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 512 !
 513 ! ABSTRACT: This subroutine unpacks Section 4 (Product Definition Section)
 514 !   starting at octet 6 of that Section.
 515 !
 516 ! PROGRAM HISTORY LOG:
 517 ! 2000-05-26  Gilbert
 518 !
 519 ! USAGE:    CALL unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,mappdslen,
 520 !    &                   coordlist,numcoord,ierr)
 521 !   INPUT ARGUMENT LIST:
 522 !     cgrib    - Character array that contains the GRIB2 message
 523 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 524 !     iofst    - Bit offset of the beginning of Section 4.
 525 !
 526 !   OUTPUT ARGUMENT LIST:
 527 !     iofst    - Bit offset of the end of Section 4, returned.
 528 !     ipdsnum  - Product Definition Template Number ( see Code Table 4.0)
 529 !     ipdstmpl - Contains the data values for the specified Product Definition
 530 !                Template ( N=ipdsnum ).  Each element of this integer
 531 !                array contains an entry (in the order specified) of Product
 532 !                Defintion Template 4.N
 533 !     mappdslen- Number of elements in ipdstmpl().  i.e. number of entries
 534 !                in Product Defintion Template 4.N  ( N=ipdsnum ).
 535 !     coordlist- Array containg floating point values intended to document
 536 !                the vertical discretisation associated to model data
 537 !                on hybrid coordinate vertical levels.  (part of Section 4)
 538 !     numcoord - number of values in array coordlist.
 539 !     ierr     - Error return code.
 540 !                0 = no error
 541 !                5 = "GRIB" message contains an undefined Product Definition
 542 !                    Template.
 543 !
 544 ! REMARKS: Uses Fortran 90 module pdstemplates.
 545 !
 546 ! ATTRIBUTES:
 547 !   LANGUAGE: Fortran 90
 548 !   MACHINE:  IBM SP
 549 !
 550 !$$$
 551
 552       use pdstemplates
 553
 554       character(len=1),intent(in) :: cgrib(lcgrib)
 555       integer,intent(in) :: lcgrib
 556       integer,intent(inout) :: iofst
 557       real,intent(out) :: coordlist(*)
 558       integer,intent(out) :: ipdsnum,ipdstmpl(*)
 559       integer,intent(out) :: ierr,numcoord
 560
 561       real(4),allocatable :: coordieee(:)
 562       integer,allocatable :: mappds(:)
 563       integer :: mappdslen
 564       logical needext
 565
 566       ierr=0
 567
 568       call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 569       iofst=iofst+32
 570       iofst=iofst+8     ! skip section number
 571       allocate(mappds(lensec))
 572
 573       call gbyte(cgrib,numcoord,iofst,16)    ! Get num of coordinate values
 574       iofst=iofst+16
 575       call gbyte(cgrib,ipdsnum,iofst,16)    ! Get Prod. Def Template num.
 576       iofst=iofst+16
 577       !   Get Product Definition Template
 578       call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret)
 579       if (iret.ne.0) then
 580         ierr=5
 581         return
 582       endif
 583       !
 584       !   Unpack each value into array ipdstmpl from the
 585       !   the appropriate number of octets, which are specified in
 586       !   corresponding entries in array mappds.
 587       !
 588       do i=1,mappdslen
 589         nbits=iabs(mappds(i))*8
 590         if ( mappds(i).ge.0 ) then
 591           call gbyte(cgrib,ipdstmpl(i),iofst,nbits)
 592         else
 593           call gbyte(cgrib,isign,iofst,1)
 594           call gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
 595           if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
 596         endif
 597         iofst=iofst+nbits
 598       enddo
 599       !
 600       !   Check to see if the Product Definition Template needs to be
 601       !   extended.
 602       !   The number of values in a specific template may vary
 603       !   depending on data specified in the "static" part of the
 604       !   template.
 605       !
 606       if ( needext ) then
 607         call extpdstemplate(ipdsnum,ipdstmpl,newmappdslen,mappds)
 608         !   Unpack the rest of the Product Definition Template
 609         do i=mappdslen+1,newmappdslen
 610           nbits=iabs(mappds(i))*8
 611           if ( mappds(i).ge.0 ) then
 612             call gbyte(cgrib,ipdstmpl(i),iofst,nbits)
 613           else
 614             call gbyte(cgrib,isign,iofst,1)
 615             call gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
 616             if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
 617           endif
 618           iofst=iofst+nbits
 619         enddo
 620         mappdslen=newmappdslen
 621       endif
 622       !
 623       !   Get Optional list of vertical coordinate values
 624       !   after the Product Definition Template, if necessary.
 625       !
 626       if ( numcoord .ne. 0 ) then
 627         allocate (coordieee(numcoord))
 628         call gbytes(cgrib,coordieee,iofst,32,0,numcoord)
 629         call rdieee(coordieee,coordlist,numcoord)
 630         deallocate (coordieee)
 631         iofst=iofst+(32*numcoord)
 632       endif
 633       if( allocated(mappds) ) deallocate(mappds)
 634       return    ! End of Section 4 processing
 635       end
 636
 637
 638       subroutine unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
 639      &                   mapdrslen,ierr)
 640 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 641 !                .      .    .                                       .
 642 ! SUBPROGRAM:    unpack5
 643 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 644 !
 645 ! ABSTRACT: This subroutine unpacks Section 5 (Data Representation Section)
 646 !   starting at octet 6 of that Section.
 647 !
 648 ! PROGRAM HISTORY LOG:
 649 ! 2000-05-26  Gilbert
 650 !
 651 ! USAGE:    CALL unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
 652 !                        mapdrslen,ierr)
 653 !   INPUT ARGUMENT LIST:
 654 !     cgrib    - Character array that contains the GRIB2 message
 655 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 656 !     iofst    - Bit offset of the beginning of Section 5.
 657 !
 658 !   OUTPUT ARGUMENT LIST:
 659 !     iofst    - Bit offset at the end of Section 5, returned.
 660 !     ndpts    - Number of data points unpacked and returned.
 661 !     idrsnum  - Data Representation Template Number ( see Code Table 5.0)
 662 !     idrstmpl - Contains the data values for the specified Data Representation
 663 !                Template ( N=idrsnum ).  Each element of this integer
 664 !                array contains an entry (in the order specified) of Data
 665 !                Representation Template 5.N
 666 !     mapdrslen- Number of elements in idrstmpl().  i.e. number of entries
 667 !                in Data Representation Template 5.N  ( N=idrsnum ).
 668 !     ierr     - Error return code.
 669 !                0 = no error
 670 !                7 = "GRIB" message contains an undefined Data
 671 !                    Representation Template.
 672 !
 673 ! REMARKS: None
 674 !
 675 ! ATTRIBUTES:
 676 !   LANGUAGE: Fortran 90
 677 !   MACHINE:  IBM SP
 678 !
 679 !$$$
 680
 681       use drstemplates
 682
 683       character(len=1),intent(in) :: cgrib(lcgrib)
 684       integer,intent(in) :: lcgrib
 685       integer,intent(inout) :: iofst
 686       integer,intent(out) :: ndpts,idrsnum,idrstmpl(*)
 687       integer,intent(out) :: ierr
 688
 689 C      integer,allocatable :: mapdrs(:)
 690       integer,allocatable :: mapdrs(:)
 691       integer :: mapdrslen
 692       logical needext
 693
 694       ierr=0
 695
 696       call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
 697       iofst=iofst+32
 698       iofst=iofst+8     ! skip section number
 699       allocate(mapdrs(lensec))
 700
 701       call gbyte(cgrib,ndpts,iofst,32)    ! Get num of data points
 702       iofst=iofst+32
 703       call gbyte(cgrib,idrsnum,iofst,16)     ! Get Data Rep Template Num.
 704       iofst=iofst+16
 705       !   Gen Data Representation Template
 706       call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret)
 707       if (iret.ne.0) then
 708         ierr=7
 709         return
 710       endif
 711       !
 712       !   Unpack each value into array ipdstmpl from the
 713       !   the appropriate number of octets, which are specified in
 714       !   corresponding entries in array mappds.
 715       !
 716       do i=1,mapdrslen
 717         nbits=iabs(mapdrs(i))*8
 718         if ( mapdrs(i).ge.0 ) then
 719           call gbyte(cgrib,idrstmpl(i),iofst,nbits)
 720         else
 721           call gbyte(cgrib,isign,iofst,1)
 722           call gbyte(cgrib,idrstmpl(i),iofst+1,nbits-1)
 723           if (isign.eq.1) idrstmpl(i)=-idrstmpl(i)
 724         endif
 725         iofst=iofst+nbits
 726       enddo
 727       !
 728       !   Check to see if the Data Representation Template needs to be
 729       !   extended.
 730       !   The number of values in a specific template may vary
 731       !   depending on data specified in the "static" part of the
 732       !   template.
 733       !
 734       if ( needext ) then
 735         call extdrstemplate(idrsnum,idrstmpl,newmapdrslen,mapdrs)
 736         !   Unpack the rest of the Data Representation Template
 737         do i=mapdrslen+1,newmapdrslen
 738           nbits=iabs(mapdrs(i))*8
 739           if ( mapdrs(i).ge.0 ) then
 740             call gbyte(cgrib,idrstmpl(i),iofst,nbits)
 741           else
 742             call gbyte(cgrib,isign,iofst,1)
 743             call gbyte(cgrib,idrstmpl(i),iofst+1,nbits-1)
 744             if (isign.eq.1) idrstmpl(i)=-idrstmpl(i)
 745           endif
 746           iofst=iofst+nbits
 747         enddo
 748         mapdrslen=newmapdrslen
 749       endif
 750       if( allocated(mapdrs) ) deallocate(mapdrs)
 751       return    ! End of Section 5 processing
 752       end
 753
 754
 755       subroutine unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
 756 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 757 !                .      .    .                                       .
 758 ! SUBPROGRAM:    unpack6
 759 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
 760 !
 761 ! ABSTRACT: This subroutine unpacks Section 6 (Bit-Map Section)
 762 !   starting at octet 6 of that Section.
 763 !
 764 ! PROGRAM HISTORY LOG:
 765 ! 2000-05-26  Gilbert
 766 !
 767 ! USAGE:    CALL unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
 768 !   INPUT ARGUMENT LIST:
 769 !     cgrib    - Character array that contains the GRIB2 message
 770 !     lcgrib   - Length (in bytes) of GRIB message array cgrib.
 771 !     iofst    - Bit offset of the beginning of Section 6.
 772 !     ngpts    - Number of grid points specified in the bit-map
 773 !
 774 !   OUTPUT ARGUMENT LIST:
 775 !     iofst    - Bit offset at the end of Section 6, returned.
 776 !     ibmap    - Bitmap indicator ( see Code Table 6.0 )
 777 !                0 = bitmap applies and is included in Section 6.
 778 !                1-253 = Predefined bitmap applies
 779 !                254 = Previously defined bitmap applies to this field
 780 !                255 = Bit map does not apply to this product.
 781 !     bmap()   - Logical*1 array containing decoded bitmap. ( if ibmap=0 )
 782 !     ierr     - Error return code.
 783 !                0 = no error
 784 !                4 = Unrecognized pre-defined bit-map.
 785 !
 786 ! REMARKS: None
 787 !
 788 ! ATTRIBUTES:
 789 !   LANGUAGE: Fortran 90
 790 !   MACHINE:  IBM SP
 791 !
 792 !$$$
 793
 794       character(len=1),intent(in) :: cgrib(lcgrib)
 795       integer,intent(in) :: lcgrib,ngpts
 796       integer,intent(inout) :: iofst
 797       integer,intent(out) :: ibmap
 798       integer,intent(out) :: ierr
 799       logical*1,intent(out) :: bmap(ngpts)
 800
 801       integer :: intbmap(ngpts)
 802
 803       ierr=0
 804
 805       iofst=iofst+32    ! skip Length of Section
 806       iofst=iofst+8     ! skip section number
 807
 808       call gbyte(cgrib,ibmap,iofst,8)    ! Get bit-map indicator
 809       iofst=iofst+8
 810
 811       if (ibmap.eq.0) then               ! Unpack bitmap
 812         call gbytes(cgrib,intbmap,iofst,1,0,ngpts)
 813         iofst=iofst+ngpts
 814         do j=1,ngpts
 815           bmap(j)=.true.
 816           if (intbmap(j).eq.0) bmap(j)=.false.
 817         enddo
 818       elseif (ibmap.eq.254) then               ! Use previous bitmap
 819         return
 820       elseif (ibmap.eq.255) then               ! No bitmap in message
 821         bmap(1:ngpts)=.true.
 822       else
 823         print *,'unpack6: Predefined bitmap ',ibmap,' not recognized.'
 824         ierr=4
 825       endif
 826
 827       return    ! End of Section 6 processing
 828       end
 829