      subroutine init(neq, y, re, rpar, ipar)

      implicit none

c --------------------------------------------------------------------
c This subroutine loads the tail of y with the data for the 
c Navier-Stokes problem. This is an interface routine. 
c --------------------------------------------------------------------

      integer neq

      integer ipar(*)

      double precision re

      double precision rpar(*)
      double precision y(neq)

      integer lc
      integer ld
      integer lf
      integer lg
      integer lvisc
      integer lwk
      integer n
      integer np2

      save

      n = ipar(1)
      np2 = ipar(2)
      lc = 1
      ld = lc + 4
      lf = ld + 4*n
      lg = lf + n*n
      lvisc = lg + np2*np2
      lwk = lvisc + 1
      call init2(n, np2, y, re, rpar(lc), rpar(ld), rpar(lf), 
     $     rpar(lg), rpar(lvisc), rpar(lwk))
      return
c --------------------------------------------------------------------
c End of subroutine init. 
c --------------------------------------------------------------------
      end

      subroutine init2(n, np2, y, re, c, d, f, g, visc, wk)

      implicit none

c --------------------------------------------------------------------
c This subroutine is called by the interface subroutine init and does 
c the work of loading the tail of y with the data for the Navier-Stokes 
c problem. 
c --------------------------------------------------------------------

      integer n
      integer np2

      double precision re
      double precision visc

      double precision c(4)
      double precision d(n,4)
      double precision f(n,n)
      double precision g(np2,np2)
      double precision wk(np2,np2)
      double precision y(n,*)

      integer i
      integer j
      integer np1

      double precision zero,       half,        one
      parameter      ( zero=0.0d0, half=0.50d0, one=1.0d0 )

      save
c
c     fill viscosity
c
      visc = 1/re
c
c     fill normal derivatives at the corners.
c     (corners are ordered: sw,se,ne,nw)
c
      c(1) = zero
      c(2) = zero
      c(3) = half
      c(4) = half
c
c     fill normal derivatives on the sides.
c     (sides are ordered s,e,n,w)
c
      do 10 j = 1, n
         d(j,1) = zero
         d(j,2) = zero
         d(j,3) = one
         d(j,4) = zero
10    continue
c
c     fill forcing term.
c
      do 30 j = 1, n
         do 20 i = 1, n
            f(i,j) = zero
20       continue
30    continue
c
c     fill dirichlet conditions on auxiliary grid.
c
      np1 = n + 1
c
      do 40 i = 1, np2
         g(i,1) = zero
         g(i,np2) = zero
40    continue
c
      do 50 j = 2, np1
         g(1,j) = zero
         g(np2,j) = zero
50    continue
c
      do 70 j = 1, n
         do 60 i = 1, n
            y(i,j) = zero
 60      continue
 70   continue
      return
c --------------------------------------------------------------------
c End of subroutine init2. 
c --------------------------------------------------------------------
      end

      subroutine fcav(neq, y, fvec, rpar, ipar, itrmf)

      implicit none

c --------------------------------------------------------------------
c This subroutine evaluates the nonlinear residual for the time independent
c Navier-Stokes problem. This is an interface routine. 
c --------------------------------------------------------------------

      integer itrmf
      integer neq

      integer ipar(*)

      double precision fvec(neq)
      double precision rpar(*)
      double precision y(neq)

      integer lc
      integer ld
      integer lf
      integer lg
      integer lvisc
      integer lwk
      integer n
      integer np2

      save

      n = ipar(1)
      np2 = ipar(2)
      lc = 1
      ld = lc + 4
      lf = ld + 4*n
      lg = lf + n*n
      lvisc = lg + np2*np2
      lwk = lvisc + 1
      call resid( n, np2, y, fvec, rpar(lvisc), rpar(lc), rpar(ld), 
     $     rpar(lf), rpar(lg), rpar(lwk))

      itrmf = 0

      return
c --------------------------------------------------------------------
c End of subroutine fcav. 
c --------------------------------------------------------------------
      end

      subroutine resid( n, np2, y, fy, visc, c, d, f, g, w )

      implicit none

c --------------------------------------------------------------------
c This subroutine is called by the interface subroutine fcav. It does 
c the work of evaluating the residual for the time independent 
c Navier-Stokes problem.
c --------------------------------------------------------------------

      integer n
      integer np2

      double precision visc

      double precision c(*)
      double precision d(n,4)
      double precision f(n,*)
      double precision fy(n,*)
      double precision g(np2,*)
      double precision w(np2,*)
      double precision y(n,*)

      integer i
      integer im1
      integer ip1
      integer j
      integer jm1
      integer jp1
      integer np1

      double precision gimjp, gijp, gipjp
      double precision gimj,  gij,  gipj
      double precision gimjm, gijm, gipjm
      double precision h
      double precision h2
      double precision t
      double precision wijp, wipjp
      double precision wimj,  wij,  wipj
      double precision wimjm, wijm
      double precision s, s1, s2, s3, s4, s5, s6

      double precision cvmgz
      external cvmgz

      save

      np1 = n + 1
      h = 1.d0/dfloat(np1)
      h2 = h*h
c
c     load grid array with interior points.
c
      do 20 j = 2, np1
         jm1 = j - 1
         do 10 i = 2, np1
            g(i,j) = y(i-1,jm1)
10       continue
20    continue
      call wload( n, np1, np2, h, g, d, c, w )
