c***********************************************************************
      SUBROUTINE PREPNAPPE(u,v,w,c,p,ls,nx,ny,nz)
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k,n0,n1s2,n2s2
      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,x0,dl,x,y,z,theta,phi,r,r1,r2
      DOUBLE PRECISION GP_DFETCH
      DOUBLE PRECISION pi
      DOUBLE PRECISION sign, lambdaz

      WRITE(*,*) "PREPARING NAPPE"
      pi = 3.1416
      r0 = GP_DFETCH('radius',1.d0,'SEVERE')
      dl = GP_DFETCH('liquid_layer',0.d0,'WARNING')
      x0 = GP_DFETCH('cyl_center',0.d0,'WARNING')

      r1 = GP_DFETCH('amplitude_of_disturbance',0.d0,'DEFAULT')
      n1s2 = GP_IFETCH('number_of_digits',0,'DEFAULT')
c***
     
      write(*,*) 'r0=',r0
      write(*,*) 'x0=',x0
      write(*,*) 'dl=',dl
      h = 1.d0/(nx-2)
      lambdaz = 2*nz*h     ! S.Z. definition
      lambdaz = n1s2*h     ! A.B. definition
      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 * cos(2 * pi * z / lambdaz + pi)
               if (x .lt. r) then
                  sign = dl - y 
               else
                  sign = dl**2 - (y**2 + (x-r)**2) 
               end if 

               if (sign .ge. 0.d0) then
                  c(i,j,k)=1.d0  
                  u(i,j,k)=0.d0
                  v(i,j,k)=0.d0
                  w(i,j,k)=0.d0
               else 
                  c(i,j,k)=0.d0
                  u(i,j,k)=0.d0
                  v(i,j,k)=0.d0
                  w(i,j,k)=0.d0
               end if

               p(i,j,k)=0.d0
            end do
         end do
      end do
c      call levelset2vof(ls,c)
c      call writeminmax('cc at end of prepjet',c,nx,ny,nz)
      return
      end


