compack.F90

Go to the documentation of this file.

 
 
subroutine cmplxpack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
 
  integer,intent(in) :: ndpts,idrsnum
  real,intent(in) :: fld(ndpts)
  character(len=1),intent(out) :: cpack(*)
  integer,intent(inout) :: idrstmpl(*)
  integer,intent(out) :: lcpack
 
  if ( idrstmpl(7) .eq. 0 ) then       ! No internal missing values
     call compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
  elseif ( idrstmpl(7).eq.1 .OR. idrstmpl(7).eq.2) then
     call misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
  else
     print *,'cmplxpack: Do not recognize Missing value option.'
     lcpack=-1
  endif
 
  return
end subroutine cmplxpack
 
subroutine compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
 
  use intmath
  implicit none
  integer,intent(in) :: ndpts,idrsnum
  real,intent(in) :: fld(ndpts)
  character(len=1),intent(out) :: cpack(*)
  integer,intent(inout) :: idrstmpl(*)
  integer,intent(out) :: lcpack
 
  real :: bscale,dscale
  integer :: j,iofst,imin,ival1,ival2,minsd,nbitsd,n
  integer :: igmax,nbitsgref,left,iwmax,i,ilmax,kk,ij
  integer :: ngwidthref,nbitsgwidth,nglenref,nglenlast
  integer :: maxorig,nbitorig,isd,ngroups,itemp,minpk
  integer :: inc,maxgrps,missopt,miss1,miss2,lg
  integer :: ibit,jbit,kbit,novref,lbitref,ier,ng,imax
  integer :: nbitsglen
  real(4) :: ref,rmin4
  real(8) :: rmin,rmax
 
  integer(4) :: iref
  integer,allocatable :: ifld(:)
  integer,allocatable :: jmin(:),jmax(:),lbit(:)
  integer,parameter :: zero=0
  integer,allocatable :: gref(:),gwidth(:),glen(:)
  integer :: glength,grpwidth
  logical :: simple_alg
 
  simple_alg = .false.
 
  bscale=2.0**real(-idrstmpl(2))
  dscale=10.0**real(idrstmpl(3))
  !
  !  Find max and min values in the data
  !
  if(ndpts>0) then
     rmax=fld(1)
     rmin=fld(1)
  else
     rmax=1.0
     rmin=1.0
  endif
  do j=2,ndpts
     if (fld(j).gt.rmax) rmax=fld(j)
     if (fld(j).lt.rmin) rmin=fld(j)
  enddo
  !
  !  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.
  !
  multival: if (rmin.ne.rmax) then
     iofst=0
     allocate(ifld(ndpts))
     allocate(gref(ndpts))
     allocate(gwidth(ndpts))
     allocate(glen(ndpts))
     !
     !  Scale original data
     !
     if (idrstmpl(2).eq.0) then        !  No binary scaling
        imin=nint(rmin*dscale)
        !imax=nint(rmax*dscale)
        rmin=real(imin)
        do j=1,ndpts
           ifld(j)=max(0,nint(fld(j)*dscale)-imin)
        enddo
     else                              !  Use binary scaling factor
        rmin=rmin*dscale
        !rmax=rmax*dscale
        do j=1,ndpts
           ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
        enddo
     endif
     !
     !  Calculate Spatial differences, if using DRS Template 5.3
     !
     alg3: if (idrsnum.eq.3) then        ! spatial differences
        if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2
        if (idrstmpl(17).eq.1) then      ! first order
           ival1=ifld(1)
           if(ival1<0) then
              print *,'G2: negative ival1',ival1
              stop 101
           endif
           do j=ndpts,2,-1
              ifld(j)=ifld(j)-ifld(j-1)
           enddo
           ifld(1)=0
        elseif (idrstmpl(17).eq.2) then      ! second order
           ival1=ifld(1)
           ival2=ifld(2)
           if(ival1<0) then
              print *,'G2: negative ival1',ival1
              stop 11
           endif
           if(ival2<0) then
              print *,'G2: negative ival2',ival2
              stop 12
           endif
           do j=ndpts,3,-1
              ifld(j)=ifld(j)-(2*ifld(j-1))+ifld(j-2)
           enddo
           ifld(1)=0
           ifld(2)=0
        endif
        !
        !  subtract min value from spatial diff field
        !
        isd=idrstmpl(17)+1
        minsd=minval(ifld(isd:ndpts))
        do j=isd,ndpts
           ifld(j)=ifld(j)-minsd
        enddo
        !
        !   find num of bits need to store minsd and add 1 extra bit
        !   to indicate sign
        !
        nbitsd=i1log2(abs(minsd))+1
        !
        !   find num of bits need to store ifld(1) ( and ifld(2)
        !   if using 2nd order differencing )
        !
        maxorig=ival1
        if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2
        nbitorig=i1log2(maxorig)+1
        if (nbitorig.gt.nbitsd) nbitsd=nbitorig
        !   increase number of bits to even multiple of 8 ( octet )
        if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8))
        !
        !  Store extra spatial differencing info into the packed
        !  data section.
        !
        if (nbitsd.ne.0) then
           !   pack first original value
           if (ival1.ge.0) then
              call g2_sbytec(cpack,ival1,iofst,nbitsd)
              iofst=iofst+nbitsd
           else
              call g2_sbytec(cpack,1,iofst,1)
              iofst=iofst+1
              call g2_sbytec(cpack,iabs(ival1),iofst,nbitsd-1)
              iofst=iofst+nbitsd-1
           endif
           if (idrstmpl(17).eq.2) then
              !  pack second original value
              if (ival2.ge.0) then
                 call g2_sbytec(cpack,ival2,iofst,nbitsd)
                 iofst=iofst+nbitsd
              else
                 call g2_sbytec(cpack,1,iofst,1)
                 iofst=iofst+1
                 call g2_sbytec(cpack,iabs(ival2),iofst,nbitsd-1)
                 iofst=iofst+nbitsd-1
              endif
           endif
           !  pack overall min of spatial differences
           if (minsd.ge.0) then
              call g2_sbytec(cpack,minsd,iofst,nbitsd)
              iofst=iofst+nbitsd
           else
              call g2_sbytec(cpack,1,iofst,1)
              iofst=iofst+1
              call g2_sbytec(cpack,iabs(minsd),iofst,nbitsd-1)
              iofst=iofst+nbitsd-1
           endif
        endif
        !print *,'SDp ',ival1,ival2,minsd,nbitsd
     endif alg3              !  end of spatial diff section
     !
     !   Determine Groups to be used.
     !
     simplealg: if ( simple_alg ) then
        !  set group length to 10 :  calculate number of groups
        !  and length of last group
        print *,'G2: use simple algorithm'
        ngroups=ndpts/10
        glen(1:ngroups)=10
        itemp=mod(ndpts,10)
        if (itemp.ne.0) then
           ngroups=ngroups+1
           glen(ngroups)=itemp
        endif
     else
        ! Use Dr. Glahn's algorithm for determining grouping.
        !
        minpk=10
        inc=1
        maxgrps=((ndpts+minpk-1)/minpk)
        allocate(jmin(maxgrps))
        allocate(jmax(maxgrps))
        allocate(lbit(maxgrps))
        missopt=0
        call pack_gp(ifld,ndpts,missopt,minpk,inc,miss1,miss2, &
             jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit, &
             kbit,novref,lbitref,ier)
        if(ier/=0) then
           ! Dr. Glahn's algorithm failed; use simple packing method instead.
