external/io_grib2/g2lib/jpcpack.F

   1       subroutine jpcpack(fld,width,height,idrstmpl,cpack,lcpack,ierr)
   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       integer,intent(out)   :: ierr
  67
  68       real(4) :: ref
  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       integer :: orig_dscalefct
  75       integer :: orig_bscalefct
  76       integer(8) :: maxdif,imin,imax
  77
  78       ierr = 0
  79
  80       ndpts=width*height
  81       bscale=2.0**real(-idrstmpl(2))
  82       dscale=10.0**real(idrstmpl(3))
  83
  84       orig_bscalefct = idrstmpl(2)
  85       orig_dscalefct = idrstmpl(3)
  86
  87 !
  88 !  Find max and min values in the data
  89 !
  90       rmax=fld(1)
  91       rmin=fld(1)
  92       do j=2,ndpts
  93         if (fld(j).gt.rmax) rmax=fld(j)
  94         if (fld(j).lt.rmin) rmin=fld(j)
  95       enddo
  96       if (idrstmpl(2).eq.0) then
  97          maxdif=nint(rmax*dscale,8)-nint(rmin*dscale,8)
  98       else
  99          maxdif=nint((rmax-rmin)*dscale*bscale,8)
 100       endif
 101 !
 102 !  If max and min values are not equal, pack up field.
 103 !  If they are equal, we have a constant field, and the reference
 104 !  value (rmin) is the value for each point in the field and
 105 !  set nbits to 0.
 106 !
 107       if (rmin.ne.rmax .AND. maxdif.ne.0) then
 108         !
 109         !  Determine which algorithm to use based on user-supplied
 110         !  binary scale factor and number of bits.
 111         !
 112
 113         if (idrstmpl(2).eq.0) then
 114            !
 115            !  No binary scaling and calculate minimum number of
 116            !  bits in which the data will fit.
 117            !
 118
 119            nbits = 25
 120
 121            do while (nbits .gt. 24)
 122               imin=nint(rmin*dscale,8)
 123               imax=nint(rmax*dscale,8)
 124               maxdif=imax-imin
 125               temp=alog(real(maxdif+1))/alog(2.0)
 126               nbits=ceiling(temp)
 127               !
 128               ! These lines assure that we do no overflow
 129               !
 130               !  Todd Hutchinson, WSI 9/23/05
 131               !
 132               if (nbits .le. 24) then
 133                  exit
 134               else
 135                  idrstmpl(3) = idrstmpl(3) - 1
 136                  dscale=10.0**real(idrstmpl(3))
 137               endif
 138            enddo
 139
 140            rmin=real(imin)
 141            !   scale data
 142            do j=1,ndpts
 143              ifld(j)=nint(fld(j)*dscale)-imin
 144            enddo
 145         else
 146            !
 147            !  Use binary scaling factor and calculate minimum number of
 148            !  bits in which the data will fit.
 149            !
 150            nbits = 25
 151
 152            do while (nbits .gt. 24)
 153               rmin=rmin*dscale
 154               rmax=rmax*dscale
 155               maxdif=nint((rmax-rmin)*bscale)
 156               temp=alog(real(maxdif+1))/alog(2.0)
 157               nbits=ceiling(temp)
 158               !
 159               ! These lines assure that we do no overflow
 160               !
 161               !  Todd Hutchinson, WSI 9/23/05
 162               !
 163               if (nbits .le. 24) then
 164                  exit
 165               else
 166                  idrstmpl(2) = idrstmpl(2) + 1
 167                  bscale=2.0**real(-idrstmpl(2))
 168               endif
 169            enddo
 170            !   scale data
 171            do j=1,ndpts
 172              ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
 173            enddo
 174         endif
 175
 176         if (idrstmpl(2) .ne. orig_bscalefct) then
 177            write(6,'(A,I2,A,I2,A)')
 178      &          ' JPCPACK: Reduced binary scale fctr from ',
 179      &          orig_bscalefct,' to ',idrstmpl(2),
 180      &          ' to prevent overflow'
 181            ierr = 12
 182         endif
 183
 184         if (idrstmpl(3) .ne. orig_dscalefct) then
 185            write(6,'(A,I2,A,I2,A)')
 186      &          ' JPCPACK: Reduced decimal scale fctr from ',
 187      &          orig_dscalefct,' to ',idrstmpl(3),
 188      &          ' to prevent overflow'
 189            ierr = 11
 190         endif
 191
 192         !
 193         !  Pack data into full octets, then do JPEG2000 encode.
 194         !  and calculate the length of the packed data in bytes
 195         !
 196         retry=0
 197         nbytes=(nbits+7)/8
 198         nsize=lcpack      ! needed for input to enc_jpeg2000
 199         allocate(ctemp(nbytes*ndpts))
 200         call g2lib_sbytes(ctemp,ifld,0,nbytes*8,0,ndpts)
 201         lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6),
 202      &                      idrstmpl(7),retry,cpack,nsize)
 203         if (lcpack.le.0) then
 204            print *,'jpcpack: ERROR Packing JPC=',lcpack
 205            if (lcpack.eq.-3) then
 206               retry=1
 207               print *,'jpcpack: Retrying....'
 208               lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6),
 209      &                         idrstmpl(7),retry,cpack,nsize)
 210               if (lcpack.le.0) then
 211                  print *,'jpcpack: Retry Failed.'
 212               else
 213                  print *,'jpcpack: Retry Successful.'
 214               endif
 215            endif
 216         endif
 217         deallocate(ctemp)
 218
 219       else
 220         nbits=0
 221         lcpack=0
 222       endif
 223
 224 !
 225 !  Fill in ref value and number of bits in Template 5.0
 226 !
 227       call mkieee(rmin,ref,1)   ! ensure reference value is IEEE format
 228 !      call g2lib_gbyte(ref,idrstmpl(1),0,32)
 229       iref=transfer(ref,iref)
 230       idrstmpl(1)=iref
 231       idrstmpl(4)=nbits
 232       idrstmpl(5)=0         ! original data were reals
 233       if (idrstmpl(6).eq.0) idrstmpl(7)=255       ! lossy not used
 234
 235       return
 236       end