c two-d driven cavity flow, MAC, R_K
c Possion solve by FFT
c boundary points at n=1 and n=nxp, dx=1/nx, staggered grid

      parameter(n=32,np=n+1)
      implicit real(a-h,o-z)
      dimension u (0:n, 0:np),
     /          u1(0:n, 0:np),
     /          u2(0:n, 0:np),
     /          u3(0:n, 0:np),
     /           f(0:n, 0:np),
     /          v (0:np,0:n),
     /          v1(0:np,0:n),
     /          v2(0:np,0:n),
     /          v3(0:np,0:n),
     /           g(0:np,0:n),
     /           p(0:n, 0:n)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)

      nx=n
      ny=nx
      nxm=nx-1
      nym=ny-1
      nxp=nx+1
      nyp=ny+1

      call setup(t,kmax,u,v,f,g)

      do 41 j=0,nyp
      do 41 i=0,nx
        u1(i,j)=u(i,j)
 41   continue

      do 42 j=0,ny
      do 42 i=0,nxp
        v1(i,j)=v(i,j)
 42   continue

      open(1,file='history.m') 
      write(1,*) 'his=['

c time marching

      write(6,*) 'step, time, energy'
 7777 do 1111 kcount=1,kmax
      if(kcount/3*3.eq.kcount) write(6,*)  kcount, t, eng
      call field(u,v,p,f,g)

      do 100 j=1,ny
      do 100 i=1,nxm
        u (i,j)=u1(i,j)+0.5*dt*(f(i,j)-(p(i+1,j)-p(i,j))/dx)
        u2(i,j)=u (i,j)
 100  continue     

      do 110 j=1,nym
      do 110 i=1,nx
        v (i,j)=v1(i,j)+0.5*dt*(g(i,j)-(p(i,j+1)-p(i,j))/dy)
        v2(i,j)=v (i,j)
 110  continue     

      call field(u,v,p,f,g)

      do 120 j=1,ny
      do 120 i=1,nxm
        u (i,j)=u1(i,j)+0.5*dt*(f(i,j)-(p(i+1,j)-p(i,j))/dx)
        u3(i,j)=u (i,j)
 120  continue     

      do 130 j=1,nym
      do 130 i=1,nx
        v (i,j)=v1(i,j)+0.5*dt*(g(i,j)-(p(i,j+1)-p(i,j))/dy)
        v3(i,j)=v (i,j)
 130  continue     

      call field(u,v,p,f,g)

      do 140 j=1,ny
      do 140 i=1,nxm
        u(i,j)=u1(i,j)+dt*(f(i,j)-(p(i+1,j)-p(i,j))/dx)
 140  continue     

      do 150 j=1,nym
      do 150 i=1,nx
        v(i,j)=v1(i,j)+dt*(g(i,j)-(p(i,j+1)-p(i,j))/dy)
 150  continue     

      call field(u,v,p,f,g)

      do 160 j=1,ny
      do 160 i=1,nxm
        u(i,j)=(-u1(i,j)+u2(i,j)+2.*u3(i,j)+u(i,j))/3.
     /         +dt/6.*(f(i,j)-(p(i+1,j)-p(i,j))/dx)
        u1(i,j)=u(i,j)
 160  continue     

      do 170 j=1,nym
      do 170 i=1,nx
        v(i,j)=(-v1(i,j)+v2(i,j)+2.*v3(i,j)+v(i,j))/3.
     /         +dt/6.*(g(i,j)-(p(i,j+1)-p(i,j))/dy)
        v1(i,j)=v(i,j)
 170  continue     

      eng=0.
      do 180 j=2,ny
      do 180 i=2,nx
        eng = eng + u(i,j)*u(i,j) + v(i,j)*v(i,j)
 180  continue     
      eng=sqrt(eng*dx*dy)

      t=t+dt
      write(1,*) t, eng
 1111 continue

      call out(u,v)

      write(6,*)'continue?  yes=1, no=0'
      read(5,*)icont
      if(icont.eq.0)go to 8888
      write(6,*)'read next tnew'
      read(5,*)tnew
      kmax=ifix(tnew/dt)
      write(6,*) 'new kmax=',kmax
      go to 7777

 9999 format(2e11.2)
 8888 rewind(7)
      write(7,*) t
      write(7,*) ((u(i,j), i=0,nx ),j=0,nyp)
      write(7,*) ((v(i,j), i=0,nxp),j=0,ny )


      write(1,*) '];'
      close(1)

      stop
      end

      subroutine setup(t,kmax,u,v,f,g)
      implicit real(a-h,o-z)
      parameter(n=32)
      dimension u (0:nx, 0:ny+1),
     /           f(0:nx, 0:ny+1),
     /          v (0:nx+1,0:ny),
     /           g(0:nx+1,0:ny)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)

      pi=3.14159265358979
      rnu=1./100.
      cfl1=0.5
      cfl2=0.75

      dx=1./float(nx)
      dy=1./float(ny)
      dt1=cfl1*dx*dx/rnu
      dt2=cfl2*dx
      if(dt1.lt.dt2) then
        dt=dt1
      else
        dt=dt2
      endif
      dtnu=dt*rnu

