c
c ---------------------------------------------------------------
c
      subroutine mesh(mx,my,x,y)
c
c generate mesh
c
      implicit   real*8 (a-h,o-z)
      dimension  x(mx,my),y(mx,my)
c
c read mesh from ft14.dat
c
      iunit     = 14
      read (iunit,*) mx,my
      do 10 i=1,mx
      do 10 j=1,my
   10 read (iunit,*) x(i,j),y(i,j)
c
c search for smallest cell
c
      amin      = 1.0E+99
      imin      = -1
      jmin      = -1
      do 15 i=1,mx-1
      do 15 j=1,my-1
      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)
      area      = xx*yy  -xy*yx
      if (area .gt. amin) go to 15
      amin      = area
      imin      = i
      jmin      = j
   15 continue
      dxmin     = sqrt(amin)
      write (6,*) ' '
      write (6,*) 'minimum mesh size dx = ',dxmin
      write (6,*) '                at i = ',imin
      write (6,*) '                   j = ',jmin
c
c save mesh for MATLAB plotting
c
      iunitx    = 12
      iunity    = 13
c
      jc        = 0
      do 60 j=1,my
      if (jc .eq. 0) then
      i1        = 1
      i2        = mx
      i3        = 1
      jc        = 1
      else
      i1        = mx
      i2        = 1
      i3        = -1
      jc        = 0
      endif
      do 60 i=i1,i2,i3
      write(iunitx,100) x(i,j)
   60 write(iunity,100) y(i,j)
c
      ic        = 0
      do 70 i=mx,1,-1
      if (ic .eq. 0) then
      j1        = 1
      j2        = my
      j3        = 1
      ic        = 1
      else
      j1        = my
      j2        = 1
      j3        = -1
      ic        = 0
      endif
      do 70 j=j1,j2,j3
      write(iunitx,100) x(i,j)
   70 write(iunity,100) y(i,j)
c
  100 format(2(2x,f12.7))
  101 format(2x,2i5)
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine ic(mx,my,x,y,g,w,p)
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),g(mx,my),w(mx,my),p(mx,my)
      common/ai/ jalpha
      common/ar/ alpha(inode,jnode)
c
      do 10 j=1,my
      do 10 i=1,mx
      g(i,j)     = 0.
      w(i,j)     = 0.
      p(i,j)     = 0.
   10 alpha(i,j) = 0.
c
      do 20 i=1,mx
      do 22 j=1,jalpha
   22 alpha(i,j) = float(j  -1)
      do 24 j=jalpha+1,my
   24 alpha(i,j) = alpha(i,jalpha)
   20 continue
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine bc(mx,my,x,y,g)
c
c set up boundary condition
c
      implicit   real*8 (a-h,o-z)
      dimension  x(mx,my),y(mx,my),g(mx,my)
      common/j/  re,dt,time,tstop
c
      pi        = 4.*atan(1.)
c
      if (time .eq. 0.) then
      write (6,*) ' '
      write (6,*) ' Flowfield impulsively started at farfield'
c
c along object
c
      do 10 i=1,mx
   10 g(i,1)    = 0.
c
c along symmetry lines
c
      do 20 j=1,my
      g(1,j)    = 0.
   20 g(mx,j)   = 0.
c
c along farfield boundary
c
      j         = my
      uinf      = 1.
      do 30 i=1,mx
      radius    = sqrt(x(i,j)**2  +y(i,j)**2)
      sine      = y(i,j)/radius
   30 g(i,j)    = uinf*radius*sine
      else
      endif
c
      if (time .ge. tstop) then
c     write (6,*) ' '
c     write (6,*) ' Flowfield impulsively stopped at farfield'
c
c along object
c
      do 60 i=1,mx
   60 g(i,1)    = 0.
c
c along symmetry lines
c
      do 70 j=1,my
      g(1,j)    = 0.
   70 g(mx,j)   = 0.
c
c along farfield boundary
c
      j         = my
      uinf      = 1.
      do 80 i=1,mx
   80 g(i,j)    = 0.
      else
      endif
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine io98 (mx,my,x,y)
c
c write output to unit 99
c mesh information for AVS
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my)
c
      iwrit     = 98
c
      write (iwrit) mx,my
      write (iwrit) ((x(i,j),i=1,mx),j=1,my)
      write (iwrit) ((y(i,j),i=1,mx),j=1,my)
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine io99 (mx,my,w,time)
c
c write output to unit 99
c write flow data for AVS
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  w(mx,my)
c
      iwrit     = 99
c
      write (iwrit) time
      write (iwrit) ((w(i,j),i=1,mx),j=1,my)
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine save (mx,my,x,y,g,w,nstep,time,dt)
c
c write output to unit 2
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),g(mx,my),w(mx,my)
      common/ar/ alpha(inode,jnode)
c
      iwrit     = 21
c
      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
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine rstart (mx,my,x,y,g,w,p,nstep1,time1,dt)
c
c write output to unit 2
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),g(mx,my),w(mx,my),p(mx,my)
      common/ar/ alpha(inode,jnode)
c
      iread     = 22
c
      read (iread) mx,my,nstep1
      read (iread) time1,dt
      read (iread) ((x(i,j),i=1,mx),j=1,my)
      read (iread) ((y(i,j),i=1,mx),j=1,my)
      read (iread) ((g(i,j),i=1,mx),j=1,my)
      read (iread) ((w(i,j),i=1,mx),j=1,my)
      read (iread) ((alpha(i,j),i=1,mx),j=1,my)
c
      do 10 i=1,mx
      do 10 j=1,my
   10 p(i,j)    = -w(i,j)
c
      return
      end

























































































