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
c***********************************************************************
      SUBROUTINE ciam_raus(u,w,c,flxm,flxp,flzm,flzp,
     %     divx,divz,nx,nz)
c***
c     Solve: dc/dt = -div(uc) + c div(u) (c changes only at cells
c     where grad(c) is not zero)
c***
      INCLUDE 'undefined.h'
      INTEGER i,k,inv,nx,nz
c      DOUBLE PRECISION u(nx,nz),w(nx,nz),c(nx,nz),cdiv(nx,nz)
      DOUBLE PRECISION u(nx,nz),w(nx,nz),c(nx,nz)
      DOUBLE PRECISION divx(nx,nz), divz(nx,nz)
      DOUBLE PRECISION flxm(nx,nz),flxp(nx,nz),flzm(nx,nz),flzp(nx,nz)
      DOUBLE PRECISION cc,u1,u2,w1,w2,normx,normz,mrx,mrz,mx,mz
      DOUBLE PRECISION mm1,mm2,alfac,alfa,are1,are2,VOL2
      DOUBLE PRECISION du1,du2,dw1,dw2 
      DOUBLE PRECISION cx,cz,cm,cd
      INTRINSIC DMAX1,DMIN1,DABS,DSQRT
      EXTERNAL VOL2
c***      
      do k=2,nz-1
         do i=2,nx-1
            u1 = u(i,k)
            u2 = u(i+1,k)
            w1 = w(i,k)
            w2 = w(i,k+1)
            if (c(i,k) .eq. 0.d0) then
               flxm(i,k) = 0.d0
               flxp(i,k) = 0.d0
               flzm(i,k) = 0.d0
               flzp(i,k) = 0.d0
               divx(i,k) = 0.d0
               divz(i,k) = 0.d0
            elseif (c(i,k) .eq. 1.d0) then
               du1 = 0.5d0*(DABS(u1) - u1)
               du2 = 0.5d0*(DABS(u2) + u2)
               dw1 = 0.5d0*(DABS(w1) - w1)
               dw2 = 0.5d0*(DABS(w2) + w2)
               flxm(i,k) = du1
               flxp(i,k) = du2
               flzm(i,k) = dw1
               flzp(i,k) = dw2
               divx(i,k) = (1.d0 - (du1+u1) + (u2-du2))
               divz(i,k) = (1.d0 - (dw1+w1) + (w2-dw2)) 

c            if (c(i,k) .eq. 0.d0 .or. c(i,k) .eq. 1.d0) then
c               du1 = 0.5d0*(DABS(u1) - u1)
c               du2 = 0.5d0*(DABS(u2) + u2)
c               dw1 = 0.5d0*(DABS(w1) - w1)
c               dw2 = 0.5d0*(DABS(w2) + w2)
c               flxm(i,k) = c(i,k)*du1
c               flxp(i,k) = c(i,k)*du2
c               flzm(i,k) = c(i,k)*dw1
c               flzp(i,k) = c(i,k)*dw2
cc               cdiv(i,k) = c(i,k)*(1.d0 - u1 - du1 + u2 - du2 - 
cc     %              w1 - dw1 + w2 - dw2)
c               divx(i,k) = c(i,k)*(1.d0 - (du1+u1) + (u2-du2))
c               divz(i,k) = c(i,k)*(1.d0 - (dw1+w1) + (w2-dw2)) 
            else 
c***  
c     find the normal to the interface, its normalized components
c     and alfa(c)
c***  
               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)
               normx = 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)
               normz = mm1 - mm2
c***  
               mrx = DABS(normx) + 1.d-50 
               mrz = DABS(normz) + 1.d-50
               mm2 = DMAX1(mrx,mrz)
               mrx = mrx/mm2
               mrz = mrz/mm2
c***
               mm1 = DMIN1(mrx,mrz)
               are1 = 0.5d0*mm1
               are2 = 1.0d0 - are1
               if (c(i,k) .le. are1) then
                  alfac = DSQRT(2.0d0*c(i,k)*mm1)
               elseif (c(i,k) .le. are2) then
                  alfac = c(i,k) + 0.5d0*mm1
               else
                  alfac = mm1 + 1.0d0 - DSQRT(2.0d0*(1.0d0-c(i,k))*mm1)
               endif
c***  
c     fluxes along x: symmetry, the new equation for the interface after
c     advection along x, fluxes to the left (flxm), to the right 
c     (flxp), inside variation (divx) and finally invert, if necessary
c***  
               inv = 0
               if (normx .lt. 0.0d0) then
                  mm1 = -u1
                  u1 = -u2
                  u2 = mm1
                  inv = 1
               endif
c***  
               mx = mrx/(1.0d0-u1+u2)
               mz = mrz
               alfa = alfac + mx*u1
