--- From the main body of triple_zyx.f -------------------------------
***
      call xopause(tbar,beta,zip,work,work(nz+2),nz,ibz,rx,ry,rz,ipro)
***
          call wavepacket(tt,ww,v3,tbar,beta,work,work(nxp+1),
     &                    work(nxp+nz+2),nx,ny,nz,n3,ibx,iby,ibz,
     &                 tgx,tgy,tgz,wp,fn,mp,ipro,mype,ncpu,ncput)
***
 --- And the subroutine -----------------------------------------------
!***********************************************************************
      subroutine wavepacket(T,W,V3,tbar,beta,x,z,emask,nx,ny,nz,n3,
     &                            ibx,iby,ibz,tgx,tgy,tgz,wp,fn,mp,
     &                                       iprob,mype,ncpu,ncput)
!-----------------------------------------------------------------------
!  Initialize a wave packet for iprob=12 and iprob=13.  The wave packet
!  will propagate into the x-o-pause created in subroutine xopause.
!
!  VARIABLES:
!    freq = omega/N (where omega and N here have dimensions of sec^-1)
!    wide = standard deviation for wave packet in z
!
!  iprob=12: 
!    wave ~ exp{-(z-z0)/2.wide^2} exp{i.(k.x+omega.t)}   -or-
!    wave ~ exp{-(z-z0)/2.wide^2} exp{i.2pi(x+freq.t)}
!  
!    since we nondimensionalize by lambda,T=2pi/N,lambda.beta1
!
!  iprob=13:
!    wave ~ slanted and limited by wide2 in wave phase plane
!
!    .              .              .
!     `^.            `^.            `^.           freq = kx/k
!        `^.            `^.            `^.
!           `^.            `^.           .|^.
!              `^.            `^.     k . |  `^.
!    .            `^.            `^.   .  |kz   `^.
!     `^.            `^.            `^.kx_|        `^.
!        `^.            `^.            `^.            `^.
!           `^.            `^.            `^.            `^.
!              `^.            `^.            `^.            `^.
!                 `^.            `^.            `^.            `^.
!                    `^.            `^.            `^.            `^.
!
!  the phase velocity is (omega/k^2)(kx,ky,kz).  so the vertical phase
!  velocity is omega.kz/k^2
!-----------------------------------------------------------------------
      ! include some parameter files
      include 'include/garten.fi'
      include 'include/pars.fi'

      ! declare arrays
      dimension T(n3),W(n3),V3(n3)
      dimension tbar(nz+ibz,2),x(nx/ncpu),z(nz+1),emask(nz+1)
      dimension tgx(6*(nx+ibx)+64),tgy(6*(ny+iby)+64),tgz(6*(nz+ibz)+64)
      dimension wp(6),mp(ncpu+1,4)

      ! define useful constants
      nxp    = nx / ncpu
      nzp    = nz / ncpu
      nxt    = nx + ibx
      nyt    = ny + iby
      nzt    = nz + ibz
      zmx1   = zmx - zmx2
      pi     = 4.0 * atan( 1.0 )
      x0     = xmx * 0.5
      z0     = zmx1-3.*wide
      ww2i   = 1./(2.*wide *wide )
      ww22i  = 1./(2.*wide2*wide2)
      fs     = freq
      fc     = sqrt(1-freq*freq)

      ! set the vertical-velocity amplitude.  don't use a minus sign
      ! because this code is goofy.
      wamp   = 2.*pi*freq

      tzr1i  = 1./(tzr-1.)
      z1     = zmx/float(nz)
      x1     = xmx/float(nx)
      if (mype.eq.ncpu) then
        nxpp = abs(ibx)
        nzpp = abs(ibz)
      else
        nxpp = nxp
        nzpp = nzp
      endif
      ioff   = mype * nxp

      ! define x
      do i=1,nxpp
        x(i) = float( ioff + i-1 )*x1
      enddo
      
      ! define z
      do k=1,nz+1
        z(k) = float(        k-1 )*z1
      enddo

      ! define mask to supress initial waves in upper region
      if (tzr.eq.1.) then
        do k=1,nzt
          emask(k) = 1.0
        enddo
      else
        do k=1,nzt
          emask(k) = (tzr-tbar(k,2))*tzr1i
        enddo
      endif

      ! initialize wave packet with envelope in z only
      wkx=2.*pi*freq
      wkz=2.*pi*sqrt(1-freq*freq)
      if (iprob.eq.12) then

        cc =-exp(-z0*z0*ww2i)
        dd = exp(-(zmx-z0)*(zmx-z0)*ww2i)
        bb = dd*((z0-zmx)*ww2i-1./zmx) + cc*(z0*ww2i-1./zmx)
        aa = dd*ww2i*(1.-z0/zmx) - cc*z0*ww2i/zmx
        
        ! specify the temperature wave field
        do i=1,nxpp
          li = nzt*nyt*(i-1)
          wkxxi = wkx*x(i)
          do j=1,nyt
            lj = nzt*(j-1) + li
            do k=1,nzt
              l = k + lj
              zdiff = (z(k)-z0)*(z(k)-z0)
              T(l)  = amp*cos(wkxxi+wkz*z(k))*emask(k)*
     &                (exp(-zdiff*ww2i)+aa*z(k)*z(k)+bb*z(k)+cc)
