c
c ---------------------------------------------------------------
c
      subroutine vortfv (mx,my,x,y,g,w,p)
c
c solve vorticity transport equation
c finite-volume method
c First-order Upwind scheme
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/bi/ mode,itmax,ires,jres,ierr,jerr,icor,jcor
      common/br/ resmax,errmax,cormax
      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
      ierr      = -1
      jerr      = -1
      errmax    = 0.
c
      nx        = mx  -1
      ny        = my  -1
c
      i         = 1
      do 10 j=1,ny
      xc(1,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
   10 yc(1,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
c
      do 20 i=2,nx
c
      do 30 j=1,ny
      xc(2,j)   = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
   30 yc(2,j)   = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
c
      do 40 j=2,ny
      flux      = 0.
c cell area
      area      = .5*abs((xc(2,j)  -xc(1,j-1))
     .                  *(yc(2,j-1)  -yc(1,j))
     .                  -(yc(2,j)  -yc(1,j-1))
     .                  *(xc(2,j-1)  -xc(1,j)))
c west face
      dx        = xc(1,j)  -xc(1,j-1)
      dy        = yc(1,j)  -yc(1,j-1)
      vxavg     = .5*(vx(i-1,j)  +vx(i-1,j-1))
      vxpos     = .5*(vxavg  +abs(vxavg))
      vxneg     = .5*(vxavg  -abs(vxavg))
      fluxf     = .5*(dwx(i-1,j)  +dwx(i-1,j-1))/re
     .            -w(i-1,j)*vxpos  -w(i,j)*vxneg
      vyavg     = .5*(vy(i-1,j)  +vy(i-1,j-1))
      vypos     = .5*(vyavg  +abs(vyavg))
      vyneg     = .5*(vyavg  -abs(vyavg))
      fluxg     = .5*(dwy(i-1,j)  +dwy(i-1,j-1))/re
     .            -w(i-1,j)*vypos  -w(i,j)*vyneg
      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)
      vxavg     = .5*(vx(i,j)  +vx(i,j-1))
      vxpos     = .5*(vxavg  +abs(vxavg))
      vxneg     = .5*(vxavg  -abs(vxavg))
      fluxf     = .5*(dwx(i,j)  +dwx(i,j-1))/re
     .            -w(i,j)*vxpos  -w(i+1,j)*vxneg
      vyavg     = .5*(vy(i,j)  +vy(i,j-1))
      vypos     = .5*(vyavg  +abs(vyavg))
      vyneg     = .5*(vyavg  -abs(vyavg))
      fluxg     = .5*(dwy(i,j)  +dwy(i,j-1))/re
     .            -w(i,j)*vypos  -w(i+1,j)*vyneg
      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)
      vxavg     = .5*(vx(i,j-1)  +vx(i-1,j-1))
      vxpos     = .5*(vxavg  +abs(vxavg))
      vxneg     = .5*(vxavg  -abs(vxavg))
      fluxf     = .5*(dwx(i,j-1)  +dwx(i-1,j-1))/re
     .            -w(i,j-1)*vxpos  -w(i,j)*vxneg
      vyavg     = .5*(vy(i,j-1)  +vy(i-1,j-1))
      vypos     = .5*(vyavg  +abs(vyavg))
      vyneg     = .5*(vyavg  -abs(vyavg))
      fluxg     = .5*(dwy(i,j-1)  +dwy(i-1,j-1))/re
     .            -w(i,j-1)*vypos  -w(i,j)*vyneg
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c north face
      dx        = xc(2,j)  -xc(1,j)
      dy        = yc(2,j)  -yc(1,j)
      vxavg     = .5*(vx(i,j)  +vx(i-1,j))
      vxpos     = .5*(vxavg  +abs(vxavg))
      vxneg     = .5*(vxavg  -abs(vxavg))
      fluxf     = .5*(dwx(i,j)  +dwx(i-1,j))/re
     .            -w(i,j)*vxpos  -w(i,j+1)*vxneg
      vyavg     = .5*(vy(i,j)  +vy(i-1,j))
      vypos     = .5*(vyavg  +abs(vyavg))
      vyneg     = .5*(vyavg  -abs(vyavg))
      fluxg     = .5*(dwy(i,j)  +dwy(i-1,j))/re
     .            -w(i,j)*vypos  -w(i,j+1)*vyneg
      flux      = flux  +(-fluxf*dy  +fluxg*dx)
c
      rhs       = flux/area
   40 wc(i,j)   = w(i,j)  +dt*rhs
c
      do 50 j=1,ny
      xc(1,j)   = xc(2,j)
   50 yc(1,j)   = yc(2,j)
c
   20 continue
c
c     do 200 j=2,my-1
c     do 200 i=2,mx-1
c     error     = wc(i,j)  -w(i,j)
c     if (abs(error) .lt. abs(errmax)) go to 200
c     ierr      = i
c     jerr      = j
c     errmax    = error
c 200 continue
c
      do 300 j=2,my-1
      do 300 i=2,mx-1
      w(i,j)    = wc(i,j)
  300 p(i,j)    = -w(i,j)
c
      return
      end 
c
c ---------------------------------------------------------------
c
      subroutine bcvort (mx,my,x,y,g,w)
c
c set up boundary condition for vorticity
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)
c
c at farfield boundary
c
      j         = my
      do 10 i=1,mx
   10 w(i,j)    = 0.
c
c along symmetry plane
c
      do 20 j=1,my
      w(1,j)    = 0.
   20 w(mx,j)   = 0.
c
c on object
c
      j         = 1
      do 30 i=2,mx-1
      dx        = x(i,j+1)  -x(i,j)
      dy        = y(i,j+1)  -y(i,j)
      dn2       = dx*dx  +dy*dy
   30 w(i,j)    = -2.*(g(i,j+1)  -g(i,j))/dn2
c
      return
      end 
c
c ---------------------------------------------------------------
c
      subroutine step (mx,my,x,y,dtinpt,dtstab)
c
c compute time step
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my)
      common/h/  vx(inode,jnode),vy(inode,jnode),
     .           dwx(inode,jnode),dwy(inode,jnode)
      common/j/  re,dt,time,tstop
c
      cfl       = 0.1
      jmax      = (my  -1)/2  +1
c
      dtstab    = 1.0E+99
c
      do 1 j=1,jmax
      do 1 i=1,mx-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
      delta     = sqrt(area)
c     tdiff     = 2.*(1./(dx*dx)  +1./(dy*dy))/re
      tdiff     = 4./(area*re)
c     tconv     = abs(vx(i,j))/dx  +abs(vy(i,j))/dy
      tconv     = (abs(vx(i,j))  +abs(vy(i,j)))/delta
      dtime     = 1./(tdiff  +tconv)
      if (dtime .lt. dtstab) dtstab = dtime
    1 continue
c
      dtstab    = cfl*dtstab
c
      if (dtstab .lt. dtinpt) then
      dt        = dtstab
      else
      dt        = dtinpt
      endif
c
      return
      end
