source: trunk/WRF.COMMON/WRFV3/external/io_grib2/g2lib/jpcpack.F @ 3567

      subroutine jpcpack(fld,width,height,idrstmpl,cpack,lcpack,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    jpcpack
!   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-12-17
!
! ABSTRACT: This subroutine packs up a data field into a JPEG2000 code stream.
!   After the data field is scaled, and the reference value is subtracted out,
!   it is treated as a grayscale image and passed to a JPEG2000 encoder.
!   It also fills in GRIB2 Data Representation Template 5.40 or 5.40000 with the
!   appropriate values.
!
! PROGRAM HISTORY LOG:
! 2002-12-17  Gilbert
! 2004-07-19  Gilbert - Added check on whether the jpeg2000 encoding was
!                       successful.  If not, try again with different encoder
!                       options.
!
! USAGE:    CALL jpcpack(fld,width,height,idrstmpl,cpack,lcpack)
!   INPUT ARGUMENT LIST:
!     fld()    - Contains the data values to pack
!     width    - number of points in the x direction
!     height   - number of points in the y direction
!     idrstmpl - Contains the array of values for Data Representation
!                Template 5.40 or 5.40000
!                (1) = Reference value - ignored on input
!                (2) = Binary Scale Factor
!                (3) = Decimal Scale Factor
!                (4) = number of bits for each data value - ignored on input
!                (5) = Original field type - currently ignored on input
!                      Data values assumed to be reals.
!                (6) = 0 - use lossless compression
!                    = 1 - use lossy compression
!                (7) = Desired compression ratio, if idrstmpl(6)=1.
!                      Set to 255, if idrstmpl(6)=0.
!     lcpack   - size of array cpack().
!
!   OUTPUT ARGUMENT LIST: 
!     idrstmpl - Contains the array of values for Data Representation
!                Template 5.0
!                (1) = Reference value - set by jpcpack routine.
!                (2) = Binary Scale Factor - unchanged from input
!                (3) = Decimal Scale Factor - unchanged from input
!                (4) = Number of bits containing each grayscale pixel value
!                (5) = Original field type - currently set = 0 on output.
!                      Data values assumed to be reals.
!                (6) = 0 - use lossless compression
!                    = 1 - use lossy compression
!                (7) = Desired compression ratio, if idrstmpl(6)=1
!     cpack    - The packed data field (character*1 array)
!     lcpack   - length of packed field in cpack().
!
! REMARKS: None
!
! ATTRIBUTES:
!   LANGUAGE: XL Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      integer,intent(in) :: width,height
      real,intent(in) :: fld(width*height)
      character(len=1),intent(out) :: cpack(*)
      integer,intent(inout) :: idrstmpl(*)
      integer,intent(inout) :: lcpack
      integer,intent(out)   :: ierr

      real(4) :: ref
      integer(4) :: iref
      integer :: ifld(width*height),retry
      integer,parameter :: zero=0
      integer :: enc_jpeg2000
      character(len=1),allocatable :: ctemp(:)
      integer :: orig_dscalefct
      integer :: orig_bscalefct
      integer(8) :: maxdif,imin,imax
      
      ierr = 0

      ndpts=width*height
      bscale=2.0**real(-idrstmpl(2))
      dscale=10.0**real(idrstmpl(3))

      orig_bscalefct = idrstmpl(2)
      orig_dscalefct = idrstmpl(3)

!
!  Find max and min values in the data
!
      rmax=fld(1)
      rmin=fld(1)
      do j=2,ndpts
        if (fld(j).gt.rmax) rmax=fld(j)
        if (fld(j).lt.rmin) rmin=fld(j)
      enddo
      if (idrstmpl(2).eq.0) then
         maxdif=nint(rmax*dscale,8)-nint(rmin*dscale,8)
      else
         maxdif=nint((rmax-rmin)*dscale*bscale,8)
      endif
!
!  If max and min values are not equal, pack up field.
!  If they are equal, we have a constant field, and the reference
!  value (rmin) is the value for each point in the field and
!  set nbits to 0.
!
      if (rmin.ne.rmax .AND. maxdif.ne.0) then
        !
        !  Determine which algorithm to use based on user-supplied 
        !  binary scale factor and number of bits.
        !

        if (idrstmpl(2).eq.0) then
           !
           !  No binary scaling and calculate minimum number of 
           !  bits in which the data will fit.
           !

           nbits = 25

           do while (nbits .gt. 24) 
              imin=nint(rmin*dscale,8)
              imax=nint(rmax*dscale,8)
              maxdif=imax-imin
              temp=alog(real(maxdif+1))/alog(2.0)
              nbits=ceiling(temp)
              !
              ! These lines assure that we do no overflow
              !
              !  Todd Hutchinson, WSI 9/23/05
              !
              if (nbits .le. 24) then 
                 exit
              else
                 idrstmpl(3) = idrstmpl(3) - 1
                 dscale=10.0**real(idrstmpl(3))
              endif
           enddo

           rmin=real(imin)
           !   scale data
           do j=1,ndpts
             ifld(j)=nint(fld(j)*dscale)-imin
           enddo
        else
           !
           !  Use binary scaling factor and calculate minimum number of 
           !  bits in which the data will fit.
           !
           nbits = 25

           do while (nbits .gt. 24) 
              rmin=rmin*dscale
              rmax=rmax*dscale
              maxdif=nint((rmax-rmin)*bscale)
              temp=alog(real(maxdif+1))/alog(2.0)
              nbits=ceiling(temp)
              !
              ! These lines assure that we do no overflow
              !
              !  Todd Hutchinson, WSI 9/23/05
              !
              if (nbits .le. 24) then 
                 exit
              else
                 idrstmpl(2) = idrstmpl(2) + 1
                 bscale=2.0**real(-idrstmpl(2))
              endif              
           enddo
           !   scale data
           do j=1,ndpts
             ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
           enddo
        endif

        if (idrstmpl(2) .ne. orig_bscalefct) then 
           write(6,'(A,I2,A,I2,A)')
     &          ' JPCPACK: Reduced binary scale fctr from ',
     &          orig_bscalefct,' to ',idrstmpl(2),
     &          ' to prevent overflow'
           ierr = 12
        endif

        if (idrstmpl(3) .ne. orig_dscalefct) then 
           write(6,'(A,I2,A,I2,A)')
     &          ' JPCPACK: Reduced decimal scale fctr from ',
     &          orig_dscalefct,' to ',idrstmpl(3), 
     &          ' to prevent overflow'
           ierr = 11
        endif

        !
        !  Pack data into full octets, then do JPEG2000 encode.
        !  and calculate the length of the packed data in bytes
        !
        retry=0
        nbytes=(nbits+7)/8
        nsize=lcpack      ! needed for input to enc_jpeg2000
        allocate(ctemp(nbytes*ndpts))
        call g2lib_sbytes(ctemp,ifld,0,nbytes*8,0,ndpts)
        lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6),
     &                      idrstmpl(7),retry,cpack,nsize)
        if (lcpack.le.0) then
           print *,'jpcpack: ERROR Packing JPC=',lcpack
           if (lcpack.eq.-3) then
              retry=1
              print *,'jpcpack: Retrying....'
              lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6),
     &                         idrstmpl(7),retry,cpack,nsize)
              if (lcpack.le.0) then
                 print *,'jpcpack: Retry Failed.'
              else
                 print *,'jpcpack: Retry Successful.'
              endif
           endif
        endif
        deallocate(ctemp)

      else
        nbits=0
        lcpack=0
      endif

!
!  Fill in ref value and number of bits in Template 5.0
!
      call mkieee(rmin,ref,1)   ! ensure reference value is IEEE format
!      call g2lib_gbyte(ref,idrstmpl(1),0,32)
      iref=transfer(ref,iref)
      idrstmpl(1)=iref
      idrstmpl(4)=nbits
      idrstmpl(5)=0         ! original data were reals
      if (idrstmpl(6).eq.0) idrstmpl(7)=255       ! lossy not used

      return
      end