c   $Id: ini_vortices.f,v 1.1 2001/03/09 12:59:00 zaleski Exp $
c   (c) S. Zaleski
c
c   $Author: zaleski $
c

      SUBROUTINE inivort(u,v,c,p,timescalesi,tau,du,nx,ny)

      include 'undefined.h'
      INTEGER nx,ny
      integer axisym
      INTEGER i,j , k,  onephase, ifetch, dfetch, jint
      DOUBLE PRECISION u(nx,ny), v(nx,ny), c(nx,ny), p(nx,ny)
      DOUBLE PRECISION y, x, h, duover2
      DOUBLE PRECISION tau,du, tpi
      DOUBLE PRECISION xc(4),yc(4),circ(4), circul , circulSI
      DOUBLE PRECISION uvort, vvort
      DOUBLE PRECISION ray2, vmax
      DOUBLE PRECISION yint, lambdaSI, timescaleSI
      DOUBLE PRECISION core, coreSI, core2
      DOUBLE PRECISION blwliq,blwgas, invliq, invgas, uliq, ugas
      double precision nnumber, mnumber, rnumber, unumber
      double precision erf
c
c Indicate with which revision we work (see initlog.f)
c
      axisym=1
      write(11,*) 
     c"$RCSfile: ini_vortices.f,v $"
      write(11,*) 
     c"$Revision: 1.1 $", "$Date: 2001/03/09 12:59:00 $"

      if(dfetch("lambda-si",lambdaSI).ne.1) i=dfetch("L-si",lambdaSI)
      coreSI = GP_dfetch("core-si",0.01,"WARNING")
      circulSI = GP_dfetch("circul-si",0.,"SEVERE")
      circul = circulSI/(lambdaSI**2/timescaleSI)
      core = coreSI/lambdaSI
      core2 = core**2
      blwliq = GP_dfetch("blwliq-si",0.1,"WARNING")/lambdaSI;
      blwgas = GP_dfetch("blwgas-si",0.1,"WARNING")/lambdaSI;
      if(blwliq.gt.0.) then 
         invliq=1./blwliq
      else
         invliq=0.
      endif
      if(blwgas.gt.0.) then 
         invgas=1./blwgas
      else
         invgas=0.
      endif
      mnumber =  GP_dfetch("viscgas-si",1.,"WARNING")
     $          /GP_dfetch("viscliq-si",1.,"WARNING")
      nnumber =  GP_dfetch("blwgas-si",0.1,"WARNING")
     $          /GP_dfetch("blwliq-si",0.1,"WARNING")
      rnumber = GP_dfetch("rhogas-si",1.,"WARNING")
     $          /GP_dfetch("rholiq-si",1.,"WARNING")
c
c  stress condition : mu_g U_g / blw_g =  mu_l U_l / blw_l
c
      unumber = nnumber/mnumber
      Ugas = du*unumber/(1.+unumber)
      Uliq = du/(1.+unumber)
      write(6,*) "U_gas = ",ugas, ",  U_liq = ", uliq

c
c Start initialisation
c
      tpi = dacos(0.d0)*4.d0
      call fill0(p,nx,ny)
      h=1.d0/(nx-2)
      duover2 = 0.5d0*du
      
cc circ est la circulation a 2 pi pres
      circ(1)=circul/tpi
      circ(2)=-circul/tpi
      circ(3)=circul/tpi
      circ(4)=-circul/tpi
cc definition de la position de l'interface liquide-gaz
      yint= -1./6.
      jint = yint/h + ny/2

cc definition des positions des centres de tourbillons
      do k=1,4
         xc(k) = 1./8. + ((k*1.)-1.)/4.
         yc(k)=-0.359375d0
      end do
      
      do i=2,nx-1
         do  j=2,ny-1
            y = (j- ny/2)*h
            x = (i - 1.5)*h

c test y. If y above yint velocity to the right, else to the left. 

            if(y.ge.yint) then
               c(i,j) = 0.d0
               u(i,j) = ugas*erf((y-yint)*invgas)
               v(i,j) = 0.d0
            else
               c(i,j) = 1.d0
               u(i,j) = uliq*erf((y-yint)*invliq)
               v(i,j) = 0.d0
            endif
            
            do k=1,4
               ray2=(x-xc(k))**2+(y-yc(k))**2

c     calcul de u due a chaque tourbillon

               if(ray2.gt.core2) then
                  uvort =-(y-yc(k))*circ(k)/ray2
                  vvort = (x-xc(k))*circ(k)/ray2
               else
                  uvort =-(y-yc(k))*circ(k)/core2
                  vvort = (x-xc(k))*circ(k)/core2
               endif

c     u and v + vortices 

               u(i,j)= u(i,j) + uvort
               v(i,j)= v(i,j) + vvort
            enddo
         enddo
      enddo
c
c    find maximum velocity on interface
c
      vmax=0.d0
      do i=2,nx-1
	  vmax = dmax1(v(i,jint),vmax)
      enddo
      vmax = vmax/du
      write(6,*) "vmax/du = ",vmax      

c     rescale all velocities
      if(axisym.eq.0) then
      do i=2,nx-1
	do j=2,ny-1
	   u(i,j) = u(i,j) * tau/h
	   v(i,j) = v(i,j) * tau/h
	   enddo
      enddo
      else  ! exchange x and y
         if(nx.ne.ny) stop
      do i=2,nx-1
	do j=2,ny-1
	   p(i,j) = v(j,i) * tau/h
	   enddo
      enddo
      do i=2,nx-1
	do j=2,ny-1
	   v(i,j) = u(j,i) * tau/h
	   enddo
      enddo
      do i=2,nx-1
	do j=2,ny-1
	   u(i,j) = p(j,i)
	   enddo
      enddo
      do i=2,nx-1
	do j=2,ny-1
	   p(i,j) = c(j,i)
	   enddo
      enddo
      do i=2,nx-1
	do j=2,ny-1
	   c(i,j) = p(i,j)
           p(i,j) = 0.d0
	   enddo
      enddo

      endif

      onephase=0
      i = ifetch("onephase",onephase)
      if(onephase.eq.1) then
         do i=1,nx
            do j=1,ny
               c(i,j) = 1.
            enddo
         enddo
      endif
      return
      end