c***  
               du1 = 0.5d0*(DABS(u1) - u1)
               mm1 = alfa + mx*du1
               flxm(i,k) = vol2(mx,mz,mm1,du1)
               
               du2 = 0.5d0*(DABS(u2) + u2)
               mm2 = alfa - mx
               flxp(i,k) = vol2(mx,mz,mm2,du2)

               mm1 = (1.d0 - (du1+u1) + (u2-du2))
               mm2 = alfa - mx*(du1+u1)
               divx(i,k) = vol2(mx,mz,mm2,mm1)
c***  
               if (inv .eq. 1) then
                  mm1 = flxm(i,k)
                  flxm(i,k) = flxp(i,k)
                  flxp(i,k) = mm1
               endif
c***  
c     fluxes along z: symmetry, the new equation for the interface after
c     advection along z, fluxes to the bottom (flzm), to the top 
c     (flzp), inside variation (divz) and finally invert, if necessary
c***  
               inv = 0
               if (normz .lt. 0.0d0) then
                  mm1 = -w1
                  w1 = -w2
                  w2 = mm1
                  inv = 1
               endif
c***  
               mz = mrz/(1.0d0-w1+w2)
               mx = mrx
               alfa = alfac + mz*w1
c***  

               dw1 = 0.5d0*(DABS(w1) - w1)
               mm1 = alfa + mz*dw1
               flzm(i,k) = vol2(mz,mx,mm1,dw1)

               dw2 = 0.5d0*(DABS(w2) + w2)
               mm2 = alfa - mz
               flzp(i,k) = vol2(mz,mx,mm2,dw2)

               mm1 = (1.d0 - (dw1+w1) + (w2-dw2)) 
               mm2 = alfa - mz*(dw1+w1)
               divz(i,k) = vol2(mz,mx,mm2,mm1)
c***  
               if (inv .eq. 1) then
                  mm1 = flzm(i,k)
                  flzm(i,k) = flzp(i,k)
                  flzp(i,k) = mm1
               endif
            endif
         enddo
      enddo
c***
c     BC on the fluxes
c***
      call makebcflux(flxm,flxp,flzm,flzp,nx,nz)
c***
c     update colour function, clip it: 0<c<1 and BC
c***  
      do k=2,nz-1
         do i=2,nx-1
            cx = flxp(i-1,k)+divx(i,k)+flxm(i+1,k)
            cz = flzp(i,k-1)+divz(i,k)+flzm(i,k+1)
            cc = cx + cz - c(i,k)
            cm = 0.5d0*(DABS(cc) + cc)
            cd = cm - 1.d0
            c(i,k) = cm - 0.5d0*(DABS(cd) + cd)
         enddo
      enddo
      call makebccf(c,nx,nz)
c***  
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION VOL2(mx,mz,alfa,b)
      include 'undefined.h'
      DOUBLE PRECISION mx,mz,alfa,alfaa,b,mb
c***      
c     0<mx or 0<mz 
c***
      mb = mx*b
      if (alfa .le. 0.0d0) then
         vol2 = 0.d0
      elseif (alfa .ge. (mb+mz) .or. b .eq. 0.0d0) then
         vol2 = b
      else
         if (mb .le. mz) then
            if (alfa .lt. mb) then
               vol2 = 0.5d0*alfa**2/(mx*mz)
            elseif (alfa .gt. mz) then
               alfaa = mb + mz - alfa
               vol2 = b - 0.5d0*alfaa**2/(mx*mz)  
            else
               vol2 = b*(alfa-0.5d0*mb)/mz
            endif
         else
            if (alfa .lt. mz) then
               vol2 = 0.5d0*alfa**2/(mx*mz)
            elseif (alfa .gt. mb) then
               alfaa = mb + mz - alfa
               vol2 = b - 0.5d0*alfaa**2/(mx*mz)  
            else
               vol2 = (alfa-0.5d0*mz)/mx
            endif
         endif
      endif
c***  
      return
      end
cc***********************************************************************
c      SUBROUTINE makebcflux(flxm,flxp,flzm,flzp,nx,nz)
cc***
c      INCLUDE 'undefined.h'
c      INTEGER i,k,nx,nz
c      DOUBLE PRECISION flxm(nx,nz),flxp(nx,nz),flzm(nx,nz),flzp(nx,nz)
cc***
c      do k=2,nz-1
c         flxp(1,k) = 0.d0
c         flxm(nx,k) = 0.d0
ccc         flxp(1,k) = flxp(nx-1,k)
ccc         flxm(nx,k) = flxm(2,k)
c      enddo 
c      do i=2,nx-1
c         flzp(i,1) = 0.d0
c         flzm(i,nz) = 0.d0
c      enddo 
cc***
c      return
c      end