1099       format('G2: fall back to simple algorithm (glahn ier=',i0,')')
           print 1099,ier
           ngroups=ndpts/10
           glen(1:ngroups)=10
           itemp=mod(ndpts,10)
           if (itemp.ne.0) then
              ngroups=ngroups+1
              glen(ngroups)=itemp
           endif
        elseif(ngroups<1) then
           ! Dr. Glahn's algorithm failed; use simple packing method instead.
           print *,'Glahn algorithm failed; use simple packing'
           ngroups=ndpts/10
           glen(1:ngroups)=10
           itemp=mod(ndpts,10)
           if (itemp.ne.0) then
              ngroups=ngroups+1
              glen(ngroups)=itemp
           endif
        else
           !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref
           do ng=1,ngroups
              glen(ng)=glen(ng)+novref
           enddo
           deallocate(jmin)
           deallocate(jmax)
           deallocate(lbit)
        endif
     endif simplealg
     !  
     !  For each group, find the group's reference value
     !  and the number of bits needed to hold the remaining values
     !
     n=1
     do ng=1,ngroups
        !    find max and min values of group
        gref(ng)=ifld(n)
        imax=ifld(n)
        j=n+1
        do lg=2,glen(ng)
           if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) 
           if (ifld(j).gt.imax) imax=ifld(j) 
           j=j+1
        enddo
        !   calc num of bits needed to hold data
        if ( gref(ng).ne.imax ) then
           gwidth(ng)=i1log2(imax-gref(ng))
        else
           gwidth(ng)=0
        endif
        !   Subtract min from data
        j=n
        do lg=1,glen(ng)
           ifld(j)=ifld(j)-gref(ng)
           j=j+1
        enddo
        !   increment fld array counter
        n=n+glen(ng)
     enddo
     !  
     !  Find max of the group references and calc num of bits needed 
     !  to pack each groups reference value, then
     !  pack up group reference values
     !
     !write(77,*)'GREFS: ',(gref(j),j=1,ngroups)
     igmax=maxval(gref(1:ngroups))
     if (igmax.ne.0) then
        nbitsgref=i1log2(igmax)
        call g2_sbytesc(cpack,gref,iofst,nbitsgref,0,ngroups)
        itemp=nbitsgref*ngroups
        iofst=iofst+itemp
        if (mod(itemp,8).ne.0) then
           left=8-mod(itemp,8)
           call g2_sbytec(cpack,zero,iofst,left)
           iofst=iofst+left
        endif
     else
        nbitsgref=0
     endif
     !
     !  Find max/min of the group widths and calc num of bits needed
     !  to pack each groups width value, then
     !  pack up group width values
     !
     !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
     iwmax=maxval(gwidth(1:ngroups))
     ngwidthref=minval(gwidth(1:ngroups))
     if (iwmax.ne.ngwidthref) then
        nbitsgwidth=i1log2(iwmax-ngwidthref)
        do i=1,ngroups
           gwidth(i)=gwidth(i)-ngwidthref
           if(gwidth(i)<0) then
              write(0,*) 'i,gw,ngw=',i,gwidth(i),ngwidthref
              stop 9
           endif
        enddo
        call g2_sbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
        itemp=nbitsgwidth*ngroups
        iofst=iofst+itemp
        !     Pad last octet with Zeros, if necessary,
        if (mod(itemp,8).ne.0) then
           left=8-mod(itemp,8)
           call g2_sbytec(cpack,zero,iofst,left)
           iofst=iofst+left
        endif
     else
        nbitsgwidth=0
        gwidth(1:ngroups)=0
     endif
     !
     !  Find max/min of the group lengths and calc num of bits needed
     !  to pack each groups length value, then
     !  pack up group length values
     !
     !write(77,*)'GLENS: ',(glen(j),j=1,ngroups)
     ilmax=maxval(glen(1:ngroups-1))
     nglenref=minval(glen(1:ngroups-1))
     nglenlast=glen(ngroups)
     if (ilmax.ne.nglenref) then
        nbitsglen=i1log2(ilmax-nglenref)
        do i=1,ngroups-1
           glen(i)=glen(i)-nglenref
           if(glen(i)<0) then
              write(0,*) 'i,glen(i) = ',i,glen(i)
              stop 23
           endif
        enddo
        call g2_sbytesc(cpack,glen,iofst,nbitsglen,0,ngroups)
        itemp=nbitsglen*ngroups
        iofst=iofst+itemp
        !     Pad last octet with Zeros, if necessary,
        if (mod(itemp,8).ne.0) then
           left=8-mod(itemp,8)
           call g2_sbytec(cpack,zero,iofst,left)
           iofst=iofst+left
        endif
     else
        nbitsglen=0
        glen(1:ngroups)=0
     endif
     !
     !  For each group, pack data values
     !
     !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
     n=1
     ij=0
     do ng=1,ngroups
        glength=glen(ng)+nglenref
        if (ng.eq.ngroups ) glength=nglenlast
        grpwidth=gwidth(ng)+ngwidthref
        !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
        if ( grpwidth.ne.0 ) then
           call g2_sbytesc(cpack,ifld(n),iofst,grpwidth,0,glength)
           iofst=iofst+(grpwidth*glength)
        endif
        do kk=1,glength
           ij=ij+1
           !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
        enddo
        n=n+glength
     enddo
     !     Pad last octet with Zeros, if necessary,
     if (mod(iofst,8).ne.0) then
        left=8-mod(iofst,8)
        call g2_sbytec(cpack,zero,iofst,left)
        iofst=iofst+left
     endif
     lcpack=iofst/8
     !
     if ( allocated(ifld) ) deallocate(ifld)
     if ( allocated(gref) ) deallocate(gref)
     if ( allocated(gwidth) ) deallocate(gwidth)
     if ( allocated(glen) ) deallocate(glen)
  else           !   Constant field ( max = min )
     lcpack=0
     nbitsgref=0
     ngroups=0
     ngwidthref=0
     nbitsgwidth=0
     nglenref=0
     nglenlast=0
     nbitsglen=0
     nbitsd=0
  endif multival
 
  !
  !  Fill in ref value and number of bits in Template 5.2
  !
  rmin4 = real(rmin, 4)
  call mkieee(rmin4,ref,1)   ! ensure reference value is IEEE format
  !    call g2_gbytec(ref,idrstmpl(1),0,32)
  iref=transfer(ref,iref)
  idrstmpl(1)=iref
  idrstmpl(4)=nbitsgref
  idrstmpl(5)=0         ! original data were reals
  idrstmpl(6)=1         ! general group splitting
  idrstmpl(7)=0         ! No internal missing values
  idrstmpl(8)=0         ! Primary missing value
  idrstmpl(9)=0         ! secondary missing value
  idrstmpl(10)=ngroups          ! Number of groups
  idrstmpl(11)=ngwidthref       ! reference for group widths
  idrstmpl(12)=nbitsgwidth      ! num bits used for group widths
  idrstmpl(13)=nglenref         ! Reference for group lengths
  idrstmpl(14)=1                ! length increment for group lengths
  idrstmpl(15)=nglenlast        ! True length of last group
  idrstmpl(16)=nbitsglen        ! num bits used for group lengths
  if (idrsnum.eq.3) then
     idrstmpl(18)=nbitsd/8      ! num bits used for extra spatial
     ! differencing values
  endif
 
