c   $Id: ini_conditions2phi.f,v 1.1 1998/10/20 13:23:03 zaleski Exp $
c   (c) S. Zaleski, CNRS, UPMC
c
c   $Author: zaleski $
c
      SUBROUTINE ini1disk(u,v,c,ccol,p,xvel,yvel,tau,nx,ny)
      include 'undefined.h'
      INTEGER nx , ny, dfetch,i,idrop,ibubble,ifetch
      DOUBLE PRECISION u(nx,ny),v(nx,ny),c(nx,ny),p(nx,ny),ccol(nx,ny)
      DOUBLE PRECISION tau,color,xc,yc,xvel,yvel,ylayer
      REAL*8 radius

      idrop=0
      i=ifetch("drop",idrop)
      ibubble=0
      i=ifetch("bubble",ibubble)

      if(ibubble+idrop.eq.2)  stop "can t have both drop and bubble"
      if(ibubble+idrop.eq.0) stop "Not a drop/bubble"

      i=dfetch("radius",radius)
      write(6,*) "radius=",radius
      i=dfetch("xcenter",xc)
      i=dfetch("ycenter",yc)
      write(6,*) "xcenter=",xc
      write(6,*) "ycenter=",yc
      write(6,*) "yvel=",yvel
      i=dfetch("layer-width",ylayer)
      color=1.d0
c      call diskpre(u,v,c,p,xvel,yvel,xc,yc,radius,tau,nx,ny)
c      call creatdisk(u,v,c,p,xvel,yvel,xc,yc,color,radius,tau,nx,ny)
      call inilayerdisk(u,v,c,ccol,p,ylayer,xvel,yvel,xc,yc,color,radius
     * ,tau,nx,ny)
c      call creatdisk(u,v,c,p,ylayer,xvel,yvel,xc,yc,color,radius,
c     *tau,nx,ny)
      return
      end

      SUBROUTINE creatdisk(u,v,c,p,xvel,yvel,xc,yc,idrop,r,tau,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny), v(nx,ny), c(nx,ny), p(nx,ny)
      DOUBLE PRECISION xvel, yvel,xc,yc,r,x,y, h, tau , color
      INTEGER i,j, idrop
      h = 1.d0/(nx-2)
      do j=1,ny
         do i=1,nx
            y = h*(j - ny/2) - yc
            x = h*(i - nx/2) - xc
            if(x*x + y*y .le. r*r) then
               c(i,j) = 1.d0*idrop
               u(i,j) = xvel
               v(i,j) = yvel
            else
               u(i,j) =0.d0
               v(i,j) =0.d0
               c(i,j) = (1 - idrop)*1.d0
            endif
         enddo
      enddo
      call makebccf(c,nx,ny) 
      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

