WPS/ungrib/src/ngl/g2/simpack.f

   1       subroutine simpack(fld,ndpts,idrstmpl,cpack,lcpack)
   2 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   3 !                .      .    .                                       .
   4 ! SUBPROGRAM:    simpack
   5 !   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2000-06-21
   6 !
   7 ! ABSTRACT: This subroutine packs up a data field using a simple
   8 !   packing algorithm as defined in the GRIB2 documention.  It
   9 !   also fills in GRIB2 Data Representation Template 5.0 with the
  10 !   appropriate values.
  11 !
  12 ! PROGRAM HISTORY LOG:
  13 ! 2000-06-21  Gilbert
  14 !
  15 ! USAGE:    CALL simpack(fld,ndpts,idrstmpl,cpack,lcpack)
  16 !   INPUT ARGUMENT LIST:
  17 !     fld()    - Contains the data values to pack
  18 !     ndpts    - The number of data values in array fld()
  19 !     idrstmpl - Contains the array of values for Data Representation
  20 !                Template 5.0
  21 !                (1) = Reference value - ignored on input
  22 !                (2) = Binary Scale Factor
  23 !                (3) = Decimal Scale Factor
  24 !                (4) = Number of bits used to pack data, if value is
  25 !                      > 0 and  <= 31.
  26 !                      If this input value is 0 or outside above range
  27 !                      then the num of bits is calculated based on given
  28 !                      data and scale factors.
  29 !                (5) = Original field type - currently ignored on input
  30 !                      Data values assumed to be reals.
  31 !
  32 !   OUTPUT ARGUMENT LIST:
  33 !     idrstmpl - Contains the array of values for Data Representation
  34 !                Template 5.0
  35 !                (1) = Reference value - set by simpack routine.
  36 !                (2) = Binary Scale Factor - unchanged from input
  37 !                (3) = Decimal Scale Factor - unchanged from input
  38 !                (4) = Number of bits used to pack data, unchanged from
  39 !                      input if value is between 0 and 31.
  40 !                      If this input value is 0 or outside above range
  41 !                      then the num of bits is calculated based on given
  42 !                      data and scale factors.
  43 !                (5) = Original field type - currently set = 0 on output.
  44 !                      Data values assumed to be reals.
  45 !     cpack    - The packed data field (character*1 array)
  46 !     lcpack   - length of packed field cpack().
  47 !
  48 ! REMARKS: None
  49 !
  50 ! ATTRIBUTES:
  51 !   LANGUAGE: XL Fortran 90
  52 !   MACHINE:  IBM SP
  53 !
  54 !$$$
  55
  56       integer,intent(in) :: ndpts
  57       real,intent(in) :: fld(ndpts)
  58       character(len=1),intent(out) :: cpack(*)
  59       integer,intent(inout) :: idrstmpl(*)
  60       integer,intent(out) :: lcpack
  61
  62       real(4) :: ref
  63       integer(4) :: iref
  64       integer :: ifld(ndpts)
  65       integer,parameter :: zero=0
  66
  67       bscale=2.0**real(-idrstmpl(2))
  68       dscale=10.0**real(idrstmpl(3))
  69       if (idrstmpl(4).le.0.OR.idrstmpl(4).gt.31) then
  70          nbits=0
  71       else
  72          nbits=idrstmpl(4)
  73       endif
  74 !
  75 !  Find max and min values in the data
  76 !
  77       rmax=fld(1)
  78       rmin=fld(1)
  79       do j=2,ndpts
  80         if (fld(j).gt.rmax) rmax=fld(j)
  81         if (fld(j).lt.rmin) rmin=fld(j)
  82       enddo
  83 !
  84 !  If max and min values are not equal, pack up field.
  85 !  If they are equal, we have a constant field, and the reference
  86 !  value (rmin) is the value for each point in the field and
  87 !  set nbits to 0.
  88 !
  89       if (rmin.ne.rmax) then
  90         !
  91         !  Determine which algorithm to use based on user-supplied
  92         !  binary scale factor and number of bits.
  93         !
  94         if (nbits.eq.0.AND.idrstmpl(2).eq.0) then
  95            !
  96            !  No binary scaling and calculate minumum number of
  97            !  bits in which the data will fit.
  98            !
  99            imin=nint(rmin*dscale)
 100            imax=nint(rmax*dscale)
 101            maxdif=imax-imin
 102            temp=alog(real(maxdif+1))/alog(2.0)
 103            nbits=ceiling(temp)
 104            rmin=real(imin)
 105            !   scale data
 106            do j=1,ndpts
 107              ifld(j)=nint(fld(j)*dscale)-imin
 108            enddo
 109         elseif (nbits.ne.0.AND.idrstmpl(2).eq.0) then
 110            !
 111            !  Use minimum number of bits specified by user and
 112            !  adjust binary scaling factor to accomodate data.
 113            !
 114            rmin=rmin*dscale
 115            rmax=rmax*dscale
 116            maxnum=(2**nbits)-1
 117            temp=alog(real(maxnum)/(rmax-rmin))/alog(2.0)
 118            idrstmpl(2)=ceiling(-1.0*temp)
 119            bscale=2.0**real(-idrstmpl(2))
 120            !   scale data
 121            do j=1,ndpts
 122              ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
 123            enddo
 124         elseif (nbits.eq.0.AND.idrstmpl(2).ne.0) then
 125            !
 126            !  Use binary scaling factor and calculate minumum number of
 127            !  bits in which the data will fit.
 128            !
 129            rmin=rmin*dscale
 130            rmax=rmax*dscale
 131            maxdif=nint((rmax-rmin)*bscale)
 132            temp=alog(real(maxdif+1))/alog(2.0)
 133            nbits=ceiling(temp)
 134            !   scale data
 135            do j=1,ndpts
 136              ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
 137            enddo
 138         elseif (nbits.ne.0.AND.idrstmpl(2).ne.0) then
 139            !
 140            !  Use binary scaling factor and use minumum number of
 141            !  bits specified by user.   Dangerous - may loose
 142            !  information if binary scale factor and nbits not set
 143            !  properly by user.
 144            !
 145            rmin=rmin*dscale
 146            !   scale data
 147            do j=1,ndpts
 148              ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
 149            enddo
 150         endif
 151         !
 152         !  Pack data, Pad last octet with Zeros, if necessary,
 153         !  and calculate the length of the packed data in bytes
 154         !
 155         call sbytes(cpack,ifld,0,nbits,0,ndpts)
 156         nbittot=nbits*ndpts
 157         left=8-mod(nbittot,8)
 158         if (left.ne.8) then
 159           call sbyte(cpack,zero,nbittot,left)    ! Pad with zeros to fill Octet
 160           nbittot=nbittot+left
 161         endif
 162         lcpack=nbittot/8
 163
 164       else
 165         nbits=0
 166         lcpack=0
 167       endif
 168
 169 !
 170 !  Fill in ref value and number of bits in Template 5.0
 171 !
 172       call mkieee(rmin,ref,1)   ! ensure reference value is IEEE format
 173       !print *,'SAGref = ',rmin,ref
 174 !      call gbyte(ref,idrstmpl(1),0,32)
 175       iref=transfer(ref,iref)
 176       idrstmpl(1)=iref
 177       idrstmpl(4)=nbits
 178       idrstmpl(5)=0         ! original data were reals
 179
 180       return
 181       end