end subroutine compack
 
subroutine comunpack(cpack,len,lensec,idrsnum,idrstmpl,ndpts, &
     fld,ier)
 
  character(len=1),intent(in) :: cpack(len)
  integer,intent(in) :: ndpts,len
  integer,intent(in) :: idrstmpl(*)
  real,intent(out) :: fld(ndpts)
 
  integer,allocatable :: ifld(:),ifldmiss(:)
  integer(4) :: ieee
  integer,allocatable :: gref(:),gwidth(:),glen(:)
  real :: ref,bscale,dscale,rmiss1,rmiss2
  !      real :: fldo(6045)
  integer :: totBit, totLen
 
  ier=0
  !print *,'IDRSTMPL: ',(idrstmpl(j),j=1,16)
  ieee = idrstmpl(1)
  call rdieee(ieee,ref,1)
  bscale = 2.0**real(idrstmpl(2))
  dscale = 10.0**real(-idrstmpl(3))
  nbitsgref = idrstmpl(4)
  itype = idrstmpl(5)
  ngroups = idrstmpl(10)
  nbitsgwidth = idrstmpl(12)
  nbitsglen = idrstmpl(16)
  if (idrsnum.eq.3) then
     nbitsd=idrstmpl(18)*8
  endif
 
  !   Constant field
 
  if (ngroups.eq.0) then
     do j=1,ndpts
        fld(j)=ref
     enddo
     return
  endif
 
  iofst=0
  allocate(ifld(ndpts),stat=is)
  !print *,'ALLOC ifld: ',is,ndpts
  allocate(gref(ngroups),stat=is)
  !print *,'ALLOC gref: ',is
  allocate(gwidth(ngroups),stat=is)
  !print *,'ALLOC gwidth: ',is
  !
  !  Get missing values, if supplied
  !
  if ( idrstmpl(7).eq.1 ) then
     if (itype.eq.0) then
        call rdieee(idrstmpl(8),rmiss1,1)
     else
        rmiss1=real(idrstmpl(8))
     endif
  elseif ( idrstmpl(7).eq.2 ) then
     if (itype.eq.0) then
        call rdieee(idrstmpl(8),rmiss1,1)
        call rdieee(idrstmpl(9),rmiss2,1)
     else
        rmiss1=real(idrstmpl(8))
        rmiss2=real(idrstmpl(9))
     endif
  endif
  !print *,'RMISSs: ',rmiss1,rmiss2,ref
  ! 
  !  Extract Spatial differencing values, if using DRS Template 5.3
  !
  if (idrsnum.eq.3) then
     if (nbitsd.ne.0) then
        call g2_gbytec(cpack,ival1,iofst,nbitsd)
        iofst=iofst+nbitsd
        if (idrstmpl(17).eq.2) then
           call g2_gbytec(cpack,ival2,iofst,nbitsd)
           iofst=iofst+nbitsd
        endif
        call g2_gbytec(cpack,isign,iofst,1)
        iofst=iofst+1
        call g2_gbytec(cpack,minsd,iofst,nbitsd-1)
        iofst=iofst+nbitsd-1
        if (isign.eq.1) minsd=-minsd
     else
        ival1=0
        ival2=0
        minsd=0
     endif
     !print *,'SDu ',ival1,ival2,minsd,nbitsd
  endif
  !
  !  Extract Each Group's reference value
  !
  !print *,'SAG1: ',nbitsgref,ngroups,iofst
  if (nbitsgref.ne.0) then
     call g2_gbytesc(cpack,gref,iofst,nbitsgref,0,ngroups)
     itemp=nbitsgref*ngroups
     iofst=iofst+(itemp)
     if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8))
  else
     gref(1:ngroups)=0
  endif
  !write(78,*)'GREFs: ',(gref(j),j=1,ngroups)
  !
  !  Extract Each Group's bit width
  !
  !print *,'SAG2: ',nbitsgwidth,ngroups,iofst,idrstmpl(11)
  if (nbitsgwidth.ne.0) then
     call g2_gbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
     itemp=nbitsgwidth*ngroups
     iofst=iofst+(itemp)
     if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8))
  else
     gwidth(1:ngroups)=0
  endif
  do j=1,ngroups
     gwidth(j)=gwidth(j)+idrstmpl(11)
  enddo
  !write(78,*)'GWIDTHs: ',(gwidth(j),j=1,ngroups)
  !
  !  Extract Each Group's length (number of values in each group)
  !
  allocate(glen(ngroups),stat=is)
  !print *,'ALLOC glen: ',is
  !print *,'SAG3: ',nbitsglen,ngroups,iofst,idrstmpl(14),idrstmpl(13)
  if (nbitsglen.ne.0) then
     call g2_gbytesc(cpack,glen,iofst,nbitsglen,0,ngroups)
     itemp=nbitsglen*ngroups
     iofst=iofst+(itemp)
     if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8))
  else
     glen(1:ngroups)=0
  endif
  do j=1,ngroups
     glen(j)=(glen(j)*idrstmpl(14))+idrstmpl(13)
  enddo
  glen(ngroups)=idrstmpl(15)
  !write(78,*)'GLENs: ',(glen(j),j=1,ngroups)
  !print *,'GLENsum: ',sum(glen)
  !
  !  Test to see if the group widths and lengths are consistent with number of
  !  values, and length of section 7.
  !
  totbit = 0
  totlen = 0
  do j=1,ngroups
     totbit = totbit + (gwidth(j)*glen(j));
     totlen = totlen + glen(j);
  enddo
  if (totlen .NE. ndpts) then
     ier=1
     return
  endif
  if ( (totbit/8) .GT. lensec) then
     ier=1
     return
  endif
  !
  !  For each group, unpack data values
  !
  if ( idrstmpl(7).eq.0 ) then        ! no missing values
     n=1
     do j=1,ngroups
        !write(78,*)'NGP ',j,gwidth(j),glen(j),gref(j)
        if (gwidth(j).ne.0) then
           call g2_gbytesc(cpack,ifld(n),iofst,gwidth(j),0,glen(j))
           do k=1,glen(j)
              ifld(n)=ifld(n)+gref(j)
              n=n+1
           enddo
        else
           ifld(n:n+glen(j)-1)=gref(j)
           n=n+glen(j)
        endif
        iofst=iofst+(gwidth(j)*glen(j))
     enddo
  elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then
     ! missing values included
     allocate(ifldmiss(ndpts))
     !ifldmiss=0
     n=1
     non=1
     do j=1,ngroups
        !print *,'SAGNGP ',j,gwidth(j),glen(j),gref(j)
        if (gwidth(j).ne.0) then
           msng1=(2**gwidth(j))-1
           msng2=msng1-1
           call g2_gbytesc(cpack,ifld(n),iofst,gwidth(j),0,glen(j))
           iofst=iofst+(gwidth(j)*glen(j))
           do k=1,glen(j)
              if (ifld(n).eq.msng1) then
                 ifldmiss(n)=1
              elseif (idrstmpl(7).eq.2.AND.ifld(n).eq.msng2) then
                 ifldmiss(n)=2
              else
                 ifldmiss(n)=0
                 ifld(non)=ifld(n)+gref(j)
                 non=non+1
              endif
              n=n+1
           enddo
        else
           msng1=(2**nbitsgref)-1
           msng2=msng1-1
           if (gref(j).eq.msng1) then
              ifldmiss(n:n+glen(j)-1)=1
              !ifld(n:n+glen(j)-1)=0
           elseif (idrstmpl(7).eq.2.AND.gref(j).eq.msng2) then
              ifldmiss(n:n+glen(j)-1)=2
              !ifld(n:n+glen(j)-1)=0
           else
              ifldmiss(n:n+glen(j)-1)=0
              ifld(non:non+glen(j)-1)=gref(j)
              non=non+glen(j)
           endif
           n=n+glen(j)
        endif
     enddo
  endif
  !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts)
 
  if ( allocated(gref) ) deallocate(gref)
  if ( allocated(gwidth) ) deallocate(gwidth)
  if ( allocated(glen) ) deallocate(glen)
  !
  !  If using spatial differences, add overall min value, and
  !  sum up recursively
  !
  if (idrsnum.eq.3) then         ! spatial differencing
     if (idrstmpl(17).eq.1) then      ! first order
        ifld(1)=ival1
        if ( idrstmpl(7).eq.0 ) then        ! no missing values
           itemp=ndpts
        else
           itemp=non-1
        endif
        do n=2,itemp
           ifld(n)=ifld(n)+minsd
           ifld(n)=ifld(n)+ifld(n-1)
        enddo
     elseif (idrstmpl(17).eq.2) then    ! second order
        ifld(1)=ival1
        ifld(2)=ival2
        if ( idrstmpl(7).eq.0 ) then        ! no missing values
           itemp=ndpts
        else
           itemp=non-1
        endif
        do n=3,itemp
           ifld(n)=ifld(n)+minsd
           ifld(n)=ifld(n)+(2*ifld(n-1))-ifld(n-2)
        enddo
     endif
     !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts)
  endif
  !
  !  Scale data back to original form
  !
  !print *,'SAGT: ',ref,bscale,dscale
  if ( idrstmpl(7).eq.0 ) then        ! no missing values
     do n=1,ndpts
        fld(n)=((real(ifld(n))*bscale)+ref)*dscale
        !write(78,*)'SAG ',n,fld(n),ifld(n),bscale,ref,1./dscale
     enddo
  elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then
     ! missing values included
     non=1
     do n=1,ndpts
        if ( ifldmiss(n).eq.0 ) then
           fld(n)=((real(ifld(non))*bscale)+ref)*dscale
           !print *,'SAG ',n,fld(n),ifld(non),bscale,ref,dscale
           non=non+1
        elseif ( ifldmiss(n).eq.1 ) then
           fld(n)=rmiss1
        elseif ( ifldmiss(n).eq.2 ) then
           fld(n)=rmiss2
        endif
     enddo
     if ( allocated(ifldmiss) ) deallocate(ifldmiss)
  endif
 
  if ( allocated(ifld) ) deallocate(ifld)
 
  !open(10,form='unformatted',recl=24180,access='direct') 
  !read(10,rec=1) (fldo(k),k=1,6045)
  !do i =1,6045
  !  print *,i,fldo(i),fld(i),fldo(i)-fld(i)
  !enddo
 
  return
