c
c ---------------------------------------------------------------
c
      subroutine chart(mx,my,x,y,g,w,res,igrid)
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),res(mx,my)
      dimension  m(jnode)
      common/ai/ jalpha
      common/ar/ alpha(inode,jnode)
      common/f/  wnode(inode,jnode),wcell(inode,jnode)
      common/h/  vx(inode,jnode),vy(inode,jnode),
     .           dwx(inode,jnode),dwy(inode,jnode)
c
      iwrit     = 2
      jint      = 1  +(my  -1)/33
c
      if (igrid .eq. 1) then
      write(iwrit,500)
      write(iwrit,*) ' x-coordinate chart (x100) '
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,my,jint)
      write(iwrit,*) ' '
      do 1 i=1,mx                                                     
      do 2 j=1,my,jint                                                
    2 m(j)      = 100.*x(i,j)
    1 write (iwrit,501) i,(m(j),j=1,my,jint)                              
      write(iwrit,500)
      write(iwrit,*) ' y-coordinate chart (x100) '
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,my,jint)
      write(iwrit,*) ' '
      do 6 i=1,mx                                                     
      do 7 j=1,my,jint                                                
    7 m(j)      = 100.*y(i,j)
    6 write (iwrit,501) i,(m(j),j=1,my,jint)                              
      else
      endif
c                                                                       
      write(iwrit,500)
      write(iwrit,*) ' streamfunction g chart (x100) '                                 
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,my,jint)
      write(iwrit,*) ' '
      do 10 i=1,mx                                                     
      do 11 j=1,my,jint                                                
   11 m(j)      = 100.*g(i,j)
   10 write (iwrit,501) i,(m(j),j=1,my,jint)                              
c                                                                       
      write(iwrit,500)
      write(iwrit,*) ' vorticity w chart (x100) '
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,my,jint)
      write(iwrit,*) ' '
      wmax      = 0.
      do 20 i=1,mx
      do 20 j=1,my
      if (abs(w(i,j)) .lt. wmax) go to 20
      wmax      = abs(w(i,j))
   20 continue
      xnorm     = 1.
      if (wmax .ne. 0.) xnorm = wmax
      do 22 i=1,mx                                                     
      do 23 j=1,my
   23 m(j)      = 100.*w(i,j)/xnorm
   22 write (iwrit,501) i,(m(j),j=1,my,jint)                              
      write(iwrit,*) ' '
      write(iwrit,*) ' max[w] = ',wmax
c
      write(iwrit,500)
      write(iwrit,*) ' u-velocity chart (x100) '
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,my,jint)
      write(iwrit,*) ' '
      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 32 j=1,my
   32 m(j)       = 100.*wnode(i,j)
   31 write (iwrit,501) i,(m(j),j=1,my,jint)                              
c                                                                       
      write(iwrit,500)
      write(iwrit,*) ' v-velocity chart (x100) '
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,my,jint)
      write(iwrit,*) ' '
      do 40 i=1,mx-1
      do 40 j=1,my-1
   40 wcell(i,j) = vy(i,j)
      iopt       = 2
      call average(mx,my,iopt)
      do 41 i=1,mx                                                   
      do 42 j=1,my
   42 m(j)       = 100.*wnode(i,j)
   41 write (iwrit,501) i,(m(j),j=1,my,jint)                              
c                                                                       
      write(iwrit,500)
      write(iwrit,*) ' Lagrangian coordinate ALPHA chart (x10) '                                 
      write(iwrit,*) ' '
      write(iwrit,502) (j,j=1,jalpha,jint)
      write(iwrit,*) ' '
      do 50 i=1,mx                                                     
      do 51 j=1,jalpha,jint                                                
   51 m(j)      = 10.*alpha(i,j)
   50 write (iwrit,501) i,(m(j),j=1,jalpha,jint)                              
c
  500 format(1h1)
  501 format(1x,i2,2x,33i4)  
  502 format(5x,33i4)  
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine circul (mx,my,x,y,w)
c
c compute circulation
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),w(mx,my)
      common/c/  press(inode),circ,fx
c
      nx        = mx  -1
      ny        = my  -1
c
      circ      = 0.
c
      do 10 j=1,ny/2
      do 10 i=1,nx
      area      = .5*abs((x(i+1,j+1)  -x(i,j))
     .                  *(y(i+1,j)  -y(i,j+1))
     .                  -(y(i+1,j+1)  -y(i,j))
     .                  *(x(i+1,j)  -x(i,j+1)))
      vort      = .25*(w(i,j)  +w(i-1,j)  +w(i,j-1)  +w(i-1,j-1))
   10 circ      = circ  +vort*area
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine force (mx,my,x,y,w)
c
c compute force acting on plate
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),w(mx,my)
      common/c/  press(inode),circ,fx
      common/j/  re,dt,time,tstop
