c
      implicit real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      common/a/  x(inode,jnode),y(inode,jnode)
      common/b/  g(inode,jnode),w(inode,jnode)
      common/ci/ jalpha
      common/c/  alpha(inode,jnode)
c
      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(98,file='ft98.dat',form='unformatted')
      open(99,file='ft99.dat',form='unformatted')
c
c read mesh from file 14
c
      iread     = 14
      read (iread,*) mx,my
      do 1 i=1,mx
      do 1 j=1,my
    1 read (iread,*) x(i,j),y(i,j)
c
c read data from file 99
c
      iread     = 99
c
c     read (iread) mx,my
c     read (iread) ((x(i,j),i=1,mx),j=1,my)
c     read (iread) ((y(i,j),i=1,mx),j=1,my)
c
      write (6,*) ' '
      write (6,*) 'Enter desired time number'
      read  (5,*) nmax
c
      do 10 n=1,nmax
      read (iread) nstep
      read (iread) time
      write (6,*) ' '
      write (6,*) 'time t =',time
      read (iread) ((g(i,j),i=1,mx),j=1,my)
      read (iread) ((w(i,j),i=1,mx),j=1,my)
   10 read (iread) ((alpha(i,j),i=1,mx),j=1,my)
c
      write (6,*) ' '
      write (6,*) 'Enter jalpha'
      read  (5,*) jalpha
      call velo(mx,my)
      call contour(mx,my)
c
      iwrit     = 98
      write (iwrit) mx,my,nstep
      write (iwrit) time,dt
      write (iwrit) ((x(i,j),i=1,mx),j=1,my)
      write (iwrit) ((y(i,j),i=1,mx),j=1,my)
      write (iwrit) ((g(i,j),i=1,mx),j=1,my)
      write (iwrit) ((w(i,j),i=1,mx),j=1,my)
      write (iwrit) ((alpha(i,j),i=1,mx),j=1,my)
c
      stop
      end
c
c ---------------------------------------------------------------
c
      subroutine contour(mx,my)
c
c write data for MATLAB 
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  work(2*inode,2*jnode)
      common/a/  x(inode,jnode),y(inode,jnode)
      common/b/  g(inode,jnode),w(inode,jnode)
      common/ci/ jalpha
      common/c/  alpha(inode,jnode)
      common/f/  wnode(inode,jnode),wcell(inode,jnode)
      common/h/  vx(inode,jnode),vy(inode,jnode)
c
      nx        = mx  -1
      ny        = my  -1
c
      istart    = 1
      iend      = 2*nx  +1
      jstart    = 1
      jend      = jalpha
c
c write data for NCAR
c
      mpx       = iend  -istart  +1
      mpy       = jend  -jstart  +1
