external/io_grib2/g2lib/addfield.F

   1       subroutine addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen,
   2      &                    coordlist,numcoord,idrsnum,idrstmpl,
   3      &                    idrstmplen,fld,ngrdpts,ibmap,bmap,ierr)
   4 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   5 !                .      .    .                                       .
   6 ! SUBPROGRAM:    addfield
   7 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-02
   8 !
   9 ! ABSTRACT: This subroutine packs up Sections 4 through 7 for a given field
  10 !   and adds them to a GRIB2 message.  They are Product Definition Section,
  11 !   Data Representation Section, Bit-Map Section and Data Section,
  12 !   respectively.
  13 !   This routine is used with routines "gribcreate", "addlocal", "addgrid",
  14 !   and "gribend" to create a complete GRIB2 message.  Subroutine
  15 !   gribcreate must be called first to initialize a new GRIB2 message.
  16 !   Also, subroutine addgrid must be called after gribcreate and
  17 !   before this routine to add the appropriate grid description to
  18 !   the GRIB2 message.   Also, a call to gribend is required to complete
  19 !   GRIB2 message after all fields have been added.
  20 !
  21 ! PROGRAM HISTORY LOG:
  22 ! 2000-05-02  Gilbert
  23 ! 2002-12-17  Gilbert  - Added support for new templates using
  24 !                        PNG and JPEG2000 algorithms/templates.
  25 ! 2004-06-22  Gilbert  - Added check to determine if packing algorithm failed.
  26 !
  27 ! USAGE:    CALL addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen,
  28 !                         coordlist,numcoord,idrsnum,idrstmpl,
  29 !                         idrstmplen,fld,ngrdpts,ibmap,bmap,ierr)
  30 !   INPUT ARGUMENT LIST:
  31 !     cgrib    - Character array to contain the GRIB2 message
  32 !     lcgrib   - Maximum length (bytes) of array cgrib.
  33 !     ipdsnum  - Product Definition Template Number ( see Code Table 4.0)
  34 !     ipdstmpl - Contains the data values for the specified Product Definition
  35 !                Template ( N=ipdsnum ).  Each element of this integer
  36 !                array contains an entry (in the order specified) of Product
  37 !                Defintion Template 4.N
  38 !   ipdstmplen - Max dimension of ipdstmpl()
  39 !     coordlist- Array containg floating point values intended to document
  40 !                the vertical discretisation associated to model data
  41 !                on hybrid coordinate vertical levels.
  42 !     numcoord - number of values in array coordlist.
  43 !     idrsnum  - Data Representation Template Number ( see Code Table 5.0 )
  44 !     idrstmpl - Contains the data values for the specified Data Representation
  45 !                Template ( N=idrsnum ).  Each element of this integer
  46 !                array contains an entry (in the order specified) of Data
  47 !                Representation Template 5.N
  48 !                Note that some values in this template (eg. reference
  49 !                values, number of bits, etc...) may be changed by the
  50 !                data packing algorithms.
  51 !                Use this to specify scaling factors and order of
  52 !                spatial differencing, if desired.
  53 !   idrstmplen - Max dimension of idrstmpl()
  54 !     fld()    - Array of data points to pack.
  55 !     ngrdpts  - Number of data points in grid.
  56 !                i.e.  size of fld and bmap.
  57 !     ibmap    - Bitmap indicator ( see Code Table 6.0 )
  58 !                0 = bitmap applies and is included in Section 6.
  59 !                1-253 = Predefined bitmap applies
  60 !                254 = Previously defined bitmap applies to this field
  61 !                255 = Bit map does not apply to this product.
  62 !     bmap()   - Logical*1 array containing bitmap to be added.
  63 !                ( if ibmap=0 or ibmap=254)
  64 !
  65 !   OUTPUT ARGUMENT LIST:
  66 !     cgrib    - Character array to contain the GRIB2 message
  67 !     ierr     - Error return code.
  68 !                0 = no error
  69 !                1 = GRIB message was not initialized.  Need to call
  70 !                    routine gribcreate first.
  71 !                2 = GRIB message already complete.  Cannot add new section.
  72 !                3 = Sum of Section byte counts does not add to total
  73 !                    byte count.
  74 !                4 = Previous Section was not 3 or 7.
  75 !                5 = Could not find requested Product Definition Template.
  76 !                6 = Section 3 (GDS) not previously defined in message
  77 !                7 = Tried to use unsupported Data Representationi Template
  78 !                8 = Specified use of a previously defined bitmap, but one
  79 !                    does not exist in the GRIB message.
  80 !                9 = GDT of one of 5.50 through 5.53 required to pack
  81 !                    using DRT 5.51
  82 !               10 = Error packing data field.
  83 !
  84 ! REMARKS: Note that the Local Use Section ( Section 2 ) can only follow
  85 !          Section 1 or Section 7 in a GRIB2 message.
  86 !
  87 ! ATTRIBUTES:
  88 !   LANGUAGE: Fortran 90
  89 !   MACHINE:  IBM SP
  90 !
  91 !$$$
  92
  93       use pdstemplates
  94       use drstemplates
  95
  96       character(len=1),intent(inout) :: cgrib(lcgrib)
  97       integer,intent(in) :: ipdsnum,ipdstmpl(*)
  98       integer,intent(in) :: idrsnum,numcoord,ipdstmplen,idrstmplen
  99       integer,intent(in) :: lcgrib,ngrdpts,ibmap
 100       real,intent(in) :: coordlist(numcoord)
 101       real,target,intent(in) :: fld(ngrdpts)
 102       integer,intent(out) :: ierr
 103       integer,intent(inout) :: idrstmpl(*)
 104       logical*1,intent(in) :: bmap(ngrdpts)
 105
 106       character(len=4),parameter :: grib='GRIB',c7777='7777'
 107       character(len=4):: ctemp
 108       character(len=1),allocatable :: cpack(:)
 109       real,pointer,dimension(:) :: pfld
 110       real(4) :: coordieee(numcoord),re00
 111       integer(4) :: ire00,allones
 112       integer :: mappds(ipdstmplen),intbmap(ngrdpts),mapdrs(idrstmplen)
 113       integer,parameter :: zero=0,one=1,four=4,five=5,six=6,seven=7
 114       integer,parameter :: minsize=50000
 115       integer iofst,ibeg,lencurr,len,mappdslen,mapdrslen,lpos3
 116       integer width,height,ndpts
 117       integer lensec3,lensec4,lensec5,lensec6,lensec7
 118       logical issec3,needext,isprevbmap
 119
 120       ierr=0
 121       do jj=0,31
 122          allones=ibset(allones,jj)
 123       enddo
 124 !
 125 !  Check to see if beginning of GRIB message exists
 126 !
 127       ctemp=cgrib(1)//cgrib(2)//cgrib(3)//cgrib(4)
 128       if ( ctemp.ne.grib ) then
 129         print *,'addfield: GRIB not found in given message.'
 130         print *,'addfield: Call to routine gribcreate required',
 131      &          ' to initialize GRIB messge.'
 132         ierr=1
 133         return
 134       endif
 135 !
 136 !  Get current length of GRIB message
 137 !
 138       call g2lib_gbyte(cgrib,lencurr,96,32)
 139 !
 140 !  Check to see if GRIB message is already complete
 141 !
 142       ctemp=cgrib(lencurr-3)//cgrib(lencurr-2)//cgrib(lencurr-1)
 143      &      //cgrib(lencurr)
 144       if ( ctemp.eq.c7777 ) then
 145         print *,'addfield: GRIB message already complete.  Cannot',
 146      &          ' add new section.'
 147         ierr=2
 148         return
 149       endif
 150 !
 151 !  Loop through all current sections of the GRIB message to
 152 !  find the last section number.
 153 !
 154       issec3=.false.
 155       isprevbmap=.false.
 156       len=16    ! length of Section 0
 157       do
 158       !    Get number and length of next section
 159         iofst=len*8
 160         call g2lib_gbyte(cgrib,ilen,iofst,32)
 161         iofst=iofst+32
 162         call g2lib_gbyte(cgrib,isecnum,iofst,8)
 163         iofst=iofst+8
 164       !  Check if previous Section 3 exists and save location of
 165       !  the section 3 in case needed later.
 166         if (isecnum.eq.3) then
 167            issec3=.true.
 168            lpos3=len+1
 169            lensec3=ilen
 170         endif
 171       !  Check if a previous defined bitmap exists
 172         if (isecnum.eq.6) then
 173           call g2lib_gbyte(cgrib,ibmprev,iofst,8)
 174           iofst=iofst+8
 175           if ((ibmprev.ge.0).and.(ibmprev.le.253)) isprevbmap=.true.
 176         endif
 177         len=len+ilen
 178       !    Exit loop if last section reached
 179         if ( len.eq.lencurr ) exit
 180       !    If byte count for each section does not match current
 181       !    total length, then there is a problem.
 182         if ( len.gt.lencurr ) then
 183           print *,'addfield: Section byte counts don''t add to total.'
 184           print *,'addfield: Sum of section byte counts = ',len
 185           print *,'addfield: Total byte count in Section 0 = ',lencurr
 186           ierr=3
 187           return
 188         endif
 189       enddo
 190 !
 191 !  Sections 4 through 7 can only be added after section 3 or 7.
 192 !
 193       if ( (isecnum.ne.3) .and. (isecnum.ne.7) ) then
 194         print *,'addfield: Sections 4-7 can only be added after',
 195      &          ' Section 3 or 7.'
 196         print *,'addfield: Section ',isecnum,' was the last found in',
 197      &          ' given GRIB message.'
 198         ierr=4
 199         return
 200 !
 201 !  Sections 4 through 7 can only be added if section 3 was previously defined.
 202 !
 203       elseif (.not.issec3) then
 204         print *,'addfield: Sections 4-7 can only be added if Section',
 205      &          ' 3 was previously included.'
 206         print *,'addfield: Section 3 was not found in',
 207      &          ' given GRIB message.'
 208         print *,'addfield: Call to routine addgrid required',
 209      &          ' to specify Grid definition.'
 210         ierr=6
 211         return
 212       endif
 213 !
 214 !  Add Section 4  - Product Definition Section
 215 !
 216       ibeg=lencurr*8        !   Calculate offset for beginning of section 4
 217       iofst=ibeg+32         !   leave space for length of section
 218       call g2lib_sbyte(cgrib,four,iofst,8)     ! Store section number ( 4 )
 219       iofst=iofst+8
 220       call g2lib_sbyte(cgrib,numcoord,iofst,16)   ! Store num of coordinate values
 221       iofst=iofst+16
 222       call g2lib_sbyte(cgrib,ipdsnum,iofst,16)    ! Store Prod Def Template num.
 223       iofst=iofst+16
 224       !
 225       !   Get Product Definition Template
 226       !
 227       call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret)
 228       if (iret.ne.0) then
 229         ierr=5
 230         return
 231       endif
 232       !
 233       !   Extend the Product Definition Template, if necessary.
 234       !   The number of values in a specific template may vary
 235       !   depending on data specified in the "static" part of the
 236       !   template.
 237       !
 238       if ( needext ) then
 239         call extpdstemplate(ipdsnum,ipdstmpl,mappdslen,mappds)
 240       endif
 241       !
 242       !   Pack up each input value in array ipdstmpl into the
 243       !   the appropriate number of octets, which are specified in
 244       !   corresponding entries in array mappds.
 245       !
 246       do i=1,mappdslen
 247         nbits=iabs(mappds(i))*8
 248         if ( (mappds(i).ge.0).or.(ipdstmpl(i).ge.0) ) then
 249           call g2lib_sbyte(cgrib,ipdstmpl(i),iofst,nbits)
 250         else
 251           call g2lib_sbyte(cgrib,one,iofst,1)
 252           call g2lib_sbyte(cgrib,iabs(ipdstmpl(i)),iofst+1,nbits-1)
 253         endif
 254         iofst=iofst+nbits
 255       enddo
 256       !
 257       !   Add Optional list of vertical coordinate values
 258       !   after the Product Definition Template, if necessary.
 259       !
 260       if ( numcoord .ne. 0 ) then
 261         call mkieee(coordlist,coordieee,numcoord)
 262         call g2lib_sbytes(cgrib,coordieee,iofst,32,0,numcoord)
 263         iofst=iofst+(32*numcoord)
 264       endif
 265       !
 266       !   Calculate length of section 4 and store it in octets
 267       !   1-4 of section 4.
 268       !
 269       lensec4=(iofst-ibeg)/8
 270       call g2lib_sbyte(cgrib,lensec4,ibeg,32)
 271 !
 272 !  Pack Data using appropriate algorithm
 273 !
 274       !
 275       !   Get Data Representation Template
 276       !
 277       call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret)
 278       if (iret.ne.0) then
 279         ierr=5
 280         return
 281       endif
 282       !
 283       !  contract data field, removing data at invalid grid points,
 284       !  if bit-map is provided with field.
 285       !
 286       if ( ibmap.eq.0 .OR. ibmap.eq.254 ) then
 287          allocate(pfld(ngrdpts))
 288          ndpts=0;
 289          do jj=1,ngrdpts
 290              intbmap(jj)=0
 291              if ( bmap(jj) ) then
 292                 intbmap(jj)=1
 293                 ndpts=ndpts+1
 294                 pfld(ndpts)=fld(jj);
 295              endif
 296          enddo
 297       else
 298          ndpts=ngrdpts;
 299          pfld=>fld;
 300       endif
 301       lcpack=0
 302       nsize=ndpts*4
 303       if (nsize .lt. minsize) nsize=minsize
 304       allocate(cpack(nsize),stat=istat)
 305       if (idrsnum.eq.0) then      !  Simple Packing
 306         call simpack(pfld,ndpts,idrstmpl,cpack,lcpack)
 307       elseif (idrsnum.eq.2.or.idrsnum.eq.3) then      !  Complex Packing
 308         call cmplxpack(pfld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
 309       elseif (idrsnum.eq.50) then      !  Sperical Harmonic Simple Packing
 310         call simpack(pfld(2),ndpts-1,idrstmpl,cpack,lcpack)
 311         call mkieee(real(pfld(1)),re00,1)  ! ensure RE(0,0) value is IEEE format
 312         !call g2lib_gbyte(re00,idrstmpl(5),0,32)
 313         ire00=transfer(re00,ire00)
 314         idrstmpl(5)=ire00
 315       elseif (idrsnum.eq.51) then      !  Sperical Harmonic Complex Packing
 316            call getpoly(cgrib(lpos3),lensec3,jj,kk,mm)
 317            if (jj.ne.0 .AND. kk.ne.0 .AND. mm.ne.0) then
 318              call specpack(pfld,ndpts,jj,kk,mm,idrstmpl,cpack,lcpack)
 319            else
 320              print *,'addfield: Cannot pack DRT 5.51.'
 321              ierr=9
 322              return
 323            endif
 324 #ifdef USE_JPEG2000
 325       elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then     !  JPEG2000 encoding
 326         if (ibmap.eq.255) then
 327            call getdim(cgrib(lpos3),lensec3,width,height,iscan)
 328            if ( width.eq.0 .OR. height.eq.0 ) then
 329               width=ndpts
 330               height=1
 331            elseif ( width.eq.allones .OR. height.eq.allones ) then
 332               width=ndpts
 333               height=1
 334            elseif ( ibits(iscan,5,1) .eq. 1) then   ! Scanning mode: bit 3
 335               itemp=width
 336               width=height
 337               height=itemp
 338            endif
 339         else
 340            width=ndpts
 341            height=1
 342         endif
 343         lcpack=nsize
 344         call jpcpack(pfld,width,height,idrstmpl,cpack,lcpack,ierr)
 345 #endif  /* USE_JPEG2000 */
 346 #ifdef USE_PNG
 347       elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then      !  PNG encoding
 348         if (ibmap.eq.255) then
 349            call getdim(cgrib(lpos3),lensec3,width,height,iscan)
 350            if ( width.eq.0 .OR. height.eq.0 ) then
 351               width=ndpts
 352               height=1
 353            elseif ( width.eq.allones .OR. height.eq.allones ) then
 354               width=ndpts
 355               height=1
 356            elseif ( ibits(iscan,5,1) .eq. 1) then   ! Scanning mode: bit 3
 357               itemp=width
 358               width=height
 359               height=itemp
 360            endif
 361         else
 362            width=ndpts
 363            height=1
 364         endif
 365         call pngpack(pfld,width,height,idrstmpl,cpack,lcpack)
 366 #endif   /* USE_PNG */
 367       else
 368         print *,'addfield: Data Representation Template 5.',idrsnum,
 369      *          ' not yet implemented.'
 370         ierr=7
 371         return
 372       endif
 373       if ( ibmap.eq.0 .OR. ibmap.eq.254 ) then
 374          deallocate(pfld)
 375       endif
 376       if ( lcpack .lt. 0 ) then
 377         if( allocated(cpack) )deallocate(cpack)
 378         ierr=10
 379         return
 380       endif
 381
 382 !
 383 !  Add Section 5  - Data Representation Section
 384 !
 385       ibeg=iofst            !   Calculate offset for beginning of section 5
 386       iofst=ibeg+32         !   leave space for length of section
 387       call g2lib_sbyte(cgrib,five,iofst,8)     ! Store section number ( 5 )
 388       iofst=iofst+8
 389       call g2lib_sbyte(cgrib,ndpts,iofst,32)    ! Store num of actual data points
 390       iofst=iofst+32
 391       call g2lib_sbyte(cgrib,idrsnum,iofst,16)    ! Store Data Repr. Template num.
 392       iofst=iofst+16
 393       !
 394       !   Pack up each input value in array idrstmpl into the
 395       !   the appropriate number of octets, which are specified in
 396       !   corresponding entries in array mapdrs.
 397       !
 398       do i=1,mapdrslen
 399         nbits=iabs(mapdrs(i))*8
 400         if ( (mapdrs(i).ge.0).or.(idrstmpl(i).ge.0) ) then
 401           call g2lib_sbyte(cgrib,idrstmpl(i),iofst,nbits)
 402         else
 403           call g2lib_sbyte(cgrib,one,iofst,1)
 404           call g2lib_sbyte(cgrib,iabs(idrstmpl(i)),iofst+1,nbits-1)
 405         endif
 406         iofst=iofst+nbits
 407       enddo
 408       !
 409       !   Calculate length of section 5 and store it in octets
 410       !   1-4 of section 5.
 411       !
 412       lensec5=(iofst-ibeg)/8
 413       call g2lib_sbyte(cgrib,lensec5,ibeg,32)
 414
 415 !
 416 !  Add Section 6  - Bit-Map Section
 417 !
 418       ibeg=iofst            !   Calculate offset for beginning of section 6
 419       iofst=ibeg+32         !   leave space for length of section
 420       call g2lib_sbyte(cgrib,six,iofst,8)     ! Store section number ( 6 )
 421       iofst=iofst+8
 422       call g2lib_sbyte(cgrib,ibmap,iofst,8)    ! Store Bit Map indicator
 423       iofst=iofst+8
 424       !
 425       !  Store bitmap, if supplied
 426       !
 427       if (ibmap.eq.0) then
 428         call g2lib_sbytes(cgrib,intbmap,iofst,1,0,ngrdpts)    ! Store BitMap
 429         iofst=iofst+ngrdpts
 430       endif
 431       !
 432       !  If specifying a previously defined bit-map, make sure
 433       !  one already exists in the current GRIB message.
 434       !
 435       if ((ibmap.eq.254).and.(.not.isprevbmap)) then
 436         print *,'addfield: Requested previously defined bitmap, ',
 437      &        ' but one does not exist in the current GRIB message.'
 438         ierr=8
 439         return
 440       endif
 441       !
 442       !   Calculate length of section 6 and store it in octets
 443       !   1-4 of section 6.  Pad to end of octect, if necessary.
 444       !
 445       left=8-mod(iofst,8)
 446       if (left.ne.8) then
 447         call g2lib_sbyte(cgrib,zero,iofst,left)     ! Pad with zeros to fill Octet
 448         iofst=iofst+left
 449       endif
 450       lensec6=(iofst-ibeg)/8
 451       call g2lib_sbyte(cgrib,lensec6,ibeg,32)
 452
 453 !
 454 !  Add Section 7  - Data Section
 455 !
 456       ibeg=iofst            !   Calculate offset for beginning of section 7
 457       iofst=ibeg+32         !   leave space for length of section
 458       call g2lib_sbyte(cgrib,seven,iofst,8)     ! Store section number ( 7 )
 459       iofst=iofst+8
 460       !      Store Packed Binary Data values, if non-constant field
 461       if (lcpack.ne.0) then
 462         ioctet=iofst/8
 463         cgrib(ioctet+1:ioctet+lcpack)=cpack(1:lcpack)
 464         iofst=iofst+(8*lcpack)
 465       endif
 466       !
 467       !   Calculate length of section 7 and store it in octets
 468       !   1-4 of section 7.
 469       !
 470       lensec7=(iofst-ibeg)/8
 471       call g2lib_sbyte(cgrib,lensec7,ibeg,32)
 472
 473       if( allocated(cpack) )deallocate(cpack)
 474 !
 475 !  Update current byte total of message in Section 0
 476 !
 477       newlen=lencurr+lensec4+lensec5+lensec6+lensec7
 478       call g2lib_sbyte(cgrib,newlen,96,32)
 479
 480       return
 481       end
 482