c***********************************************************************
      SUBROUTINE prepjet(u,v,w,c,dist,nx,ny,nz,PARINT,PARDPR,COMM2D)
C
C     Name: prepjet
C     Author: Anthony Leboissetier (10/01)
C     Objective: initialization of the u,v,w & c if restart = 0
C     Called by: initialize
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INCLUDE 'mpif.h'
      INTEGER nx,ny,nz,i,j,k,n0,PARINT(100),PPRANK,NPROCS,IERR00,COMM2D
      DOUBLE PRECISION c(nx,ny,nz), u(nx,ny,nz), v(nx,ny,nz),PARDPR(100)
      DOUBLE PRECISION w(nx,ny,nz), dist(nx,ny,nz)
      DOUBLE PRECISION u0,h,r0,x0,y0,z0,zz,yy,xx,outlet
      INTEGER GP_IFETCH,l,STAT00(MPI_STATUS_SIZE)
      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***
      PPRANK = PARINT(55)
      NPROCS = PARINT(56)
      r0 = PARDPR(10)
      u0 = PARDPR(11)
      write(*,*) 'r0=',r0
      write(*,*) 'u0=',u0
      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.
      IF(PPRANK.EQ.0) THEN
      write (*,*) n0,h,y0,z0
      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.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).gt.r0.and.i.eq.1) then
                 u(1,j,k) = 0.d0
                 c(1,j,k) = 0.d0
               endif
               if (dist(2,j,k).le.r0.and.i.eq.2) then
                 u(2,j,k) = u0
                 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
      CALL MPI_SEND(outlet,1,MPI_DOUBLE_PRECISION,NPROCS-1,100,COMM2D,
     &              IERR00)
c      WRITE(*,*)'0 envoie',outlet,' à',NPROCS-1
      ENDIF
      IF(PPRANK.EQ.(NPROCS-1)) THEN
      CALL MPI_RECV(outlet,1,MPI_DOUBLE_PRECISION,0,100,COMM2D,STAT00,
     &              IERR00)
c      WRITE(*,*)NPROCS-1,' reçoit',outlet,' de 0'
        do k=1,nz
          do j=1,ny
              u(nx,j,k) = outlet / (ny)**2
          enddo
        enddo
      ENDIF
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