c
c write streamfunction
c
      icont     = 15
      write (icont,*) mpx,mpy
      ii        = 2*mx
      do 10 i=1,mx
      ii        = ii  -1
      do 10 j=1,my
      work(i,j) = x(i,j)
   10 work(ii,j)= x(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
      ii        = 2*mx
      do 12 i=1,mx
      ii        = ii  -1
      do 12 j=1,my
      work(i,j) =  y(i,j)
   12 work(ii,j)= -y(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
      ii        = 2*mx
      do 14 i=1,mx
      ii        = ii  -1
      do 14 j=1,my
      work(i,j) =  g(i,j)
   14 work(ii,j)= -g(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
c
c write vorticity
c
      icont     = 16
      write (icont,*) mpx,mpy
      ii        = 2*mx
      do 20 i=1,mx
      ii        = ii  -1
      do 20 j=1,my
      work(i,j) = x(i,j)
   20 work(ii,j)= x(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
      ii        = 2*mx
      do 22 i=1,mx
      ii        = ii  -1
      do 22 j=1,my
      work(i,j) =  y(i,j)
   22 work(ii,j)= -y(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
      ii        = 2*mx
      do 24 i=1,mx
      ii        = ii  -1
      do 24 j=1,my
      work(i,j) =  w(i,j)
   24 work(ii,j)= -w(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
c
c write velocity
c
      icont     = 17
      write (icont,*) mpx,mpy
      ii        = 2*mx
      do 40 i=1,mx
      ii        = ii  -1
      do 40 j=1,my
      work(i,j) = x(i,j)
   40 work(ii,j)= x(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
c
      ii        = 2*mx
      do 42 i=1,mx
      ii        = ii  -1
      do 42 j=1,my
      work(i,j) =  y(i,j)
   42 work(ii,j)= -y(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
c
      do 30 i=1,mx-1
      do 30 j=1,my-1
   30 wcell(i,j) = vx(i,j)
      iopt       = 1
      call average(mx,my,iopt)
      do 31 i=1,mx
      do 31 j=1,my
   31 work(i,j)  = wnode(i,j)**2
      do 32 i=1,mx-1
      do 32 j=1,my-1
   32 wcell(i,j) = vy(i,j)
      iopt       = 2
      call average(mx,my,iopt)
      do 33 i=1,mx
      do 33 j=1,my
   33 wcell(i,j) = sqrt(work(i,j)  +wnode(i,j)**2)
c
      ii        = 2*mx
      do 44 i=1,mx
      ii        = ii  -1
      do 44 j=1,my
      work(i,j) =  wcell(i,j)
   44 work(ii,j)=  wcell(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
c
c write Lagrangian coordinate
c
      icont     = 18
      write (icont,*) mpx,mpy
      ii        = 2*mx
      do 50 i=1,mx
      ii        = ii  -1
      do 50 j=1,my
      work(i,j) = x(i,j)
   50 work(ii,j)= x(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
      ii        = 2*mx
      do 52 i=1,mx
      ii        = ii  -1
      do 52 j=1,my
      work(i,j) =  y(i,j)
   52 work(ii,j)= -y(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
      ii        = 2*mx
      do 54 i=1,mx
      ii        = ii  -1
      do 54 j=1,my
      work(i,j) =  alpha(i,j)
   54 work(ii,j)=  alpha(i,j)
      write (icont,*) ((work(i,j),j=jstart,jend),i=istart,iend)
c
  100 format(65e12.4)
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine velo (mx,my)
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      common/a/  x(inode,jnode),y(inode,jnode)
      common/b/  g(inode,jnode),w(inode,jnode)
      common/h/  vx(inode,jnode),vy(inode,jnode)
c
      nx        = mx  -1
      ny        = my  -1
c
      do 10 i=1,nx
      do 10 j=1,ny
      xx        = x(i+1,j+1)  -x(i,j+1)  +x(i+1,j)  -x(i,j)
      xy        = x(i+1,j+1)  +x(i,j+1)  -x(i+1,j)  -x(i,j)
      yx        = y(i+1,j+1)  -y(i,j+1)  +y(i+1,j)  -y(i,j)
      yy        = y(i+1,j+1)  +y(i,j+1)  -y(i+1,j)  -y(i,j)
      det       = 1./(xx*yy  -xy*yx)
      gx        = g(i+1,j+1)  -g(i,j+1)  +g(i+1,j)  -g(i,j)
      gy        = g(i+1,j+1)  +g(i,j+1)  -g(i+1,j)  -g(i,j)
      vy(i,j)   = -(yy*gx  -yx*gy)*det
   10 vx(i,j)   =  (xx*gy  -xy*gx)*det
c
      return
      end 
c
c ----------------------------------------------------------------
c
      subroutine average(mx,my,iopt)
c
c compute nodal value using cell-averaged value
c iopt = 1 ---> x-velocity
c iopt = 2 ---> y-velocity
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      common/f/  wnode(inode,jnode),wcell(inode,jnode)
c
      nx            = mx  -1
      ny            = my  -1
c
c inner points
c
      do 1 i=2,nx
      do 1 j=2,ny
    1 wnode(i,j)    = .25*(wcell(i,j)  +wcell(i,j-1)
     .                    +wcell(i-1,j)  +wcell(i-1,j-1))
c
c along objects
c
      j         = 1
      do 2 i=2,nx
    2 wnode(i,j)    = 0.
      wnode(1,j)    = 0.
      wnode(mx,j)   = 0.
c
c along farfield
c
      j         = my
      do 3 i=2,nx
    3 wnode(i,j)    = .5*(wcell(i,j-1)  +wcell(i-1,j-1))
      if (iopt .eq. 1) then
      wnode(1,j)    = wnode(2,j)
      wnode(mx,j)   = wnode(mx-1,j)
      else
      wnode(1,j)    = 0.
      wnode(mx,j)   = 0.
      endif
c
c along symmetric plane
c
      if (iopt .eq. 1) then
      do 4 j=2,ny
      i             = 1
      wnode(i,j)    = .5*(wcell(i,j)  +wcell(i,j-1))
      i             = mx
    4 wnode(i,j)    = .5*(wcell(i-1,j)  +wcell(i-1,j-1))
      else
      do 5 j=2,ny
      i             = 1
      wnode(i,j)    = 0.
      i             = mx
    5 wnode(i,j)    = 0.
      endif
c
      return
      end
