program square
c**********************************************************************
c
c                          MAIN PROGRAM                               
c
c     solve the navier-stokes equations for two incompressible fluids 
c     in a rectangular region with a uniform gird. The two fluids are
c     separated by an interface
c
c     parameters:                                                     
c
c     tend:                            end of calculation             
c     rey:                             Reynolds number                
c     cfl:                             cfl number
c     (nx,ny):                         grid size
c
c
c**********************************************************************
      integer nx,ny,nxy,karea,nxm,nym,nm,nme,mxy,mxym,
     +        nrew,npit,kpit,nout,i,nstep,kk,kunit,krei
      logical piter
      include 'gridsz'
      parameter (nrew=20,npit=2,karea=5,krei=10)
      parameter (nm=nxm*nym, nme=(nxm+1)*(nym+1))
      double precision uf(nm),vf(nm),px(nm),py(nm),
     +                 pf(0:nxm+1,0:nym+1)
      double precision us(nm),vs(nm),
     +                 fx(nm),fy(nm),wp(nm),wy(nm),
     +                 ux(nme),uy(nme),vx(nme),vy(nme)
      double precision tout(10)
      double precision pi,tend,rey,cfl,h,hinv,hhinv,dtin,
     +                 dt,t,ek,dth,ab,rin,ru,rb,grav
      common/reynod/rey
      common/meshsz/h,hinv
      common/hfract/hhinv
      common/timest/dt
      common/timspc/dth
      common/valupi/pi
      common/gravity/grav
      data kpit,kunit/0,20/, t,ru,rb,grav/0.D0,1.D0,0.5D0,10.D0/
      data piter/.true./
c     *********************************************************

      read (5,*) nx,ny
      read (5,*) tend,rey,cfl
      read (5,*) nout
      read (5,*) (tout(i),i=1,nout)

      pi = 4.D0*datan(1.D0)
      nxy = nx*ny
      mxy = 2*nxy
      hinv = dble(ny)
      hhinv = 0.5D0*hinv
      h  = 1.D0/hinv
      dtin = 0.5D0*h
      nstep = idint((tend+1.D-10)/dtin)
      rin = 0.1D0

c     initialize the flow

      call stflow(uf,vf,px,py,nx,ny)
      call inifrt(pf,rin,nx,ny)
      call densit(pf,ru,rb,nx,ny)
      call output(pf,uf,vf,nx,ny,kunit)

c*******************************************************************
c
c---                 -- time evolution --
c   
c      do 5 kk=1,nstep+npit
      kk = 0
 5    continue
      call chtmst(uf,vf,ek,dtin,cfl,nx,ny)
      dth = dt*hinv
      ab = 0.5D0*dt*hinv*hinv/rey
      call edgevl(uf,vf,us,vs,ux,vx,uy,vy,nx,ny)
      call macpro(ux,vx,uy,vy,nx,ny,mxy)
      call timrhs(uf,vf,px,py,us,vs,ux,vx,uy,vy,fx,fy,ab,nx,ny)
      if( .not. piter ) call convec(uf,vf,ux,vy,pf,nx,ny)
      if( mod(kk,krei) .eq. 0 ) call reinit(pf,nx,ny)
c---  obtain the intermediate velocity u star
      call conjug(fx,us,uf,wp,wy,ab,nxy,nx,ny)
      call conjug(fy,vs,vf,wp,wy,ab,nxy,nx,ny)
c---  project to the incompressible velocity space
      call projec(uf,vf,us,vs,px,py,nx,ny,mxy,piter)
c
      if( piter ) then
        kpit = kpit + 1
        if( kpit .eq. npit) piter = .false.
      else
        t = t + dt
        kk = kk + 1
      end if
      print 110, dt,t,ek
c     .........................................................
      do 150 i=1,nout
      if ( t-tout(i) .ge. 0. .and. t-tout(i) .lt. dt) then
          kunit = kunit+1
c          call vortic(uf,vf,wp,nx,ny)
          call output(pf,uf,vf,nx,ny,kunit)
c          call stream(uf,vf,vx,nx,ny)
      endif
 150  continue

      if ( t .le. tend) goto 5
c 5    continue

c     *************************************************

 110  format(' dt=',f10.7,'   t=',f8.5,' ek=',f12.8 )

      end


      subroutine densit(pf,ru,rb,nx,ny)
      integer i,j,nxm,nym,mxym,nx,ny
      include 'gridsz'
      double precision pf(0:nx+1,0:ny+1)
      double precision roa,ru,rb
      common/density/roa(0:nxm+1,0:nym+1)
      do 10 j=0,ny+1
      do 10 i=0,nx+1
         if (pf(i,j) .ge. 0.D0) then
            roa(i,j) = rb
         else
            roa(i,j) = ru
         endif
 10   continue
      end