subroutine addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen, 1,43
& coordlist,numcoord,idrsnum,idrstmpl,
& idrstmplen,fld,ngrdpts,ibmap,bmap,ierr)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . . .
! SUBPROGRAM: addfield
! PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-05-02
!
! ABSTRACT: This subroutine packs up Sections 4 through 7 for a given field
! and adds them to a GRIB2 message. They are Product Definition Section,
! Data Representation Section, Bit-Map Section and Data Section,
! respectively.
! This routine is used with routines "gribcreate", "addlocal", "addgrid",
! and "gribend" to create a complete GRIB2 message. Subroutine
! gribcreate must be called first to initialize a new GRIB2 message.
! Also, subroutine addgrid must be called after gribcreate and
! before this routine to add the appropriate grid description to
! the GRIB2 message. Also, a call to gribend is required to complete
! GRIB2 message after all fields have been added.
!
! PROGRAM HISTORY LOG:
! 2000-05-02 Gilbert
! 2002-12-17 Gilbert - Added support for new templates using
! PNG and JPEG2000 algorithms/templates.
! 2004-06-22 Gilbert - Added check to determine if packing algorithm failed.
!
! USAGE: CALL addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen,
! coordlist,numcoord,idrsnum,idrstmpl,
! idrstmplen,fld,ngrdpts,ibmap,bmap,ierr)
! INPUT ARGUMENT LIST:
! cgrib - Character array to contain the GRIB2 message
! lcgrib - Maximum length (bytes) of array cgrib.
! ipdsnum - Product Definition Template Number ( see Code Table 4.0)
! ipdstmpl - Contains the data values for the specified Product Definition
! Template ( N=ipdsnum ). Each element of this integer
! array contains an entry (in the order specified) of Product
! Defintion Template 4.N
! ipdstmplen - Max dimension of ipdstmpl()
! coordlist- Array containg floating point values intended to document
! the vertical discretisation associated to model data
! on hybrid coordinate vertical levels.
! numcoord - number of values in array coordlist.
! idrsnum - Data Representation Template Number ( see Code Table 5.0 )
! idrstmpl - Contains the data values for the specified Data Representation
! Template ( N=idrsnum ). Each element of this integer
! array contains an entry (in the order specified) of Data
! Representation Template 5.N
! Note that some values in this template (eg. reference
! values, number of bits, etc...) may be changed by the
! data packing algorithms.
! Use this to specify scaling factors and order of
! spatial differencing, if desired.
! idrstmplen - Max dimension of idrstmpl()
! fld() - Array of data points to pack.
! ngrdpts - Number of data points in grid.
! i.e. size of fld and bmap.
! ibmap - Bitmap indicator ( see Code Table 6.0 )
! 0 = bitmap applies and is included in Section 6.
! 1-253 = Predefined bitmap applies
! 254 = Previously defined bitmap applies to this field
! 255 = Bit map does not apply to this product.
! bmap() - Logical*1 array containing bitmap to be added.
! ( if ibmap=0 or ibmap=254)
!
! OUTPUT ARGUMENT LIST:
! cgrib - Character array to contain the GRIB2 message
! ierr - Error return code.
! 0 = no error
! 1 = GRIB message was not initialized. Need to call
! routine gribcreate first.
! 2 = GRIB message already complete. Cannot add new section.
! 3 = Sum of Section byte counts does not add to total
! byte count.
! 4 = Previous Section was not 3 or 7.
! 5 = Could not find requested Product Definition Template.
! 6 = Section 3 (GDS) not previously defined in message
! 7 = Tried to use unsupported Data Representationi Template
! 8 = Specified use of a previously defined bitmap, but one
! does not exist in the GRIB message.
! 9 = GDT of one of 5.50 through 5.53 required to pack
! using DRT 5.51
! 10 = Error packing data field.
!
! REMARKS: Note that the Local Use Section ( Section 2 ) can only follow
! Section 1 or Section 7 in a GRIB2 message.
!
! ATTRIBUTES:
! LANGUAGE: Fortran 90
! MACHINE: IBM SP
!
!$$$
use pdstemplates
use drstemplates
character(len=1),intent(inout) :: cgrib(lcgrib)
integer,intent(in) :: ipdsnum,ipdstmpl(*)
integer,intent(in) :: idrsnum,numcoord,ipdstmplen,idrstmplen
integer,intent(in) :: lcgrib,ngrdpts,ibmap
real,intent(in) :: coordlist(numcoord)
real,target,intent(in) :: fld(ngrdpts)
integer,intent(out) :: ierr
integer,intent(inout) :: idrstmpl(*)
logical*1,intent(in) :: bmap(ngrdpts)
character(len=4),parameter :: grib='GRIB',c7777='7777'
character(len=4):: ctemp
character(len=1),allocatable :: cpack(:)
real,pointer,dimension(:) :: pfld
real(4) :: coordieee(numcoord),re00
integer(4) :: ire00,allones
integer :: mappds(ipdstmplen),intbmap(ngrdpts),mapdrs(idrstmplen)
integer,parameter :: zero=0,one=1,four=4,five=5,six=6,seven=7
integer iofst,ibeg,lencurr,len,mappdslen,mapdrslen,lpos3
integer width,height,ndpts
integer lensec3,lensec4,lensec5,lensec6,lensec7
logical issec3,needext,isprevbmap
ierr=0
do jj=0,31
allones=ibset(allones,jj)
enddo
!
! Check to see if beginning of GRIB message exists
!
ctemp=cgrib(1)//cgrib(2)//cgrib(3)//cgrib(4)
if ( ctemp.ne.grib ) then
print *,'addfield: GRIB not found in given message.'
print *,'addfield: Call to routine gribcreate required',
& ' to initialize GRIB messge.'
ierr=1
return
endif
!
! Get current length of GRIB message
!
call gbyte(cgrib,lencurr,96,32)
!
! Check to see if GRIB message is already complete
!
ctemp=cgrib(lencurr-3)//cgrib(lencurr-2)//cgrib(lencurr-1)
& //cgrib(lencurr)
if ( ctemp.eq.c7777 ) then
print *,'addfield: GRIB message already complete. Cannot',
& ' add new section.'
ierr=2
return
endif
!
! Loop through all current sections of the GRIB message to
! find the last section number.
!
issec3=.false.
isprevbmap=.false.
len=16 ! length of Section 0
do
! Get number and length of next section
iofst=len*8
call gbyte(cgrib,ilen,iofst,32)
iofst=iofst+32
call gbyte(cgrib,isecnum,iofst,8)
iofst=iofst+8
! Check if previous Section 3 exists and save location of
! the section 3 in case needed later.
if (isecnum.eq.3) then
issec3=.true.
lpos3=len+1
lensec3=ilen
endif
! Check if a previous defined bitmap exists
if (isecnum.eq.6) then
call gbyte(cgrib,ibmprev,iofst,8)
iofst=iofst+8
if ((ibmprev.ge.0).and.(ibmprev.le.253)) isprevbmap=.true.
endif
len=len+ilen
! Exit loop if last section reached
if ( len.eq.lencurr ) exit
! If byte count for each section does not match current
! total length, then there is a problem.
if ( len.gt.lencurr ) then
print *,'addfield: Section byte counts don''t add to total.'
print *,'addfield: Sum of section byte counts = ',len
print *,'addfield: Total byte count in Section 0 = ',lencurr
ierr=3
return
endif
enddo
!
! Sections 4 through 7 can only be added after section 3 or 7.
!
if ( (isecnum.ne.3) .and. (isecnum.ne.7) ) then
print *,'addfield: Sections 4-7 can only be added after',
& ' Section 3 or 7.'
print *,'addfield: Section ',isecnum,' was the last found in',
& ' given GRIB message.'
ierr=4
return
!
! Sections 4 through 7 can only be added if section 3 was previously defined.
!
elseif (.not.issec3) then
print *,'addfield: Sections 4-7 can only be added if Section',
& ' 3 was previously included.'
print *,'addfield: Section 3 was not found in',
& ' given GRIB message.'
print *,'addfield: Call to routine addgrid required',
& ' to specify Grid definition.'
ierr=6
return
endif
!
! Add Section 4 - Product Definition Section
!
ibeg=lencurr*8 ! Calculate offset for beginning of section 4
iofst=ibeg+32 ! leave space for length of section
call sbyte(cgrib,four,iofst,8) ! Store section number ( 4 )
iofst=iofst+8
call sbyte(cgrib,numcoord,iofst,16) ! Store num of coordinate values
iofst=iofst+16
call sbyte(cgrib,ipdsnum,iofst,16) ! Store Prod Def Template num.
iofst=iofst+16
!
! Get Product Definition Template
!
call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret)
if (iret.ne.0) then
ierr=5
return
endif
!
! Extend the Product Definition Template, if necessary.
! The number of values in a specific template may vary
! depending on data specified in the "static" part of the
! template.
!
if ( needext ) then
call extpdstemplate(ipdsnum,ipdstmpl,mappdslen,mappds)
endif
!
! Pack up each input value in array ipdstmpl into the
! the appropriate number of octets, which are specified in
! corresponding entries in array mappds.
!
do i=1,mappdslen
nbits=iabs(mappds(i))*8
if ( (mappds(i).ge.0).or.(ipdstmpl(i).ge.0) ) then
call sbyte(cgrib,ipdstmpl(i),iofst,nbits)
else
call sbyte(cgrib,one,iofst,1)
call sbyte(cgrib,iabs(ipdstmpl(i)),iofst+1,nbits-1)
endif
iofst=iofst+nbits
enddo
!
! Add Optional list of vertical coordinate values
! after the Product Definition Template, if necessary.
!
if ( numcoord .ne. 0 ) then
call mkieee(coordlist,coordieee,numcoord)
call sbytes(cgrib,coordieee,iofst,32,0,numcoord)
iofst=iofst+(32*numcoord)
endif
!
! Calculate length of section 4 and store it in octets
! 1-4 of section 4.
!
lensec4=(iofst-ibeg)/8
call sbyte(cgrib,lensec4,ibeg,32)
!
! Pack Data using appropriate algorithm
!
!
! Get Data Representation Template
!
call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret)
if (iret.ne.0) then
ierr=5
return
endif
!
! contract data field, removing data at invalid grid points,
! if bit-map is provided with field.
!
if ( ibmap.eq.0 .OR. ibmap.eq.254 ) then
allocate(pfld(ngrdpts))
ndpts=0;
do jj=1,ngrdpts
intbmap(jj)=0
if ( bmap(jj) ) then
intbmap(jj)=1
ndpts=ndpts+1
pfld(ndpts)=fld(jj);
endif
enddo
else
ndpts=ngrdpts;
pfld=>fld;
endif
lcpack=0
allocate(cpack(ndpts*4),stat=istat)
if (idrsnum.eq.0) then ! Simple Packing
call simpack(pfld,ndpts,idrstmpl,cpack,lcpack)
elseif (idrsnum.eq.2.or.idrsnum.eq.3) then ! Complex Packing
call cmplxpack(pfld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
elseif (idrsnum.eq.50) then ! Sperical Harmonic Simple Packing
call simpack(pfld(2),ndpts-1,idrstmpl,cpack,lcpack)
call mkieee(real(pfld(1)),re00,1) ! ensure RE(0,0) value is IEEE format
!call gbyte(re00,idrstmpl(5),0,32)
ire00=transfer(re00,ire00)
idrstmpl(5)=ire00
elseif (idrsnum.eq.51) then ! Sperical Harmonic Complex Packing
call getpoly(cgrib(lpos3),lensec3,jj,kk,mm)
if (jj.ne.0 .AND. kk.ne.0 .AND. mm.ne.0) then
call specpack(pfld,ndpts,jj,kk,mm,idrstmpl,cpack,lcpack)
else
print *,'addfield: Cannot pack DRT 5.51.'
ierr=9
return
endif
#ifdef USE_JPEG2000
elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then ! JPEG2000 encoding
if (ibmap.eq.255) then
call getdim(cgrib(lpos3),lensec3,width,height,iscan)
if ( width.eq.0 .OR. height.eq.0 ) then
width=ndpts
height=1
elseif ( width.eq.allones .OR. height.eq.allones ) then
width=ndpts
height=1
elseif ( ibits(iscan,5,1) .eq. 1) then ! Scanning mode: bit 3
itemp=width
width=height
height=itemp
endif
else
width=ndpts
height=1
endif
call jpcpack(pfld,width,height,idrstmpl,cpack,lcpack)
#endif /* USE_JPEG2000 */
#ifdef USE_PNG
elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then ! PNG encoding
if (ibmap.eq.255) then
call getdim(cgrib(lpos3),lensec3,width,height,iscan)
if ( width.eq.0 .OR. height.eq.0 ) then
width=ndpts
height=1
elseif ( width.eq.allones .OR. height.eq.allones ) then
width=ndpts
height=1
elseif ( ibits(iscan,5,1) .eq. 1) then ! Scanning mode: bit 3
itemp=width
width=height
height=itemp
endif
else
width=ndpts
height=1
endif
call pngpack(pfld,width,height,idrstmpl,cpack,lcpack)
#endif /* USE_PNG */
else
print *,'addfield: Data Representation Template 5.',idrsnum,
* ' not yet implemented.'
ierr=7
return
endif
if ( ibmap.eq.0 .OR. ibmap.eq.254 ) then
deallocate(pfld)
endif
if ( lcpack .lt. 0 ) then
if( allocated(cpack) )deallocate(cpack)
ierr=10
return
endif
!
! Add Section 5 - Data Representation Section
!
ibeg=iofst ! Calculate offset for beginning of section 5
iofst=ibeg+32 ! leave space for length of section
call sbyte(cgrib,five,iofst,8) ! Store section number ( 5 )
iofst=iofst+8
call sbyte(cgrib,ndpts,iofst,32) ! Store num of actual data points
iofst=iofst+32
call sbyte(cgrib,idrsnum,iofst,16) ! Store Data Repr. Template num.
iofst=iofst+16
!
! Pack up each input value in array idrstmpl into the
! the appropriate number of octets, which are specified in
! corresponding entries in array mapdrs.
!
do i=1,mapdrslen
nbits=iabs(mapdrs(i))*8
if ( (mapdrs(i).ge.0).or.(idrstmpl(i).ge.0) ) then
call sbyte(cgrib,idrstmpl(i),iofst,nbits)
else
call sbyte(cgrib,one,iofst,1)
call sbyte(cgrib,iabs(idrstmpl(i)),iofst+1,nbits-1)
endif
iofst=iofst+nbits
enddo
!
! Calculate length of section 5 and store it in octets
! 1-4 of section 5.
!
lensec5=(iofst-ibeg)/8
call sbyte(cgrib,lensec5,ibeg,32)
!
! Add Section 6 - Bit-Map Section
!
ibeg=iofst ! Calculate offset for beginning of section 6
iofst=ibeg+32 ! leave space for length of section
call sbyte(cgrib,six,iofst,8) ! Store section number ( 6 )
iofst=iofst+8
call sbyte(cgrib,ibmap,iofst,8) ! Store Bit Map indicator
iofst=iofst+8
!
! Store bitmap, if supplied
!
if (ibmap.eq.0) then
call sbytes(cgrib,intbmap,iofst,1,0,ngrdpts) ! Store BitMap
iofst=iofst+ngrdpts
endif
!
! If specifying a previously defined bit-map, make sure
! one already exists in the current GRIB message.
!
if ((ibmap.eq.254).and.(.not.isprevbmap)) then
print *,'addfield: Requested previously defined bitmap, ',
& ' but one does not exist in the current GRIB message.'
ierr=8
return
endif
!
! Calculate length of section 6 and store it in octets
! 1-4 of section 6. Pad to end of octect, if necessary.
!
left=8-mod(iofst,8)
if (left.ne.8) then
call sbyte(cgrib,zero,iofst,left) ! Pad with zeros to fill Octet
iofst=iofst+left
endif
lensec6=(iofst-ibeg)/8
call sbyte(cgrib,lensec6,ibeg,32)
!
! Add Section 7 - Data Section
!
ibeg=iofst ! Calculate offset for beginning of section 7
iofst=ibeg+32 ! leave space for length of section
call sbyte(cgrib,seven,iofst,8) ! Store section number ( 7 )
iofst=iofst+8
! Store Packed Binary Data values, if non-constant field
if (lcpack.ne.0) then
ioctet=iofst/8
cgrib(ioctet+1:ioctet+lcpack)=cpack(1:lcpack)
iofst=iofst+(8*lcpack)
endif
!
! Calculate length of section 7 and store it in octets
! 1-4 of section 7.
!
lensec7=(iofst-ibeg)/8
call sbyte(cgrib,lensec7,ibeg,32)
if( allocated(cpack) )deallocate(cpack)
!
! Update current byte total of message in Section 0
!
newlen=lencurr+lensec4+lensec5+lensec6+lensec7
call sbyte(cgrib,newlen,96,32)
return
end