c***********************************************************************
      SUBROUTINE prepjet(u,v,w,c,dist,interf,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k,n0
      DOUBLE PRECISION c(nx,ny,nz), u(nx,ny,nz), v(nx,ny,nz)
      DOUBLE PRECISION w(nx,ny,nz), dist(nx,ny,nz),interf(nx,ny,nz)
      DOUBLE PRECISION u0,h,r0,x0,y0,z0,zz,yy,xx,ecl,outlet
      INTEGER GP_IFETCH
      DOUBLE PRECISION GP_DFETCH
      INTRINSIC DSQRT,DFLOAT,INT,DABS,DEXP
      EXTERNAL GP_DFETCH,GP_IFETCH

      write(*,*) '--- prepjet ---'
c***
c     center of the 1/8 of the sphere
c***
      r0 = GP_DFETCH('radius',1.d0,'SEVERE')
      u0 = GP_DFETCH('liquid_velocity',1.d0,'SEVERE')
      ecl = GP_DFETCH('boundary_layer_width',1.d0,'SEVERE')
      write(*,*) 'nx=',nx
      write(*,*) 'ny=',ny
      write(*,*) 'nz=',nz
      h = 1.d0/(nx-2)
      n0 = INT(r0/h) + 1
      x0 = (2+n0)*h  
      y0 = (ny-2)*h/2+h
      z0 = (nz-2)*h/2+h
c***
c     get distances of the vertices from the center of the drop,
c     set everything to zero
c***
      outlet=0.
      write(*,*) "b1"
      do k=1,nz
         zz = k*h
         do j=1,ny
            yy = j*h
            do i=1,nx
               xx = i*h
               dist(i,j,k) = DSQRT((yy-y0)*(yy-y0) + (zz-z0)*(zz-z0))
               u(i,j,k) = 0.d0
               c(i,j,k) = 0.d0
               if (dist(i,j,k).le.((4./5.)*r0).and.i.eq.1) then
                 u(1,j,k) = u0
                 c(1,j,k) = 1.d0
                 outlet = outlet + u0
               elseif (dist(i,j,k).le.r0.and.i.eq.1) then
                 u(1,j,k) = 5.*u0*(1. - dist(i,j,k)/r0)
                 c(1,j,k) = 1.d0
                 outlet = outlet + 5.*u0*(1. - dist(i,j,k)/r0)
               elseif(dist(i,j,k).gt.r0.and.i.eq.1) then
                 u(1,j,k) = 0.d0
                 c(1,j,k) = 0.d0
               endif
               if (dist(i,j,k).le.((4./5.)*r0).and.i.eq.2) then
                 u(2,j,k) = u0
                 c(2,j,k) = 1.d0
               elseif (dist(i,j,k).le.r0.and.i.eq.2) then
                 u(2,j,k) = 5.*u0*(1. - dist(i,j,k)/r0)
                 c(2,j,k) = 1.d0
               elseif(dist(i,j,k).gt.r0.and.i.eq.2) then
                 u(2,j,k) = 0.d0
                 c(2,j,k) = 0.d0
               endif
               v(i,j,k) = 0.d0
               w(i,j,k) = 0.d0
            enddo 
         enddo 
      enddo 
      do k=1,nz
        do j=1,ny
            u(nx,j,k) = outlet / (ny)**2
        enddo
      enddo      


c***
c     get the c fraction in each cell, in a very crude way,
c     just not to have a brutal change from 0 to 1
c***
 10   format (3(1x,i3),2(1x,1pe10.3))
c***
c     put pressure field in dist
c***
      do k=1,nz
         do j=1,ny
            do i=1,nx
               dist(i,j,k) = 0.d0
            enddo 
         enddo 
      enddo 
c***
      return
      end

c***********************************************************************
      SUBROUTINE PREPDISTDROP(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 u0,r0,x0,dl,x,y,z,theta,phi,r,r1,r2
      DOUBLE PRECISION GP_DFETCH
      DOUBLE PRECISION FINDN2S
      DOUBLE PRECISION y0

      r0 = GP_DFETCH('radius',1.d0,'SEVERE')
      u0 = GP_DFETCH('liquid_velocity',1.d0,'SEVERE')
      dl = GP_DFETCH('liquid_layer',0.d0,'WARNING')
      x0 = GP_DFETCH('drop_center',0.d0,'WARNING')
      y0 = GP_DFETCH('drop_center_y',0.d0,'WARNING')

      r1 = GP_DFETCH('amplitude_of_disturbance',0.d0,'DEFAULT')
      n1s2 = GP_IFETCH('number_of_digits',0,'DEFAULT')
      r2 = GP_DFETCH('amplitude_of_disturbance_2',0.d0,'DEFAULT')
      n2s2 = GP_IFETCH('number_of_digits_2',0,'DEFAULT')
c***
     
      write(*,*) 'r0=',r0
      write(*,*) 'u0=',u0
      write(*,*) 'x0=',x0
      write(*,*) 'y0=',y0
      write(*,*) 'dl=',dl

      h = 1.d0/(nx-2)
      do j=1,ny
         y=(j-1.5)*h - y0
         do k=1,nz
            z=(k-1.5)*h
            phi = atan(y/z)
            do i=1,nx
               x=(i-1.5)*h+PPRANK
               theta = atan((x-x0)/(sqrt(y**2+z**2)))
               r = r0 + r1 * cos(2 * n1s2 * phi) * cos(theta)
     &                + r2 * cos(2 * n2s2 * phi) * cos(theta)

               ls(i,j,k) = r**2 - ((x-x0)**2 + y**2 + z**2)
               if (ls(i,j,k).gt.0.d0) then
                  u(i,j,k)= -u0
                  v(i,j,k)=0.d0
                  w(i,j,k)=0.d0
                  p(i,j,k)=0.d0
               else if (x .lt. dl) then
                  u(i,j,k)=0.d0
                  v(i,j,k)=0.d0
                  w(i,j,k)=0.d0
                  p(i,j,k)=0.d0
               else
                  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 if
            end do
         end do
      end do
      i = INT(x0/h +1.5)
      write(*,*) 'at center before ls = '
     $        ,ls(i,2,2), " i = ",i
      call levelset2vof(ls,c)
      write (*,*) "Norm2 p initial = ",h**2*FINDN2S(p,nx,ny,nz)
      write(*,*) 'at center cc =',c(i,2,2),' ls = '
     $        ,ls(i,2,2)
      do j=1,ny
         do k=1,nz
            do i=1,nx
               x=(i-1.5)*h
               if(x .lt. dl) then
                  c(i,j,k)=1.d0
               end if
            enddo
         enddo
      enddo
      call writeminmax('cc at end of prepjet',c,nx,ny,nz)
      return
      end