
c
      implicit real*8 (a-h,o-z)
      parameter  (imesh = 70000)
      parameter  (inode=257,jnode=129)
      dimension  x(imesh),y(imesh),g(imesh),w(imesh),
     .           p(imesh),res(imesh)
      common/ai/ jalpha
      common/ar/ alpha(inode,jnode)
      common/bi/ mode,itmax,ires,jres,ierr,jerr,icor,jcor
      common/br/ resmax,errmax,cormax
      common/c/  press(inode),circ,fx
      common/j/  re,dt,time,tstop
      common/h/  vx(inode,jnode),vy(inode,jnode),
     .           dwx(inode,jnode),dwy(inode,jnode)
      character*60 title
c
      open( 1,file='ft01.dat')
      open( 2,file='ft02.dat')
      open( 3,file='ft03.dat')
      open( 4,file='ft04.dat')
      open( 7,file='ft07.dat')
      open( 8,file='ft08.dat')
      open( 9,file='ft09.dat')
      open(10,file='ft10.dat')
      open(11,file='ft11.dat')
      open(12,file='x.dat'   )
      open(13,file='y.dat'   )
      open(14,file='ft14.dat')
      open(15,file='ft15.dat')
      open(16,file='ft16.dat')
      open(17,file='ft17.dat')
      open(18,file='ft18.dat')
      open(19,file='ft19.dat')
      open(20,file='ft20.dat')
      open(21,file='ft21.dat')
      open(22,file='ft22.dat')
      open(98,file='ft98.dat')
      open(99,file='ft99.dat')
c
      iread     = 1
      iwrit     = 2
      iconv     = 11
c
      dum       = 0.
      nmesh     = fmesh
      itmax     = fitmax
      read (iread,501) title
      write(iwrit,502) title
      read (iread,505) re,dtinpt,tstop,dtio
      write(iwrit,504) re,dtinpt,tstop,dtio
      dt        = dtinpt
c
      mx        = nx  +1
      my        = ny  +1
      kf        = 1
c
c calculation with no restart
c
      if (istart .eq. 0) then
      call mesh(mx,my,x(kf),y(kf))
      call ic(mx,my,x(kf),y(kf),g(kf),w(kf),p(kf))
      call bc(mx,my,x(kf),y(kf),g(kf))
      nstep1    = 0
      time1     = -dt
      timeio    = 1.01*dtio
      else
      endif
c
c restart calculation
c
      if (istart .eq. 1) then
      call rstart(mx,my,x(kf),y(kf),g(kf),w(kf),p(kf),nstep1,time1,dt)
      timeio    = 0.
      else
      endif
c
c set up data for multigrid calculation
c
      if (nmesh. ne. 1) then
      do 10 n=2,nmesh
      mmx       = (mx  -1)/2  +1
      mmy       = (my  -1)/2  +1
      kc        = kf  +mx*my
      call inject(mx,my,x(kf),mmx,mmy,x(kc))
      call inject(mx,my,y(kf),mmx,mmy,y(kc))
      mx        = mmx
      my        = mmy
   10 kf        = kc
      else
      endif
c
      if (nmesh. ne. 1) then
      do 15 n=2,nmesh
      mx        = 2*(mx  -1)  +1
      my        = 2*(my  -1)  +1
   15 kf        = kf  -mx*my
      else
      endif
c
c *********************************
c Streamfunction/Vorticity solution
c *********************************
c
      write(6    ,*  ) ' ' 
      write(6    ,*  ) '     ---------------------------'
      write(6    ,*  ) '         CONVERGENCE HISTORY'
      nx        = mx  -1
      ny        = my  -1
      write(6    ,550)
      write(6    ,552) nx,ny,nmesh
      write(6    ,*  ) '     ---------------------------'
      write(6    ,553)
      write(iwrit,*  ) ' '
      write(iwrit,*  ) '     ---------------------------'
      write(iwrit,*  ) '         CONVERGENCE HISTORY'
      write(iwrit,550)
      write(iwrit,552) nx,ny,nmesh
      write(iwrit,*  ) '     ---------------------------'
      write(iwrit,553)
  550 format(5x,8h      nx,8h      ny,8h  n mesh) 
  552 format(5x,3i8)
  553 format(5h    n,12h    ratio dt,12h        time,
     .               12h          fx,12h        circ,
     .       5h  cyc,12h   max resid,4h   i,4h   j,
     .               12h    max corr,4h   i,4h   j)
c
      if (mstep .lt. 0) go to 998
      if (mstep .eq. 0) go to 997
c
      nstep     = nstep1
      time      = time1
c
    1 nstep     = nstep  +1
      time      = time  +dt
      timeio    = timeio  +dt
c
c set up boundary condition
c
      call bc(mx,my,x(kf),y(kf),g(kf))
