specunpack.F90

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 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
  endif
end subroutine specunpack