end subroutine comunpack
 
subroutine misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
 
  use intmath
  integer,intent(in) :: ndpts,idrsnum
  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,allocatable :: ifld(:),ifldmiss(:),jfld(:)
  integer,allocatable :: jmin(:),jmax(:),lbit(:)
  integer,parameter :: zero=0
  integer,allocatable :: gref(:),gwidth(:),glen(:)
  integer :: glength,grpwidth
  logical :: simple_alg
  logical :: have_rmin
 
  simple_alg = .false.
  have_rmin = .false.
  bscale=2.0**real(-idrstmpl(2))
  dscale=10.0**real(idrstmpl(3))
  missopt=idrstmpl(7)
  if ( missopt.ne.1 .AND. missopt.ne.2 ) then
     print *,'misspack: Unrecognized option.'
     lcpack=-1
     return
  else     !  Get missing values
     call rdieee(idrstmpl(8),rmissp,1)
     if (missopt.eq.2) call rdieee(idrstmpl(9),rmisss,1)
  endif
 
  !  Find min value of non-missing values in the data,
  !  AND set up missing value mapping of the field.
  allocate(ifldmiss(ndpts))
 
  if ( missopt .eq. 1 ) then        ! Primary missing value only
     do j=1,ndpts
        if (fld(j).eq.rmissp) then
           ifldmiss(j)=1
        else
           ifldmiss(j)=0
           if(have_rmin) then
              if (fld(j).lt.rmin) rmin=fld(j)
           else
              rmin=fld(j)
              have_rmin=.true.
           endif
        endif
     enddo
     if(.not.have_rmin) rmin=rmissp
  endif
  if ( missopt .eq. 2 ) then        ! Primary and secondary missing values
     do j=1,ndpts
        if (fld(j).eq.rmissp) then
           ifldmiss(j)=1
        elseif (fld(j).eq.rmisss) then
           ifldmiss(j)=2
        else
           ifldmiss(j)=0
           if(have_rmin) then
              if (fld(j).lt.rmin) rmin=fld(j)
           else
              rmin=fld(j)
              have_rmin=.true.
           endif
        endif
        if(.not.have_rmin) rmin=rmissp
     enddo
  endif
 
  !  Allocate work arrays:
  !  Note: -ifldmiss(j),j=1,ndpts is a map of original field indicating 
  !         which of the original data values
  !         are primary missing (1), sencondary missing (2) or non-missing (0).
  !        -jfld(j),j=1,nonmiss is a subarray of just the non-missing values from
  !         the original field.
  iofst=0
  allocate(ifld(ndpts))
  allocate(jfld(ndpts))
  allocate(gref(ndpts))
  allocate(gwidth(ndpts))
  allocate(glen(ndpts))
  !
  !  Scale original data
  !
  nonmiss=0
  if (idrstmpl(2).eq.0) then        !  No binary scaling
     imin=nint(rmin*dscale)
     !imax=nint(rmax*dscale)
     rmin=real(imin)
     do j=1,ndpts
        if (ifldmiss(j).eq.0) then
           nonmiss=nonmiss+1
           jfld(nonmiss)=max(0,nint(fld(j)*dscale)-imin)
        endif
     enddo
  else                              !  Use binary scaling factor
     rmin=rmin*dscale
     do j=1,ndpts
        if (ifldmiss(j).eq.0) then
           nonmiss=nonmiss+1
           jfld(nonmiss)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
        endif
     enddo
  endif
  !
  !  Calculate Spatial differences, if using DRS Template 5.3
  !
  if (idrsnum.eq.3) then        ! spatial differences
     if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2
     if (idrstmpl(17).eq.1) then      ! first order
        if(nonmiss<1) then
           ival1=1.0
        else
           ival1=jfld(1)
        endif
        do j=nonmiss,2,-1
           jfld(j)=jfld(j)-jfld(j-1)
        enddo
        if(nonmiss>0)             jfld(1)=0
     elseif (idrstmpl(17).eq.2) then      ! second order
        if(nonmiss==1) then
           ival1=jfld(1)
           ival2=jfld(1)
        elseif(nonmiss<1) then
           ival1=1.0
           ival2=1.0
        else
           ival1=jfld(1)
           ival2=jfld(2)
        endif
        do j=nonmiss,3,-1
           jfld(j)=jfld(j)-(2*jfld(j-1))+jfld(j-2)
        enddo
        if(nonmiss>=1) jfld(1)=0
        if(nonmiss>=2) jfld(2)=0
     endif
 
     !  subtract min value from spatial diff field
     isd=idrstmpl(17)+1
     minsd=minval(jfld(isd:nonmiss))
     do j=isd,nonmiss
        jfld(j)=jfld(j)-minsd
     enddo
 
     !   find num of bits need to store minsd and add 1 extra bit
     !   to indicate sign
     temp=i1log2(abs(minsd))
     nbitsd=ceiling(temp)+1
 
     !   find num of bits need to store ifld(1) ( and ifld(2)
     !   if using 2nd order differencing )
     maxorig=ival1
     if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2
     temp=i1log2(maxorig)
     nbitorig=ceiling(temp)+1
     if (nbitorig.gt.nbitsd) nbitsd=nbitorig
     !   increase number of bits to even multiple of 8 ( octet )
     if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8))
 
     !  Store extra spatial differencing info into the packed
     !  data section.
     if (nbitsd.ne.0) then
        !   pack first original value
        if (ival1.ge.0) then
           call g2_sbytec(cpack,ival1,iofst,nbitsd)
           iofst=iofst+nbitsd
        else
           call g2_sbytec(cpack,1,iofst,1)
           iofst=iofst+1
           call g2_sbytec(cpack,iabs(ival1),iofst,nbitsd-1)
           iofst=iofst+nbitsd-1
        endif
        if (idrstmpl(17).eq.2) then
           !  pack second original value
           if (ival2.ge.0) then
              call g2_sbytec(cpack,ival2,iofst,nbitsd)
              iofst=iofst+nbitsd
           else
              call g2_sbytec(cpack,1,iofst,1)
              iofst=iofst+1
              call g2_sbytec(cpack,iabs(ival2),iofst,nbitsd-1)
              iofst=iofst+nbitsd-1
           endif
        endif
        !  pack overall min of spatial differences
        if (minsd.ge.0) then
           call g2_sbytec(cpack,minsd,iofst,nbitsd)
           iofst=iofst+nbitsd
        else
           call g2_sbytec(cpack,1,iofst,1)
           iofst=iofst+1
           call g2_sbytec(cpack,iabs(minsd),iofst,nbitsd-1)
           iofst=iofst+nbitsd-1
        endif
     endif
     !print *,'SDp ',ival1,ival2,minsd,nbitsd
  endif     !  end of spatial diff section
 
  !  Expand non-missing data values to original grid.
  miss1=minval(jfld(1:nonmiss))-1
  miss2=miss1-1
  n=0
  do j=1,ndpts
     if ( ifldmiss(j).eq.0 ) then
        n=n+1
        ifld(j)=jfld(n)
     elseif ( ifldmiss(j).eq.1 ) then
        ifld(j)=miss1
     elseif ( ifldmiss(j).eq.2 ) then
        ifld(j)=miss2
     endif
  enddo
  if(ndpts<2) simple_alg=.true.
 
  !   Determine Groups to be used.
  if ( simple_alg ) then
     !  set group length to 10 :  calculate number of groups
     !  and length of last group
     ngroups=ndpts/10
     glen(1:ngroups)=10
     itemp=mod(ndpts,10)
     if (itemp.ne.0) then
        ngroups=ngroups+1
        glen(ngroups)=itemp
     endif
  else
     ! Use Dr. Glahn's algorithm for determining grouping.
     !
     minpk=10 
     inc=1
     maxgrps=(ndpts/minpk)+1
     allocate(jmin(maxgrps))
     allocate(jmax(maxgrps))
     allocate(lbit(maxgrps))
     call pack_gp(ifld,ndpts,missopt,minpk,inc,miss1,miss2, &
          jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit, &
          kbit,novref,lbitref,ier)
     !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref
     do ng=1,ngroups
        glen(ng)=glen(ng)+novref
     enddo
     deallocate(jmin)
     deallocate(jmax)
     deallocate(lbit)
  endif
 
  !  For each group, find the group's reference value (min)
  !  and the number of bits needed to hold the remaining values
  n=1
  do ng=1,ngroups
     !  how many of each type?
     num0=count(ifldmiss(n:n+glen(ng)-1) .EQ. 0)
     num1=count(ifldmiss(n:n+glen(ng)-1) .EQ. 1)
     num2=count(ifldmiss(n:n+glen(ng)-1) .EQ. 2)
     if ( num0.eq.0 ) then      ! all missing values
        if ( num1.eq.0 ) then       ! all secondary missing
           gref(ng)=-2
           gwidth(ng)=0
        elseif ( num2.eq.0 ) then       ! all primary missing
           gref(ng)=-1
           gwidth(ng)=0
        else                           ! both primary and secondary
           gref(ng)=0
           gwidth(ng)=1
        endif
     else                       ! contains some non-missing data
        !    find max and min values of group
        gref(ng)=huge(n)
        imax=-1*huge(n)
        j=n
        do lg=1,glen(ng)
           if ( ifldmiss(j).eq.0 ) then
              if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) 
              if (ifld(j).gt.imax) imax=ifld(j) 
           endif
           j=j+1
        enddo
        if (missopt.eq.1) imax=imax+1
        if (missopt.eq.2) imax=imax+2
        !   calc num of bits needed to hold data
        if ( gref(ng).ne.imax ) then
           temp=i1log2(imax-gref(ng))
           gwidth(ng)=ceiling(temp)
        else
           gwidth(ng)=0
        endif
     endif
     !   Subtract min from data
     j=n
     mtemp=2**gwidth(ng)
     do lg=1,glen(ng)
        if (ifldmiss(j).eq.0) then       ! non-missing
           ifld(j)=ifld(j)-gref(ng)
        elseif (ifldmiss(j).eq.1) then    ! primary missing
           ifld(j)=mtemp-1
        elseif (ifldmiss(j).eq.2) then    ! secondary missing
           ifld(j)=mtemp-2
        endif
        j=j+1
     enddo
     !   increment fld array counter
     n=n+glen(ng)
  enddo
 
  !  Find max of the group references and calc num of bits needed 
  !  to pack each groups reference value, then
  !  pack up group reference values
  !write(77,*)'GREFS: ',(gref(j),j=1,ngroups)
  igmax=maxval(gref(1:ngroups))
  if (missopt.eq.1) igmax=igmax+1
  if (missopt.eq.2) igmax=igmax+2
  if (igmax.ne.0) then
     temp=i1log2(igmax)
     nbitsgref=ceiling(temp)
     ! restet the ref values of any "missing only" groups.
     mtemp=2**nbitsgref
     do j=1,ngroups
        if (gref(j).eq.-1) gref(j)=mtemp-1
        if (gref(j).eq.-2) gref(j)=mtemp-2
     enddo
     call g2_sbytesc(cpack,gref,iofst,nbitsgref,0,ngroups)
     itemp=nbitsgref*ngroups
     iofst=iofst+itemp
     !         Pad last octet with Zeros, if necessary,
     if (mod(itemp,8).ne.0) then
        left=8-mod(itemp,8)
        call g2_sbytec(cpack,zero,iofst,left)
        iofst=iofst+left
     endif
  else
     nbitsgref=0
  endif
 
  !  Find max/min of the group widths and calc num of bits needed
  !  to pack each groups width value, then
  !  pack up group width values
  !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
  iwmax=maxval(gwidth(1:ngroups))
  ngwidthref=minval(gwidth(1:ngroups))
  if (iwmax.ne.ngwidthref) then
     temp=i1log2(iwmax-ngwidthref)
     nbitsgwidth=ceiling(temp)
     do i=1,ngroups
        gwidth(i)=gwidth(i)-ngwidthref
     enddo
     call g2_sbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
     itemp=nbitsgwidth*ngroups
     iofst=iofst+itemp
     !         Pad last octet with Zeros, if necessary,
     if (mod(itemp,8).ne.0) then
        left=8-mod(itemp,8)
        call g2_sbytec(cpack,zero,iofst,left)
        iofst=iofst+left
     endif
  else
     nbitsgwidth=0
     gwidth(1:ngroups)=0
  endif
 
  !  Find max/min of the group lengths and calc num of bits needed
  !  to pack each groups length value, then
  !  pack up group length values
  !write(77,*)'GLENS: ',(glen(j),j=1,ngroups)
  ilmax=maxval(glen(1:ngroups-1))
  nglenref=minval(glen(1:ngroups-1))
  if(ngroups>0) then
     nglenlast=glen(ngroups)
  else
     nglenlast=0
  endif
  if (ilmax.ne.nglenref) then
     temp=i1log2(ilmax-nglenref)
     nbitsglen=ceiling(temp)
     do i=1,ngroups-1
        glen(i)=glen(i)-nglenref
     enddo
     call g2_sbytesc(cpack,glen,iofst,nbitsglen,0,ngroups)
     itemp=nbitsglen*ngroups
     iofst=iofst+itemp
     !         Pad last octet with Zeros, if necessary,
     if (mod(itemp,8).ne.0) then
        left=8-mod(itemp,8)
        call g2_sbytec(cpack,zero,iofst,left)
        iofst=iofst+left
     endif
  else
     nbitsglen=0
     glen(1:ngroups)=0
  endif
 
  !  For each group, pack data values
  !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
  n=1
  ij=0
  do ng=1,ngroups
     glength=glen(ng)+nglenref
     if (ng.eq.ngroups ) glength=nglenlast
     grpwidth=gwidth(ng)+ngwidthref
     !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
     if ( grpwidth.ne.0 ) then
        call g2_sbytesc(cpack,ifld(n),iofst,grpwidth,0,glength)
        iofst=iofst+(grpwidth*glength)
     endif
     do kk=1,glength
        ij=ij+1
        !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
     enddo
     n=n+glength
  enddo
  !         Pad last octet with Zeros, if necessary,
  if (mod(iofst,8).ne.0) then
     left=8-mod(iofst,8)
     call g2_sbytec(cpack,zero,iofst,left)
     iofst=iofst+left
  endif
  lcpack=iofst/8
  !
  if ( allocated(ifld) ) deallocate(ifld)
  if ( allocated(jfld) ) deallocate(jfld)
  if ( allocated(ifldmiss) ) deallocate(ifldmiss)
  if ( allocated(gref) ) deallocate(gref)
  if ( allocated(gwidth) ) deallocate(gwidth)
  if ( allocated(glen) ) deallocate(glen)
 
  !     Fill in ref value and number of bits in Template 5.2
  rmin4 = real(rmin, 4)
  call mkieee(rmin4,ref,1)   ! ensure reference value is IEEE format
  iref=transfer(ref,iref)
  idrstmpl(1)=iref
  idrstmpl(4)=nbitsgref
  idrstmpl(5)=0         ! original data were reals
  idrstmpl(6)=1         ! general group splitting
  idrstmpl(10)=ngroups          ! Number of groups
  idrstmpl(11)=ngwidthref       ! reference for group widths
  idrstmpl(12)=nbitsgwidth      ! num bits used for group widths
  idrstmpl(13)=nglenref         ! Reference for group lengths
  idrstmpl(14)=1                ! length increment for group lengths
  idrstmpl(15)=nglenlast        ! True length of last group
  idrstmpl(16)=nbitsglen        ! num bits used for group lengths
  if (idrsnum.eq.3) then
     idrstmpl(18)=nbitsd/8      ! num bits used for extra spatial
     ! differencing values
  endif
 
end subroutine misspack