g2sim.F90

Go to the documentation of this file.

 
 
subroutine simpack(fld,ndpts,idrstmpl,cpack,lcpack)
  use intmath
  implicit none
 
  integer,intent(in) :: ndpts
  real,intent(in) :: fld(ndpts)
  character(len=1),intent(out) :: cpack(*)
  integer,intent(inout) :: idrstmpl(*)
  integer,intent(out) :: lcpack
 
  real(4) :: ref,rmin4
 
  integer(4) :: iref
  integer :: ifld(ndpts)
  integer,parameter :: zero=0
  integer :: nbittot, nbits, maxnum, maxdif, left
  integer :: imax, imin, j
  real :: rmax, rmin, temp, dscale, bscale
 
  bscale=2.0**real(-idrstmpl(2))
  dscale=10.0**real(idrstmpl(3))
  if (idrstmpl(4).le.0.OR.idrstmpl(4).gt.31) then
     nbits=0
  else
     nbits=idrstmpl(4)
  endif
 
  !     Find max and min values in the data
  if(ndpts<1) then
     rmin=0
     rmax=0
  else
     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
  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) then
 
     !  Determine which algorithm to use based on user-supplied 
     !  binary scale factor and number of bits.
     if (nbits.eq.0.AND.idrstmpl(2).eq.0) then
 
        !  No binary scaling and calculate minumum number of 
        !  bits in which the data will fit.
        imin=nint(rmin*dscale)
        imax=nint(rmax*dscale)
        maxdif=imax-imin
        temp=i1log2(maxdif+1)
        nbits=ceiling(temp)
        rmin=real(imin)
        !   scale data
        do j=1,ndpts
           ifld(j)=nint(fld(j)*dscale)-imin
        enddo
     elseif (nbits.ne.0.AND.idrstmpl(2).eq.0) then
 
        !  Use minimum number of bits specified by user and
        !  adjust binary scaling factor to accomodate data.
        rmin=rmin*dscale
        rmax=rmax*dscale
        maxnum=(2**nbits)-1
        temp=ilog2(nint(real(maxnum)/(rmax-rmin)))
        idrstmpl(2)=ceiling(-1.0*temp)
        bscale=2.0**real(-idrstmpl(2))
        !   scale data
        do j=1,ndpts
           ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
        enddo
     elseif (nbits.eq.0.AND.idrstmpl(2).ne.0) then
 
        !  Use binary scaling factor and calculate minumum number of 
        !  bits in which the data will fit.
        rmin=rmin*dscale
        rmax=rmax*dscale
        maxdif=nint((rmax-rmin)*bscale)
        temp=i1log2(maxdif)
        nbits=ceiling(temp)
        !   scale data
        do j=1,ndpts
           ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
        enddo
     elseif (nbits.ne.0.AND.idrstmpl(2).ne.0) then
 
        !  Use binary scaling factor and use minumum number of 
        !  bits specified by user.   Dangerous - may loose
        !  information if binary scale factor and nbits not set
        !  properly by user.
        rmin=rmin*dscale
        !   scale data
        do j=1,ndpts
           ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
        enddo
     endif
 
     !  Pack data, Pad last octet with Zeros, if necessary,
     !  and calculate the length of the packed data in bytes
     call g2_sbytesc(cpack,ifld,0,nbits,0,ndpts)
     nbittot=nbits*ndpts
     left=8-mod(nbittot,8)
     if (left.ne.8) then
        call g2_sbytec(cpack,zero,nbittot,left)    ! Pad with zeros to fill Octet
        nbittot=nbittot+left
     endif
     lcpack=nbittot/8
 
  else
     !print *,'nbits 0'
     nbits=0
     lcpack=0
  endif
 
 
  !     Fill in ref value and number of bits in Template 5.0
  rmin4 = real(rmin, 4)
  call mkieee(rmin4,ref,1)   ! ensure reference value is IEEE format
  iref=transfer(ref,iref)
  idrstmpl(1)=iref
  idrstmpl(4)=nbits
  idrstmpl(5)=0         ! original data were reals
 
end subroutine simpack
 
subroutine simunpack(cpack, len, idrstmpl, ndpts, fld)
  implicit none
  
  character(len=1), intent(in) :: cpack(len)
  integer, intent(in) :: ndpts, len
  integer, intent(in) :: idrstmpl(*)
  real, intent(out) :: fld(ndpts)
 
  integer :: ifld(ndpts)
  integer(4) :: ieee
  real :: ref, bscale, dscale
  integer :: itype, j, nbits
 
  ieee = idrstmpl(1)
  call rdieee(ieee, ref, 1)
  bscale = 2.0**real(idrstmpl(2))
  dscale = 10.0**real(-idrstmpl(3))
  nbits = idrstmpl(4)
  itype = idrstmpl(5)
 
  ! If nbits equals 0,  we have a constant field where the reference value
  ! is the data value at each gridpoint.
  if (nbits .ne. 0) then
     call g2_gbytesc(cpack, ifld, 0, nbits, 0, ndpts)
     do j=1, ndpts
        fld(j) = ((real(ifld(j)) * bscale) + ref) * dscale
     enddo
  else
     !print *, 'unpack ref ', ref
     !print *, 'unpack ndpts ', ndpts
     do j=1, ndpts
        fld(j) = ref
     enddo
  endif
end subroutine simunpack