c***********************************************************************
      SUBROUTINE PREPRAYJET(u,v,w,c,p,ls,nx,ny,nz)
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k,n1
      INTEGER GP_IFETCH
      DOUBLE PRECISION c(nx,ny,nz), u(nx,ny,nz), v(nx,ny,nz)
      DOUBLE PRECISION w(nx,ny,nz),p(nx,ny,nz),h,ls(nx,ny,nz)
      DOUBLE PRECISION r0,x,y,z,r,r1,kz
      DOUBLE PRECISION GP_DFETCH
      DOUBLE PRECISION pi


      WRITE(*,*) "PREPARING RAYLEIGH JET"
      pi=4.d0*datan(1.d0)

      r0 = GP_DFETCH('radius',1.d0,'SEVERE')
      r1 = GP_DFETCH('amplitude_of_disturbance',0.d0,'DEFAULT')
      n1 = GP_IFETCH('number_of_digits',1,'DEFAULT')
c***

      write(*,*) 'r0=',r0
      h = 1.d0/(nx-2)

c  Compute z wavenumber. for 1 wave it is 2 pi/ Lz 
c  using symmetry it is  pi / Lz and Lz = nz * h 

      kz = pi * n1 / (h * nz)
      do j=1,ny
         y=(j-1.5)*h
         do i=1,nx
            x=(i-1.5)*h+PPRANK
            do k=1,nz
               z=(k-1.5)*h
               r = r0 - r1 * (sin ( 2 * kz * z ) + 0.1 * cos ( kz * z))
               ls(i,j,k) = r**2 - (y**2 + x**2) 
               u(i,j,k)=0.d0
               v(i,j,k)=0.d0
               w(i,j,k)=0.d0
               p(i,j,k)=0.d0
            end do
         end do
      end do
      call levelset2vof(ls,c)
c      call writeminmax('cc at end of prepjet',c,nx,ny,nz)
      return
      end


