c
c ---------------------------------------------------------------
c
      subroutine track (mx,my,x,y,g)
c
c track material lines using D(alpha)/Dt = 0
c finite-volume method
c McCormack scheme
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),g(mx,my)
      common/ai/ jalpha
      common/ar/ alpha(inode,jnode)
      common/e/  resid(jnode),xc(2,jnode),yc(2,jnode),
     .           gxw(jnode),gyw(jnode),gxe(jnode),gye(jnode),
     .           gxns(jnode),gyns(jnode)
      common/f/  wp(inode,jnode),wc(inode,jnode)
      common/h/  vx(inode,jnode),vy(inode,jnode),
     .           dwx(inode,jnode),dwy(inode,jnode)
      common/j/  re,dt,time,tstop
c
      nx        = mx  -1
      imid      = nx/2  +1
      mj        = jalpha
      nj        = mj  -1
c
      j         = 1
      do 1 i=1,mx
      wp(i,j)   = alpha(i,j)
    1 wc(i,j)   = alpha(i,j)
c
c **************
c predictor step
c **************
c
c update alpha inside domain
c
      i         = 1
      do 100 j=1,mj
      xc(1,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
  100 yc(1,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
c
      do 110 i=2,nx
c
      do 120 j=1,mj
      xc(2,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
  120 yc(2,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
c
      do 130 j=2,nj
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j+1)  -xc(1,j))
     .                  *(yc(2,j)  -yc(1,j+1))
     .                  -(yc(2,j+1)  -yc(1,j))
     .                  *(xc(2,j)  -xc(1,j+1)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      fluxf     = -.5*alpha(i-1,j)*(vx(i-1,j)  +vx(i-1,j-1))
      fluxg     = -.5*alpha(i-1,j)*(vy(i-1,j)  +vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c east face
      dx        = xc(2,j-1)  -xc(2,j)
      dy        = yc(2,j-1)  -yc(2,j)
      fluxf     = -.5*alpha(i,j)*(vx(i,j)  +vx(i,j-1))
      fluxg     = -.5*alpha(i,j)*(vy(i,j)  +vy(i,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c south face
      dx        = xc(1,j-1)  -xc(2,j-1)
      dy        = yc(1,j-1)  -yc(2,j-1)
      fluxf     = -.5*alpha(i,j-1)*(vx(i,j-1)  +vx(i-1,j-1))
      fluxg     = -.5*alpha(i,j-1)*(vy(i,j-1)  +vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      fluxf     = -.5*alpha(i,j)*(vx(i,j)  +vx(i-1,j))
      fluxg     = -.5*alpha(i,j)*(vy(i,j)  +vy(i-1,j))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
      rhs       = flux/area
  130 wp(i,j)   = alpha(i,j)  +dt*rhs
c
      do 150 j=1,mj
      xc(1,j)   = xc(2,j)
  150 yc(1,j)   = yc(2,j)
c
  110 continue
c
c update alpha along symmetry plane upstream of object
c
      i         = 1
c
      do 200 j=1,mj
      xc(2,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
      xc(1,j)   = xc(2,j)
      yc(2,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
  200 yc(1,j)   = -yc(2,j)
c
      do 230 j=2,nj
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j+1)  -xc(1,j))
     .                  *(yc(2,j)  -yc(1,j+1))
     .                  -(yc(2,j+1)  -yc(1,j))
     .                  *(xc(2,j)  -xc(1,j+1)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      fluxf     = -.5*alpha(i+1,j)*(vx(i,j)  +vx(i,j-1))
      fluxg     = -.5*alpha(i+1,j)*(-vy(i,j)  -vy(i,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c east face
      dx        = xc(2,j-1)  -xc(2,j)
      dy        = yc(2,j-1)  -yc(2,j)
      fluxf     = -.5*alpha(i,j)*(vx(i,j)  +vx(i,j-1))
      fluxg     = -.5*alpha(i,j)*(vy(i,j)  +vy(i,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c south face
      dx        = xc(1,j-1)  -xc(2,j-1)
      dy        = yc(1,j-1)  -yc(2,j-1)
      fluxf     = -alpha(i,j-1)*vx(i,j-1)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      fluxf     = -alpha(i,j)*vx(i,j)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
      rhs       = flux/area
  230 wp(i,j)   = alpha(i,j)  +dt*rhs
c
c update alpha along symmetry plane downstream of object
c
      i         = mx
c
      do 300 j=1,mj
      xc(1,j)   = .25*(x(i,j+1)  +x(i-1,j+1)  +x(i,j)  +x(i-1,j))
      xc(2,j)   = xc(1,j)
      yc(1,j)   = .25*(y(i,j+1)  +y(i-1,j+1)  +y(i,j)  +y(i-1,j))
  300 yc(2,j)   = -yc(1,j)
c
      do 330 j=2,nj
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j+1)  -xc(1,j))
     .                  *(yc(2,j)  -yc(1,j+1))
     .                  -(yc(2,j+1)  -yc(1,j))
     .                  *(xc(2,j)  -xc(1,j+1)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      fluxf     = -.5*alpha(i-1,j)*(vx(i-1,j)  +vx(i-1,j-1))
      fluxg     = -.5*alpha(i-1,j)*(vy(i-1,j)  +vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c east face
      dx        = xc(2,j-1)  -xc(2,j)
      dy        = yc(2,j-1)  -yc(2,j)
      fluxf     = -.5*alpha(i,j)*(vx(i-1,j)  +vx(i-1,j-1))
      fluxg     = -.5*alpha(i,j)*(-vy(i-1,j)  -vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c south face
      dx        = xc(1,j-1)  -xc(2,j-1)
      dy        = yc(1,j-1)  -yc(2,j-1)
      fluxf     = -alpha(i,j-1)*vx(i-1,j-1)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      fluxf     = -alpha(i,j)*vx(i-1,j)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
      rhs       = flux/area
  330 wp(i,j)   = alpha(i,j)  +dt*rhs
c
c extrapolation at farfield location 
c
      j         = mj
      do 400 i=1,mx
  400 wp(i,j)   = (4.*wp(i,j-1)  -wp(i,j-2))/3.
c
c **************
c corrector step
c **************
c
c update alpha inside domain
c
      i         = 1
      do 600 j=1,mj
      xc(1,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
  600 yc(1,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
c
      do 610 i=2,nx
c
      do 620 j=1,mj
      xc(2,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
  620 yc(2,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
c
      do 630 j=2,nj
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j+1)  -xc(1,j))
     .                  *(yc(2,j)  -yc(1,j+1))
     .                  -(yc(2,j+1)  -yc(1,j))
     .                  *(xc(2,j)  -xc(1,j+1)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      fluxf     = -.5*wp(i,j)*(vx(i-1,j)  +vx(i-1,j-1))
      fluxg     = -.5*wp(i,j)*(vy(i-1,j)  +vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c east face
      dx        = xc(2,j-1)  -xc(2,j)
      dy        = yc(2,j-1)  -yc(2,j)
      fluxf     = -.5*wp(i+1,j)*(vx(i,j)  +vx(i,j-1))
      fluxg     = -.5*wp(i+1,j)*(vy(i,j)  +vy(i,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c south face
      dx        = xc(1,j-1)  -xc(2,j-1)
      dy        = yc(1,j-1)  -yc(2,j-1)
      fluxf     = -.5*wp(i,j)*(vx(i,j-1)  +vx(i-1,j-1))
      fluxg     = -.5*wp(i,j)*(vy(i,j-1)  +vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      fluxf     = -.5*wp(i,j+1)*(vx(i,j)  +vx(i-1,j))
      fluxg     = -.5*wp(i,j+1)*(vy(i,j)  +vy(i-1,j))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
      rhs       = flux/area
  630 wc(i,j)   = .5*(alpha(i,j)  +wp(i,j)  +dt*rhs)
c
      do 650 j=1,mj
      xc(1,j)   = xc(2,j)
  650 yc(1,j)   = yc(2,j)
c
  610 continue
c
c update alpha along symmetry plane upstream of object
c
      i         = 1
c
      do 700 j=1,mj
      xc(2,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
      xc(1,j)   = xc(2,j)
      yc(2,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
  700 yc(1,j)   = -yc(2,j)
c
      do 730 j=2,nj
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j+1)  -xc(1,j))
     .                  *(yc(2,j)  -yc(1,j+1))
     .                  -(yc(2,j+1)  -yc(1,j))
     .                  *(xc(2,j)  -xc(1,j+1)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      fluxf     = -.5*wp(i,j)*(vx(i,j)  +vx(i,j-1))
      fluxg     = -.5*wp(i,j)*(-vy(i,j)  -vy(i,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c east face
      dx        = xc(2,j-1)  -xc(2,j)
      dy        = yc(2,j-1)  -yc(2,j)
      fluxf     = -.5*wp(i+1,j)*(vx(i,j)  +vx(i,j-1))
      fluxg     = -.5*wp(i+1,j)*(vy(i,j)  +vy(i,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c south face
      dx        = xc(1,j-1)  -xc(2,j-1)
      dy        = yc(1,j-1)  -yc(2,j-1)
      fluxf     = -wp(i,j)*vx(i,j-1)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      fluxf     = -wp(i,j+1)*vx(i,j)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
      rhs       = flux/area
  730 wc(i,j)   = .5*(alpha(i,j)  +wp(i,j)  +dt*rhs)
c
c update alpha along symmetry plane downstream of object
c
      i         = mx
c
      do 800 j=1,mj
      xc(1,j)   = .25*(x(i,j+1)  +x(i-1,j+1)  +x(i,j)  +x(i-1,j))
      xc(2,j)   = xc(1,j)
      yc(1,j)   = .25*(y(i,j+1)  +y(i-1,j+1)  +y(i,j)  +y(i-1,j))
  800 yc(2,j)   = -yc(1,j)
c
      do 830 j=2,nj
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j+1)  -xc(1,j))
     .                  *(yc(2,j)  -yc(1,j+1))
     .                  -(yc(2,j+1)  -yc(1,j))
     .                  *(xc(2,j)  -xc(1,j+1)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      fluxf     = -.5*wp(i,j)*(vx(i-1,j)  +vx(i-1,j-1))
      fluxg     = -.5*wp(i,j)*(vy(i-1,j)  +vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c east face
      dx        = xc(2,j-1)  -xc(2,j)
      dy        = yc(2,j-1)  -yc(2,j)
      fluxf     = -.5*wp(i-1,j)*(vx(i-1,j)  +vx(i-1,j-1))
      fluxg     = -.5*wp(i-1,j)*(-vy(i-1,j)  -vy(i-1,j-1))
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c south face
      dx        = xc(1,j-1)  -xc(2,j-1)
      dy        = yc(1,j-1)  -yc(2,j-1)
      fluxf     = -wp(i,j)*vx(i-1,j-1)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      fluxf     = -wp(i,j+1)*vx(i-1,j)
      fluxg     = 0.
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
      rhs       = flux/area
  830 wc(i,j)   = .5*(alpha(i,j)  +wp(i,j)  +dt*rhs)
c
c extrapolation at farfield location 
c
      j         = mj
      do 900 i=1,mx
  900 wc(i,j)   = (4.*wc(i,j-1)  -wc(i,j-2))/3.
c
c ************
c update alpha 
c ************
c
      do 1000 i=1,mx
      do 1000 j=1,mj
 1000 alpha(i,j) = wc(i,j)
c
      return
      end