c
c write mesh for AVS
c
      call io98 (mx,my,x(kf),y(kf))
c
c -----------------------------------------------------------------
c relaxation loop for streamfunction with multigrid (statement 100)
c -----------------------------------------------------------------
c
c compute inital residual
c
      mode       = 100
      call relxfv(mx,my,x(kf),y(kf),g(kf),p(kf),res(kf))
      write(6,508)     resmax,ires,jres
c
      ncyc      = 0
  100 ncyc      = ncyc  +1
c
c multigrid process - fine mesh to coarse mesh
c
      mode      = 0
      n         = 0
  110 n         = n  +1
      call relxfv(mx,my,x(kf),y(kf),g(kf),p(kf),res(kf))
      if (n .eq. nmesh) go to 120
      mode      = 1
      mmx       = (mx  -1)/2  +1
      mmy       = (my  -1)/2  +1
      kc        = kf  +mx*my
      call inject(mx,my,g(kf),mmx,mmy,g(kc))
      call inject(mx,my,res(kf),mmx,mmy,res(kc))
      kf        = kc
      mx        = mmx
      my        = mmy
      itmax     = itmax  +1
      go to 110
c
c multigrid process - coarse mesh to fine mesh
c
  120 n         = n  -1
      if (n .lt. 1) go to 130
      mode      = -1
      mmx       = mx
      mmy       = my
      kc        = kf
      mx        = 2*(mx  -1)  +1
      my        = 2*(my  -1)  +1
      kf        = kf  -mx*my
      call add(mx,my,g(kf),mmx,mmy,g(kc))
      itmax     = itmax  -1
      call relxfv(mx,my,x(kf),y(kf),g(kf),p(kf),res(kf))
      go to 120
c
c monitor convergence
c
  130 mode       = 100
      call relxfv(mx,my,x(kf),y(kf),g(kf),p(kf),res(kf))
      if (abs(resmax) .lt. 0.000001) go to 140
c
      if (ncyc .lt. mcyc) go to 100
c
c ----------------------------------
c solve vorticity transport equation
c ----------------------------------
c
  140 continue
      call   velo (mx,my,x(kf),y(kf),g(kf),w(kf))
      call   step (mx,my,x(kf),y(kf),dtinpt,dtstab)
      call bcvort (mx,my,x(kf),y(kf),g(kf),w(kf))
      call vortfv (mx,my,x(kf),y(kf),g(kf),w(kf),p(kf))
c
      call  track (mx,my,x(kf),y(kf),g(kf))
      call circul (mx,my,x(kf),y(kf),w(kf))
      call  force (mx,my,x(kf),y(kf),w(kf))
      cd        = -fx
      gamma     = -circ
      write (19,511) time,cd,gamma
c
      dtrat     = dtinpt/dtstab
      write(6,507)     nstep,dtrat,time,fx,circ,
     .                 ncyc,resmax,ires,jres,cormax,icor,jcor
      write(iwrit,507) nstep,dtrat,time,fx,circ,
     .                 ncyc,resmax,ires,jres,cormax,icor,jcor
c
c     if (timeio .gt. dtio) then
c     write (6,*) ' '
c     write (6,*) 'write data at time t =',time
c
c write flow data for AVS
c
      call io99 (mx,my,w(kf),time)
c     timeio     = 0.
c     else
c     endif
c
      if (nstep .lt. mstep) go to 1
c
c generate output
c
  997 continue
      write (6,*) ' '
      write (6,*) 'program saving data for restart calculation'
      call save(mx,my,x(kf),y(kf),g(kf),w(kf),nstep,time,dt)
      write (6,*) ' done'
c
      write (6,*) ' '
      write (6,*) 'program writing data in ft02.dat'
      call chart(mx,my,x(kf),y(kf),g(kf),w(kf),res(kf),igrid)
      write (6,*) ' done'
c
      write (6,*) ' '
      write (6,*) 'program writing data for MATLAB plot'
      call data(mx,my,x(kf),y(kf),w(kf))
      write (6,*) ' done'
c
      write (6,*) ' '
      write (6,*) 'program saving data for contour plot'
      call contour(mx,my,x(kf),y(kf),g(kf),w(kf),icont)
      write (6,*) ' done'
c
  501 format(1a60)
  502 format(1x,1a60)
  503 format(8f10.5)
  504 format(1x,8f10.0)
  505 format(8f10.3)
  506 format(1x,8f10.3)
  507 format(i5,4e12.3,i5,2(e12.3,2i4))
  508 format(58x,e12.3,2i4)
  509 format(6x,i6,2(e12.3,2i4))
  510 format(6x,i6,e12.3,2i4)
  511 format(3(2x,e12.5))
c
  998 stop
      end
















































