!
!  don't add in the mean part for T yet.
!    &                  - tbar(k,1)
!                       ! minus sign because this code is goofy!
            enddo
          enddo
        enddo
        do i=nxpp+1,nxp
          li = nzt*nyt*(i-1)
          do j=1,nyt
            lj = nzt*(j-1) + li
            do k=1,nzt
              l = k + lj
              T(l) = 0.0
            enddo
          enddo
        enddo

      ! initialize wave packet with envelope in z and x
      elseif (iprob.eq.13) then
        do i=1,nxpp
          li = nzt*nyt*(i-1)
          xdiff = x(i)-x0
          wkxxi = wkx*x(i)
          do j=1,nyt
            lj = nzt*(j-1) + li
            do k=1,nzt
              zdiff = z(k)-z0
              l = k + lj
              T(l) = amp
     &                   * exp( -(xdiff*fc-zdiff*fs)**2 * ww22i )
     &                   * exp( -(xdiff*fs+zdiff*fc)**2 * ww2i  )
     &                   * cos( wkxxi + wkz*z(k) )
     &                   * emask(k)
! 
!    &                 - tbar(k,1)
!                      ! minus sign because this code is goofy!
            enddo
          enddo
        enddo
        do i=nxpp+1,nxp
          li = nzt*nyt*(i-1)
          do j=1,nyt
            lj = nzt*(j-1) + li
            do k=1,nzt
              l = k + lj
              T(l) = 0.0
            enddo
          enddo
        enddo