c
      imid         = (mx  -1)/2  +1
      press(imid)  = 0.
      j            = 1
c
c pressure acting on left side of object
c
      i            = imid
      dx           = x(i,j+1)  -x(i,j)
      dy           = y(i,j+1)  -y(i,j)
      dn           = sqrt(dx*dx  +dy*dy)
      dwdn1        = (w(i,j+1)  -w(i,j))/dn
      do 10 i=imid-1,1,-1
      dx           = x(i,j+1)  -x(i,j)
      dy           = y(i,j+1)  -y(i,j)
      dn           = sqrt(dx*dx  +dy*dy)
      dwdn2        = (w(i,j+1)  -w(i,j))/dn
      dwdn         = .5*(dwdn1  +dwdn2)
      dx           = x(i+1,j)  -x(i,j)
      dy           = y(i+1,j)  -y(i,j)
      ds           = sqrt(dx*dx  +dy*dy)
      press(i)     = press(i+1)  +dwdn*ds/re
   10 dwdn1        = dwdn2
c
c pressure acting on right side of object
c
      i            = imid
      dx           = x(i,j+1)  -x(i,j)
      dy           = y(i,j+1)  -y(i,j)
      dn           = sqrt(dx*dx  +dy*dy)
      dwdn1        = (w(i,j+1)  -w(i,j))/dn
      do 20 i=imid,mx-1
      dx           = x(i,j+1)  -x(i,j)
      dy           = y(i,j+1)  -y(i,j)
      dn           = sqrt(dx*dx  +dy*dy)
      dwdn2        = (w(i,j+1)  -w(i,j))/dn
      dwdn         = .5*(dwdn1  +dwdn2)
      dx           = x(i+1,j)  -x(i,j)
      dy           = y(i+1,j)  -y(i,j)
      ds           = sqrt(dx*dx  +dy*dy)
      press(i+1)   = press(i)  -dwdn*ds/re
   20 dwdn1        = dwdn2
c
c net drag force acting on plate
c
      fx           = 0.
      do 30 i=1,mx-1
      dx           = x(i+1,j)  -x(i,j)
      dy           = y(i+1,j)  -y(i,j)
      ds           = sqrt(dx*dx  +dy*dy)
      slop         = dy/dx
      xnx          = -slop/sqrt(1.  +slop*slop)
      pavg         = .5*(press(i+1)  +press(i))
   30 fx           = fx  -pavg*xnx*ds
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine data(mx,my,x,y,w)
c
c write data for MATLAB 
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),w(mx,my)
      common/c/  press(inode),circ,fx
      common/h/  vx(inode,jnode),vy(inode,jnode),
     .           dwx(inode,jnode),dwy(inode,jnode)
c
      pi        = 4.*atan(1.)
      imid      = (mx  -1)/2  +1
c
c write data in ft03.dat
c
      iwrit     = 3
      j         = 1
      s         = 0.
      i         = 1
      write (iwrit,502) s,w(i,j),press(i)
      do 10 i=2,mx                                                 
      dx        = x(i,j)  -x(i-1,j)
      dy        = y(i,j)  -y(i-1,j)
      ds        = sqrt(dx*dx  +dy*dy)
      s         = s  +ds
   10 write (iwrit,502) s,w(i,j),press(i)
c
c write data in ft04.dat
c
      iwrit     = 4
      i         = mx
      uavg      = 0.
      write(iwrit,502) x(i,1),uavg
      do 30 j=2,my-1
      uavg      = .5*(vx(i-1,j)  +vx(i-1,j-1))
   30 write(iwrit,502) x(i,j),uavg
c
  502 format(1x,3e15.5)
c
      return
      end
c
c ---------------------------------------------------------------
c
      subroutine contour(mx,my,x,y,g,w,icont)
c
c write data for MATLAB 
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)
      dimension  work(2*inode,2*jnode)
      common/ai/ jalpha
      common/ar/ alpha(inode,jnode)
      common/f/  wnode(inode,jnode),wcell(inode,jnode)
      common/h/  vx(inode,jnode),vy(inode,jnode),
     .           dwx(inode,jnode),dwy(inode,jnode)
c
      if (icont .lt. 0) return
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