c  override C with RS's ellipse
c      call doellipse(c,nx,ny)
      return
      end

      SUBROUTINE inilayerdisk(u,v,c,ccol,p,yl,xvel,yvel,xc,yc,color,r,
     * tau,nx,ny)

      include 'undefined.h'
      INTEGER nx,ny,ifetch,flagbubble,dfetch
      DOUBLE PRECISION u(nx,ny),v(nx,ny),c(nx,ny),p(nx,ny),ccol(nx,ny)
      DOUBLE PRECISION xvel, yvel,xc,yc,r,yl,x,y, h, tau , color
      INTEGER i,j,iw
      DOUBLE PRECISION width,height

      h = 1.d0/(nx-2)
              write(6,*) "yvel/iuni=",yvel
      if(ifetch("bubble",flagbubble).ne.1) then
         color =color
      else 
         color=1.0d0-color
      endif
      write(6,*) "color=",color
      do j=1,ny
         do i=1,nx
               c(i,j) =1.d0-color
               u(i,j) = 0.d0
               v(i,j) = 0.d0
         enddo
      enddo
      do j=1,ny
         do i=1,nx
            y=h*((j*1.)-1.5)
            if (y .le. yl) then
	       c(i,j) = color
               ccol(i,j) = 1.d0-color
               u(i,j) = 0.d0
               v(i,j) = 0.d0
           endif
         enddo
      enddo
      
      if(dfetch("width",width).ne.1) then 
        do j=1,ny
         do i=1,nx
            y = h*(j - ny/2) - yc
            x = h*(i - nx/2) - xc
            if (x*x + y*y .le. r*r) then
               c(i,j) = color
               ccol(i,j) = color
               u(i,j) = 0.0d0
               v(i,j) = yvel
            endif
         enddo
        enddo

      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
       i=dfetch("height",height)
       width=width*h
       height=height*h
       do j=1,ny
         do i=1,nx
            y = h*(j - ny/2) - yc
            x = h*(i - nx/2) - xc
            if (x*x + y*y .le. r*r) then
               c(i,j) = color
               ccol(i,j) = color
               u(i,j) = 0.0d0
               v(i,j) = yvel
            endif
            y=h*j
            if(y.ge. yl+width/2.) then
               if(y.le.yl+width) then
                  if(x .le.height) then
                    c(i,j) = color
                    ccol(i,j) = color
                    u(i,j) = 0.0d0
                    v(i,j) = yvel
                  endif
               endif
            endif
            if(y.ge. yl) then
               if(y.le.yl+width/2.) then
                  if(x .le.height) then
                     c(i,j) = color
                     ccol(i,j) = 1.d0-color
                    u(i,j) = 0.0d0
                    v(i,j) = 0.0d0
                  endif
               endif
            endif 
         enddo
      enddo
         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
      endif
      return
      end


      SUBROUTINE inisquare(u,v,c,p,tau,xvel,yvel,xc,yc,r,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny), v(nx,ny), c(nx,ny), p(nx,ny)
      DOUBLE PRECISION xvel, yvel,xc,yc,r,x,y,tau,h
      INTEGER i,j
      h = 1.d0/(nx - 2)
      do j=1,ny
         do i=1,nx
            y = float(j -nx/2)/(nx-2) - yc
            x = float(i - ny/2)/(nx-2) - xc
            c(i,j) = 0.d0
            u(i,j) = xvel*tau/h
            v(i,j) = yvel*tau/h
            p(i,j) = 0.d0
            if( dmax1(dabs(x),dabs(y)) .le. r) then
               c(i,j) = 1.d0
            endif
         enddo
      enddo
      return
      end


      SUBROUTINE inirotrec(u,v,c,p,tau,xvel,yvel,xc,yc,r,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny, ifetch, dfetch, invert,tmax
      DOUBLE PRECISION u(nx,ny), v(nx,ny), c(nx,ny), p(nx,ny)
      DOUBLE PRECISION xvel, yvel,xc,yc,r,x,y,tpi,tau,h,arot
      INTEGER i,j
      h = 1.d0/(nx - 2)
      i = ifetch("tmaxrot",tmax)
      arot=1.0
      i = dfetch("anglerot",arot)
      tpi = arot*dacos(0.d0)*4.d0*(nx-2)/tmax
      invert=1
      i = ifetch("invert",invert)
      write(6,*) "tpi=",tpi
      do j=1,ny
         do i=1,nx
            y = float(j -nx/2)/(nx-2) - yc
            x = float(i - ny/2)/(nx-2) - xc
            c(i,j) = REAL(invert - 1)
            u(i,j) = -tpi*y
            v(i,j) =  tpi*x
            p(i,j) = 0.d0
            if(dabs(x) .le. 1.66666666666d0*r
     $                .and.dabs(y).lt.r) then
               c(i,j) = REAL(invert)
            endif
         enddo
      enddo
      return
      end
c
c
      SUBROUTINE inidisk(u,v,c,p,xvel,yvel,xc,yc,r,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny), v(nx,ny), c(nx,ny), p(nx,ny)
      DOUBLE PRECISION xvel, yvel,xc,yc,r,x,y
      INTEGER i,j
cld
c      do j=1,ny
c         do i=1,nx
c            y = dfloat(j -ny/2)/(ny-2) - yc
c            x = dfloat(i -nx/2)/(nx-2) - xc
c            c(i,j) = 1.0d0
c            if(x*x + y*y .le. r*r) then
c               c(i,j) = 0.0d0
c               u(i,j) = xvel
c               v(i,j) = yvel
c            endif
c            if (j.gt.2*ny/3) c(i,j) = 0.0d0
c         enddo
c      enddo
cld
      do j=1,ny
         do i=1,nx
            y = dfloat(j) - yc
            x = dfloat(i) - xc
            c(i,j) = 1.0d0
            if(x*x + y*y .le. r*r) then
               c(i,j) = 0.0d0
               u(i,j) = xvel
               v(i,j) = yvel
            endif
cld            if (j.gt.2*ny/3) c(i,j) = 0.0d0
            if (j.gt.ny/3) c(i,j) = 0.0d0
	 enddo
      enddo
      return
      end


      SUBROUTINE inizero(u,v,c,p,tau,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny), v(nx,ny), c(nx,ny), p(nx,ny)
      DOUBLE PRECISION xc,yc,x,tau,r1,h
      INTEGER i,j
      r1=0.3056d0
      h=1.d0/(nx-2)
      do j=1,ny
         do i=1,nx
            c(i,j)= 0.d0
            u(i,j)= 0.d0
            v(i,j)= 0.d0
            p(i,j)= 0.d0
         enddo
      enddo
      return
      end