specunpack.f

Go to the documentation of this file.

 
 
 
      subroutine specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld)
 
      character(len=1),intent(in) :: cpack(len)
      integer,intent(in) :: ndpts,len,JJ,KK,MM
      integer,intent(in) :: idrstmpl(*)
      real,intent(out) :: fld(ndpts)
 
      integer :: ifld(ndpts),Ts
      integer(4) :: ieee
      real :: ref,bscale,dscale,unpk(ndpts)
      real,allocatable :: pscale(:)
 
      ieee = idrstmpl(1)
      call rdieee(ieee,ref,1)
      bscale = 2.0**real(idrstmpl(2))
      dscale = 10.0**real(-idrstmpl(3))
      nbits = idrstmpl(4)
      js=idrstmpl(6)
      ks=idrstmpl(7)
      ms=idrstmpl(8)
      ts=idrstmpl(9)
 
      if (idrstmpl(10).eq.1) then           ! unpacked floats are 32-bit IEEE
         !call g2_gbytesc(cpack,ifld,0,32,0,Ts)
         call rdieee(cpack,unpk,ts)          ! read IEEE unpacked floats
         iofst=32*ts
         call g2_gbytesc(cpack,ifld,iofst,nbits,0,ndpts-ts)  ! unpack scaled data
!
!   Calculate Laplacian scaling factors for each possible wave number.
!
         allocate(pscale(jj+mm))
         tscale=real(idrstmpl(5))*1e-6
         do n=js,jj+mm
            pscale(n)=real(n*(n+1))**(-tscale)
         enddo
!
!   Assemble spectral coeffs back to original order.
!
         inc=1
         incu=1
         incp=1
         do m=0,mm
            nm=jj      ! triangular or trapezoidal
            if ( kk .eq. jj+mm ) nm=jj+m          ! rhombodial
            ns=js      ! triangular or trapezoidal
            if ( ks .eq. js+ms ) ns=js+m          ! rhombodial
            do n=m,nm
               if (n.le.ns .AND. m.le.ms) then    ! grab unpacked value
                  fld(inc)=unpk(incu)         ! real part
                  fld(inc+1)=unpk(incu+1)     ! imaginary part
                  inc=inc+2
                  incu=incu+2
               else                         ! Calc coeff from packed value
                  fld(inc)=((real(ifld(incp))*bscale)+ref)*
     &                      dscale*pscale(n)           ! real part
                  fld(inc+1)=((real(ifld(incp+1))*bscale)+ref)*
     &                      dscale*pscale(n)           ! imaginary part
                  inc=inc+2
                  incp=incp+2
               endif
            enddo
         enddo
 
         deallocate(pscale)
 
      else
         print *,'specunpack: Cannot handle 64 or 128-bit floats.'
         fld=0.0
         return
      endif
 
      return
      end