hwrf-baseline-20120103-1354/ungrib/src/ngl/g2/jpcpack.F

   1       subroutine jpcpack(fld,width,height,idrstmpl,cpack,lcpack)
   2 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   3 !                .      .    .                                       .
   4 ! SUBPROGRAM:    jpcpack
   5 !   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-12-17
   6 !
   7 ! ABSTRACT: This subroutine packs up a data field into a JPEG2000 code stream.
   8 !   After the data field is scaled, and the reference value is subtracted out,
   9 !   it is treated as a grayscale image and passed to a JPEG2000 encoder.
  10 !   It also fills in GRIB2 Data Representation Template 5.40 or 5.40000 with the
  11 !   appropriate values.
  12 !
  13 ! PROGRAM HISTORY LOG:
  14 ! 2002-12-17  Gilbert
  15 ! 2004-07-19  Gilbert - Added check on whether the jpeg2000 encoding was
  16 !                       successful.  If not, try again with different encoder
  17 !                       options.
  18 !
  19 ! USAGE:    CALL jpcpack(fld,width,height,idrstmpl,cpack,lcpack)
  20 !   INPUT ARGUMENT LIST:
  21 !     fld()    - Contains the data values to pack
  22 !     width    - number of points in the x direction
  23 !     height   - number of points in the y direction
  24 !     idrstmpl - Contains the array of values for Data Representation
  25 !                Template 5.40 or 5.40000
  26 !                (1) = Reference value - ignored on input
  27 !                (2) = Binary Scale Factor
  28 !                (3) = Decimal Scale Factor
  29 !                (4) = number of bits for each data value - ignored on input
  30 !                (5) = Original field type - currently ignored on input
  31 !                      Data values assumed to be reals.
  32 !                (6) = 0 - use lossless compression
  33 !                    = 1 - use lossy compression
  34 !                (7) = Desired compression ratio, if idrstmpl(6)=1.
  35 !                      Set to 255, if idrstmpl(6)=0.
  36 !     lcpack   - size of array cpack().
  37 !
  38 !   OUTPUT ARGUMENT LIST:
  39 !     idrstmpl - Contains the array of values for Data Representation
  40 !                Template 5.0
  41 !                (1) = Reference value - set by jpcpack routine.
  42 !                (2) = Binary Scale Factor - unchanged from input
  43 !                (3) = Decimal Scale Factor - unchanged from input
  44 !                (4) = Number of bits containing each grayscale pixel value
  45 !                (5) = Original field type - currently set = 0 on output.
  46 !                      Data values assumed to be reals.
  47 !                (6) = 0 - use lossless compression
  48 !                    = 1 - use lossy compression
  49 !                (7) = Desired compression ratio, if idrstmpl(6)=1
  50 !     cpack    - The packed data field (character*1 array)
  51 !     lcpack   - length of packed field in cpack().
  52 !
  53 ! REMARKS: None
  54 !
  55 ! ATTRIBUTES:
  56 !   LANGUAGE: XL Fortran 90
  57 !   MACHINE:  IBM SP
  58 !
  59 !$$$
  60
  61       integer,intent(in) :: width,height
  62       real,intent(in) :: fld(width*height)
  63       character(len=1),intent(out) :: cpack(*)
  64       integer,intent(inout) :: idrstmpl(*)
  65       integer,intent(inout) :: lcpack
  66
  67       real(4) :: ref,rmin4
  68       real(8) :: rmin,rmax
  69       integer(4) :: iref
  70       integer :: ifld(width*height),retry
  71       integer,parameter :: zero=0
  72       integer :: enc_jpeg2000
  73       character(len=1),allocatable :: ctemp(:)
  74 #ifdef USE_JPEG2000
  75
  76       ndpts=width*height
  77       bscale=2.0**real(-idrstmpl(2))
  78       dscale=10.0**real(idrstmpl(3))
  79 !
  80 !  Find max and min values in the data
  81 !
  82       rmax=fld(1)
  83       rmin=fld(1)
  84       do j=2,ndpts
  85         if (fld(j).gt.rmax) rmax=fld(j)
  86         if (fld(j).lt.rmin) rmin=fld(j)
  87       enddo
  88       if (idrstmpl(2).eq.0) then
  89          maxdif=nint(rmax*dscale)-nint(rmin*dscale)
  90       else
  91          maxdif=nint((rmax-rmin)*dscale*bscale)
  92       endif
  93 !
  94 !  If max and min values are not equal, pack up field.
  95 !  If they are equal, we have a constant field, and the reference
  96 !  value (rmin) is the value for each point in the field and
  97 !  set nbits to 0.
  98 !
  99       if (rmin.ne.rmax .AND. maxdif.ne.0) then
 100         !
 101         !  Determine which algorithm to use based on user-supplied
 102         !  binary scale factor and number of bits.
 103         !
 104         if (idrstmpl(2).eq.0) then
 105            !
 106            !  No binary scaling and calculate minimum number of
 107            !  bits in which the data will fit.
 108            !
 109            imin=nint(rmin*dscale)
 110            imax=nint(rmax*dscale)
 111            maxdif=imax-imin
 112            temp=alog(real(maxdif+1))/alog(2.0)
 113            nbits=ceiling(temp)
 114            rmin=real(imin)
 115            !   scale data
 116            do j=1,ndpts
 117              ifld(j)=nint(fld(j)*dscale)-imin
 118            enddo
 119         else
 120            !
 121            !  Use binary scaling factor and calculate minimum number of
 122            !  bits in which the data will fit.
 123            !
 124            rmin=rmin*dscale
 125            rmax=rmax*dscale
 126            maxdif=nint((rmax-rmin)*bscale)
 127            temp=alog(real(maxdif+1))/alog(2.0)
 128            nbits=ceiling(temp)
 129            !   scale data
 130            do j=1,ndpts
 131              ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
 132            enddo
 133         endif
 134         !
 135         !  Pack data into full octets, then do JPEG2000 encode.
 136         !  and calculate the length of the packed data in bytes
 137         !
 138         retry=0
 139         nbytes=(nbits+7)/8
 140         nsize=lcpack      ! needed for input to enc_jpeg2000
 141         allocate(ctemp(nbytes*ndpts))
 142         call sbytes(ctemp,ifld,0,nbytes*8,0,ndpts)
 143         lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6),
 144      &                      idrstmpl(7),retry,cpack,nsize)
 145         if (lcpack.le.0) then
 146            print *,'jpcpack: ERROR Packing JPC=',lcpack
 147            if (lcpack.eq.-3) then
 148               retry=1
 149               print *,'jpcpack: Retrying....'
 150               lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6),
 151      &                         idrstmpl(7),retry,cpack,nsize)
 152               if (lcpack.le.0) then
 153                  print *,'jpcpack: Retry Failed.'
 154               else
 155                  print *,'jpcpack: Retry Successful.'
 156               endif
 157            endif
 158         endif
 159         deallocate(ctemp)
 160
 161       else
 162         nbits=0
 163         lcpack=0
 164       endif
 165
 166 !
 167 !  Fill in ref value and number of bits in Template 5.0
 168 !
 169       rmin4 = rmin
 170       call mkieee(rmin4,ref,1)   ! ensure reference value is IEEE format
 171 !      call gbyte(ref,idrstmpl(1),0,32)
 172       iref=transfer(ref,iref)
 173       idrstmpl(1)=iref
 174       idrstmpl(4)=nbits
 175       idrstmpl(5)=0         ! original data were reals
 176       if (idrstmpl(6).eq.0) idrstmpl(7)=255       ! lossy not used
 177
 178 #endif /* USE_JPEG2000 */
 179       return
 180       end