c initialize

      call vcosqi(nx,wsave)

      do 10 i=0,nxm
        si=pi*float(i)/float(nx)/2.
        s(i)=-4.*sin(si)*sin(si)/dx/dx
 10    continue

      do 15 i=0,nx
        xi=dx*float(i)
c        ud(i)= 16.*xi*xi*(1.-xi)*(1.-xi)
        ud(i)= 1.
 15   continue

      do 20 j=0,nyp
      yj=dy*float(j)-dy/2.
      do 20 i=0,nx
c        u(i,j)=-ud(i)*(2.*yj-3.*yj*yj)
        u(i,j)=0.
        f(i,j)=0.
 20   continue

      do 30 i=0,nxp
      xi=dx*float(i)-dx/2.
      do 30 j=0,ny
        yj=dy*float(j)
c        v(i,j)= 16.*(2.*xi-6.*xi*xi+4.*xi*xi*xi)*(yj*yj-yj*yj*yj)
        v(i,j)=0.
        g(i,j)=0.
 30   continue

      t=0
      write(6,*) 'irest'
      read(5,*) irest
      if(irest.eq.1) then
        read(7,*) t
        read(7,*) ((u(i,j), i=0,nx ),j=0,nyp)
        read(7,*) ((v(i,j), i=0,nxp),j=0,ny )
      endif

      write(6,*) 'input tnew'
      read(5,*) tnew
      kmax=ifix(tnew/dt)
      write(6,*) 'kmax=', kmax

      return
      end

      subroutine out(u,v)
      implicit real(a-h,o-z)
      parameter(n=32,np=n+1)
      dimension u(0:nx, 0:ny+1),
     /          v(0:nx+1,0:ny),
     /        phi(0:np,0:np)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)

      call bd(u,v)

      i=0
      do 200 j=0,ny
        phi(i,j)=0.
 200  continue     

      do 210 i=1,nx
      do 210 j=0,ny
        phi(i,j)=phi(i-1,j)-dx*v(i,j)
 210  continue     

      open(8,file='out.m')
      write(8,*) 'out=['
      do 220 i=0,nx
      do 220 j=0,ny
        omi=(v(i+1,j)-v(i,j))/dx -(u(i,j+1)-u(i,j))/dy
        write(8,*) phi(i,j), omi
 220  continue     
      write(8,*) '];'
      close(8)

      return
      end

      subroutine bd(u,v)
      implicit real(a-h,o-z)
      parameter(n=32)
      dimension u (0:nx, 0:ny+1),
     /          v (0:nx+1,0:ny)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)

      i=0
      do 151 j=0,ny
        u(i,j)=0.
        v(i,j)=-v(i+1,j)
 151  continue     

      i=nx
      do 152 j=0,ny
        u(i,  j)=0.
        v(i+1,j)=-v(i,j)
 152  continue     

      j=0
      do 153 i=0,nx
        u(i,j)=-u(i,j+1)
        v(i,j)=0.
 153  continue     

      j=ny
      do 154 i=0,nx
        u(i,j+1)=-u(i,j) + 2.*ud(i)
        v(i,j  )=0.
 154  continue     
      return
      end

      subroutine field(u,v,p,f,g)
      implicit real(a-h,o-z)
      parameter(n=32)
      dimension u (0:nx, 0:ny+1),
     /           f(0:nx, 0:ny+1),
     /          v (0:nx+1,0:ny),
     /           g(0:nx+1,0:ny),
     /           p(0:n, 0:n)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)

      call bd(u,v)