c
c     residual at interior points.
c
      do 40 j = 2, np1
         jp1 = j + 1
         jm1 = j - 1
         do 30 i = 2, np1
            ip1 = i + 1
            im1 = i - 1
            gij = g(i,j)
            gijp = g(i,jp1)
            gijm = g(i,jm1)
            gimj = g(im1,j)
            gimjp = g(im1,jp1)
            gimjm = g(im1,jm1)
            gipj = g(ip1,j)
            gipjp = g(ip1,jp1)
            gipjm = g(ip1,jm1)
c
c           linear term
c
            t = 20.d0*gij - 8.d0*( gimj + gipj + gijm + gijp )
     *          + 2.d0*( gimjp + gipjp + gimjm + gipjm )
            if ( i .eq. 2 ) then
               t = t + gij + 2.d0*h*d(jm1,4)
            else
               t = t + g(i-2,j)
            end if
            if ( i .eq. np1 ) then
               t = t + gij + 2.d0*h*d(jm1,2)
            else 
               t = t + g(i+2,j)
            end if
            if ( j .eq. 2 ) then
               t = t + gij + 2.d0*h*d(im1,1)
            else
               t = t + g(i,j-2)
            end if
            if ( j .eq. np1 ) then
               t = t + gij + 2.d0*h*d(im1,3)
            else 
               t = t + g(i,j+2)
            end if
            t = t*visc/h2
c
c           nonlinear term.
c
            wij = w(i,j)
            wipj = w(ip1,j)
            wimj = w(im1,j)
            wijp = w(i,jp1)
            wipjp = w(ip1,jp1)
            wijm = w(i,jm1)
            wimjm = w(im1,jm1)
            s1 = ( wipj - wij )*( gipjp - gipj )
     *           -  ( wipjp - wipj )*( gipj - gij )
            s2 = ( wipjp - wijp )*( gijp - gij )
     *           -  ( wijp - wij )*( gipjp - gijp )
            s3 = ( wij - wimj )*( gijp - gij )
     *           -  ( wijp - wij )*( gij - gimj )
            s4 = ( wij - wimj )*( gimj - gimjm )
     *           -  ( wimj - wimjm )*( gij - gimj )
            s5 = ( wijm - wimjm )*( gij - gijm )
     *           -  ( wij - wijm )*( gijm - gimjm )
            s6 = ( wipj - wij )*( gij - gijm )
     *           -  ( wij - wijm )*( gipj - gij )
            s = ( s1 + s2 + s3 + s4 + s5 + s6 )/(6.d0)
            fy(im1,jm1) = -t - s + h2*f(im1,jm1)
30       continue
40    continue
c
c     include influence of boundary conditions for the poisson
c     problem to be solved next.
c
      do 50 i = 1, n
         ip1 = i + 1
         fy(i,1) = fy(i,1) + g(ip1,1)
         fy(1,i) = fy(1,i) + g(1,ip1)
         fy(n,i) = fy(n,i) + g(np2,ip1)
         fy(i,n) = fy(i,n) + g(ip1,np2)
 50   continue
      return
c --------------------------------------------------------------------
c End of subroutine resid. 
c --------------------------------------------------------------------
      end

      subroutine wload( n, np1, np2, h, g, d, c, w )

      implicit none

c --------------------------------------------------------------------
c This subroutine is called by resid to calculate the 5-spot difference 
c scheme at the i,j grid point, including the boundary.
c --------------------------------------------------------------------

      integer n
      integer np1
      integer np2

      double precision h

      double precision c(*)
      double precision d(n,4)
      double precision g(np2,*)
      double precision w(np2,*)

      double precision two,       three,       four,       six 
      parameter      ( two=2.0d0, three=3.0d0, four=4.0d0, six=6.0d0 )

      integer i
      integer j

      double precision h2

      save
c
c     calculate w = g   + g  .
c                    xx    yy
c
      h2 = h*h
      i = 1
        j = 1
          w(i,j) = (three/two)*(two*g(i,j)-g(i+1,j)-g(i,j+1))/h2
     *             - (three/h)*c(1)
        do 10 j = 2, np1
          w(i,j) = (four*g(i,j)-two*g(i+1,j)-g(i,j-1)-g(i,j+1) )/h2
     *             - (two/h)*d(j-1,4)
 10       continue
        j = np2
          w(i,j) = three*(two*g(i,j)-g(i+1,j)-g(i,j-1))/h2
     *             - (three/h)*c(4)
      do 30 i = 2, np1
        do 20 j = 2, np1
          w(i,j) = (four*g(i,j)-g(i-1,j)-g(i+1,j)-g(i,j-1)-g(i,j+1))/h2
 20       continue
 30     continue
      i = np2
        j = 1
          w(i,j) = three*(two*g(i,j)-g(i-1,j)-g(i,j+1) )/h2
     *             - (six/h)*c(2)
        do 40 j = 2, np1
          w(i,j) = (four*g(i,j)-two*g(i-1,j)-g(i,j-1)-g(i,j+1))/h2
     *             - (two/h)*d(j-1,2)
 40       continue
        j = np2
          w(i,j) = (three/two)*(two*g(i,j)-g(i,j-1)-g(i-1,j))/h2
     *             - (three/h)*c(3)
      j = 1
        do 50 i = 2, np1
          w(i,j) = (four*g(i,j)-two*g(i,j+1)-g(i+1,j)-g(i-1,j))/h2
     *             - (two/h)*d(i-1,1)
 50       continue
      j = np2
        do 60 i = 2, np1
          w(i,j) = (four*g(i,j)-two*g(i,j-1)-g(i-1,j)-g(i+1,j))/h2
     *             - (two/h)*d(i-1,3)
 60       continue
      return
c --------------------------------------------------------------------
c End of subroutine wload. 
c --------------------------------------------------------------------
      end
