c
c   S. Zaleski - 25/02/2006
c   converts a level set field into a VOF field
c
c   On entry:
c
c
c    cc: the volume fraction
c    ls: the level set function
c
c   On output: 
c
c    the same
c
c
      subroutine levelset2vof(ls,cc)
      include "undefined.h"
      include "dimension.h"

      double precision cc(nx,ny,nz),ls(nx,ny,nz)
      double precision zero, one, normL1
      double precision mm1,mm2,mm3,mx,my,mz,alpha
      double precision fl3d
      integer i,j,k

      zero=0.d0
      one=1.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
c***
c     (1) normal vector: mx,my,mz; (2) mx,my,mz>0. and mx+my+mz = 1.;
c     (3) shift alpha to origin    (4) get volume from alpha;  
c***

               mm1 = ls(i-1,j-1,k-1)+ls(i-1,j-1,k+1)+ls(i-1,j+1,k-1)
     %              +ls(i-1,j+1,k+1)+2.0d0*(ls(i-1,j-1,k)+ls(i-1,j+1,k)
     %              +ls(i-1,j,k-1)+ls(i-1,j,k+1))+4.0d0*ls(i-1,j,k)
               mm2 = ls(i+1,j-1,k-1)+ls(i+1,j-1,k+1)+ls(i+1,j+1,k-1)
     %              +ls(i+1,j+1,k+1)+2.0d0*(ls(i+1,j-1,k)+ls(i+1,j+1,k)
     %              +ls(i+1,j,k-1)+ls(i+1,j,k+1))+4.0d0*ls(i+1,j,k)
               mx = (mm1 - mm2)/32.d0
               
               mm1 = ls(i-1,j-1,k-1)+ls(i-1,j-1,k+1)+ls(i+1,j-1,k-1)
     %              +ls(i+1,j-1,k+1)+2.0d0*(ls(i-1,j-1,k)+ls(i+1,j-1,k)
     %              +ls(i,j-1,k-1)+ls(i,j-1,k+1))+4.0d0*ls(i,j-1,k)
               mm2 = ls(i-1,j+1,k-1)+ls(i-1,j+1,k+1)+ls(i+1,j+1,k-1)
     %              +ls(i+1,j+1,k+1)+2.0d0*(ls(i-1,j+1,k)+ls(i+1,j+1,k)
     %              +ls(i,j+1,k-1)+ls(i,j+1,k+1))+4.0d0*ls(i,j+1,k)
               my = (mm1 - mm2)/32.d0
               
               mm1 = ls(i-1,j-1,k-1)+ls(i-1,j+1,k-1)+ls(i+1,j-1,k-1)
     %              +ls(i+1,j+1,k-1)+2.0d0*(ls(i-1,j,k-1)+ls(i+1,j,k-1)
     %              +ls(i,j-1,k-1)+ls(i,j+1,k-1))+4.0d0*ls(i,j,k-1)
               mm2 = ls(i-1,j-1,k+1)+ls(i-1,j+1,k+1)+ls(i+1,j-1,k+1)
     %              +ls(i+1,j+1,k+1)+2.0d0*(ls(i-1,j,k+1)+ls(i+1,j,k+1)
     %              +ls(i,j-1,k+1)+ls(i,j+1,k+1))+4.0d0*ls(i,j,k+1)
               mz = (mm1 - mm2)/32.D0
c     *(2)*  
               mx = DABS(MX)
               MY = DABS(MY)
               MZ = DABS(MZ)
               normL1 = mx+my+mz
               mx = mx/normL1
               my = my/normL1
               mz = mz/normL1

               alpha = ls(i,j,k)/normL1
c     *(3)*  
               alpha = alpha + 0.5d0
c     *(4)*  
               if(alpha.ge.1.d0) then 
                  cc(i,j,k) = 1.d0
               else if (alpha.le.0.d0) then
                  cc(i,j,k) = 0.d0
               else 
                  cc(i,j,k) = FL3D(mx,my,mz,alpha,zero,one)
c                  write(*,*) cc(i,j,k)
               end if
            enddo
         enddo
      enddo 
      return
      end