c
c ---------------------------------------------------------------
c
      subroutine relxfv (mx,my,x,y,g,p,res)
c
c Gauss-Seidel on residual equation
c finite-volume scheme
c
      implicit   real*8 (a-h,o-z)
      parameter  (inode=257,jnode=129)
      dimension  x(mx,my),y(mx,my),g(mx,my),p(mx,my),res(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/g/  cor(inode,jnode)
      dimension  c(2,jnode)
c
      omega     = 1.
c
c compute source term p(i,j)
c
      if (mode .eq. 1) then
      i         = 1
      call resrelx(mx,my,x,y,g,i)
      do 100 i=2,mx-1
      call resrelx(mx,my,x,y,g,i)
      do 100 j=2,my-1
  100 p(i,j)    = res(i,j)  +resid(j)
      else
      endif
c
c solve residual equation
c
      if (mode .le. 1) then
      do 200 it=1,itmax
      do 210 j=2,my-1
  210 c(1,j)    = 0.
      i         = 1
      call resrelx(mx,my,x,y,g,i)
      do 220 i=2,mx-1
      call resrelx(mx,my,x,y,g,i)
      do 230 j=2,my-1
      xx        = .5*(x(i+1,j)  -x(i-1,j))
      xy        = .5*(x(i,j+1)  -x(i,j-1))
      yx        = .5*(y(i+1,j)  -y(i-1,j))
      yy        = .5*(y(i,j+1)  -y(i,j-1))
      det       = 1./(xx*yy  -xy*yx)
      b11       = yy*yy
      b12       = yx*yx
      b13       = -2.*yx*yy
      b21       = xy*xy
      b22       = xx*xx
      b23       = -2.*xx*xy
      c1        = (b11  +b21)*det*det
      c2        = (b12  +b22)*det*det
      c3        = (b13  +b23)*det*det
      rhs       = p(i,j)  -resid(j)
  230 c(2,j)    = omega*(-rhs  +c1*c(1,j)  +c2*c(2,j-1)
     .                   -.25*c3*(c(1,j+1)  -c(1,j-1)))/
     .                  (2.*(c1  +c2))
      if (i .eq. 2) go to 250
      do 240 j=2,my-1
  240 g(i-1,j)  = g(i-1,j)  +c(1,j)
  250 do 260 j=2,my-1
  260 c(1,j)    = c(2,j)
  220 continue
      do 270 j=2,my-1
  270 g(mx-1,j) = g(mx-1,j)  +c(1,j)
  200 continue
      else
      endif
c
c store residual for multigrid use (fine to coarse or injection step)
c
      if (mode .ge. 0) then
      i         = 1
      call resrelx(mx,my,x,y,g,i)
      do 300 i=2,mx-1
      call resrelx(mx,my,x,y,g,i)
      do 300 j=2,my-1
  300 res(i,j)  = p(i,j)  -resid(j)
      else
      endif
c
c locate maximum residual, error, and correction
c
      if (mode .eq. 100) then
      ires      = -1
      jres      = -1
      resmax    = 0.
      do 400 i=2,mx-1
      do 400 j=2,my-1
      if (abs(res(i,j)) .lt. abs(resmax)) go to 400
      ires      = i
      jres      = j
      resmax    = res(i,j)
  400 continue
      icor      = -1
      jcor      = -1
      cormax    = 0.
      do 410 i=2,mx-1
      do 410 j=2,my-1
      if (abs(cor(i,j)) .lt. abs(cormax)) go to 410
      icor      = i
      jcor      = j
      cormax    = cor(i,j)
  410 continue
      else
      endif
c
      return
      end 
c
c ---------------------------------------------------------------
c
      subroutine resrelx (mx,my,x,y,g,i)
c
c compute residual
c
      parameter  (jnode=129)
      implicit   real*8 (a-h,o-z)
      dimension  x(mx,my),y(mx,my),g(mx,my)
      common/e/  resid(jnode),xc(2,jnode),yc(2,jnode),
     .           gxw(jnode),gyw(jnode),gxe(jnode),gye(jnode),
     .           gxns(jnode),gyns(jnode)
c
      if (i .ne. 1) go to 100
      do 10 j=1,my-1
      xc(1,j+1) = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
   10 yc(1,j+1) = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
      do 20 j=2,my-1
      xx        = x(i+1,j)  -x(i,j)
      xy        = .25*(x(i+1,j+1)  -x(i+1,j-1)  +x(i,j+1)  -x(i,j-1))
      yx        = y(i+1,j)  -y(i,j)
      yy        = .25*(y(i+1,j+1)  -y(i+1,j-1)  +y(i,j+1)  -y(i,j-1))
      det       = 1./(xx*yy  -xy*yx)
      dgx       = g(i+1,j)  -g(i,j)
      dgy       = .25*(g(i+1,j+1)  -g(i+1,j-1)  +g(i,j+1)  -g(i,j-1))
      gxw(j)    = (yy*dgx  -yx*dgy)*det
   20 gyw(j)    = (xx*dgy  -xy*dgx)*det
      return
c
  100 continue
      do 30 j=1,my-1
      xc(2,j+1) = .25*(x(i+1,j+1)  +x(i,j+1)  +x(i+1,j)  +x(i,j))
   30 yc(2,j+1) = .25*(y(i+1,j+1)  +y(i,j+1)  +y(i+1,j)  +y(i,j))
      do 40 j=2,my-1
      xx        = x(i+1,j)  -x(i,j)
      xy        = .25*(x(i+1,j+1)  -x(i+1,j-1)  +x(i,j+1)  -x(i,j-1))
      yx        = y(i+1,j)  -y(i,j)
      yy        = .25*(y(i+1,j+1)  -y(i+1,j-1)  +y(i,j+1)  -y(i,j-1))
      det       = 1./(xx*yy  -xy*yx)
      dgx       = g(i+1,j)  -g(i,j)
      dgy       = .25*(g(i+1,j+1)  -g(i+1,j-1)  +g(i,j+1)  -g(i,j-1))
      gxe(j)    = (yy*dgx  -yx*dgy)*det
   40 gye(j)    = (xx*dgy  -xy*dgx)*det
      do 50 j=1,my-1
      xx        = .25*(x(i+1,j+1)  -x(i-1,j+1)  +x(i+1,j)  -x(i-1,j))
      xy        = x(i,j+1)  -x(i,j)
      yx        = .25*(y(i+1,j+1)  -y(i-1,j+1)  +y(i+1,j)  -y(i-1,j))
      yy        = y(i,j+1)  -y(i,j)
      det       = 1./(xx*yy  -xy*yx)
      dgx       = .25*(g(i+1,j+1)  -g(i-1,j+1)  +g(i+1,j)  -g(i-1,j))
      dgy       = g(i,j+1)  -g(i,j)
      gxns(j+1) = (yy*dgx  -yx*dgy)*det
   50 gyns(j+1) = (xx*dgy  -xy*dgx)*det
c
      do 60 j=2,my-1
      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+1)  -xc(1,j)
      dy        = yc(1,j+1)  -yc(1,j)
      flux      = flux  +(-gxw(j)*dy  +gyw(j)*dx)
c east face
      dx        = xc(2,j)  -xc(2,j+1)
      dy        = yc(2,j)  -yc(2,j+1)
      flux      = flux  +(-gxe(j)*dy  +gye(j)*dx)
c south face
      dx        = xc(1,j)  -xc(2,j)
      dy        = yc(1,j)  -yc(2,j)
      flux      = flux  +(-gxns(j)*dy  +gyns(j)*dx)
c north face
      dx        = xc(2,j+1)  -xc(1,j+1)
      dy        = yc(2,j+1)  -yc(1,j+1)
      flux      = flux  +(-gxns(j+1)*dy  +gyns(j+1)*dx)
   60 resid(j)  = flux/area
c
      do 70 j=1,my-1
      xc(1,j+1) = xc(2,j+1)
   70 yc(1,j+1) = yc(2,j+1)
      do 80 j=2,my-1
      gxw(j)   = gxe(j)
   80 gyw(j)   = gye(j)
c
      return
      end
