c$$$    SURFER fluid interface simulation program
c$$$    Copyright (C) 2001 Stephane Zaleski and others
c$$$
c$$$    This library is free software; you can redistribute it and/or
c$$$    modify it under the terms of the GNU Lesser General Public
c$$$    License as published by the Free Software Foundation; either
c$$$    version 2.1 of the License, or (at your option) any later version.
c$$$
c$$$    This library is distributed in the hope that it will be useful,
c$$$    but WITHOUT ANY WARRANTY; without even the implied warranty of
c$$$    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
c$$$    Lesser General Public License for more details.
c$$$
c$$$    You should have received a copy of the GNU Lesser General Public
c$$$    License along with this library; if not, write to the Free Software
c$$$    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
c$$$
c$$$    Stephane Zaleski zaleski@lmm.jussieu.fr
c
      SUBROUTINE x_swp(u,c,vof1,vof2,vof3,nx,nz)

      INCLUDE 'undefined.h'
      INTEGER i,k,inv,nx,nz
      DOUBLE PRECISION u(nx,nz),c(nx,nz),mx,mz,alfa
      DOUBLE PRECISION s1,s2,mm1,mm2,V1,V2,vol2
      DOUBLE PRECISION vof1(nx,nz),vof2(nx,nz),vof3(nx,nz)
      
       
      
      do k=2,nz-1
            do i=2,nx-1
                   s1=u(i,k)
                   s2=u(i+1,k)

               if(c(i,k).eq.0.0d0)then
                     vof1(i,k)=0.0d0
                     vof2(i,k)=0.0d0
                     vof3(i,k)=0.0d0
               else if(c(i,k).eq.1.0d0)then
                     vof1(i,k)=dmax1(-s1,0.0d0)
                     vof2(i,k)=1.0d0-dmax1(s1,0.0d0)+dmin1(s2,0.0d0)
                     vof3(i,k)=dmax1(s2,0.0d0)
               else

c   normal vector of the interface
                     mm1=C(i-1,k-1)+2.0d0*C(i-1,k)+C(i-1,k+1)
                     mm2=C(i+1,k-1)+2.0d0*C(i+1,k)+C(i+1,k+1)
                     mx=mm1-mm2
                     mm1=C(i-1,k-1)+2.0d0*C(i,k-1)+C(i+1,k-1)
                     mm2=C(i-1,k+1)+2.0d0*C(i,k+1)+C(i+1,k+1)
                     mz=mm1-mm2

                   inv=0
                   if(mx.lt.0.0d0)then
                       mm1=-s1
                       s1=-s2
                       s2=mm1
                       mx=-mx
                       inv=1
                   endif
                   mx=mx+1.d-50 
                   mz=dabs(mz)+1.d-50
                   mm2=dmax1(mx,mz)
                   mx=mx/mm2
                   mz=mz/mm2

c computation of alfa to determine the interface equation

                   mm1=dmin1(mx,mz)
                   V1=0.5d0*mm1
                   V2=1.0d0-V1
                   if(c(i,k).le.V1)then
                         alfa=dsqrt(2.0d0*c(i,k)*mm1)
                   else if(c(i,k).le.V2)then
                         alfa=c(i,k)+0.5*mm1
                   else
                         alfa=mm1+1.0d0-dsqrt(2.0d0*(1.0d0-c(i,k))*mm1)
                   endif
c the new interface equation after advection

                   mx=mx/(1.0d0-s1+s2)
                   alfa=alfa+mx*s1
 
c calcul vof

                   mm1=dmax1(-s1,0.0d0)
                   mm2=alfa+mx*mm1
                   vof1(i,k)=vol2(mx,mz,mm2,mm1)

                   mm1=1.0d0-dmax1(s1,0.0d0)+dmin1(s2,0.0d0)
                   mm2=alfa-mx*dmax1(s1,0.0d0)
                   vof2(i,k)=vol2(mx,mz,mm2,mm1)

                   mm1=dmax1(s2,0.0d0)
                   mm2=alfa-mx
                   vof3(i,k)=vol2(mx,mz,mm2,mm1)
c symmetry
                   if(inv.eq.1) then
                        mm1=vof1(i,k)
                        vof1(i,k)=vof3(i,k)
                        vof3(i,k)=mm1
                   endif
               endif

            enddo
      enddo