!       do i=1,nxpp
!         li = nzt*nyt*(i-1)
!         xdiff = x(i)-x0
!         wkxxi = wkx*x(i)
!         do j=1,nyt
!           lj = nzt*(j-1) + li
!           do k=1,nzt
!             zdiff = z(k)-z0
!             l = k + lj
!             W(l) =wamp 
!    &                   * exp( -(xdiff*fc-zdiff*fs)**2 * ww22i )
!    &                   * exp( -(xdiff*fs+zdiff*fc)**2 * ww2i  )
!    &                   * sin( wkxxi + wkz*z(k) )
!    &                   * emask(k)
!           enddo
!         enddo
!       enddo
!       do i=nxpp+1,nxp
!         li = nzt*nyt*(i-1)
!         do j=1,nyt
!           lj = nzt*(j-1) + li
!           do k=1,nzt
!             l = k + lj
!             W(l) = 0.0
!           enddo
!         enddo
!       enddo

      ! unexpected value for iprob
      else
        write(6,*)"  * * * * * * * * * * * * * * * * * * * * *"
        write(6,*)"  *                                       *"
        write(6,*)"  *    subroutine wavepacket received     *"
        write(6,*)"  *    (iprob.ne.12).and.(iprob.ne.13)    *"
        write(6,*)"  *                                       *"
        write(6,*)"  *             BAILING OUT!              *"
        write(6,*)"  *                                       *"
        write(6,*)"  * * * * * * * * * * * * * * * * * * * * *"
        stop
      endif

      ! initialize Vort_3, W to zero
      do l=1,n3
        V3(l) = 0.0
      enddo
      do l=1,n3
        W(l) = 0.
      enddo

      ! Compute Fourier transform of T
      call transf(T,V3,tgx,tgy,tgz,wp,fn,nx,ny,nz, ibx, iby,-ibz,
     &                                   1,mype,ncpu,ncput,n3,mp)

      ! make sure ny=1
      if (ny.ne.1) then
        if (mype.eq.0) then
          write(6,*)"* * * * * * * * * * * * * * * * * * * *"
          write(6,*)"*                                     *"
          write(6,*)"*   wavepacket() called with ny.ne.1  *"
          write(6,*)"*                                     *"
          write(6,*)"*            bailing out!             *"
          write(6,*)"*                                     *"
          write(6,*)"* * * * * * * * * * * * * * * * * * * *"
        endif
        stop
      endif

      !
      ! Correct T spectrum so that only waves moving in (-x,-z)
      ! direction are permitted.  This requires that the following
      ! functional form be satisfied:
      !
      !      A[cos(ix)cos(kz)-sin(ix)sin(kz)] 
      !    + B[cos(ix)sin(kz)+sin(ix)cos(kz)]
      ! 
      ! hence, the coefficient for sin(ix)sin(kz) should be the
      ! negative of the coefficient for cos(ix)cos(kz), and the 
      ! coefficients for cos(ix)sin(kz) and sin(ix)cos(kz) should
      ! be the same.
      !
      ! for this you are going to have to understand the oddities
      ! of fftpack-data storage.  
      !
      ! for kx (first index), storage is
      !(cos0,cos1,-sin1,cos2,-sin2,cos3,-sin3,...,cosN-1,-sinN-1,cosN)
      !
      ! for kz (last index), storage is
      !(cos0,cosN,cos1,-sin1,cos2,-sin2,cos3,-sin3,...,cosN-1,-sinN-1)
      !
      ! Also initialize W according to T_t = - T_bar_z * W
      !
      ! If T=cos(ix+kz+ft), T_t=-f.sin(ix+kz+ft)=-T_bar_z*W
      !    T =      cos(ix)cos(kz)-sin(ix)sin(kz)
      !    W =  B { sin(ix)cos(kz)+cos(ix)sin(kz) } ; B=f/T_bar_z
      !
      ! If T=sin(ix+kz+ft), T_t= f.cos(ix+kz+ft)=-T_bar_z*W
      !    T =      sin(ix)cos(kz)+cos(ix)sin(kz)
      !    W = -B { cos(ix)cos(kz)-sin(ix)sin(kz) }
      !

      ! if mype=0, skip cos0 and cosN modes
      if (mype.eq.0) then
        kmin = 3
      else
        kmin = 1
      endif

      ! correct T and initialize W
      do k=kmin,nzpp,2
        j=1
          do i=2,nx,2

            lcc = i   + nx*(j-1) + nx*ny*(k-1)
            lcs = i   + nx*(j-1) + nx*ny*( k )
            lsc = i+1 + nx*(j-1) + nx*ny*(k-1)
            lss = i+1 + nx*(j-1) + nx*ny*( k )

            ! (-sin(ix))(-sin(kz)) = sin(ix)sin(kz)
            ! sin(ix)sin(kz) = -cos(ix)cos(kz)
            T(lss) = -T(lcc)
            W(lsc) =  T(lcc) * wamp
            W(lcs) =  T(lcc) * wamp

            ! each term has a -sin
            ! cos(ix)sin(kz) =  sin(ix)cos(kz)
            T(lcs) =  T(lsc)
            W(lcc) = -T(lsc) * wamp
            W(lss) =  T(lsc) * wamp

          enddo
      enddo

      ! Transform back to physical space
      call transf(T,V3,tgx,tgy,tgz,wp,fn,nx,ny,nz, ibx, iby,-ibz,
     &                                  -1,mype,ncpu,ncput,n3,mp)
      call transf(W,V3,tgx,tgy,tgz,wp,fn,nx,ny,nz, ibx, iby,-ibz,
     &                                  -1,mype,ncpu,ncput,n3,mp)

!     do i=1,nxpp
!       do j=1,nyt
!         do k=1,nzt
!           zdiff = (z(k)-z0)*(z(k)-z0)
!           W(k,j,i) =-wamp*exp(-zdiff*ww2i)*sin(wkx*x(i)+wkz*z(k))
!    &                                                   * emask(k)
!         enddo
!       enddo
!     enddo
!     do i=nxpp+1,nxp
!       do j=1,nyt
!         do k=1,nzt
!           W(k,j,i) = 0.0
!         enddo
!       enddo
!     enddo
      do i=1,nxp
        li = nzt*nyt*(i-1)
        do j=1,nyt
          lj = nzt*(j-1) + li
          do k=1,nzt
            l = k + lj
            V3(l) = 0.0
          enddo
        enddo
      enddo

      ! add T_bar to the temperature field
      do i=1,nxpp
        li = nzt*nyt*(i-1)
        do j=1,nyt
          lj = nzt*(j-1) + li
          do k=1,nzt
            l = k + lj
            T(l) =  T(l) - tbar(k,1)
                         ! minus sign because this code is goofy!
          enddo
        enddo
      enddo
!-----------------------------------------------------------------------
      return
      end