c     viscosity + convection terms

      do 61 j=1,ny
      do 61 i=1,nxm
        f(i,j)=rnu*(u(i+1,j)+u(i,j+1)
     /             +u(i-1,j)+u(i,j-1)-4.*u(i,j))/dx/dx
     /     -u(i,j)*(u(i+1,j)-u(i-1,j))/(2.*dx)
     /        -.25*(v(i,j-1)+v(i,j)+v(i+1,j-1)+v(i+1,j))
     /            *(u(i,j+1)-u(i,j-1))/(2.*dy)
 61   continue     

      do 62 j=1,nym
      do 62 i=1,nx
        g(i,j)=rnu*(v(i+1,j)+v(i,j+1)
     /             +v(i-1,j)+v(i,j-1)-4.*v(i,j))/dx/dx
     /        -.25*(u(i-1,j)+u(i-1,j+1)+u(i,j)+u(i,j+1))
     /            *(v(i+1,j)-v(i-1,j))/(2.*dx)
     /     -v(i,j)*(v(i,j+1)-v(i,j-1))/(2.*dy)
 62   continue     
 
c pressure

      do 70 j=1,ny
      do 70 i=1,nx
        p(i,j)=(f(i,j)-f(i-1,j))/dx 
     /       + (g(i,j)-g(i,j-1))/dy
 70   continue

      call poisson(p)
   
      return
      end

      subroutine poisson(u)
      implicit real(a-h,o-z)
      parameter(n=32)
      dimension u(0:nx,0:ny)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)
      
      call cos2d(u,-1)

      u(1,1)=0.

      i=1
      do 20 j=2,ny
        u(i,j)=u(i,j)/(s(i-1)+s(j-1))
 20   continue

      do 30 j=1,ny
      do 30 i=2,nx
        u(i,j)=u(i,j)/(s(i-1)+s(j-1))
 30   continue
          
 1111 call cos2d(u,1)

      return
      end

      subroutine cos2d(u,isign)
      implicit real(a-h,o-z)
      parameter(n=32)
      dimension u(0:nx,0:ny),
     /          aux(n,n),
     /          bux(n,n)
      common /geo/ nx,ny,nxm,nym,nxp,nyp,dx,dy,dt
      common /dat/ rnu,ud(0:200)
      common /fft/ pi,s(0:200),wsave(2000)

      if(isign.eq.1) goto 1111
      if(isign.eq.-1) goto 2222

 1111 do 20 i=1,nx
      do 20 j=1,ny
        aux(i,j)=u(j,i)
 20   continue
      call vcosqf(nx,ny,aux,bux,nx,wsave)
      do 30 i=1,nx
      do 30 j=1,ny
        u(i,j)=aux(j,i)
 30   continue

      do 40 j=1,ny
      do 40 i=1,nx
        aux(i,j)=u(i,j)
 40   continue
      call vcosqf(nx,ny,aux,bux,ny,wsave)
      do 50 j=1,ny
      do 50 i=1,nx
        u(i,j)=aux(i,j)
 50   continue
      goto 3333

 2222 do 60 i=1,nx
      do 60 j=1,ny
        aux(i,j)=u(j,i)
 60   continue
      call vcosqb(nx,ny,aux,bux,nx,wsave)
      do 70 i=1,nx
      do 70 j=1,ny
        u(i,j)=aux(j,i)
 70   continue

      do 80 j=1,ny
      do 80 i=1,nx
        aux(i,j)=u(i,j)
 80   continue
      call vcosqb(nx,ny,aux,bux,nx,wsave)
      do 90 j=1,ny
      do 90 i=1,nx
        u(i,j)=aux(i,j)
 90   continue
      goto 3333

 3333 return
      end