c  periodic boundary condition 

      call makebcxvof(vof1,vof3,nx,nz)

      do k=2,nz-1
            do i=2,nx-1
              c(i,k)=vof3(i-1,k)+vof2(i,k)+vof1(i+1,k)
            enddo
      enddo

      do k=2,nz-1
            do i=2,nx-1
              c(i,k)=dmin1(1.0d0,dmax1(0.0d0,c(i,k)))
            enddo
      enddo
      call makebccf(c,nx,nz)

      return
      end
      SUBROUTINE z_swp(v,c,vof1,vof2,vof3,nx,nz)

      INCLUDE 'undefined.h'
      INTEGER i,k,inv,nx,nz
      DOUBLE PRECISION v(nx,nz),c(nx,nz),mx,mz,alfa
      DOUBLE PRECISION s1,s2,mm1,mm2,V1,V2,vol2
      DOUBLE PRECISION vof1(nx,nz),vof2(nx,nz),vof3(nx,nz)


      do k=2,nz-1
            do i=2,nx-1
                   s1=v(i,k)
c                   s1=0.4d0
                   s2=v(i,k+1)
c                   s2=0.4d0

               if(c(i,k).eq.0.0d0)then
                     vof1(i,k)=0.0d0
                     vof2(i,k)=0.0d0
                     vof3(i,k)=0.0d0
               else if(c(i,k).eq.1.0d0)then
                     vof1(i,k)=dmax1(-s1,0.0d0)
                     vof2(i,k)=1.0d0-dmax1(s1,0.0d0)+dmin1(s2,0.0d0)
                     vof3(i,k)=dmax1(s2,0.0d0)
               else

c   le normal of the interface
                     mm1=C(i-1,k-1)+2.0d0*C(i-1,k)+C(i-1,k+1)
                     mm2=C(i+1,k-1)+2.0d0*C(i+1,k)+C(i+1,k+1)
                     mx=mm1-mm2
                     mm1=C(i-1,k-1)+2.0d0*C(i,k-1)+C(i+1,k-1)
                     mm2=C(i-1,k+1)+2.0d0*C(i,k+1)+C(i+1,k+1)
                     mz=mm1-mm2
c symetry
                   inv=0
                   if(mz.lt.0.0d0)then
                       mm1=-s1
                       s1=-s2
                       s2=mm1
                       mz=-mz
                       inv=1
                   endif
                   mx=dabs(mx)+1.d-50
                   mz=mz+1.d-50
                   mm2=dmax1(mx,mz)
                   mx=mx/mm2
                   mz=mz/mm2

c calcul of alfa to determine the equation of the interface                    

                   mm1=dmin1(mx,mz)
                   V1=0.5d0*mm1
                   V2=1.0d0-V1
                   if(c(i,k).le.V1)then
                         alfa=dsqrt(2.0d0*c(i,k)*mm1)
                   else if(c(i,k).le.V2)then
                         alfa=c(i,k)+0.5*mm1
                   else 
                         alfa=mm1+1.0d0-dsqrt(2.0d0*(1.0d0-c(i,k))*mm1)
                   endif

c the new equation for the interface after advection

                   mz=mz/(1.0d0-s1+s2)
                   alfa=alfa+mz*s1
 
c calcul vof
                   mm1=dmax1(-s1,0.0d0)
                   mm2=alfa+mz*mm1
                   vof1(i,k)=vol2(mz,mx,mm2,mm1)

                   mm1=1.0d0-dmax1(s1,0.0d0)+dmin1(s2,0.0d0)
                   mm2=alfa-mz*dmax1(s1,0.0d0)
                   vof2(i,k)=vol2(mz,mx,mm2,mm1)

                   mm1=dmax1(s2,0.0d0)
                   mm2=alfa-mz
                   vof3(i,k)=vol2(mz,mx,mm2,mm1)

c symetry
                   if(inv.eq.1) then
                        mm1=vof1(i,k)
                        vof1(i,k)=vof3(i,k)
                        vof3(i,k)=mm1
                   endif
               endif

            enddo
      enddo

c  free-slip boudary condition 
 
      call makebczvof(vof1,vof3,nx,nz) 
      
      do k=2,nz-1
            do i=2,nx-1
              c(i,k)=vof3(i,k-1)+vof2(i,k)+vof1(i,k+1)
            enddo
      enddo
 
      do k=2,nz-1
            do i=2,nx-1
              c(i,k)=dmin1(1.0d0,dmax1(0.0d0,c(i,k)))
            enddo
      enddo
      call makebccf(c,nx,nz)

      return
      end