specpack.f

Go to the documentation of this file.

 
 
 
      subroutine specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
 
      real,intent(in) :: fld(ndpts)
      integer,intent(in) :: ndpts,JJ,KK,MM
      integer,intent(inout) :: idrstmpl(*)
      character(len=1),intent(out) :: cpack(*)
      integer,intent(out) :: lcpack
 
      integer :: ifld(ndpts),Ts,tmplsim(5)
      real :: bscale,dscale,unpk(ndpts),tfld(ndpts)
      real,allocatable :: pscale(:)
 
      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)
 
!
!   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
!
!   Separate spectral coeffs into two lists; one to contain unpacked
!   values within the sub-spectrum Js, Ks, Ms, and the other with values 
!   outside of the sub-spectrum to be packed.
!
      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    ! save unpacked value
               unpk(incu)=fld(inc)         ! real part
               unpk(incu+1)=fld(inc+1)     ! imaginary part
               inc=inc+2
               incu=incu+2
            else                         ! Save value to be packed and scale
                                         ! Laplacian scale factor
               tfld(incp)=fld(inc)*pscale(n)         ! real part
               tfld(incp+1)=fld(inc+1)*pscale(n)     ! imaginary part
               inc=inc+2
               incp=incp+2
            endif
         enddo
      enddo
 
      deallocate(pscale)
 
      incu=incu-1
      if (incu .ne. ts) then
         print *,'specpack: Incorrect number of unpacked values ',
     &           'given:',ts     
         print *,'specpack: Resetting idrstmpl(9) to ',incu
         ts=incu
      endif
!
!  Add unpacked values to the packed data array in 32-bit IEEE format
!
      call mkieee(unpk,cpack,ts)
      ipos=4*ts
!
!  Scale and pack the rest of the coefficients
! 
      tmplsim(2)=idrstmpl(2)
      tmplsim(3)=idrstmpl(3)
      tmplsim(4)=idrstmpl(4)
      call simpack(tfld,ndpts-ts,tmplsim,cpack(ipos+1),lcpack)
      lcpack=lcpack+ipos
!
!  Fill in Template 5.51
!
      idrstmpl(1)=tmplsim(1)
      idrstmpl(2)=tmplsim(2)
      idrstmpl(3)=tmplsim(3)
      idrstmpl(4)=tmplsim(4)
      idrstmpl(9)=ts
      idrstmpl(10)=1         ! Unpacked spectral data is 32-bit IEEE
 
      return
      end