      program appbt
c
c***driver for the performance evaluation of the solver for
c   five coupled, nonlinear partial differential equations.
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
c***open file for input data
c
      open (unit=5,file='inputbt.data',status='old',
     *      access='sequential',form='formatted')
      rewind 5
c
c***read the unit number for output data
c
      read (5,*)
      read (5,*)
      read (5,*) iout
c
c***flag that controls printing of the progress of iterations
c
      read (5,*)
      read (5,*)
      read (5,*) ipr, inorm
c
c***set the maximum number of pseudo-time steps to be taken
c
      read (5,*)
      read (5,*)
      read (5,*) itmax
c
c***set the magnitude of the time step
c
      read (5,*)
      read (5,*)
      read (5,*) dt
c
c***set the value of over-relaxation factor for SSOR iterations
c
      read (5,*)
      read (5,*)
      read (5,*) omega
c
c***set the steady-state residual tolerance levels
c
      read (5,*)
      read (5,*)
      read (5,*) tolrsd(1),tolrsd(2),
     $           tolrsd(3),tolrsd(4),tolrsd(5)
c
c***read problem specification parameters
c
c
c***specify the number of grid points in xi, eta and zeta directions
c
      read (5,*)
      read (5,*)
      read (5,*) nx, ny, nz
c
c***open the file for output data
c
      if ( iout .eq. 7 ) then
c
         open (unit=7,file='output.data',status='unknown',
     $        access='sequential',form='formatted')
         rewind 7
c
      end if
c
      if ( ( nx .lt. 5 ) .or.
     $     ( ny .lt. 5 ) .or.
     $     ( nz .lt. 5 ) ) then
c
         write (*,2001)
 2001    format (5x,'PROBLEM SIZE IS TOO SMALL - ',
     $        /5x,'SET EACH OF NX, NY AND NZ AT LEAST EQUAL TO 5')
         stop
c
      end if
c
      if ( ( nx .gt. isiz1 ) .or.
     $     ( ny .gt. isiz2 ) .or.
     $     ( nz .gt. isiz3 ) ) then
c
         write (*,2002)
 2002    format (5x,'PROBLEM SIZE IS TOO LARGE - ',
     $        /5x,'NX, NY AND NZ SHOULD BE LESS THAN OR EQUAL TO ',
     $        /5x,'ISIZ1, ISIZ2 AND ISIZ3 RESPECTIVELY')
c
      end if
c
      dxi = 1.0d+00 / ( nx - 1 )
      deta = 1.0d+00 / ( ny - 1 )
      dzeta = 1.0d+00 / ( nz - 1 )
c
      tx1 = 1.0d+00 / ( dxi * dxi )
      tx2 = 1.0d+00 / ( 2.0d+00 * dxi )
      tx3 = 1.0d+00 / dxi
c
      ty1 = 1.0d+00 / ( deta * deta )
      ty2 = 1.0d+00 / ( 2.0d+00 * deta )
      ty3 = 1.0d+00 / deta
c
      tz1 = 1.0d+00 / ( dzeta * dzeta )
      tz2 = 1.0d+00 / ( 2.0d+00 * dzeta )
      tz3 = 1.0d+00 / dzeta
c
      ii1 = 2
      ii2 = nx - 1
      ji1 = 2
      ji2 = ny - 2
      ki1 = 3
      ki2 = nz - 1
c
c***diffusion coefficients
c
      dx1 = 0.75d+00
      dx2 = dx1
      dx3 = dx1
      dx4 = dx1
      dx5 = dx1
c
      dy1 = 0.75d+00
      dy2 = dy1
      dy3 = dy1
      dy4 = dy1
      dy5 = dy1
c
      dz1 = 1.00d+00
      dz2 = dz1
      dz3 = dz1
      dz4 = dz1 
      dz5 = dz1
c
c***fourth difference dissipation 
c
      dssp = ( max (dx1, dy1, dz1 ) ) / 4.0d+00
c
c
c***coefficients of the exact solution to the first pde
c
      ce(1,1) = 2.0d+00
      ce(1,2) = 0.0d+00
      ce(1,3) = 0.0d+00
      ce(1,4) = 4.0d+00
      ce(1,5) = 5.0d+00
      ce(1,6) = 3.0d+00
      ce(1,7) = 5.0d-01
      ce(1,8) = 2.0d-02
      ce(1,9) = 1.0d-02
      ce(1,10) = 3.0d-02
      ce(1,11) = 5.0d-01
      ce(1,12) = 4.0d-01
      ce(1,13) = 3.0d-01
c
c***coefficients of the exact solution to the second pde
c
      ce(2,1) = 1.0d+00
      ce(2,2) = 0.0d+00
      ce(2,3) = 0.0d+00
      ce(2,4) = 0.0d+00
      ce(2,5) = 1.0d+00
      ce(2,6) = 2.0d+00
      ce(2,7) = 3.0d+00
      ce(2,8) = 1.0d-02
      ce(2,9) = 3.0d-02
      ce(2,10) = 2.0d-02
      ce(2,11) = 4.0d-01
      ce(2,12) = 3.0d-01
      ce(2,13) = 5.0d-01
c
c***coefficients of the exact solution to the third pde
c
      ce(3,1) = 2.0d+00
      ce(3,2) = 2.0d+00
      ce(3,3) = 0.0d+00
      ce(3,4) = 0.0d+00
      ce(3,5) = 0.0d+00
      ce(3,6) = 2.0d+00
      ce(3,7) = 3.0d+00
      ce(3,8) = 4.0d-02
      ce(3,9) = 3.0d-02
      ce(3,10) = 5.0d-02
      ce(3,11) = 3.0d-01
      ce(3,12) = 5.0d-01
      ce(3,13) = 4.0d-01
c
c***coefficients of the exact solution to the fourth pde
c
      ce(4,1) = 2.0d+00
      ce(4,2) = 2.0d+00
      ce(4,3) = 0.0d+00
      ce(4,4) = 0.0d+00
      ce(4,5) = 0.0d+00
      ce(4,6) = 2.0d+00
      ce(4,7) = 3.0d+00
      ce(4,8) = 3.0d-02
      ce(4,9) = 5.0d-02
      ce(4,10) = 4.0d-02
      ce(4,11) = 2.0d-01
      ce(4,12) = 1.0d-01
      ce(4,13) = 3.0d-01
c
c***coefficients of the exact solution to the fifth pde
c
      ce(5,1) = 5.0d+00
      ce(5,2) = 4.0d+00
      ce(5,3) = 3.0d+00
      ce(5,4) = 2.0d+00
      ce(5,5) = 1.0d-01
      ce(5,6) = 4.0d-01
      ce(5,7) = 3.0d-01
      ce(5,8) = 5.0d-02
      ce(5,9) = 4.0d-02
      ce(5,10) = 3.0d-02
      ce(5,11) = 1.0d-01
      ce(5,12) = 3.0d-01
      ce(5,13) = 2.0d-01
c
c***set the boundary values for dependent variables
c
      call setbv
c
c***set the initial values for dependent variables
c
      call setiv
c
c***compute the forcing term based on prescribed exact solution
c
      call erhs
c
c***perform Block approximate factorization iterations
c
      call badi
c
c***compute the solution error
c
      call error
c
c***compute the surface integral
c
      call pintgr
c
c***verification test
c
      call verify ( rsdnm, errnm, frc )
c
c***print the CPU time
c
      write (iout,1001) ttotal
 1001 format (//5x,'Total CPU time = ',1pe12.4,' Sec.')
c
      stop
      end
c
c
c
c
c
      subroutine setbv
c
c***set the boundary values of dependent variables
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
c***set the dependent variable values along the top and bottom faces
c
      do j = 1, ny
         do i = 1, nx
c
            call exact ( i, j, 1, u( 1, i, j, 1 ) )
            call exact ( i, j, nz, u( 1, i, j, nz ) )
c
         end do
      end do
c
c***set the dependent variable values along north and south faces
c
      do k = 1, nz
         do i = 1, nx
c
            call exact ( i, 1, k, u( 1, i, 1, k ) )
            call exact ( i, ny, k, u( 1, i, ny, k ) )
c
         end do
      end do
c
c***set the dependent variable values along east and west faces
c
      do k = 1, nz
         do j = 1, ny
c
            call exact ( 1, j, k, u( 1, 1, j, k ) )
            call exact ( nx, j, k, u( 1, nx, j, k ) )
c
         end do
      end do
c
      return
      end
c
c
c
c
c
      subroutine setiv
c
c***set the initial values of independent variables based on tri-linear
c   interpolation of boundary values in the computational space.
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      do k = 2, nz-1
c
         zeta = ( dfloat (k-1) ) / (nz-1)
c
         do j = 2, ny-1
c
            eta = ( dfloat (j-1) ) / (ny-1)
c
            do i = 2, nx-1
c
               xi = ( dfloat (i-1) ) / (nx-1)
c
               do m = 1, 5
c
                  pxi =   ( 1.0d+00 - xi ) * u(m,1,j,k)
     $                              + xi   * u(m,nx,j,k)
c
                  peta =  ( 1.0d+00 - eta ) * u(m,i,1,k)
     $                              + eta   * u(m,i,ny,k)
c
                  pzeta = ( 1.0d+00 - zeta ) * u(m,i,j,1)
     $                              + zeta   * u(m,i,j,nz)
c
                  u( m, i, j, k ) = pxi + peta + pzeta
     $                 - pxi * peta - peta * pzeta - pzeta * pxi
     $                 + pxi * peta * pzeta
c
               end do
c
            end do
c
         end do
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine badi
c
c***to perform pseudo-time stepping iterations
c   for five coupled, nonlinear pde's.
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension idmax(5), jdmax(5), kdmax(5),
     $          imax(5), jmax(5), kmax(5),
     $          delunm(5)
      parameter ( one = 1.0d+00 )
c
      lnorm = 2
c
c***compute the steady-state residuals
c     
      call rhs
c     
      if ( lnorm .eq. 1 ) then
c
         call maxnorm ( isiz1, isiz2, isiz3,
     $                  nx, ny, nz,
     $                  imax, jmax, kmax,
     $                  rsd, rsdnm )
c
         if ( ipr .eq. 1 ) then
c
            write (iout,*) '          Initial Residual norms'
            write (iout,*)
            write (iout,1003) ( rsdnm(m),
     $           imax(m), jmax(m), kmax(m), m = 1, 5 )
c
         end if
c     
      else if ( lnorm .eq. 2 ) then
c
         call l2norm ( isiz1, isiz2, isiz3,
     $                 nx, ny, nz,
     $                 rsd, rsdnm )
c     
         if ( ipr .eq. 1 ) then
c
            write (iout,*) '          Initial Residual norms'
            write (iout,*)
            write (iout,1007) ( rsdnm(m), m = 1, 5 )
c
         end if
c     
      end if
c     
c***begin pseudo-time stepping iterations
c
      tstart = timer()
c
      do istep = 1, itmax
c
         if ( ( mod ( istep, inorm ) .eq. 0 ) .and.
     $        ( ipr .eq. 1 ) ) then
c
            write ( iout, 1001 ) istep
c
         end if
c
         do k = 1, nz
            do j = 1, ny
               do i = 1, nx
                  do m = 1, 5
c
                     rsd(m,i,j,k) = dt * rsd(m,i,j,k)
c
                  end do
               end do
            end do
         end do
c
c***perform 3-factor, Block ADI iterations
c
c***perform the xi-direction sweep
c
         call jacx
c
         call btridx ( isiz1, isiz2, isiz3,
     $                 nx, ny, nz,
     $                 a,b,c,
     $                 rsd )
c
c***perform the eta-direction sweep
c
         call jacy
c
         call btridy ( isiz1, isiz2, isiz3,
     $                 nx, ny, nz,
     $                 a,b,c,
     $                 rsd )
c
c***perform the zeta-direction sweep
c
         call jacz
c
         call btridz ( isiz1, isiz2, isiz3,
     $                 nx, ny, nz,
     $                 a,b,c,
     $                 rsd )
c
c***update the variables
c
         do k = 1, nz
            do j = 1, ny
               do i = 1, nx
                  do m = 1, 5
c
                     u( m, i, j, k ) =  u( m, i, j, k )
     $                                + rsd( m, i, j, k )
c
                  end do
               end do
            end do
         end do
c
c***compute the max-norms of pseudo-time iteration corrections
c
         if ( mod ( istep, inorm ) .eq. 0 ) then
c
            if ( lnorm .eq. 1 ) then
c
               call maxnorm ( isiz1, isiz2, isiz3,
     $                        nx, ny, nz,
     $                        idmax, jdmax, kdmax,
     $                        rsd, delunm )
c
               if ( ipr .eq. 1 ) then
c
                  write (iout,1002) ( delunm(m),
     $                 idmax(m), jdmax(m), kdmax(m), m = 1, 5 )
c
               end if
c
            else if ( lnorm .eq. 2 ) then
c
               call l2norm ( isiz1, isiz2, isiz3,
     $                       nx, ny, nz,
     $                      rsd, delunm )
c
               if ( ipr .eq. 1 ) then
c
                  write (iout,1006) ( delunm(m), m = 1, 5 )
c
               end if
c
            end if
c
         end if
c
c***compute the steady-state residuals
c
         call rhs
c
         if ( ( mod ( istep, inorm ) .eq. 0 ) .or.
     $        ( istep .eq. itmax ) ) then
c
            if ( lnorm .eq. 1 ) then
c
               call maxnorm ( isiz1, isiz2, isiz3,
     $                        nx, ny, nz,
     $                        imax, jmax, kmax,
     $                        rsd, rsdnm )
c
               if ( ipr .eq. 1 ) then
c
                  write (iout,1003) ( rsdnm(m),
     $                 imax(m), jmax(m), kmax(m), m = 1, 5 )
c
               end if
c     
            else if ( lnorm .eq. 2 ) then
c
               call l2norm ( isiz1, isiz2, isiz3,
     $              nx, ny, nz,
     $              rsd, rsdnm )
c
               if ( ipr .eq. 1 ) then
c
                  write (iout,1007) ( rsdnm(m), m = 1, 5 )
c
               end if
c
            end if
c
         end if
c
c***check the pseudo-time iteration residuals against the tolerance levels
c
         if ( ( rsdnm(1) .lt. tolrsd(1) ) .and.
     $        ( rsdnm(2) .lt. tolrsd(2) ) .and.
     $        ( rsdnm(3) .lt. tolrsd(3) ) .and.
     $        ( rsdnm(4) .lt. tolrsd(4) ) .and.
     $        ( rsdnm(5) .lt. tolrsd(5) ) ) then
c
            write (iout,1004) istep
            return
c
         end if
c
      end do
c
      tend = timer( )
c
      ttotal = tend - tstart
c
      return
c     
 1001 format (1x/5x,'pseudo-time Block ADI iteration no.=',i4/)
c
 1002 format (1x/1x,'max-norm of Block ADI-iteration correction ',
     $ 'for first pde  = ',1pe12.5/,
     $ 58x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Block ADI-iteration correction ',
     $ 'for second pde = ',1pe12.5/,
     $ 58x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Block ADI-iteration correction ',
     $ 'for third pde  = ',1pe12.5/,
     $ 58x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Block ADI-iteration correction ',
     $ 'for fourth pde = ',1pe12.5/,
     $ 58x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Block ADI-iteration correction ',
     $ 'for fifth pde  = ',1pe12.5/,
     $ 58x,'(',i4,',',i4,',',i4,')' )
c
 1003 format (1x/1x,'max-norm of steady-state residual for ',
     $ 'first pde  = ',1pe12.5/,
     $ 49x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of steady-state residual for ',
     $ 'second pde = ',1pe12.5/,
     $ 49x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of steady-state residual for ',
     $ 'third pde  = ',1pe12.5/,
     $ 49x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of steady-state residual for ',
     $ 'fourth pde = ',1pe12.5/,
     $ 49x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of steady-state residual for ',
     $ 'fifth pde  = ',1pe12.5/,
     $ 49x,'(',i4,',',i4,',',i4,')' )
c
 1004 format (1x/1x,'convergence was achieved after ',i4,
     $   ' pseudo-time steps' )
c
 1006 format (1x/1x,'RMS-norm of Block ADI-iteration correction ',
     $ 'for first pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of Block ADI-iteration correction ',
     $ 'for second pde = ',1pe12.5/,
     $ 1x,'RMS-norm of Block ADI-iteration correction ',
     $ 'for third pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of Block ADI-iteration correction ',
     $ 'for fourth pde = ',1pe12.5/,
     $ 1x,'RMS-norm of Block ADI-iteration correction ',
     $ 'for fifth pde  = ',1pe12.5)
c
 1007 format (1x/1x,'RMS-norm of steady-state residual for ',
     $ 'first pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of steady-state residual for ',
     $ 'second pde = ',1pe12.5/,
     $ 1x,'RMS-norm of steady-state residual for ',
     $ 'third pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of steady-state residual for ',
     $ 'fourth pde = ',1pe12.5/,
     $ 1x,'RMS-norm of steady-state residual for ',
     $ 'fifth pde  = ',1pe12.5)
c
      end
c
c
c
c
c
      subroutine maxnorm ( ldx, ldy, ldz,
     $                     nx, ny, nz,
     $                     imax, jmax, kmax, v, vnm )
c
c***to compute the max-norm of vector v.
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      implicit real*8 (a-h,o-z)
      dimension v(5,ldx,ldy,*),
     $          vnm(*), imax(*), jmax(*), kmax(*)
c
      do m = 1, 5
c
         vnm(m) = -1.0d+10
c
      end do
c
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
c
               do m = 1, 5
c
                  t1 = abs ( v( m, i, j, k ) )
c
                  if ( vnm(m) .lt. t1 ) then
                     vnm(m) = t1
                     imax(m) = i
                     jmax(m) = j
                     kmax(m) = k
                  end if
c
               end do
c
            end do
c
         end do
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine l2norm ( ldx, ldy, ldz,
     $                    nx, ny, nz,
     $                    v, sum )
c
c***to compute the l2-norm of vector v.
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      implicit real*8 (a-h,o-z)
      dimension v(5,ldx,ldy,*),
     $          sum(*)
c
      do m = 1, 5
c
         sum(m) = 0.0d+00
c
      end do
c
      do k = 2, nz-1
         do j = 2, ny-1
            do i = 2, nx-1
c
               do m = 1, 5
c
                  sum(m) = sum(m) + v(m,i,j,k) * v(m,i,j,k)
c
               end do
c
            end do
         end do
      end do
c
      do m = 1, 5
c
         sum(m) = sqrt ( sum(m) / ( (nx-2)*(ny-2)*(nz-2) ) )
c
      end do
c
      return
      end
c
c
c
c
c
c
      subroutine btridx ( isiz1, isiz2, isiz3,
     $                    nx, ny, nz,
     $                    a, b, c, f )
c
c***solution of multiple, independent systems of block tridiagonal systems 
c   using Gaussian elimination (without pivoting) algorithm
c   (block size = 5 X 5)
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      implicit real*8 (a-h,o-z)
c
      parameter ( zero = 0.0d+00, one = 1.0d+00 )
c
      dimension a(5,5,isiz1,isiz2,*),
     $          b(5,5,isiz1,isiz2,*),
     $          c(5,5,isiz1,isiz2,*),
     $          f(5,isiz1,isiz2,*)
c
      dimension tm(5,5),tv(5)
c
c***forward elimination phase
c
      do k = 2, nz-1
c
         do j = 2, ny-1
c
            tmp = one / b(1,1,1,j,k)
c
            b(2,1,1,j,k) = b(2,1,1,j,k) * tmp
            b(3,1,1,j,k) = b(3,1,1,j,k) * tmp
            b(4,1,1,j,k) = b(4,1,1,j,k) * tmp
            b(5,1,1,j,k) = b(5,1,1,j,k) * tmp
c
            b(2,2,1,j,k) =  b(2,2,1,j,k)
     $                    - b(2,1,1,j,k) * b(1,2,1,j,k)
c
            tmp = one / b(2,2,1,j,k)
c
            b(3,2,1,j,k) = (  b(3,2,1,j,k)
     $                      - b(3,1,1,j,k) * b(1,2,1,j,k) ) * tmp
            b(4,2,1,j,k) = (  b(4,2,1,j,k)
     $                      - b(4,1,1,j,k) * b(1,2,1,j,k) ) * tmp
            b(5,2,1,j,k) = (  b(5,2,1,j,k)
     $                      - b(5,1,1,j,k) * b(1,2,1,j,k) ) * tmp
c
            b(2,3,1,j,k) =  b(2,3,1,j,k)
     $                    - b(2,1,1,j,k) * b(1,3,1,j,k)
            b(3,3,1,j,k) =  b(3,3,1,j,k)
     $                    - b(3,1,1,j,k) * b(1,3,1,j,k)
     $                    - b(3,2,1,j,k) * b(2,3,1,j,k)
c
            tmp = one / b(3,3,1,j,k)
c
            b(4,3,1,j,k) = (  b(4,3,1,j,k)
     $                      - b(4,1,1,j,k) * b(1,3,1,j,k)
     $                      - b(4,2,1,j,k) * b(2,3,1,j,k) ) * tmp
            b(5,3,1,j,k) = (  b(5,3,1,j,k)
     $                      - b(5,1,1,j,k) * b(1,3,1,j,k)
     $                      - b(5,2,1,j,k) * b(2,3,1,j,k) ) * tmp
c
            b(2,4,1,j,k) =  b(2,4,1,j,k)
     $                    - b(2,1,1,j,k) * b(1,4,1,j,k)
            b(3,4,1,j,k) =  b(3,4,1,j,k)
     $                    - b(3,1,1,j,k) * b(1,4,1,j,k)
     $                    - b(3,2,1,j,k) * b(2,4,1,j,k)
            b(4,4,1,j,k) =  b(4,4,1,j,k)
     $                    - b(4,1,1,j,k) * b(1,4,1,j,k)
     $                    - b(4,2,1,j,k) * b(2,4,1,j,k)
     $                    - b(4,3,1,j,k) * b(3,4,1,j,k)
c
            b(5,4,1,j,k) = (  b(5,4,1,j,k)
     $                      - b(5,1,1,j,k) * b(1,4,1,j,k)
     $                      - b(5,2,1,j,k) * b(2,4,1,j,k)
     $                      - b(5,3,1,j,k) * b(3,4,1,j,k) )
     $                     / b(4,4,1,j,k)
c
            b(2,5,1,j,k) =  b(2,5,1,j,k)
     $                    - b(2,1,1,j,k) * b(1,5,1,j,k)
            b(3,5,1,j,k) =  b(3,5,1,j,k)
     $                    - b(3,1,1,j,k) * b(1,5,1,j,k)
     $                    - b(3,2,1,j,k) * b(2,5,1,j,k)
            b(4,5,1,j,k) =  b(4,5,1,j,k)
     $                    - b(4,1,1,j,k) * b(1,5,1,j,k)
     $                    - b(4,2,1,j,k) * b(2,5,1,j,k)
     $                    - b(4,3,1,j,k) * b(3,5,1,j,k)
c
            b(5,5,1,j,k) =  b(5,5,1,j,k)
     $                    - b(5,1,1,j,k) * b(1,5,1,j,k)
     $                    - b(5,2,1,j,k) * b(2,5,1,j,k)
     $                    - b(5,3,1,j,k) * b(3,5,1,j,k)
     $                    - b(5,4,1,j,k) * b(4,5,1,j,k)
c
            do i = 2, nx
c
               do m = 1, 5
c
                  tm(1,m) = c(1,m,i-1,j,k)
                  tm(2,m) = c(2,m,i-1,j,k) - b(2,1,i-1,j,k) * tm(1,m)
                  tm(3,m) = c(3,m,i-1,j,k) - b(3,1,i-1,j,k) * tm(1,m)
     $                                     - b(3,2,i-1,j,k) * tm(2,m)
                  tm(4,m) = c(4,m,i-1,j,k) - b(4,1,i-1,j,k) * tm(1,m)
     $                                     - b(4,2,i-1,j,k) * tm(2,m)
     $                                     - b(4,3,i-1,j,k) * tm(3,m)
                  tm(5,m) = c(5,m,i-1,j,k) - b(5,1,i-1,j,k) * tm(1,m)
     $                                     - b(5,2,i-1,j,k) * tm(2,m)
     $                                     - b(5,3,i-1,j,k) * tm(3,m)
     $                                     - b(5,4,i-1,j,k) * tm(4,m)
c
               end do
c
               do m = 1, 5
c
                  tm(5,m) = tm(5,m) / b(5,5,i-1,j,k)
                  tm(4,m) = (  tm(4,m)
     $                       - b(4,5,i-1,j,k) * tm(5,m) )
     $                      / b(4,4,i-1,j,k)
                  tm(3,m) = (  tm(3,m)
     $                       - b(3,4,i-1,j,k) * tm(4,m)
     $                       - b(3,5,i-1,j,k) * tm(5,m) )
     $                      / b(3,3,i-1,j,k)
                  tm(2,m) = (  tm(2,m)
     $                       - b(2,3,i-1,j,k) * tm(3,m)
     $                       - b(2,4,i-1,j,k) * tm(4,m)
     $                       - b(2,5,i-1,j,k) * tm(5,m) )
     $                      / b(2,2,i-1,j,k)
                  tm(1,m) = (  tm(1,m)
     $                       - b(1,2,i-1,j,k) * tm(2,m)
     $                       - b(1,3,i-1,j,k) * tm(3,m)
     $                       - b(1,4,i-1,j,k) * tm(4,m)
     $                       - b(1,5,i-1,j,k) * tm(5,m) )
     $                      / b(1,1,i-1,j,k)
c
               end do
c
               do m = 1, 5
c
                  b(1,m,i,j,k) =  b(1,m,i,j,k)
     $                          - a(1,1,i,j,k) * tm(1,m)
     $                          - a(1,2,i,j,k) * tm(2,m)
     $                          - a(1,3,i,j,k) * tm(3,m)
     $                          - a(1,4,i,j,k) * tm(4,m)
     $                          - a(1,5,i,j,k) * tm(5,m)
c
                  b(2,m,i,j,k) =  b(2,m,i,j,k)
     $                          - a(2,1,i,j,k) * tm(1,m)
     $                          - a(2,2,i,j,k) * tm(2,m)
     $                          - a(2,3,i,j,k) * tm(3,m)
     $                          - a(2,4,i,j,k) * tm(4,m)
     $                          - a(2,5,i,j,k) * tm(5,m)
c
                  b(3,m,i,j,k) =  b(3,m,i,j,k)
     $                          - a(3,1,i,j,k) * tm(1,m)
     $                          - a(3,2,i,j,k) * tm(2,m)
     $                          - a(3,3,i,j,k) * tm(3,m)
     $                          - a(3,4,i,j,k) * tm(4,m)
     $                          - a(3,5,i,j,k) * tm(5,m)
c
                  b(4,m,i,j,k) =  b(4,m,i,j,k)
     $                          - a(4,1,i,j,k) * tm(1,m)
     $                          - a(4,2,i,j,k) * tm(2,m)
     $                          - a(4,3,i,j,k) * tm(3,m)
     $                          - a(4,4,i,j,k) * tm(4,m)
     $                          - a(4,5,i,j,k) * tm(5,m)
c
                  b(5,m,i,j,k) =  b(5,m,i,j,k)
     $                          - a(5,1,i,j,k) * tm(1,m)
     $                          - a(5,2,i,j,k) * tm(2,m)
     $                          - a(5,3,i,j,k) * tm(3,m)
     $                          - a(5,4,i,j,k) * tm(4,m)
     $                          - a(5,5,i,j,k) * tm(5,m)
c
               end do
c
               tv(1) = f(1,i-1,j,k)
               tv(2) = f(2,i-1,j,k) - b(2,1,i-1,j,k) * tv(1)
               tv(3) = f(3,i-1,j,k) - b(3,1,i-1,j,k) * tv(1)
     $                              - b(3,2,i-1,j,k) * tv(2)
               tv(4) = f(4,i-1,j,k) - b(4,1,i-1,j,k) * tv(1)
     $                              - b(4,2,i-1,j,k) * tv(2)
     $                              - b(4,3,i-1,j,k) * tv(3)
               tv(5) = f(5,i-1,j,k) - b(5,1,i-1,j,k) * tv(1)
     $                              - b(5,2,i-1,j,k) * tv(2)
     $                              - b(5,3,i-1,j,k) * tv(3)
     $                              - b(5,4,i-1,j,k) * tv(4)
c
               tv(5) = tv(5) / b(5,5,i-1,j,k)
               tv(4) = (  tv(4)
     $                  - b(4,5,i-1,j,k) * tv(5) )
     $                 / b(4,4,i-1,j,k)
               tv(3) = (  tv(3)
     $                  - b(3,4,i-1,j,k) * tv(4)
     $                  - b(3,5,i-1,j,k) * tv(5) )
     $                 / b(3,3,i-1,j,k)
               tv(2) = (  tv(2)
     $                  - b(2,3,i-1,j,k) * tv(3)
     $                  - b(2,4,i-1,j,k) * tv(4)
     $                  - b(2,5,i-1,j,k) * tv(5) )
     $                 / b(2,2,i-1,j,k)
               tv(1) = (  tv(1)
     $                  - b(1,2,i-1,j,k) * tv(2)
     $                  - b(1,3,i-1,j,k) * tv(3)
     $                  - b(1,4,i-1,j,k) * tv(4)
     $                  - b(1,5,i-1,j,k) * tv(5) )
     $                 / b(1,1,i-1,j,k)
c
               do m = 1, 5
c
                  f(m,i,j,k) =  f(m,i,j,k)
     $                        - a(m,1,i,j,k) * tv(1)
     $                        - a(m,2,i,j,k) * tv(2)
     $                        - a(m,3,i,j,k) * tv(3)
     $                        - a(m,4,i,j,k) * tv(4)
     $                        - a(m,5,i,j,k) * tv(5)
c
               end do
c
               tmp = one / b(1,1,i,j,k)
c
               b(2,1,i,j,k) = b(2,1,i,j,k) * tmp
               b(3,1,i,j,k) = b(3,1,i,j,k) * tmp
               b(4,1,i,j,k) = b(4,1,i,j,k) * tmp
               b(5,1,i,j,k) = b(5,1,i,j,k) * tmp
c
               b(2,2,i,j,k) =  b(2,2,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,2,i,j,k)
c
               tmp = one / b(2,2,i,j,k)
c
               b(3,2,i,j,k) = (  b(3,2,i,j,k)
     $                         - b(3,1,i,j,k) * b(1,2,i,j,k) ) * tmp
               b(4,2,i,j,k) = (  b(4,2,i,j,k)
     $                         - b(4,1,i,j,k) * b(1,2,i,j,k) ) * tmp
               b(5,2,i,j,k) = (  b(5,2,i,j,k)
     $                         - b(5,1,i,j,k) * b(1,2,i,j,k) ) * tmp
c
               b(2,3,i,j,k) =  b(2,3,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,3,i,j,k)
               b(3,3,i,j,k) =  b(3,3,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,3,i,j,k)
c
               tmp = one / b(3,3,i,j,k)
c
               b(4,3,i,j,k) = (  b(4,3,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,3,i,j,k) ) * tmp
               b(5,3,i,j,k) = (  b(5,3,i,j,k)
     $                       - b(5,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(5,2,i,j,k) * b(2,3,i,j,k) ) * tmp
c
               b(2,4,i,j,k) =  b(2,4,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,4,i,j,k)
               b(3,4,i,j,k) =  b(3,4,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,4,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,4,i,j,k)
               b(4,4,i,j,k) =  b(4,4,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,4,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,4,i,j,k)
     $                       - b(4,3,i,j,k) * b(3,4,i,j,k)
c
               b(5,4,i,j,k) = (  b(5,4,i,j,k)
     $                         - b(5,1,i,j,k) * b(1,4,i,j,k)
     $                         - b(5,2,i,j,k) * b(2,4,i,j,k)
     $                         - b(5,3,i,j,k) * b(3,4,i,j,k) )
     $                        / b(4,4,i,j,k)
c
               b(2,5,i,j,k) =  b(2,5,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,5,i,j,k)
               b(3,5,i,j,k) =  b(3,5,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,5,i,j,k)
               b(4,5,i,j,k) =  b(4,5,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,5,i,j,k)
     $                       - b(4,3,i,j,k) * b(3,5,i,j,k)
c
               b(5,5,i,j,k) =  b(5,5,i,j,k)
     $                       - b(5,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(5,2,i,j,k) * b(2,5,i,j,k)
     $                       - b(5,3,i,j,k) * b(3,5,i,j,k)
     $                       - b(5,4,i,j,k) * b(4,5,i,j,k)
c
            end do
c     
         end do
c
      end do
c
c***back-substitution phase
c
      do k = 2, nz-1
c
         do j = 2, ny-1
c
            f(2,nx,j,k) =  f(2,nx,j,k)
     $                   - b(2,1,nx,j,k) * f(1,nx,j,k)  
            f(3,nx,j,k) =  f(3,nx,j,k)
     $                   - b(3,1,nx,j,k) * f(1,nx,j,k)
     $                   - b(3,2,nx,j,k) * f(2,nx,j,k)
            f(4,nx,j,k) =  f(4,nx,j,k)
     $                   - b(4,1,nx,j,k) * f(1,nx,j,k)
     $                   - b(4,2,nx,j,k) * f(2,nx,j,k)
     $                   - b(4,3,nx,j,k) * f(3,nx,j,k)
            f(5,nx,j,k) =  f(5,nx,j,k)
     $                   - b(5,1,nx,j,k) * f(1,nx,j,k)
     $                   - b(5,2,nx,j,k) * f(2,nx,j,k)
     $                   - b(5,3,nx,j,k) * f(3,nx,j,k)
     $                   - b(5,4,nx,j,k) * f(4,nx,j,k)
c
            f(5,nx,j,k) = f(5,nx,j,k) / b(5,5,nx,j,k)
            f(4,nx,j,k) = (  f(4,nx,j,k)
     $                     - b(4,5,nx,j,k) * f(5,nx,j,k) )
     $                    / b(4,4,nx,j,k)
            f(3,nx,j,k) = (  f(3,nx,j,k)
     $                     - b(3,4,nx,j,k) * f(4,nx,j,k)
     $                     - b(3,5,nx,j,k) * f(5,nx,j,k) )
     $                    / b(3,3,nx,j,k)
            f(2,nx,j,k) = (  f(2,nx,j,k)
     $                     - b(2,3,nx,j,k) * f(3,nx,j,k)
     $                     - b(2,4,nx,j,k) * f(4,nx,j,k)
     $                     - b(2,5,nx,j,k) * f(5,nx,j,k) )
     $                    / b(2,2,nx,j,k)
            f(1,nx,j,k) = (  f(1,nx,j,k)
     $                     - b(1,2,nx,j,k) * f(2,nx,j,k)
     $                     - b(1,3,nx,j,k) * f(3,nx,j,k)
     $                     - b(1,4,nx,j,k) * f(4,nx,j,k)
     $                     - b(1,5,nx,j,k) * f(5,nx,j,k) )
     $                    / b(1,1,nx,j,k)
c
            do i = nx-1, 1, -1
c
               do m = 1, 5
c
                  f(m,i,j,k) =  f(m,i,j,k)
     $                        - c(m,1,i,j,k) * f(1,i+1,j,k)
     $                        - c(m,2,i,j,k) * f(2,i+1,j,k)
     $                        - c(m,3,i,j,k) * f(3,i+1,j,k)
     $                        - c(m,4,i,j,k) * f(4,i+1,j,k)
     $                        - c(m,5,i,j,k) * f(5,i+1,j,k)
c
               end do
c
               f(2,i,j,k) =  f(2,i,j,k)
     $                     - b(2,1,i,j,k) * f(1,i,j,k)  
               f(3,i,j,k) =  f(3,i,j,k)
     $                     - b(3,1,i,j,k) * f(1,i,j,k)
     $                     - b(3,2,i,j,k) * f(2,i,j,k)
               f(4,i,j,k) =  f(4,i,j,k)
     $                     - b(4,1,i,j,k) * f(1,i,j,k)
     $                     - b(4,2,i,j,k) * f(2,i,j,k)
     $                     - b(4,3,i,j,k) * f(3,i,j,k)
               f(5,i,j,k) =  f(5,i,j,k)
     $                     - b(5,1,i,j,k) * f(1,i,j,k)
     $                     - b(5,2,i,j,k) * f(2,i,j,k)
     $                     - b(5,3,i,j,k) * f(3,i,j,k)
     $                     - b(5,4,i,j,k) * f(4,i,j,k)
c
               f(5,i,j,k) = f(5,i,j,k) / b(5,5,i,j,k)
               f(4,i,j,k) = (  f(4,i,j,k)
     $                       - b(4,5,i,j,k) * f(5,i,j,k) )
     $                      / b(4,4,i,j,k)
               f(3,i,j,k) = (  f(3,i,j,k)
     $                       - b(3,4,i,j,k) * f(4,i,j,k)
     $                       - b(3,5,i,j,k) * f(5,i,j,k) )
     $                      / b(3,3,i,j,k)
               f(2,i,j,k) = (  f(2,i,j,k)
     $                       - b(2,3,i,j,k) * f(3,i,j,k)
     $                       - b(2,4,i,j,k) * f(4,i,j,k)
     $                       - b(2,5,i,j,k) * f(5,i,j,k) )
     $                      / b(2,2,i,j,k)
               f(1,i,j,k) = (  f(1,i,j,k)
     $                       - b(1,2,i,j,k) * f(2,i,j,k)
     $                       - b(1,3,i,j,k) * f(3,i,j,k)
     $                       - b(1,4,i,j,k) * f(4,i,j,k)
     $                       - b(1,5,i,j,k) * f(5,i,j,k) )
     $                      / b(1,1,i,j,k)
c               
            end do
c
         end do
c
      end do
c
      return
      end         
c
c
c
c
c
c
      subroutine btridy ( isiz1, isiz2, isiz3,
     $                    nx, ny, nz,
     $                    a, b, c, f )
c
c***solution of multiple, independent systems of block tridiagonal systems 
c   using Gaussian elimination (without pivoting) algorithm
c   (block size = 5 X 5)
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      implicit real*8 (a-h,o-z)
c
      parameter ( zero = 0.0d+00, one = 1.0d+00 )
c
      dimension a(5,5,isiz1,isiz2,*),
     $          b(5,5,isiz1,isiz2,*),
     $          c(5,5,isiz1,isiz2,*),
     $          f(5,isiz1,isiz2,*)
c
      dimension tm(5,5), tv(5)
c
c***forward elimination phase
c
      do k = 2, nz-1
c
         do i = 2, nx-1
c
            tmp = one / b(1,1,i,1,k)
c
            b(2,1,i,1,k) = b(2,1,i,1,k) * tmp
            b(3,1,i,1,k) = b(3,1,i,1,k) * tmp
            b(4,1,i,1,k) = b(4,1,i,1,k) * tmp
            b(5,1,i,1,k) = b(5,1,i,1,k) * tmp
c
            b(2,2,i,1,k) =  b(2,2,i,1,k)
     $                    - b(2,1,i,1,k) * b(1,2,i,1,k)
c
            tmp = one / b(2,2,i,1,k)
c
            b(3,2,i,1,k) = (  b(3,2,i,1,k)
     $                      - b(3,1,i,1,k) * b(1,2,i,1,k) ) * tmp
            b(4,2,i,1,k) = (  b(4,2,i,1,k)
     $                      - b(4,1,i,1,k) * b(1,2,i,1,k) ) * tmp
            b(5,2,i,1,k) = (  b(5,2,i,1,k)
     $                      - b(5,1,i,1,k) * b(1,2,i,1,k) ) * tmp
c
            b(2,3,i,1,k) =  b(2,3,i,1,k)
     $                    - b(2,1,i,1,k) * b(1,3,i,1,k)
            b(3,3,i,1,k) =  b(3,3,i,1,k)
     $                    - b(3,1,i,1,k) * b(1,3,i,1,k)
     $                    - b(3,2,i,1,k) * b(2,3,i,1,k)
c
            tmp = one / b(3,3,i,1,k)
c
            b(4,3,i,1,k) = (  b(4,3,i,1,k)
     $                      - b(4,1,i,1,k) * b(1,3,i,1,k)
     $                      - b(4,2,i,1,k) * b(2,3,i,1,k) ) * tmp
            b(5,3,i,1,k) = (  b(5,3,i,1,k)
     $                      - b(5,1,i,1,k) * b(1,3,i,1,k)
     $                      - b(5,2,i,1,k) * b(2,3,i,1,k) ) * tmp
c
            b(2,4,i,1,k) =  b(2,4,i,1,k)
     $                    - b(2,1,i,1,k) * b(1,4,i,1,k)
            b(3,4,i,1,k) =  b(3,4,i,1,k)
     $                    - b(3,1,i,1,k) * b(1,4,i,1,k)
     $                    - b(3,2,i,1,k) * b(2,4,i,1,k)
            b(4,4,i,1,k) =  b(4,4,i,1,k)
     $                    - b(4,1,i,1,k) * b(1,4,i,1,k)
     $                    - b(4,2,i,1,k) * b(2,4,i,1,k)
     $                    - b(4,3,i,1,k) * b(3,4,i,1,k)
c
            b(5,4,i,1,k) = (  b(5,4,i,1,k)
     $                      - b(5,1,i,1,k) * b(1,4,i,1,k)
     $                      - b(5,2,i,1,k) * b(2,4,i,1,k)
     $                      - b(5,3,i,1,k) * b(3,4,i,1,k) )
     $                     / b(4,4,i,1,k)
c
            b(2,5,i,1,k) =  b(2,5,i,1,k)
     $                    - b(2,1,i,1,k) * b(1,5,i,1,k)
            b(3,5,i,1,k) =  b(3,5,i,1,k)
     $                    - b(3,1,i,1,k) * b(1,5,i,1,k)
     $                    - b(3,2,i,1,k) * b(2,5,i,1,k)
            b(4,5,i,1,k) =  b(4,5,i,1,k)
     $                    - b(4,1,i,1,k) * b(1,5,i,1,k)
     $                    - b(4,2,i,1,k) * b(2,5,i,1,k)
     $                    - b(4,3,i,1,k) * b(3,5,i,1,k)
c
            b(5,5,i,1,k) =  b(5,5,i,1,k)
     $                    - b(5,1,i,1,k) * b(1,5,i,1,k)
     $                    - b(5,2,i,1,k) * b(2,5,i,1,k)
     $                    - b(5,3,i,1,k) * b(3,5,i,1,k)
     $                    - b(5,4,i,1,k) * b(4,5,i,1,k)
c
            do j = 2, ny
c
               do m = 1, 5
c
                  tm(1,m) = c(1,m,i,j-1,k)
                  tm(2,m) = c(2,m,i,j-1,k) - b(2,1,i,j-1,k) * tm(1,m)
                  tm(3,m) = c(3,m,i,j-1,k) - b(3,1,i,j-1,k) * tm(1,m)
     $                                     - b(3,2,i,j-1,k) * tm(2,m)
                  tm(4,m) = c(4,m,i,j-1,k) - b(4,1,i,j-1,k) * tm(1,m)
     $                                     - b(4,2,i,j-1,k) * tm(2,m)
     $                                     - b(4,3,i,j-1,k) * tm(3,m)
                  tm(5,m) = c(5,m,i,j-1,k) - b(5,1,i,j-1,k) * tm(1,m)
     $                                     - b(5,2,i,j-1,k) * tm(2,m)
     $                                     - b(5,3,i,j-1,k) * tm(3,m)
     $                                     - b(5,4,i,j-1,k) * tm(4,m)
c
               end do
c
               do m = 1, 5
c
                  tm(5,m) = tm(5,m) / b(5,5,i,j-1,k)
                  tm(4,m) = (  tm(4,m)
     $                       - b(4,5,i,j-1,k) * tm(5,m) )
     $                      / b(4,4,i,j-1,k)
                  tm(3,m) = (  tm(3,m)
     $                       - b(3,4,i,j-1,k) * tm(4,m)
     $                       - b(3,5,i,j-1,k) * tm(5,m) )
     $                      / b(3,3,i,j-1,k)
                  tm(2,m) = (  tm(2,m)
     $                       - b(2,3,i,j-1,k) * tm(3,m)
     $                       - b(2,4,i,j-1,k) * tm(4,m)
     $                       - b(2,5,i,j-1,k) * tm(5,m) )
     $                      / b(2,2,i,j-1,k)
                  tm(1,m) = (  tm(1,m)
     $                       - b(1,2,i,j-1,k) * tm(2,m)
     $                       - b(1,3,i,j-1,k) * tm(3,m)
     $                       - b(1,4,i,j-1,k) * tm(4,m)
     $                       - b(1,5,i,j-1,k) * tm(5,m) )
     $                      / b(1,1,i,j-1,k)
c
               end do
c
               do m = 1, 5
c
                  b(1,m,i,j,k) =  b(1,m,i,j,k)
     $                          - a(1,1,i,j,k) * tm(1,m)
     $                          - a(1,2,i,j,k) * tm(2,m)
     $                          - a(1,3,i,j,k) * tm(3,m)
     $                          - a(1,4,i,j,k) * tm(4,m)
     $                          - a(1,5,i,j,k) * tm(5,m)
c
                  b(2,m,i,j,k) =  b(2,m,i,j,k)
     $                          - a(2,1,i,j,k) * tm(1,m)
     $                          - a(2,2,i,j,k) * tm(2,m)
     $                          - a(2,3,i,j,k) * tm(3,m)
     $                          - a(2,4,i,j,k) * tm(4,m)
     $                          - a(2,5,i,j,k) * tm(5,m)
c
                  b(3,m,i,j,k) =  b(3,m,i,j,k)
     $                          - a(3,1,i,j,k) * tm(1,m)
     $                          - a(3,2,i,j,k) * tm(2,m)
     $                          - a(3,3,i,j,k) * tm(3,m)
     $                          - a(3,4,i,j,k) * tm(4,m)
     $                          - a(3,5,i,j,k) * tm(5,m)
c
                  b(4,m,i,j,k) =  b(4,m,i,j,k)
     $                          - a(4,1,i,j,k) * tm(1,m)
     $                          - a(4,2,i,j,k) * tm(2,m)
     $                          - a(4,3,i,j,k) * tm(3,m)
     $                          - a(4,4,i,j,k) * tm(4,m)
     $                          - a(4,5,i,j,k) * tm(5,m)
c
                  b(5,m,i,j,k) =  b(5,m,i,j,k)
     $                          - a(5,1,i,j,k) * tm(1,m)
     $                          - a(5,2,i,j,k) * tm(2,m)
     $                          - a(5,3,i,j,k) * tm(3,m)
     $                          - a(5,4,i,j,k) * tm(4,m)
     $                          - a(5,5,i,j,k) * tm(5,m)
c
               end do
c
               tv(1) = f(1,i,j-1,k)
               tv(2) = f(2,i,j-1,k) - b(2,1,i,j-1,k) * tv(1)
               tv(3) = f(3,i,j-1,k) - b(3,1,i,j-1,k) * tv(1)
     $                              - b(3,2,i,j-1,k) * tv(2)
               tv(4) = f(4,i,j-1,k) - b(4,1,i,j-1,k) * tv(1)
     $                              - b(4,2,i,j-1,k) * tv(2)
     $                              - b(4,3,i,j-1,k) * tv(3)
               tv(5) = f(5,i,j-1,k) - b(5,1,i,j-1,k) * tv(1)
     $                              - b(5,2,i,j-1,k) * tv(2)
     $                              - b(5,3,i,j-1,k) * tv(3)
     $                              - b(5,4,i,j-1,k) * tv(4)
c
               tv(5) = tv(5) / b(5,5,i,j-1,k)
               tv(4) = (  tv(4)
     $                  - b(4,5,i,j-1,k) * tv(5) )
     $                 / b(4,4,i,j-1,k)
               tv(3) = (  tv(3)
     $                  - b(3,4,i,j-1,k) * tv(4)
     $                  - b(3,5,i,j-1,k) * tv(5) )
     $                 / b(3,3,i,j-1,k)
               tv(2) = (  tv(2)
     $                  - b(2,3,i,j-1,k) * tv(3)
     $                  - b(2,4,i,j-1,k) * tv(4)
     $                  - b(2,5,i,j-1,k) * tv(5) )
     $                 / b(2,2,i,j-1,k)
               tv(1) = (  tv(1)
     $                  - b(1,2,i,j-1,k) * tv(2)
     $                  - b(1,3,i,j-1,k) * tv(3)
     $                  - b(1,4,i,j-1,k) * tv(4)
     $                  - b(1,5,i,j-1,k) * tv(5) )
     $                 / b(1,1,i,j-1,k)
c
               do m = 1, 5
c
                  f(m,i,j,k) =  f(m,i,j,k)
     $                        - a(m,1,i,j,k) * tv(1)
     $                        - a(m,2,i,j,k) * tv(2)
     $                        - a(m,3,i,j,k) * tv(3)
     $                        - a(m,4,i,j,k) * tv(4)
     $                        - a(m,5,i,j,k) * tv(5)
c
               end do
c
               tmp = one / b(1,1,i,j,k)
c
               b(2,1,i,j,k) = b(2,1,i,j,k) * tmp
               b(3,1,i,j,k) = b(3,1,i,j,k) * tmp
               b(4,1,i,j,k) = b(4,1,i,j,k) * tmp
               b(5,1,i,j,k) = b(5,1,i,j,k) * tmp
c
               b(2,2,i,j,k) =  b(2,2,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,2,i,j,k)
c
               tmp = one / b(2,2,i,j,k)
c
               b(3,2,i,j,k) = (  b(3,2,i,j,k)
     $                         - b(3,1,i,j,k) * b(1,2,i,j,k) ) * tmp
               b(4,2,i,j,k) = (  b(4,2,i,j,k)
     $                         - b(4,1,i,j,k) * b(1,2,i,j,k) ) * tmp
               b(5,2,i,j,k) = (  b(5,2,i,j,k)
     $                         - b(5,1,i,j,k) * b(1,2,i,j,k) ) * tmp
c
               b(2,3,i,j,k) =  b(2,3,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,3,i,j,k)
               b(3,3,i,j,k) =  b(3,3,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,3,i,j,k)
c
               tmp = one / b(3,3,i,j,k)
c
               b(4,3,i,j,k) = (  b(4,3,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,3,i,j,k) ) * tmp
               b(5,3,i,j,k) = (  b(5,3,i,j,k)
     $                       - b(5,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(5,2,i,j,k) * b(2,3,i,j,k) ) * tmp
c
               b(2,4,i,j,k) =  b(2,4,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,4,i,j,k)
               b(3,4,i,j,k) =  b(3,4,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,4,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,4,i,j,k)
               b(4,4,i,j,k) =  b(4,4,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,4,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,4,i,j,k)
     $                       - b(4,3,i,j,k) * b(3,4,i,j,k)
c
               b(5,4,i,j,k) = (  b(5,4,i,j,k)
     $                         - b(5,1,i,j,k) * b(1,4,i,j,k)
     $                         - b(5,2,i,j,k) * b(2,4,i,j,k)
     $                         - b(5,3,i,j,k) * b(3,4,i,j,k) )
     $                        / b(4,4,i,j,k)
c
               b(2,5,i,j,k) =  b(2,5,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,5,i,j,k)
               b(3,5,i,j,k) =  b(3,5,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,5,i,j,k)
               b(4,5,i,j,k) =  b(4,5,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,5,i,j,k)
     $                       - b(4,3,i,j,k) * b(3,5,i,j,k)
c
               b(5,5,i,j,k) =  b(5,5,i,j,k)
     $                       - b(5,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(5,2,i,j,k) * b(2,5,i,j,k)
     $                       - b(5,3,i,j,k) * b(3,5,i,j,k)
     $                       - b(5,4,i,j,k) * b(4,5,i,j,k)
c
            end do
c     
         end do
c
      end do
c
c***back-substitution phase
c
      do k = 2, nz-1
c
         do i = 2, nx-1
c
            f(2,i,ny,k) =  f(2,i,ny,k)
     $                   - b(2,1,i,ny,k) * f(1,i,ny,k)  
            f(3,i,ny ,k) =  f(3,i,ny,k)
     $                   - b(3,1,i,ny,k) * f(1,i,ny,k)
     $                   - b(3,2,i,ny,k) * f(2,i,ny,k)
            f(4,i,ny,k) =  f(4,i,ny,k)
     $                   - b(4,1,i,ny,k) * f(1,i,ny,k)
     $                   - b(4,2,i,ny,k) * f(2,i,ny,k)
     $                   - b(4,3,i,ny,k) * f(3,i,ny,k)
            f(5,i,ny,k) =  f(5,i,ny,k)
     $                   - b(5,1,i,ny,k) * f(1,i,ny,k)
     $                   - b(5,2,i,ny,k) * f(2,i,ny,k)
     $                   - b(5,3,i,ny,k) * f(3,i,ny,k)
     $                   - b(5,4,i,ny,k) * f(4,i,ny,k)
c
            f(5,i,ny,k) = f(5,i,ny,k) / b(5,5,i,ny,k)
            f(4,i,ny,k) = (  f(4,i,ny,k)
     $                     - b(4,5,i,ny,k) * f(5,i,ny,k) )
     $                    / b(4,4,i,ny,k)
            f(3,i,ny,k) = (  f(3,i,ny,k)
     $                     - b(3,4,i,ny,k) * f(4,i,ny,k)
     $                     - b(3,5,i,ny,k) * f(5,i,ny,k) )
     $                    / b(3,3,i,ny,k)
            f(2,i,ny,k) = (  f(2,i,ny,k)
     $                     - b(2,3,i,ny,k) * f(3,i,ny,k)
     $                     - b(2,4,i,ny,k) * f(4,i,ny,k)
     $                     - b(2,5,i,ny,k) * f(5,i,ny,k) )
     $                    / b(2,2,i,ny,k)
            f(1,i,ny,k) = (  f(1,i,ny,k)
     $                     - b(1,2,i,ny,k) * f(2,i,ny,k)
     $                     - b(1,3,i,ny,k) * f(3,i,ny,k)
     $                     - b(1,4,i,ny,k) * f(4,i,ny,k)
     $                     - b(1,5,i,ny,k) * f(5,i,ny,k) )
     $                    / b(1,1,i,ny,k)
c
            do j = ny-1, 1, -1
c
               do m = 1, 5
c
                  f(m,i,j,k) =  f(m,i,j,k)
     $                        - c(m,1,i,j,k) * f(1,i,j+1,k)
     $                        - c(m,2,i,j,k) * f(2,i,j+1,k)
     $                        - c(m,3,i,j,k) * f(3,i,j+1,k)
     $                        - c(m,4,i,j,k) * f(4,i,j+1,k)
     $                        - c(m,5,i,j,k) * f(5,i,j+1,k)
c
               end do
c
               f(2,i,j,k) =  f(2,i,j,k)
     $                     - b(2,1,i,j,k) * f(1,i,j,k)  
               f(3,i,j,k) =  f(3,i,j,k)
     $                     - b(3,1,i,j,k) * f(1,i,j,k)
     $                     - b(3,2,i,j,k) * f(2,i,j,k)
               f(4,i,j,k) =  f(4,i,j,k)
     $                     - b(4,1,i,j,k) * f(1,i,j,k)
     $                     - b(4,2,i,j,k) * f(2,i,j,k)
     $                     - b(4,3,i,j,k) * f(3,i,j,k)
               f(5,i,j,k) =  f(5,i,j,k)
     $                     - b(5,1,i,j,k) * f(1,i,j,k)
     $                     - b(5,2,i,j,k) * f(2,i,j,k)
     $                     - b(5,3,i,j,k) * f(3,i,j,k)
     $                     - b(5,4,i,j,k) * f(4,i,j,k)
c
               f(5,i,j,k) = f(5,i,j,k) / b(5,5,i,j,k)
               f(4,i,j,k) = (  f(4,i,j,k)
     $                       - b(4,5,i,j,k) * f(5,i,j,k) )
     $                      / b(4,4,i,j,k)
               f(3,i,j,k) = (  f(3,i,j,k)
     $                       - b(3,4,i,j,k) * f(4,i,j,k)
     $                       - b(3,5,i,j,k) * f(5,i,j,k) )
     $                      / b(3,3,i,j,k)
               f(2,i,j,k) = (  f(2,i,j,k)
     $                       - b(2,3,i,j,k) * f(3,i,j,k)
     $                       - b(2,4,i,j,k) * f(4,i,j,k)
     $                       - b(2,5,i,j,k) * f(5,i,j,k) )
     $                      / b(2,2,i,j,k)
               f(1,i,j,k) = (  f(1,i,j,k)
     $                       - b(1,2,i,j,k) * f(2,i,j,k)
     $                       - b(1,3,i,j,k) * f(3,i,j,k)
     $                       - b(1,4,i,j,k) * f(4,i,j,k)
     $                       - b(1,5,i,j,k) * f(5,i,j,k) )
     $                      / b(1,1,i,j,k)
c               
            end do
c
         end do
c
      end do
c
      return
      end         
c
c
c
c
c
c
      subroutine btridz ( isiz1, isiz2, isiz3,
     $                    nx, ny, nz,
     $                    a, b, c, f )
c
c***solution of multiple, independent systems of block tridiagonal systems 
c   using Gaussian elimination (without pivoting) algorithm
c   (block size = 5 X 5)
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      implicit real*8 (a-h,o-z)
c
      parameter ( zero = 0.0d+00, one = 1.0d+00 )
c
      dimension a(5,5,isiz1,isiz2,*),
     $          b(5,5,isiz1,isiz2,*),
     $          c(5,5,isiz1,isiz2,*),
     $          f(5,isiz1,isiz2,*)
c
      dimension tm(5,5), tv(5)
c
c***forward elimination phase
c
      do j = 2, ny-1
c
         do i = 2, nx-1
c
            tmp = one / b(1,1,i,j,1)
c
            b(2,1,i,j,1) = b(2,1,i,j,1) * tmp
            b(3,1,i,j,1) = b(3,1,i,j,1) * tmp
            b(4,1,i,j,1) = b(4,1,i,j,1) * tmp
            b(5,1,i,j,1) = b(5,1,i,j,1) * tmp
c
            b(2,2,i,j,1) =  b(2,2,i,j,1)
     $                    - b(2,1,i,j,1) * b(1,2,i,j,1)
c
            tmp = one / b(2,2,i,j,1)
c
            b(3,2,i,j,1) = (  b(3,2,i,j,1)
     $                      - b(3,1,i,j,1) * b(1,2,i,j,1) ) * tmp
            b(4,2,i,j,1) = (  b(4,2,i,j,1)
     $                      - b(4,1,i,j,1) * b(1,2,i,j,1) ) * tmp
            b(5,2,i,j,1) = (  b(5,2,i,j,1)
     $                      - b(5,1,i,j,1) * b(1,2,i,j,1) ) * tmp
c
            b(2,3,i,j,1) =  b(2,3,i,j,1)
     $                    - b(2,1,i,j,1) * b(1,3,i,j,1)
            b(3,3,i,j,1) =  b(3,3,i,j,1)
     $                    - b(3,1,i,j,1) * b(1,3,i,j,1)
     $                    - b(3,2,i,j,1) * b(2,3,i,j,1)
c
            tmp = one / b(3,3,i,j,1)
c
            b(4,3,i,j,1) = (  b(4,3,i,j,1)
     $                      - b(4,1,i,j,1) * b(1,3,i,j,1)
     $                      - b(4,2,i,j,1) * b(2,3,i,j,1) ) * tmp
            b(5,3,i,j,1) = (  b(5,3,i,j,1)
     $                      - b(5,1,i,j,1) * b(1,3,i,j,1)
     $                      - b(5,2,i,j,1) * b(2,3,i,j,1) ) * tmp
c
            b(2,4,i,j,1) =  b(2,4,i,j,1)
     $                    - b(2,1,i,j,1) * b(1,4,i,j,1)
            b(3,4,i,j,1) =  b(3,4,i,j,1)
     $                    - b(3,1,i,j,1) * b(1,4,i,j,1)
     $                    - b(3,2,i,j,1) * b(2,4,i,j,1)
            b(4,4,i,j,1) =  b(4,4,i,j,1)
     $                    - b(4,1,i,j,1) * b(1,4,i,j,1)
     $                    - b(4,2,i,j,1) * b(2,4,i,j,1)
     $                    - b(4,3,i,j,1) * b(3,4,i,j,1)
c
            b(5,4,i,j,1) = (  b(5,4,i,j,1)
     $                      - b(5,1,i,j,1) * b(1,4,i,j,1)
     $                      - b(5,2,i,j,1) * b(2,4,i,j,1)
     $                      - b(5,3,i,j,1) * b(3,4,i,j,1) )
     $                     / b(4,4,i,j,1)
c
            b(2,5,i,j,1) =  b(2,5,i,j,1)
     $                    - b(2,1,i,j,1) * b(1,5,i,j,1)
            b(3,5,i,j,1) =  b(3,5,i,j,1)
     $                    - b(3,1,i,j,1) * b(1,5,i,j,1)
     $                    - b(3,2,i,j,1) * b(2,5,i,j,1)
            b(4,5,i,j,1) =  b(4,5,i,j,1)
     $                    - b(4,1,i,j,1) * b(1,5,i,j,1)
     $                    - b(4,2,i,j,1) * b(2,5,i,j,1)
     $                    - b(4,3,i,j,1) * b(3,5,i,j,1)
c
            b(5,5,i,j,1) =  b(5,5,i,j,1)
     $                    - b(5,1,i,j,1) * b(1,5,i,j,1)
     $                    - b(5,2,i,j,1) * b(2,5,i,j,1)
     $                    - b(5,3,i,j,1) * b(3,5,i,j,1)
     $                    - b(5,4,i,j,1) * b(4,5,i,j,1)
c
            do k = 2, nz
c
               do m = 1, 5
c
                  tm(1,m) = c(1,m,i,j,k-1)
                  tm(2,m) = c(2,m,i,j,k-1) - b(2,1,i,j,k-1) * tm(1,m)
                  tm(3,m) = c(3,m,i,j,k-1) - b(3,1,i,j,k-1) * tm(1,m)
     $                                     - b(3,2,i,j,k-1) * tm(2,m)
                  tm(4,m) = c(4,m,i,j,k-1) - b(4,1,i,j,k-1) * tm(1,m)
     $                                     - b(4,2,i,j,k-1) * tm(2,m)
     $                                     - b(4,3,i,j,k-1) * tm(3,m)
                  tm(5,m) = c(5,m,i,j,k-1) - b(5,1,i,j,k-1) * tm(1,m)
     $                                     - b(5,2,i,j,k-1) * tm(2,m)
     $                                     - b(5,3,i,j,k-1) * tm(3,m)
     $                                     - b(5,4,i,j,k-1) * tm(4,m)
c
               end do
c
               do m = 1, 5
c
                  tm(5,m) = tm(5,m) / b(5,5,i,j,k-1)
                  tm(4,m) = (  tm(4,m)
     $                       - b(4,5,i,j,k-1) * tm(5,m) )
     $                      / b(4,4,i,j,k-1)
                  tm(3,m) = (  tm(3,m)
     $                       - b(3,4,i,j,k-1) * tm(4,m)
     $                       - b(3,5,i,j,k-1) * tm(5,m) )
     $                      / b(3,3,i,j,k-1)
                  tm(2,m) = (  tm(2,m)
     $                       - b(2,3,i,j,k-1) * tm(3,m)
     $                       - b(2,4,i,j,k-1) * tm(4,m)
     $                       - b(2,5,i,j,k-1) * tm(5,m) )
     $                      / b(2,2,i,j,k-1)
                  tm(1,m) = (  tm(1,m)
     $                       - b(1,2,i,j,k-1) * tm(2,m)
     $                       - b(1,3,i,j,k-1) * tm(3,m)
     $                       - b(1,4,i,j,k-1) * tm(4,m)
     $                       - b(1,5,i,j,k-1) * tm(5,m) )
     $                      / b(1,1,i,j,k-1)
c
               end do
c
               do m = 1, 5
c
                  b(1,m,i,j,k) =  b(1,m,i,j,k)
     $                          - a(1,1,i,j,k) * tm(1,m)
     $                          - a(1,2,i,j,k) * tm(2,m)
     $                          - a(1,3,i,j,k) * tm(3,m)
     $                          - a(1,4,i,j,k) * tm(4,m)
     $                          - a(1,5,i,j,k) * tm(5,m)
c
                  b(2,m,i,j,k) =  b(2,m,i,j,k)
     $                          - a(2,1,i,j,k) * tm(1,m)
     $                          - a(2,2,i,j,k) * tm(2,m)
     $                          - a(2,3,i,j,k) * tm(3,m)
     $                          - a(2,4,i,j,k) * tm(4,m)
     $                          - a(2,5,i,j,k) * tm(5,m)
c
                  b(3,m,i,j,k) =  b(3,m,i,j,k)
     $                          - a(3,1,i,j,k) * tm(1,m)
     $                          - a(3,2,i,j,k) * tm(2,m)
     $                          - a(3,3,i,j,k) * tm(3,m)
     $                          - a(3,4,i,j,k) * tm(4,m)
     $                          - a(3,5,i,j,k) * tm(5,m)
c
                  b(4,m,i,j,k) =  b(4,m,i,j,k)
     $                          - a(4,1,i,j,k) * tm(1,m)
     $                          - a(4,2,i,j,k) * tm(2,m)
     $                          - a(4,3,i,j,k) * tm(3,m)
     $                          - a(4,4,i,j,k) * tm(4,m)
     $                          - a(4,5,i,j,k) * tm(5,m)
c
                  b(5,m,i,j,k) =  b(5,m,i,j,k)
     $                          - a(5,1,i,j,k) * tm(1,m)
     $                          - a(5,2,i,j,k) * tm(2,m)
     $                          - a(5,3,i,j,k) * tm(3,m)
     $                          - a(5,4,i,j,k) * tm(4,m)
     $                          - a(5,5,i,j,k) * tm(5,m)
c
               end do
c
               tv(1) = f(1,i,j,k-1)
               tv(2) = f(2,i,j,k-1) - b(2,1,i,j,k-1) * tv(1)
               tv(3) = f(3,i,j,k-1) - b(3,1,i,j,k-1) * tv(1)
     $                              - b(3,2,i,j,k-1) * tv(2)
               tv(4) = f(4,i,j,k-1) - b(4,1,i,j,k-1) * tv(1)
     $                              - b(4,2,i,j,k-1) * tv(2)
     $                              - b(4,3,i,j,k-1) * tv(3)
               tv(5) = f(5,i,j,k-1) - b(5,1,i,j,k-1) * tv(1)
     $                              - b(5,2,i,j,k-1) * tv(2)
     $                              - b(5,3,i,j,k-1) * tv(3)
     $                              - b(5,4,i,j,k-1) * tv(4)
c
               tv(5) = tv(5) / b(5,5,i,j,k-1)
               tv(4) = (  tv(4)
     $                  - b(4,5,i,j,k-1) * tv(5) )
     $                 / b(4,4,i,j,k-1)
               tv(3) = (  tv(3)
     $                  - b(3,4,i,j,k-1) * tv(4)
     $                  - b(3,5,i,j,k-1) * tv(5) )
     $                 / b(3,3,i,j,k-1)
               tv(2) = (  tv(2)
     $                  - b(2,3,i,j,k-1) * tv(3)
     $                  - b(2,4,i,j,k-1) * tv(4)
     $                  - b(2,5,i,j,k-1) * tv(5) )
     $                 / b(2,2,i,j,k-1)
               tv(1) = (  tv(1)
     $                  - b(1,2,i,j,k-1) * tv(2)
     $                  - b(1,3,i,j,k-1) * tv(3)
     $                  - b(1,4,i,j,k-1) * tv(4)
     $                  - b(1,5,i,j,k-1) * tv(5) )
     $                 / b(1,1,i,j,k-1)
c
               do m = 1, 5
c
                  f(m,i,j,k) =  f(m,i,j,k)
     $                        - a(m,1,i,j,k) * tv(1)
     $                        - a(m,2,i,j,k) * tv(2)
     $                        - a(m,3,i,j,k) * tv(3)
     $                        - a(m,4,i,j,k) * tv(4)
     $                        - a(m,5,i,j,k) * tv(5)
c
               end do
c
               tmp = one / b(1,1,i,j,k)
c
               b(2,1,i,j,k) = b(2,1,i,j,k) * tmp
               b(3,1,i,j,k) = b(3,1,i,j,k) * tmp
               b(4,1,i,j,k) = b(4,1,i,j,k) * tmp
               b(5,1,i,j,k) = b(5,1,i,j,k) * tmp
c
               b(2,2,i,j,k) =  b(2,2,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,2,i,j,k)
c
               tmp = one / b(2,2,i,j,k)
c
               b(3,2,i,j,k) = (  b(3,2,i,j,k)
     $                         - b(3,1,i,j,k) * b(1,2,i,j,k) ) * tmp
               b(4,2,i,j,k) = (  b(4,2,i,j,k)
     $                         - b(4,1,i,j,k) * b(1,2,i,j,k) ) * tmp
               b(5,2,i,j,k) = (  b(5,2,i,j,k)
     $                         - b(5,1,i,j,k) * b(1,2,i,j,k) ) * tmp
c
               b(2,3,i,j,k) =  b(2,3,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,3,i,j,k)
               b(3,3,i,j,k) =  b(3,3,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,3,i,j,k)
c
               tmp = one / b(3,3,i,j,k)
c
               b(4,3,i,j,k) = (  b(4,3,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,3,i,j,k) ) * tmp
               b(5,3,i,j,k) = (  b(5,3,i,j,k)
     $                       - b(5,1,i,j,k) * b(1,3,i,j,k)
     $                       - b(5,2,i,j,k) * b(2,3,i,j,k) ) * tmp
c
               b(2,4,i,j,k) =  b(2,4,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,4,i,j,k)
               b(3,4,i,j,k) =  b(3,4,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,4,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,4,i,j,k)
               b(4,4,i,j,k) =  b(4,4,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,4,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,4,i,j,k)
     $                       - b(4,3,i,j,k) * b(3,4,i,j,k)
c
               b(5,4,i,j,k) = (  b(5,4,i,j,k)
     $                         - b(5,1,i,j,k) * b(1,4,i,j,k)
     $                         - b(5,2,i,j,k) * b(2,4,i,j,k)
     $                         - b(5,3,i,j,k) * b(3,4,i,j,k) )
     $                        / b(4,4,i,j,k)
c
               b(2,5,i,j,k) =  b(2,5,i,j,k)
     $                       - b(2,1,i,j,k) * b(1,5,i,j,k)
               b(3,5,i,j,k) =  b(3,5,i,j,k)
     $                       - b(3,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(3,2,i,j,k) * b(2,5,i,j,k)
               b(4,5,i,j,k) =  b(4,5,i,j,k)
     $                       - b(4,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(4,2,i,j,k) * b(2,5,i,j,k)
     $                       - b(4,3,i,j,k) * b(3,5,i,j,k)
c
               b(5,5,i,j,k) =  b(5,5,i,j,k)
     $                       - b(5,1,i,j,k) * b(1,5,i,j,k)
     $                       - b(5,2,i,j,k) * b(2,5,i,j,k)
     $                       - b(5,3,i,j,k) * b(3,5,i,j,k)
     $                       - b(5,4,i,j,k) * b(4,5,i,j,k)
c
            end do
c     
         end do
c
      end do
c
c***back-substitution phase
c
      do j = 2, ny-1
c
         do i = 2, nx-1
c
            f(2,i,j,nz) =  f(2,i,j,nz)
     $                   - b(2,1,i,j,nz) * f(1,i,j,nz)  
            f(3,i,j,nz) =  f(3,i,j,nz)
     $                   - b(3,1,i,j,nz) * f(1,i,j,nz)
     $                   - b(3,2,i,j,nz) * f(2,i,j,nz)
            f(4,i,j,nz) =  f(4,i,j,nz)
     $                   - b(4,1,i,j,nz) * f(1,i,j,nz)
     $                   - b(4,2,i,j,nz) * f(2,i,j,nz)
     $                   - b(4,3,i,j,nz) * f(3,i,j,nz)
            f(5,i,j,nz) =  f(5,i,j,nz)
     $                   - b(5,1,i,j,nz) * f(1,i,j,nz)
     $                   - b(5,2,i,j,nz) * f(2,i,j,nz)
     $                   - b(5,3,i,j,nz) * f(3,i,j,nz)
     $                   - b(5,4,i,j,nz) * f(4,i,j,nz)
c
            f(5,i,j,nz) = f(5,i,j,nz) / b(5,5,i,j,nz)
            f(4,i,j,nz) = (  f(4,i,j,nz)
     $                     - b(4,5,i,j,nz) * f(5,i,j,nz) )
     $                    / b(4,4,i,j,nz)
            f(3,i,j,nz) = (  f(3,i,j,nz)
     $                     - b(3,4,i,j,nz) * f(4,i,j,nz)
     $                     - b(3,5,i,j,nz) * f(5,i,j,nz) )
     $                    / b(3,3,i,j,nz)
            f(2,i,j,nz) = (  f(2,i,j,nz)
     $                     - b(2,3,i,j,nz) * f(3,i,j,nz)
     $                     - b(2,4,i,j,nz) * f(4,i,j,nz)
     $                     - b(2,5,i,j,nz) * f(5,i,j,nz) )
     $                    / b(2,2,i,j,nz)
            f(1,i,j,nz) = (  f(1,i,j,nz)
     $                     - b(1,2,i,j,nz) * f(2,i,j,nz)
     $                     - b(1,3,i,j,nz) * f(3,i,j,nz)
     $                     - b(1,4,i,j,nz) * f(4,i,j,nz)
     $                     - b(1,5,i,j,nz) * f(5,i,j,nz) )
     $                    / b(1,1,i,j,nz)
c
            do k = nz-1, 1, -1
c
               do m = 1, 5
c
                  f(m,i,j,k) =  f(m,i,j,k)
     $                        - c(m,1,i,j,k) * f(1,i,j,k+1)
     $                        - c(m,2,i,j,k) * f(2,i,j,k+1)
     $                        - c(m,3,i,j,k) * f(3,i,j,k+1)
     $                        - c(m,4,i,j,k) * f(4,i,j,k+1)
     $                        - c(m,5,i,j,k) * f(5,i,j,k+1)
c
               end do
c
               f(2,i,j,k) =  f(2,i,j,k)
     $                     - b(2,1,i,j,k) * f(1,i,j,k)  
               f(3,i,j,k) =  f(3,i,j,k)
     $                     - b(3,1,i,j,k) * f(1,i,j,k)
     $                     - b(3,2,i,j,k) * f(2,i,j,k)
               f(4,i,j,k) =  f(4,i,j,k)
     $                     - b(4,1,i,j,k) * f(1,i,j,k)
     $                     - b(4,2,i,j,k) * f(2,i,j,k)
     $                     - b(4,3,i,j,k) * f(3,i,j,k)
               f(5,i,j,k) =  f(5,i,j,k)
     $                     - b(5,1,i,j,k) * f(1,i,j,k)
     $                     - b(5,2,i,j,k) * f(2,i,j,k)
     $                     - b(5,3,i,j,k) * f(3,i,j,k)
     $                     - b(5,4,i,j,k) * f(4,i,j,k)
c
               f(5,i,j,k) = f(5,i,j,k) / b(5,5,i,j,k)
               f(4,i,j,k) = (  f(4,i,j,k)
     $                       - b(4,5,i,j,k) * f(5,i,j,k) )
     $                      / b(4,4,i,j,k)
               f(3,i,j,k) = (  f(3,i,j,k)
     $                       - b(3,4,i,j,k) * f(4,i,j,k)
     $                       - b(3,5,i,j,k) * f(5,i,j,k) )
     $                      / b(3,3,i,j,k)
               f(2,i,j,k) = (  f(2,i,j,k)
     $                       - b(2,3,i,j,k) * f(3,i,j,k)
     $                       - b(2,4,i,j,k) * f(4,i,j,k)
     $                       - b(2,5,i,j,k) * f(5,i,j,k) )
     $                      / b(2,2,i,j,k)
               f(1,i,j,k) = (  f(1,i,j,k)
     $                       - b(1,2,i,j,k) * f(2,i,j,k)
     $                       - b(1,3,i,j,k) * f(3,i,j,k)
     $                       - b(1,4,i,j,k) * f(4,i,j,k)
     $                       - b(1,5,i,j,k) * f(5,i,j,k) )
     $                      / b(1,1,i,j,k)
c               
            end do
c
         end do
c
      end do
c
      return
      end         
c
c
c
c
c
      subroutine erhs
c
c***compute the right hand side based on exact solution
c
c Author: Sisira Weeratunga
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension flux(5,isiz1), ue(5,isiz1)
c
      dsspm = dssp
c
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  frct( m, i, j, k ) = 0.0d+00
               end do
            end do
         end do
      end do
c
c***xi-direction flux differences
c
      do k = 2, nz - 1
c
         zeta = ( dfloat(k-1) ) / ( nz - 1 )
c
         do j = 2, ny - 1
c
            eta = ( dfloat(j-1) ) / ( ny - 1 )
c
            do i = 1, nx
c
               xi = ( dfloat(i-1) ) / ( nx - 1 )
c
               do m = 1, 5
c
                  ue(m,i) =  ce(m,1)
     $                 + ce(m,2) * xi
     $                 + ce(m,3) * eta
     $                 + ce(m,4) * zeta
     $                 + ce(m,5) * xi * xi
     $                 + ce(m,6) * eta * eta
     $                 + ce(m,7) * zeta * zeta
     $                 + ce(m,8) * xi * xi * xi
     $                 + ce(m,9) * eta * eta * eta
     $                 + ce(m,10) * zeta * zeta * zeta
     $                 + ce(m,11) * xi * xi * xi * xi
     $                 + ce(m,12) * eta * eta * eta * eta
     $                 + ce(m,13) * zeta * zeta * zeta * zeta
c
               end do
c
               flux(1,i) = ue(2,i)
c
               u21 = ue(2,i) / ue(1,i)
c
               q = 0.50d+00 * (  ue(2,i) * ue(2,i)
     $                         + ue(3,i) * ue(3,i)
     $                         + ue(4,i) * ue(4,i) )
     $                      / ue(1,i)
c
               flux(2,i) = ue(2,i) * u21 + c2 * ( ue(5,i) - q )
c
               flux(3,i) = ue(3,i) * u21
c
               flux(4,i) = ue(4,i) * u21
c
               flux(5,i) = ( c1 * ue(5,i) - c2 * q ) * u21
c    
            end do
c
            do i = 2, nx - 1
c
               do m = 1, 5
c
                  frct(m,i,j,k) =  frct(m,i,j,k)
     $                       - tx2 * ( flux(m,i+1) - flux(m,i-1) )
c
               end do
c
            end do
c
            do i = 2, nx
c
               tmp = 1.0d+00 / ue(1,i)
c
               u21i = tmp * ue(2,i)
               u31i = tmp * ue(3,i)
               u41i = tmp * ue(4,i)
               u51i = tmp * ue(5,i)
c
               tmp = 1.0d+00 / ue(1,i-1)
c
               u21im1 = tmp * ue(2,i-1)
               u31im1 = tmp * ue(3,i-1)
               u41im1 = tmp * ue(4,i-1)
               u51im1 = tmp * ue(5,i-1)
c
               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * ( u21i - u21im1 )
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     $                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * tx3 * ( u21i**2 - u21im1**2 )
     $              + c1 * c5 * tx3 * ( u51i - u51im1 )
c
            end do
c
            do i = 2, nx - 1
c
               frct(1,i,j,k) = frct(1,i,j,k)
     $              + dx1 * tx1 * (            ue(1,i-1)
     $                             - 2.0d+00 * ue(1,i)
     $                             +           ue(1,i+1) )
c
               frct(2,i,j,k) = frct(2,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     $              + dx2 * tx1 * (            ue(2,i-1)
     $                             - 2.0d+00 * ue(2,i)
     $                             +           ue(2,i+1) )
c
               frct(3,i,j,k) = frct(3,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     $              + dx3 * tx1 * (            ue(3,i-1)
     $                             - 2.0d+00 * ue(3,i)
     $                             +           ue(3,i+1) )
c
               frct(4,i,j,k) = frct(4,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     $              + dx4 * tx1 * (            ue(4,i-1)
     $                             - 2.0d+00 * ue(4,i)
     $                             +           ue(4,i+1) )
c
               frct(5,i,j,k) = frct(5,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     $              + dx5 * tx1 * (            ue(5,i-1)
     $                             - 2.0d+00 * ue(5,i)
     $                             +           ue(5,i+1) )
c
            end do
c
c***Fourth-order dissipation
c
            do m = 1, 5
c
               frct(m,2,j,k) = frct(m,2,j,k)
     $           - dsspm * ( + 5.0d+00 * ue(m,2)
     $                       - 4.0d+00 * ue(m,3)
     $                       +           ue(m,4) )
c
               frct(m,3,j,k) = frct(m,3,j,k)
     $           - dsspm * ( - 4.0d+00 * ue(m,2)
     $                       + 6.0d+00 * ue(m,3)
     $                       - 4.0d+00 * ue(m,4)
     $                       +           ue(m,5) )
c
            end do
c
            do i = 4, nx - 3
c
               do m = 1, 5
c
                  frct(m,i,j,k) = frct(m,i,j,k)
     $              - dsspm * (            ue(m,i-2)
     $                         - 4.0d+00 * ue(m,i-1)
     $                         + 6.0d+00 * ue(m,i)
     $                         - 4.0d+00 * ue(m,i+1)
     $                         +           ue(m,i+2) )
c
               end do
c
            end do
c
            do m = 1, 5
c
               frct(m,nx-2,j,k) = frct(m,nx-2,j,k)
     $           - dsspm * (             ue(m,nx-4)
     $                       - 4.0d+00 * ue(m,nx-3)
     $                       + 6.0d+00 * ue(m,nx-2)
     $                       - 4.0d+00 * ue(m,nx-1)  )
c
               frct(m,nx-1,j,k) = frct(m,nx-1,j,k)
     $           - dsspm * (             ue(m,nx-3)
     $                       - 4.0d+00 * ue(m,nx-2)
     $                       + 5.0d+00 * ue(m,nx-1) )
c
            end do
c
         end do
c
      end do
c
c***eta-direction flux differences
c
      do k = 2, nz - 1
c
         zeta = ( dfloat(k-1) ) / ( nz - 1 )
c
         do i = 2, nx - 1
c
            xi = ( dfloat(i-1) ) / ( nx - 1 )
c 
            do j = 1, ny
c
               eta = ( dfloat(j-1) ) / ( ny - 1 )
c
               do m = 1, 5
c
                  ue(m,j) =  ce(m,1)
     $                 + ce(m,2) * xi
     $                 + ce(m,3) * eta
     $                 + ce(m,4) * zeta
     $                 + ce(m,5) * xi * xi
     $                 + ce(m,6) * eta * eta
     $                 + ce(m,7) * zeta * zeta
     $                 + ce(m,8) * xi * xi * xi
     $                 + ce(m,9) * eta * eta * eta
     $                 + ce(m,10) * zeta * zeta * zeta
     $                 + ce(m,11) * xi * xi * xi * xi
     $                 + ce(m,12) * eta * eta * eta * eta
     $                 + ce(m,13) * zeta * zeta * zeta * zeta
c
               end do
c
               flux(1,j) = ue(3,j)
c
               u31 = ue(3,j) / ue(1,j)
c
               q = 0.50d+00 * (  ue(2,j) * ue(2,j)
     $                         + ue(3,j) * ue(3,j)
     $                         + ue(4,j) * ue(4,j) )
     $                      / ue(1,j)
c
               flux(2,j) = ue(2,j) * u31 
c
               flux(3,j) = ue(3,j) * u31 + c2 * ( ue(5,j) - q )
c
               flux(4,j) = ue(4,j) * u31
c
               flux(5,j) = ( c1 * ue(5,j) - c2 * q ) * u31
c    
            end do
c
            do j = 2, ny - 1
c
               do m = 1, 5
c
                  frct(m,i,j,k) =  frct(m,i,j,k)
     $                       - ty2 * ( flux(m,j+1) - flux(m,j-1) )
c
               end do
c
            end do
c
            do j = 2, ny
c
               tmp = 1.0d+00 / ue(1,j)
c
               u21j = tmp * ue(2,j)
               u31j = tmp * ue(3,j)
               u41j = tmp * ue(4,j)
               u51j = tmp * ue(5,j)
c
               tmp = 1.0d+00 / ue(1,j-1)
c
               u21jm1 = tmp * ue(2,j-1)
               u31jm1 = tmp * ue(3,j-1)
               u41jm1 = tmp * ue(4,j-1)
               u51jm1 = tmp * ue(5,j-1)
c
               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * ( u31j - u31jm1 )
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     $                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * ty3 * ( u31j**2 - u31jm1**2 )
     $              + c1 * c5 * ty3 * ( u51j - u51jm1 )
c
            end do
c
            do j = 2, ny - 1
c
               frct(1,i,j,k) = frct(1,i,j,k)
     $              + dy1 * ty1 * (            ue(1,j-1)
     $                             - 2.0d+00 * ue(1,j)
     $                             +           ue(1,j+1) )
c
               frct(2,i,j,k) = frct(2,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     $              + dy2 * ty1 * (            ue(2,j-1)
     $                             - 2.0d+00 * ue(2,j)
     $                             +           ue(2,j+1) )
c
               frct(3,i,j,k) = frct(3,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     $              + dy3 * ty1 * (            ue(3,j-1)
     $                             - 2.0d+00 * ue(3,j)
     $                             +           ue(3,j+1) )
c
               frct(4,i,j,k) = frct(4,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     $              + dy4 * ty1 * (            ue(4,j-1)
     $                             - 2.0d+00 * ue(4,j)
     $                             +           ue(4,j+1) )
c
               frct(5,i,j,k) = frct(5,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     $              + dy5 * ty1 * (            ue(5,j-1)
     $                             - 2.0d+00 * ue(5,j)
     $                             +           ue(5,j+1) )
c
            end do
c
c***fourth-order dissipation
c
            do m = 1, 5
c
               frct(m,i,2,k) = frct(m,i,2,k)
     $           - dsspm * ( + 5.0d+00 * ue(m,2)
     $                       - 4.0d+00 * ue(m,3)
     $                       +           ue(m,4) )
c
               frct(m,i,3,k) = frct(m,i,3,k)
     $           - dsspm * ( - 4.0d+00 * ue(m,2)
     $                       + 6.0d+00 * ue(m,3)
     $                       - 4.0d+00 * ue(m,4)
     $                       +           ue(m,5) )
c
            end do
c
            do j = 4, ny - 3
c
               do m = 1, 5
c
                  frct(m,i,j,k) = frct(m,i,j,k)
     $              - dsspm * (            ue(m,j-2)
     $                        - 4.0d+00 * ue(m,j-1)
     $                        + 6.0d+00 * ue(m,j)
     $                        - 4.0d+00 * ue(m,j+1)
     $                        +           ue(m,j+2) )
c
               end do
c
            end do
c
            do m = 1, 5
c
               frct(m,i,ny-2,k) = frct(m,i,ny-2,k)
     $           - dsspm * (             ue(m,ny-4)
     $                       - 4.0d+00 * ue(m,ny-3)
     $                       + 6.0d+00 * ue(m,ny-2)
     $                       - 4.0d+00 * ue(m,ny-1)  )
c
               frct(m,i,ny-1,k) = frct(m,i,ny-1,k)
     $           - dsspm * (             ue(m,ny-3)
     $                       - 4.0d+00 * ue(m,ny-2)
     $                       + 5.0d+00 * ue(m,ny-1)  )
c
            end do
c
         end do
c
      end do
c
c***zeta-direction flux differences
c
      do j = 2, ny - 1
c
         eta = ( dfloat(j-1) ) / ( ny - 1 )
c
         do i = 2, nx - 1
c
            xi = ( dfloat(i-1) ) / ( nx - 1 )
c
            do k = 1, nz
c
               zeta = ( dfloat(k-1) ) / ( nz - 1 )
c
               do m = 1, 5
c
                  ue(m,k) =  ce(m,1)
     $                 + ce(m,2) * xi
     $                 + ce(m,3) * eta
     $                 + ce(m,4) * zeta
     $                 + ce(m,5) * xi * xi
     $                 + ce(m,6) * eta * eta
     $                 + ce(m,7) * zeta * zeta
     $                 + ce(m,8) * xi * xi * xi
     $                 + ce(m,9) * eta * eta * eta
     $                 + ce(m,10) * zeta * zeta * zeta
     $                 + ce(m,11) * xi * xi * xi * xi
     $                 + ce(m,12) * eta * eta * eta * eta
     $                 + ce(m,13) * zeta * zeta * zeta * zeta
c
               end do
c
               flux(1,k) = ue(4,k)
c
               u41 = ue(4,k) / ue(1,k)
c
               q = 0.50d+00 * (  ue(2,k) * ue(2,k)
     $                         + ue(3,k) * ue(3,k)
     $                         + ue(4,k) * ue(4,k) )
     $                      / ue(1,k)
c
               flux(2,k) = ue(2,k) * u41 
c
               flux(3,k) = ue(3,k) * u41 
c
               flux(4,k) = ue(4,k) * u41 + c2 * ( ue(5,k) - q )
c
               flux(5,k) = ( c1 * ue(5,k) - c2 * q ) * u41
c    
            end do
c
            do k = 2, nz - 1
c
               do m = 1, 5
c
                  frct(m,i,j,k) =  frct(m,i,j,k)
     $                       - tz2 * ( flux(m,k+1) - flux(m,k-1) )
c
               end do
c
            end do
c
            do k = 2, nz
c
               tmp = 1.0d+00 / ue(1,k)
c
               u21k = tmp * ue(2,k)
               u31k = tmp * ue(3,k)
               u41k = tmp * ue(4,k)
               u51k = tmp * ue(5,k)
c
               tmp = 1.0d+00 / ue(1,k-1)
c
               u21km1 = tmp * ue(2,k-1)
               u31km1 = tmp * ue(3,k-1)
               u41km1 = tmp * ue(4,k-1)
               u51km1 = tmp * ue(5,k-1)
c
               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * ( u41k - u41km1 )
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     $                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * tz3 * ( u41k**2 - u41km1**2 )
     $              + c1 * c5 * tz3 * ( u51k - u51km1 )
c
            end do
c
            do k = 2, nz - 1
c
               frct(1,i,j,k) = frct(1,i,j,k)
     $              + dz1 * tz1 * (            ue(1,k+1)
     $                             - 2.0d+00 * ue(1,k)
     $                             +           ue(1,k-1) )
c
               frct(2,i,j,k) = frct(2,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     $              + dz2 * tz1 * (            ue(2,k+1)
     $                             - 2.0d+00 * ue(2,k)
     $                             +           ue(2,k-1) )
c
               frct(3,i,j,k) = frct(3,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     $              + dz3 * tz1 * (            ue(3,k+1)
     $                             - 2.0d+00 * ue(3,k)
     $                             +           ue(3,k-1) )
c
               frct(4,i,j,k) = frct(4,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     $              + dz4 * tz1 * (            ue(4,k+1)
     $                             - 2.0d+00 * ue(4,k)
     $                             +           ue(4,k-1) )
c
               frct(5,i,j,k) = frct(5,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     $              + dz5 * tz1 * (            ue(5,k+1)
     $                             - 2.0d+00 * ue(5,k)
     $                             +           ue(5,k-1) )
c
            end do
c
c***fourth-order dissipation
c
            do m = 1, 5
c
               frct(m,i,j,2) = frct(m,i,j,2)
     $           - dsspm * ( + 5.0d+00 * ue(m,2)
     $                       - 4.0d+00 * ue(m,3)
     $                       +           ue(m,4) )
c
               frct(m,i,j,3) = frct(m,i,j,3)
     $           - dsspm * (- 4.0d+00 * ue(m,2)
     $                      + 6.0d+00 * ue(m,3)
     $                      - 4.0d+00 * ue(m,4)
     $                      +           ue(m,5) )
c
            end do
c
            do k = 4, nz - 3
c
               do m = 1, 5
c
                  frct(m,i,j,k) = frct(m,i,j,k)
     $              - dsspm * (           ue(m,k-2)
     $                        - 4.0d+00 * ue(m,k-1)
     $                        + 6.0d+00 * ue(m,k)
     $                        - 4.0d+00 * ue(m,k+1)
     $                        +           ue(m,k+2) )
c
               end do
c
            end do
c
            do m = 1, 5
c
               frct(m,i,j,nz-2) = frct(m,i,j,nz-2)
     $           - dsspm * (            ue(m,nz-4)
     $                      - 4.0d+00 * ue(m,nz-3)
     $                      + 6.0d+00 * ue(m,nz-2)
     $                      - 4.0d+00 * ue(m,nz-1)  )
c
               frct(m,i,j,nz-1) = frct(m,i,j,nz-1)
     $           - dsspm * (             ue(m,nz-3)
     $                       - 4.0d+00 * ue(m,nz-2)
     $                       + 5.0d+00 * ue(m,nz-1)  )
c
            end do
c
         end do
c
      end do
      return
      end
c
c
c
c
c
      subroutine exact ( i, j, k, u000ijk )
c
c***compute the exact solution at (i,j,k)
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension u000ijk(*)
c
      xi  = ( dfloat ( i - 1 ) ) / ( nx - 1 )
      eta  = ( dfloat ( j - 1 ) ) / ( ny - 1 )
      zeta = ( dfloat ( k - 1 ) ) / ( nz - 1 )
c
      do m = 1, 5
c
         u000ijk(m) =  ce(m,1)
     $        + ce(m,2) * xi
     $        + ce(m,3) * eta
     $        + ce(m,4) * zeta
     $        + ce(m,5) * xi * xi
     $        + ce(m,6) * eta * eta
     $        + ce(m,7) * zeta * zeta
     $        + ce(m,8) * xi * xi * xi
     $        + ce(m,9) * eta * eta * eta
     $        + ce(m,10) * zeta * zeta * zeta
     $        + ce(m,11) * xi * xi * xi * xi
     $        + ce(m,12) * eta * eta * eta * eta
     $        + ce(m,13) * zeta * zeta * zeta * zeta
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine error
c
c***compute the solution error
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension imax(5), jmax(5), kmax(5),
     $          u000ijk(5), errmax(5)
c
      lnorm = 2
c
      if ( lnorm .eq. 1 ) then
c
         do m = 1, 5
            errmax(m) = - 1.0d+20
         end do
c
         do k = 2, nz-1
            do j = 2, ny-1
               do i = 2, nx-1
c
                  call exact ( i, j, k, u000ijk )
c
                  do m = 1, 5
c
                     tmp = abs ( u000ijk(m) - u(m,i,j,k) )
c
                     if ( tmp .gt. errmax(m) ) then
c
                        errmax(m) = tmp
                        imax(m) = i
                        jmax(m) = j
                        kmax(m) = k
c
                     end if
c
                  end do
c
               end do
            end do
         end do
c
         write (iout,1001) ( errmax(m),
     $        imax(m), jmax(m), kmax(m), m = 1, 5 )
c
      else if ( lnorm .eq. 2 ) then
c
         do m = 1, 5
            errnm(m) = 0.0d+00
         end do
c
         do k = 2, nz-1
            do j = 2, ny-1
               do i = 2, nx-1
c
                  call exact ( i, j, k, u000ijk )
c
                  do m = 1, 5
c
                     tmp = ( u000ijk(m) - u(m,i,j,k) )
c
                     errnm(m) = errnm(m) + tmp ** 2
c
                  end do
c
               end do
            end do
         end do
c
         do m = 1, 5
            errnm(m) = sqrt ( errnm(m) / ( (nx-2)*(ny-2)*(nz-2) ) )
         end do
c
         write (iout,1002) ( errnm(m), m = 1, 5 )
c
      end if
c
 1001 format (/5x,'max. error in soln. to first pde  =',1pe12.4/,
     $5x,'and its location                  = (',i4','i4,','i4,' )'/,
     $/5x,'max. error in soln. to second pde =',1pe12.4/,
     $5x,'and its location                  = (',i4','i4,','i4,' )'/,
     $/5x,'max. error in soln. to third pde  =',1pe12.4/,
     $5x,'and its location                  = (',i4','i4,','i4,' )'/,
     $/5x,'max. error in soln. to fourth pde =',1pe12.4/,
     $5x,'and its location                  = (',i4','i4,','i4,' )'/,
     $/5x,'max. error in soln. to fifth pde  =',1pe12.4/,
     $5x,'and its location                  = (',i4','i4,','i4,' )' )
c
 1002 format (1x/1x,'RMS-norm of error in soln. to ',
     $ 'first pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of error in soln. to ',
     $ 'second pde = ',1pe12.5/,
     $ 1x,'RMS-norm of error in soln. to ',
     $ 'third pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of error in soln. to ',
     $ 'fourth pde = ',1pe12.5/,
     $ 1x,'RMS-norm of error in soln. to ',
     $ 'fifth pde  = ',1pe12.5)
c
      return
      end
c
c
c
c
c
      subroutine jacx
c
c***form the block tridiagonal system for xi-direction sweep
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      real*8 njac
c
      dimension fjac(5,5,isiz1), njac(5,5,isiz1)
c
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c
      do k = 2, nz - 1
c
         do j = 2, ny - 1
c
c***compute the xi-direction flux jacobian
c
            do i = 1, nx
c
               fjac(1,1,i) = 0.0d+00
               fjac(1,2,i) = 1.0d+00
               fjac(1,3,i) = 0.0d+00
               fjac(1,4,i) = 0.0d+00
               fjac(1,5,i) = 0.0d+00
c
               fjac(2,1,i) = - ( u(2,i,j,k) / u(1,i,j,k) ) ** 2
     $           + c2 * 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
     $                              + u(3,i,j,k) * u(3,i,j,k)
     $                              + u(4,i,j,k) * u(4,i,j,k) )
     $                           / ( u(1,i,j,k) * u(1,i,j,k) )
               fjac(2,2,i) = ( 2.0d+00 - c2 )
     $                     * ( u(2,i,j,k) / u(1,i,j,k) )
               fjac(2,3,i) = - c2 * ( u(3,i,j,k) / u(1,i,j,k) )
               fjac(2,4,i) = - c2 * ( u(4,i,j,k) / u(1,i,j,k) )
               fjac(2,5,i) = c2
c
               fjac(3,1,i) = - ( u(2,i,j,k)*u(3,i,j,k) )
     $                       / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(3,2,i) = u(3,i,j,k) / u(1,i,j,k)
               fjac(3,3,i) = u(2,i,j,k) / u(1,i,j,k)
               fjac(3,4,i) = 0.0d+00
               fjac(3,5,i) = 0.0d+00
c
               fjac(4,1,i) = - ( u(2,i,j,k)*u(4,i,j,k) )
     $                       / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(4,2,i) = u(4,i,j,k) / u(1,i,j,k)
               fjac(4,3,i) = 0.0d+00
               fjac(4,4,i) = u(2,i,j,k) / u(1,i,j,k)
               fjac(4,5,i) = 0.0d+00
c
               fjac(5,1,i) = ( c2 * (  u(2,i,j,k) * u(2,i,j,k)
     $                               + u(3,i,j,k) * u(3,i,j,k)
     $                               + u(4,i,j,k) * u(4,i,j,k) )
     $                            / ( u(1,i,j,k) * u(1,i,j,k) )
     $                         - c1 * ( u(5,i,j,k) / u(1,i,j,k) ) )
     $                        * ( u(2,i,j,k) / u(1,i,j,k) )
               fjac(5,2,i) = c1 * ( u(5,i,j,k) / u(1,i,j,k) )
     $                      - 0.50d+00 * c2
     $                      * ( (  3.0d+00*u(2,i,j,k)*u(2,i,j,k)
     $                           + u(3,i,j,k)*u(3,i,j,k)
     $                           + u(4,i,j,k)*u(4,i,j,k) )
     $                          / ( u(1,i,j,k)*u(1,i,j,k) ) )
               fjac(5,3,i) = - c2 * ( u(3,i,j,k)*u(2,i,j,k) )
     $              / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(5,4,i) = - c2 * ( u(4,i,j,k)*u(2,i,j,k) )
     $              / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(5,5,i) = c1 * ( u(2,i,j,k) / u(1,i,j,k) )
c
               tmp1 = 1.0d+00 / u(1,i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2
c
               njac(1,1,i) = 0.0d+00
               njac(1,2,i) = 0.0d+00
               njac(1,3,i) = 0.0d+00
               njac(1,4,i) = 0.0d+00
               njac(1,5,i) = 0.0d+00
c
               njac(2,1,i) = - r43 * c34 * tmp2 * u(2,i,j,k)
               njac(2,2,i) =   r43 * c34 * tmp1
               njac(2,3,i) =   0.0d+00
               njac(2,4,i) =   0.0d+00
               njac(2,5,i) =   0.0d+00
c     
               njac(3,1,i) = - c34 * tmp2 * u(3,i,j,k)
               njac(3,2,i) =   0.0d+00
               njac(3,3,i) =   c34 * tmp1
               njac(3,4,i) =   0.0d+00
               njac(3,5,i) =   0.0d+00
c
               njac(4,1,i) = - c34 * tmp2 * u(4,i,j,k)
               njac(4,2,i) =   0.0d+00 
               njac(4,3,i) =   0.0d+00
               njac(4,4,i) =   c34 * tmp1
               njac(4,5,i) =   0.0d+00
c
               njac(5,1,i) = - ( r43 * c34
     $                          - c1345 ) * tmp3 * (u(2,i,j,k)**2)
     $              - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k)**2)
     $              - ( c34 - c1345 ) * tmp3 * (u(4,i,j,k)**2)
     $              - c1345 * tmp2 * u(5,i,j,k)
c
               njac(5,2,i) = ( r43 * c34
     $                        - c1345 ) * tmp2 * u(2,i,j,k)
               njac(5,3,i) = ( c34 - c1345 ) * tmp2 * u(3,i,j,k)
               njac(5,4,i) = ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
               njac(5,5,i) = ( c1345 ) * tmp1
c
            end do
c
c***dirichlet boundary conditions
c
            do l = 1, 5
               do m = 1, 5
c     
                  a(m,l,1,j,k) = 0.0d+00
                  b(m,l,1,j,k) = 0.0d+00
                  c(m,l,1,j,k) = 0.0d+00
c
               end do
               b(l,l,1,j,k) = 1.0d+00
            end do
c
            do i = 2, nx - 1
c
               tmp1 = dt * tx1
               tmp2 = dt * tx2
c
               a(1,1,i,j,k) = - tmp2 * fjac(1,1,i-1)
     $                        - tmp1 * njac(1,1,i-1)
     $                        - tmp1 * dx1 
               a(1,2,i,j,k) = - tmp2 * fjac(1,2,i-1)
     $                        - tmp1 * njac(1,2,i-1)
               a(1,3,i,j,k) = - tmp2 * fjac(1,3,i-1)
     $                        - tmp1 * njac(1,3,i-1)
               a(1,4,i,j,k) = - tmp2 * fjac(1,4,i-1)
     $                        - tmp1 * njac(1,4,i-1)
               a(1,5,i,j,k) = - tmp2 * fjac(1,5,i-1)
     $                        - tmp1 * njac(1,5,i-1)
c
               a(2,1,i,j,k) = - tmp2 * fjac(2,1,i-1)
     $                        - tmp1 * njac(2,1,i-1)
               a(2,2,i,j,k) = - tmp2 * fjac(2,2,i-1)
     $                        - tmp1 * njac(2,2,i-1)
     $                        - tmp1 * dx2
               a(2,3,i,j,k) = - tmp2 * fjac(2,3,i-1)
     $                        - tmp1 * njac(2,3,i-1)
               a(2,4,i,j,k) = - tmp2 * fjac(2,4,i-1)
     $                        - tmp1 * njac(2,4,i-1)
               a(2,5,i,j,k) = - tmp2 * fjac(2,5,i-1)
     $                        - tmp1 * njac(2,5,i-1)
c
               a(3,1,i,j,k) = - tmp2 * fjac(3,1,i-1)
     $                        - tmp1 * njac(3,1,i-1)
               a(3,2,i,j,k) = - tmp2 * fjac(3,2,i-1)
     $                        - tmp1 * njac(3,2,i-1)
               a(3,3,i,j,k) = - tmp2 * fjac(3,3,i-1)
     $                        - tmp1 * njac(3,3,i-1)
     $                        - tmp1 * dx3 
               a(3,4,i,j,k) = - tmp2 * fjac(3,4,i-1)
     $                        - tmp1 * njac(3,4,i-1)
               a(3,5,i,j,k) = - tmp2 * fjac(3,5,i-1)
     $                        - tmp1 * njac(3,5,i-1)
c
               a(4,1,i,j,k) = - tmp2 * fjac(4,1,i-1)
     $                        - tmp1 * njac(4,1,i-1)
               a(4,2,i,j,k) = - tmp2 * fjac(4,2,i-1)
     $                        - tmp1 * njac(4,2,i-1)
               a(4,3,i,j,k) = - tmp2 * fjac(4,3,i-1)
     $                        - tmp1 * njac(4,3,i-1)
               a(4,4,i,j,k) = - tmp2 * fjac(4,4,i-1)
     $                        - tmp1 * njac(4,4,i-1)
     $                        - tmp1 * dx4
               a(4,5,i,j,k) = - tmp2 * fjac(4,5,i-1)
     $                        - tmp1 * njac(4,5,i-1)
c
               a(5,1,i,j,k) = - tmp2 * fjac(5,1,i-1)
     $                        - tmp1 * njac(5,1,i-1)
               a(5,2,i,j,k) = - tmp2 * fjac(5,2,i-1)
     $                        - tmp1 * njac(5,2,i-1)
               a(5,3,i,j,k) = - tmp2 * fjac(5,3,i-1)
     $                        - tmp1 * njac(5,3,i-1)
               a(5,4,i,j,k) = - tmp2 * fjac(5,4,i-1)
     $                        - tmp1 * njac(5,4,i-1)
               a(5,5,i,j,k) = - tmp2 * fjac(5,5,i-1)
     $                        - tmp1 * njac(5,5,i-1)
     $                        - tmp1 * dx5
c     
               b(1,1,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * njac(1,1,i)
     $                      + tmp1 * 2.0d+00 * dx1
               b(1,2,i,j,k) = tmp1 * 2.0d+00 * njac(1,2,i)
               b(1,3,i,j,k) = tmp1 * 2.0d+00 * njac(1,3,i)
               b(1,4,i,j,k) = tmp1 * 2.0d+00 * njac(1,4,i)
               b(1,5,i,j,k) = tmp1 * 2.0d+00 * njac(1,5,i)
c     
               b(2,1,i,j,k) = tmp1 * 2.0d+00 * njac(2,1,i)
               b(2,2,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * njac(2,2,i)
     $                      + tmp1 * 2.0d+00 * dx2
               b(2,3,i,j,k) = tmp1 * 2.0d+00 * njac(2,3,i)
               b(2,4,i,j,k) = tmp1 * 2.0d+00 * njac(2,4,i)
               b(2,5,i,j,k) = tmp1 * 2.0d+00 * njac(2,5,i)
c     
               b(3,1,i,j,k) = tmp1 * 2.0d+00 * njac(3,1,i)
               b(3,2,i,j,k) = tmp1 * 2.0d+00 * njac(3,2,i)
               b(3,3,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * njac(3,3,i)
     $                      + tmp1 * 2.0d+00 * dx3
               b(3,4,i,j,k) = tmp1 * 2.0d+00 * njac(3,4,i)
               b(3,5,i,j,k) = tmp1 * 2.0d+00 * njac(3,5,i)
c     
               b(4,1,i,j,k) = tmp1 * 2.0d+00 * njac(4,1,i)
               b(4,2,i,j,k) = tmp1 * 2.0d+00 * njac(4,2,i)
               b(4,3,i,j,k) = tmp1 * 2.0d+00 * njac(4,3,i)
               b(4,4,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * njac(4,4,i)
     $                      + tmp1 * 2.0d+00 * dx4
               b(4,5,i,j,k) = tmp1 * 2.0d+00 * njac(4,5,i)
c
               b(5,1,i,j,k) = tmp1 * 2.0d+00 * njac(5,1,i)
               b(5,2,i,j,k) = tmp1 * 2.0d+00 * njac(5,2,i)
               b(5,3,i,j,k) = tmp1 * 2.0d+00 * njac(5,3,i)
               b(5,4,i,j,k) = tmp1 * 2.0d+00 * njac(5,4,i)
               b(5,5,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * njac(5,5,i)
     $                      + tmp1 * 2.0d+00 * dx5
c     
               c(1,1,i,j,k) =  tmp2 * fjac(1,1,i+1)
     $                       - tmp1 * njac(1,1,i+1)
     $                       - tmp1 * dx1
               c(1,2,i,j,k) =  tmp2 * fjac(1,2,i+1)
     $                       - tmp1 * njac(1,2,i+1)
               c(1,3,i,j,k) =  tmp2 * fjac(1,3,i+1)
     $                       - tmp1 * njac(1,3,i+1)
               c(1,4,i,j,k) =  tmp2 * fjac(1,4,i+1)
     $                       - tmp1 * njac(1,4,i+1)
               c(1,5,i,j,k) =  tmp2 * fjac(1,5,i+1)
     $                       - tmp1 * njac(1,5,i+1)
c     
               c(2,1,i,j,k) =  tmp2 * fjac(2,1,i+1)
     $                       - tmp1 * njac(2,1,i+1)
               c(2,2,i,j,k) =  tmp2 * fjac(2,2,i+1)
     $                       - tmp1 * njac(2,2,i+1)
     $                       - tmp1 * dx2
               c(2,3,i,j,k) =  tmp2 * fjac(2,3,i+1)
     $                       - tmp1 * njac(2,3,i+1)
               c(2,4,i,j,k) =  tmp2 * fjac(2,4,i+1)
     $                       - tmp1 * njac(2,4,i+1)
               c(2,5,i,j,k) =  tmp2 * fjac(2,5,i+1)
     $                       - tmp1 * njac(2,5,i+1)
c     
               c(3,1,i,j,k) =  tmp2 * fjac(3,1,i+1)
     $                       - tmp1 * njac(3,1,i+1)
               c(3,2,i,j,k) =  tmp2 * fjac(3,2,i+1)
     $                       - tmp1 * njac(3,2,i+1)
               c(3,3,i,j,k) =  tmp2 * fjac(3,3,i+1)
     $                       - tmp1 * njac(3,3,i+1)
     $                       - tmp1 * dx3
               c(3,4,i,j,k) =  tmp2 * fjac(3,4,i+1)
     $                       - tmp1 * njac(3,4,i+1)
               c(3,5,i,j,k) =  tmp2 * fjac(3,5,i+1)
     $                       - tmp1 * njac(3,5,i+1)
c     
               c(4,1,i,j,k) =  tmp2 * fjac(4,1,i+1)
     $                       - tmp1 * njac(4,1,i+1)
               c(4,2,i,j,k) =  tmp2 * fjac(4,2,i+1)
     $                       - tmp1 * njac(4,2,i+1)
               c(4,3,i,j,k) =  tmp2 * fjac(4,3,i+1)
     $                       - tmp1 * njac(4,3,i+1)
               c(4,4,i,j,k) =  tmp2 * fjac(4,4,i+1)
     $                       - tmp1 * njac(4,4,i+1)
     $                       - tmp1 * dx4
               c(4,5,i,j,k) =  tmp2 * fjac(4,5,i+1)
     $                       - tmp1 * njac(4,5,i+1)
c     
               c(5,1,i,j,k) =  tmp2 * fjac(5,1,i+1)
     $                       - tmp1 * njac(5,1,i+1)
               c(5,2,i,j,k) =  tmp2 * fjac(5,2,i+1)
     $                       - tmp1 * njac(5,2,i+1)
               c(5,3,i,j,k) =  tmp2 * fjac(5,3,i+1)
     $                       - tmp1 * njac(5,3,i+1)
               c(5,4,i,j,k) =  tmp2 * fjac(5,4,i+1)
     $                       - tmp1 * njac(5,4,i+1)
               c(5,5,i,j,k) =  tmp2 * fjac(5,5,i+1)
     $                       - tmp1 * njac(5,5,i+1)
     $                       - tmp1 * dx5
c     
            end do
c
c***dirichlet boundary conditions
c
            do l = 1, 5
               do m = 1, 5
c     
                  a(m,l,nx,j,k) = 0.0d+00
                  b(m,l,nx,j,k) = 0.0d+00
                  c(m,l,nx,j,k) = 0.0d+00
c
               end do
               b(l,l,nx,j,k) = 1.0d+00
            end do
c
         end do
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine jacy
c
c***form the block tridiagonal system for eta-direction sweep
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension fjac(5,5,isiz2), qjac(5,5,isiz2)
c
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c
      do k = 2, nz - 1
c
         do i = 2, nx - 1
c
c***compute the eta-direction flux jacobians
c
            do j = 1, ny
c
               fjac(1,1,j) = 0.0d+00
               fjac(1,2,j) = 0.0d+00
               fjac(1,3,j) = 1.0d+00
               fjac(1,4,j) = 0.0d+00
               fjac(1,5,j) = 0.0d+00
c
               fjac(2,1,j) = - ( u(2,i,j,k)*u(3,i,j,k) )
     $                       / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(2,2,j) = u(3,i,j,k) / u(1,i,j,k)
               fjac(2,3,j) = u(2,i,j,k) / u(1,i,j,k)
               fjac(2,4,j) = 0.0d+00
               fjac(2,5,j) = 0.0d+00
c
               fjac(3,1,j) = - ( u(3,i,j,k) / u(1,i,j,k) ) ** 2
     $           + 0.50d+00 * c2 * ( (  u(2,i,j,k) * u(2,i,j,k)
     $                                + u(3,i,j,k) * u(3,i,j,k)
     $                                + u(4,i,j,k) * u(4,i,j,k) )
     $                              / ( u(1,i,j,k) * u(1,i,j,k) ) )
               fjac(3,2,j) = - c2* ( u(2,i,j,k) / u(1,i,j,k) )
               fjac(3,3,j) = ( 2.0d+00 - c2 )
     $                     * ( u(3,i,j,k) / u(1,i,j,k) )
               fjac(3,4,j) = - c2 * ( u(4,i,j,k) / u(1,i,j,k) )
               fjac(3,5,j) = c2
c
               fjac(4,1,j) = - ( u(3,i,j,k)*u(4,i,j,k) )
     $                       / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(4,2,j) = 0.0d+00
               fjac(4,3,j) = u(4,i,j,k) / u(1,i,j,k)
               fjac(4,4,j) = u(3,i,j,k) / u(1,i,j,k)
               fjac(4,5,j) = 0.0d+00
c
               fjac(5,1,j) = ( c2 * (  u(2,i,j,k) * u(2,i,j,k)
     $                               + u(3,i,j,k) * u(3,i,j,k)
     $                               + u(4,i,j,k) * u(4,i,j,k) )
     $                            / ( u(1,i,j,k) * u(1,i,j,k) )
     $                         - c1 * ( u(5,i,j,k) / u(1,i,j,k) ) )
     $                        * ( u(3,i,j,k) / u(1,i,j,k) )
               fjac(5,2,j) = - c2 * ( u(2,i,j,k)*u(3,i,j,k) )
     $              / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(5,3,j) = c1 * ( u(5,i,j,k) / u(1,i,j,k) )
     $                      - 0.50d+00 * c2 
     $                      * ( (  u(2,i,j,k)*u(2,i,j,k)
     $                           + 3.0d+00 * u(3,i,j,k)*u(3,i,j,k)
     $                           + u(4,i,j,k)*u(4,i,j,k) )
     $                        / ( u(1,i,j,k)*u(1,i,j,k) ) )
               fjac(5,4,j) = - c2 * ( u(3,i,j,k)*u(4,i,j,k) )
     $              / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(5,5,j) = c1 * ( u(3,i,j,k) / u(1,i,j,k) )
c
               tmp1 = 1.0d+00 / u(1,i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2
c
               qjac(1,1,j) = 0.0d+00
               qjac(1,2,j) = 0.0d+00
               qjac(1,3,j) = 0.0d+00
               qjac(1,4,j) = 0.0d+00
               qjac(1,5,j) = 0.0d+00
c
               qjac(2,1,j) = - c34 * tmp2 * u(2,i,j,k)
               qjac(2,2,j) =   c34 * tmp1
               qjac(2,3,j) =   0.0d+00
               qjac(2,4,j) =   0.0d+00
               qjac(2,5,j) =   0.0d+00
c     
               qjac(3,1,j) = - r43 * c34 * tmp2 * u(3,i,j,k)
               qjac(3,2,j) =   0.0d+00
               qjac(3,3,j) =   r43 * c34 * tmp1
               qjac(3,4,j) =   0.0d+00
               qjac(3,5,j) =   0.0d+00
c
               qjac(4,1,j) = - c34 * tmp2 * u(4,i,j,k)
               qjac(4,2,j) =   0.0d+00
               qjac(4,3,j) =   0.0d+00
               qjac(4,4,j) =   c34 * tmp1
               qjac(4,5,j) =   0.0d+00
c
               qjac(5,1,j) = - (  c34
     $                          - c1345 ) * tmp3 * (u(2,i,j,k)**2)
     $              - ( r43 * c34
     $                 - c1345 ) * tmp3 * (u(3,i,j,k)**2)
     $              - ( c34 - c1345 ) * tmp3 * (u(4,i,j,k)**2)
     $              - c1345 * tmp2 * u(5,i,j,k)
c
               qjac(5,2,j) = (  c34 - c1345 ) * tmp2 * u(2,i,j,k)
               qjac(5,3,j) = ( r43 * c34
     $                        - c1345 ) * tmp2 * u(3,i,j,k)
               qjac(5,4,j) = ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
               qjac(5,5,j) = ( c1345 ) * tmp1
c
            end do
c
            do l = 1, 5
               do m = 1, 5
c                  
                  a(m,l,i,1,k) = 0.0d+00
                  b(m,l,i,1,k) = 0.0d+00
                  c(m,l,i,1,k) = 0.0d+00
c
               end do
               b(l,l,i,1,k) = 1.0d+00
            end do
c
            do j = 2, ny - 1
c
               tmp1 = dt * ty1
               tmp2 = dt * ty2
c
               a(1,1,i,j,k) = - tmp2 * fjac(1,1,j-1)
     $                        - tmp1 * qjac(1,1,j-1)
     $                        - tmp1 * dy1 
               a(1,2,i,j,k) = - tmp2 * fjac(1,2,j-1)
     $                        - tmp1 * qjac(1,2,j-1)
               a(1,3,i,j,k) = - tmp2 * fjac(1,3,j-1)
     $                        - tmp1 * qjac(1,3,j-1)
               a(1,4,i,j,k) = - tmp2 * fjac(1,4,j-1)
     $                        - tmp1 * qjac(1,4,j-1)
               a(1,5,i,j,k) = - tmp2 * fjac(1,5,j-1)
     $                        - tmp1 * qjac(1,5,j-1)
c
               a(2,1,i,j,k) = - tmp2 * fjac(2,1,j-1)
     $                        - tmp1 * qjac(2,1,j-1)
               a(2,2,i,j,k) = - tmp2 * fjac(2,2,j-1)
     $                        - tmp1 * qjac(2,2,j-1)
     $                        - tmp1 * dy2
               a(2,3,i,j,k) = - tmp2 * fjac(2,3,j-1)
     $                        - tmp1 * qjac(2,3,j-1)
               a(2,4,i,j,k) = - tmp2 * fjac(2,4,j-1)
     $                        - tmp1 * qjac(2,4,j-1)
               a(2,5,i,j,k) = - tmp2 * fjac(2,5,j-1)
     $                        - tmp1 * qjac(2,5,j-1)
c
               a(3,1,i,j,k) = - tmp2 * fjac(3,1,j-1)
     $                        - tmp1 * qjac(3,1,j-1)
               a(3,2,i,j,k) = - tmp2 * fjac(3,2,j-1)
     $                        - tmp1 * qjac(3,2,j-1)
               a(3,3,i,j,k) = - tmp2 * fjac(3,3,j-1)
     $                        - tmp1 * qjac(3,3,j-1)
     $                        - tmp1 * dy3 
               a(3,4,i,j,k) = - tmp2 * fjac(3,4,j-1)
     $                        - tmp1 * qjac(3,4,j-1)
               a(3,5,i,j,k) = - tmp2 * fjac(3,5,j-1)
     $                        - tmp1 * qjac(3,5,j-1)
c
               a(4,1,i,j,k) = - tmp2 * fjac(4,1,j-1)
     $                        - tmp1 * qjac(4,1,j-1)
               a(4,2,i,j,k) = - tmp2 * fjac(4,2,j-1)
     $                        - tmp1 * qjac(4,2,j-1)
               a(4,3,i,j,k) = - tmp2 * fjac(4,3,j-1)
     $                        - tmp1 * qjac(4,3,j-1)
               a(4,4,i,j,k) = - tmp2 * fjac(4,4,j-1)
     $                        - tmp1 * qjac(4,4,j-1)
     $                        - tmp1 * dy4
               a(4,5,i,j,k) = - tmp2 * fjac(4,5,j-1)
     $                        - tmp1 * qjac(4,5,j-1)
c
               a(5,1,i,j,k) = - tmp2 * fjac(5,1,j-1)
     $                        - tmp1 * qjac(5,1,j-1)
               a(5,2,i,j,k) = - tmp2 * fjac(5,2,j-1)
     $                        - tmp1 * qjac(5,2,j-1)
               a(5,3,i,j,k) = - tmp2 * fjac(5,3,j-1)
     $                        - tmp1 * qjac(5,3,j-1)
               a(5,4,i,j,k) = - tmp2 * fjac(5,4,j-1)
     $                        - tmp1 * qjac(5,4,j-1)
               a(5,5,i,j,k) = - tmp2 * fjac(5,5,j-1)
     $                        - tmp1 * qjac(5,5,j-1)
     $                        - tmp1 * dy5
c     
               b(1,1,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * qjac(1,1,j)
     $                      + tmp1 * 2.0d+00 * dy1
               b(1,2,i,j,k) = tmp1 * 2.0d+00 * qjac(1,2,j)
               b(1,3,i,j,k) = tmp1 * 2.0d+00 * qjac(1,3,j)
               b(1,4,i,j,k) = tmp1 * 2.0d+00 * qjac(1,4,j)
               b(1,5,i,j,k) = tmp1 * 2.0d+00 * qjac(1,5,j)
c     
               b(2,1,i,j,k) = tmp1 * 2.0d+00 * qjac(2,1,j)
               b(2,2,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * qjac(2,2,j)
     $                      + tmp1 * 2.0d+00 * dy2
               b(2,3,i,j,k) = tmp1 * 2.0d+00 * qjac(2,3,j)
               b(2,4,i,j,k) = tmp1 * 2.0d+00 * qjac(2,4,j)
               b(2,5,i,j,k) = tmp1 * 2.0d+00 * qjac(2,5,j)
c     
               b(3,1,i,j,k) = tmp1 * 2.0d+00 * qjac(3,1,j)
               b(3,2,i,j,k) = tmp1 * 2.0d+00 * qjac(3,2,j)
               b(3,3,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * qjac(3,3,j)
     $                      + tmp1 * 2.0d+00 * dy3
               b(3,4,i,j,k) = tmp1 * 2.0d+00 * qjac(3,4,j)
               b(3,5,i,j,k) = tmp1 * 2.0d+00 * qjac(3,5,j)
c     
               b(4,1,i,j,k) = tmp1 * 2.0d+00 * qjac(4,1,j)
               b(4,2,i,j,k) = tmp1 * 2.0d+00 * qjac(4,2,j)
               b(4,3,i,j,k) = tmp1 * 2.0d+00 * qjac(4,3,j)
               b(4,4,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * qjac(4,4,j)
     $                      + tmp1 * 2.0d+00 * dy4
               b(4,5,i,j,k) = tmp1 * 2.0d+00 * qjac(4,5,j)
c
               b(5,1,i,j,k) = tmp1 * 2.0d+00 * qjac(5,1,j)
               b(5,2,i,j,k) = tmp1 * 2.0d+00 * qjac(5,2,j)
               b(5,3,i,j,k) = tmp1 * 2.0d+00 * qjac(5,3,j)
               b(5,4,i,j,k) = tmp1 * 2.0d+00 * qjac(5,4,j)
               b(5,5,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * qjac(5,5,j) 
     $                      + tmp1 * 2.0d+00 * dy5
c     
               c(1,1,i,j,k) =  tmp2 * fjac(1,1,j+1)
     $                       - tmp1 * qjac(1,1,j+1)
     $                       - tmp1 * dy1
               c(1,2,i,j,k) =  tmp2 * fjac(1,2,j+1)
     $                       - tmp1 * qjac(1,2,j+1)
               c(1,3,i,j,k) =  tmp2 * fjac(1,3,j+1)
     $                       - tmp1 * qjac(1,3,j+1)
               c(1,4,i,j,k) =  tmp2 * fjac(1,4,j+1)
     $                       - tmp1 * qjac(1,4,j+1)
               c(1,5,i,j,k) =  tmp2 * fjac(1,5,j+1)
     $                       - tmp1 * qjac(1,5,j+1)
c     
               c(2,1,i,j,k) =  tmp2 * fjac(2,1,j+1)
     $                       - tmp1 * qjac(2,1,j+1)
               c(2,2,i,j,k) =  tmp2 * fjac(2,2,j+1)
     $                       - tmp1 * qjac(2,2,j+1)
     $                       - tmp1 * dy2
               c(2,3,i,j,k) =  tmp2 * fjac(2,3,j+1)
     $                       - tmp1 * qjac(2,3,j+1)
               c(2,4,i,j,k) =  tmp2 * fjac(2,4,j+1)
     $                       - tmp1 * qjac(2,4,j+1)
               c(2,5,i,j,k) =  tmp2 * fjac(2,5,j+1)
     $                       - tmp1 * qjac(2,5,j+1)
c     
               c(3,1,i,j,k) =  tmp2 * fjac(3,1,j+1)
     $                       - tmp1 * qjac(3,1,j+1)
               c(3,2,i,j,k) =  tmp2 * fjac(3,2,j+1)
     $                       - tmp1 * qjac(3,2,j+1)
               c(3,3,i,j,k) =  tmp2 * fjac(3,3,j+1)
     $                       - tmp1 * qjac(3,3,j+1)
     $                       - tmp1 * dy3
               c(3,4,i,j,k) =  tmp2 * fjac(3,4,j+1)
     $                       - tmp1 * qjac(3,4,j+1)
               c(3,5,i,j,k) =  tmp2 * fjac(3,5,j+1)
     $                       - tmp1 * qjac(3,5,j+1)
c     
               c(4,1,i,j,k) =  tmp2 * fjac(4,1,j+1)
     $                       - tmp1 * qjac(4,1,j+1)
               c(4,2,i,j,k) =  tmp2 * fjac(4,2,j+1)
     $                       - tmp1 * qjac(4,2,j+1)
               c(4,3,i,j,k) =  tmp2 * fjac(4,3,j+1)
     $                       - tmp1 * qjac(4,3,j+1)
               c(4,4,i,j,k) =  tmp2 * fjac(4,4,j+1)
     $                       - tmp1 * qjac(4,4,j+1)
     $                       - tmp1 * dy4
               c(4,5,i,j,k) =  tmp2 * fjac(4,5,j+1)
     $                       - tmp1 * qjac(4,5,j+1)
c     
               c(5,1,i,j,k) =  tmp2 * fjac(5,1,j+1)
     $                       - tmp1 * qjac(5,1,j+1)
               c(5,2,i,j,k) =  tmp2 * fjac(5,2,j+1)
     $                       - tmp1 * qjac(5,2,j+1)
               c(5,3,i,j,k) =  tmp2 * fjac(5,3,j+1)
     $                       - tmp1 * qjac(5,3,j+1)
               c(5,4,i,j,k) =  tmp2 * fjac(5,4,j+1)
     $                       - tmp1 * qjac(5,4,j+1)
               c(5,5,i,j,k) =  tmp2 * fjac(5,5,j+1)
     $                       - tmp1 * qjac(5,5,j+1)
     $                       - tmp1 * dy5
c     
            end do
c
            do l = 1, 5
               do m = 1, 5
c                  
                  a(m,l,i,ny,k) = 0.0d+00
                  b(m,l,i,ny,k) = 0.0d+00
                  c(m,l,i,ny,k) = 0.0d+00
c
               end do
               b(l,l,i,ny,k) = 1.0d+00
            end do
c
         end do
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine jacz
c
c***form the block tridiagonal system for zeta-direction sweep
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension fjac(5,5,isiz3), sjac(5,5,isiz3)
c
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c
      do j = 2, ny - 1
c
         do i = 2, nx - 1
c
c***compute the zeta-direction flux jacobian
c
            do k = 1, nz
c
               fjac(1,1,k) = 0.0d+00
               fjac(1,2,k) = 0.0d+00
               fjac(1,3,k) = 0.0d+00
               fjac(1,4,k) = 1.0d+00
               fjac(1,5,k) = 0.0d+00
c
               fjac(2,1,k) = - ( u(2,i,j,k)*u(4,i,j,k) )
     $                       / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(2,2,k) = u(4,i,j,k) / u(1,i,j,k)
               fjac(2,3,k) = 0.0d+00
               fjac(2,4,k) = u(2,i,j,k) / u(1,i,j,k)
               fjac(2,5,k) = 0.0d+00
c
               fjac(3,1,k) = - ( u(3,i,j,k)*u(4,i,j,k) )
     $                       / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(3,2,k) = 0.0d+00
               fjac(3,3,k) = u(4,i,j,k) / u(1,i,j,k)
               fjac(3,4,k) = u(3,i,j,k) / u(1,i,j,k)
               fjac(3,5,k) = 0.0d+00
c
               fjac(4,1,k) = - ( u(4,i,j,k) / u(1,i,j,k) ) ** 2
     $           + 0.50d+00 * c2 * ( (  u(2,i,j,k) * u(2,i,j,k)
     $                                + u(3,i,j,k) * u(3,i,j,k)
     $                                + u(4,i,j,k) * u(4,i,j,k) )
     $                             / ( u(1,i,j,k) * u(1,i,j,k) ) )
               fjac(4,2,k) = - c2 * ( u(2,i,j,k) / u(1,i,j,k) )
               fjac(4,3,k) = - c2 * ( u(3,i,j,k) / u(1,i,j,k) )
               fjac(4,4,k) = ( 2.0d+00 - c2 )
     $                     * ( u(4,i,j,k) / u(1,i,j,k) )
               fjac(4,5,k) = c2
c
               fjac(5,1,k) = ( c2 * (  u(2,i,j,k) * u(2,i,j,k)
     $                               + u(3,i,j,k) * u(3,i,j,k)
     $                               + u(4,i,j,k) * u(4,i,j,k) )
     $                            / ( u(1,i,j,k) * u(1,i,j,k) )
     $                         - c1 * ( u(5,i,j,k) / u(1,i,j,k) ) )
     $                        * ( u(4,i,j,k) / u(1,i,j,k) )
               fjac(5,2,k) = - c2 * ( u(2,i,j,k)*u(4,i,j,k) )
     $              / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(5,3,k) = - c2 * ( u(3,i,j,k)*u(4,i,j,k) )
     $              / ( u(1,i,j,k)*u(1,i,j,k) )
               fjac(5,4,k) = c1 * ( u(5,i,j,k) / u(1,i,j,k) )
     $                      - 0.50d+00 * c2
     $                      * ( (  u(2,i,j,k)*u(2,i,j,k)
     $                           + u(3,i,j,k)*u(3,i,j,k)
     $                           + 3.0d+00*u(4,i,j,k)*u(4,i,j,k) )
     $                      / ( u(1,i,j,k)*u(1,i,j,k) ) )
               fjac(5,5,k) = c1 * ( u(4,i,j,k) / u(1,i,j,k) )
c
               tmp1 = 1.0d+00 / u(1,i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2
c
               sjac(1,1,k) = 0.0d+00
               sjac(1,2,k) = 0.0d+00
               sjac(1,3,k) = 0.0d+00
               sjac(1,4,k) = 0.0d+00
               sjac(1,5,k) = 0.0d+00
c
               sjac(2,1,k) = - c34 * tmp2 * u(2,i,j,k)
               sjac(2,2,k) =   c34 * tmp1
               sjac(2,3,k) =   0.0d+00
               sjac(2,4,k) =   0.0d+00
               sjac(2,5,k) =   0.0d+00
c     
               sjac(3,1,k) = - c34 * tmp2 * u(3,i,j,k)
               sjac(3,2,k) =   0.0d+00
               sjac(3,3,k) =   c34 * tmp1
               sjac(3,4,k) =   0.0d+00
               sjac(3,5,k) =   0.0d+00
c
               sjac(4,1,k) = - r43 * c34 * tmp2 * u(4,i,j,k)
               sjac(4,2,k) =   0.0d+00
               sjac(4,3,k) =   0.0d+00
               sjac(4,4,k) =   r43 * c3 * c4 * tmp1
               sjac(4,5,k) =   0.0d+00
c
               sjac(5,1,k) = - (  c34
     $                          - c1345 ) * tmp3 * (u(2,i,j,k)**2)
     $              - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k)**2)
     $              - ( r43 * c34
     $                 - c1345 ) * tmp3 * (u(4,i,j,k)**2)
     $              - c1345 * tmp2 * u(5,i,j,k)
c
               sjac(5,2,k) = (  c34 - c1345 ) * tmp2 * u(2,i,j,k)
               sjac(5,3,k) = (  c34 - c1345 ) * tmp2 * u(3,i,j,k)
               sjac(5,4,k) = ( r43 * c34
     $                       - c1345 ) * tmp2 * u(4,i,j,k)
               sjac(5,5,k) = ( c1345 )* tmp1
c
            end do
c
            do l = 1, 5
               do m = 1, 5
c                  
                  a(m,l,i,j,1) = 0.0d+00
                  b(m,l,i,j,1) = 0.0d+00
                  c(m,l,i,j,1) = 0.0d+00
c
               end do
               b(l,l,i,j,1) = 1.0d+00
            end do
c
            do k = 2, nz - 1
c
               tmp1 = dt * tz1
               tmp2 = dt * tz2
c
               a(1,1,i,j,k) = - tmp2 * fjac(1,1,k-1)
     $                        - tmp1 * sjac(1,1,k-1)
     $                        - tmp1 * dz1 
               a(1,2,i,j,k) = - tmp2 * fjac(1,2,k-1)
     $                        - tmp1 * sjac(1,2,k-1)
               a(1,3,i,j,k) = - tmp2 * fjac(1,3,k-1)
     $                        - tmp1 * sjac(1,3,k-1)
               a(1,4,i,j,k) = - tmp2 * fjac(1,4,k-1)
     $                        - tmp1 * sjac(1,4,k-1)
               a(1,5,i,j,k) = - tmp2 * fjac(1,5,k-1)
     $                        - tmp1 * sjac(1,5,k-1)
c
               a(2,1,i,j,k) = - tmp2 * fjac(2,1,k-1)
     $                        - tmp1 * sjac(2,1,k-1)
               a(2,2,i,j,k) = - tmp2 * fjac(2,2,k-1)
     $                        - tmp1 * sjac(2,2,k-1)
     $                        - tmp1 * dz2
               a(2,3,i,j,k) = - tmp2 * fjac(2,3,k-1)
     $                        - tmp1 * sjac(2,3,k-1)
               a(2,4,i,j,k) = - tmp2 * fjac(2,4,k-1)
     $                        - tmp1 * sjac(2,4,k-1)
               a(2,5,i,j,k) = - tmp2 * fjac(2,5,k-1)
     $                        - tmp1 * sjac(2,5,k-1)
c
               a(3,1,i,j,k) = - tmp2 * fjac(3,1,k-1)
     $                        - tmp1 * sjac(3,1,k-1)
               a(3,2,i,j,k) = - tmp2 * fjac(3,2,k-1)
     $                        - tmp1 * sjac(3,2,k-1)
               a(3,3,i,j,k) = - tmp2 * fjac(3,3,k-1)
     $                        - tmp1 * sjac(3,3,k-1)
     $                        - tmp1 * dz3 
               a(3,4,i,j,k) = - tmp2 * fjac(3,4,k-1)
     $                        - tmp1 * sjac(3,4,k-1)
               a(3,5,i,j,k) = - tmp2 * fjac(3,5,k-1)
     $                        - tmp1 * sjac(3,5,k-1)
c
               a(4,1,i,j,k) = - tmp2 * fjac(4,1,k-1)
     $                        - tmp1 * sjac(4,1,k-1)
               a(4,2,i,j,k) = - tmp2 * fjac(4,2,k-1)
     $                        - tmp1 * sjac(4,2,k-1)
               a(4,3,i,j,k) = - tmp2 * fjac(4,3,k-1)
     $                        - tmp1 * sjac(4,3,k-1)
               a(4,4,i,j,k) = - tmp2 * fjac(4,4,k-1)
     $                        - tmp1 * sjac(4,4,k-1)
     $                        - tmp1 * dz4
               a(4,5,i,j,k) = - tmp2 * fjac(4,5,k-1)
     $                        - tmp1 * sjac(4,5,k-1)
c
               a(5,1,i,j,k) = - tmp2 * fjac(5,1,k-1)
     $                        - tmp1 * sjac(5,1,k-1)
               a(5,2,i,j,k) = - tmp2 * fjac(5,2,k-1)
     $                        - tmp1 * sjac(5,2,k-1)
               a(5,3,i,j,k) = - tmp2 * fjac(5,3,k-1)
     $                        - tmp1 * sjac(5,3,k-1)
               a(5,4,i,j,k) = - tmp2 * fjac(5,4,k-1)
     $                        - tmp1 * sjac(5,4,k-1)
               a(5,5,i,j,k) = - tmp2 * fjac(5,5,k-1)
     $                        - tmp1 * sjac(5,5,k-1)
     $                        - tmp1 * dz5
c     
               b(1,1,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * sjac(1,1,k)
     $                      + tmp1 * 2.0d+00 * dz1
               b(1,2,i,j,k) = tmp1 * 2.0d+00 * sjac(1,2,k)
               b(1,3,i,j,k) = tmp1 * 2.0d+00 * sjac(1,3,k)
               b(1,4,i,j,k) = tmp1 * 2.0d+00 * sjac(1,4,k)
               b(1,5,i,j,k) = tmp1 * 2.0d+00 * sjac(1,5,k)
c     
               b(2,1,i,j,k) = tmp1 * 2.0d+00 * sjac(2,1,k)
               b(2,2,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * sjac(2,2,k)
     $                      + tmp1 * 2.0d+00 * dz2
               b(2,3,i,j,k) = tmp1 * 2.0d+00 * sjac(2,3,k)
               b(2,4,i,j,k) = tmp1 * 2.0d+00 * sjac(2,4,k)
               b(2,5,i,j,k) = tmp1 * 2.0d+00 * sjac(2,5,k)
c     
               b(3,1,i,j,k) = tmp1 * 2.0d+00 * sjac(3,1,k)
               b(3,2,i,j,k) = tmp1 * 2.0d+00 * sjac(3,2,k)
               b(3,3,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * sjac(3,3,k)
     $                      + tmp1 * 2.0d+00 * dz3
               b(3,4,i,j,k) = tmp1 * 2.0d+00 * sjac(3,4,k)
               b(3,5,i,j,k) = tmp1 * 2.0d+00 * sjac(3,5,k)
c     
               b(4,1,i,j,k) = tmp1 * 2.0d+00 * sjac(4,1,k)
               b(4,2,i,j,k) = tmp1 * 2.0d+00 * sjac(4,2,k)
               b(4,3,i,j,k) = tmp1 * 2.0d+00 * sjac(4,3,k)
               b(4,4,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * sjac(4,4,k)
     $                      + tmp1 * 2.0d+00 * dz4
               b(4,5,i,j,k) = tmp1 * 2.0d+00 * sjac(4,5,k)
c
               b(5,1,i,j,k) = tmp1 * 2.0d+00 * sjac(5,1,k)
               b(5,2,i,j,k) = tmp1 * 2.0d+00 * sjac(5,2,k)
               b(5,3,i,j,k) = tmp1 * 2.0d+00 * sjac(5,3,k)
               b(5,4,i,j,k) = tmp1 * 2.0d+00 * sjac(5,4,k)
               b(5,5,i,j,k) = 1.0d+00
     $                      + tmp1 * 2.0d+00 * sjac(5,5,k) 
     $                      + tmp1 * 2.0d+00 * dz5
c     
               c(1,1,i,j,k) =  tmp2 * fjac(1,1,k+1)
     $                       - tmp1 * sjac(1,1,k+1)
     $                       - tmp1 * dz1
               c(1,2,i,j,k) =  tmp2 * fjac(1,2,k+1)
     $                       - tmp1 * sjac(1,2,k+1)
               c(1,3,i,j,k) =  tmp2 * fjac(1,3,k+1)
     $                       - tmp1 * sjac(1,3,k+1)
               c(1,4,i,j,k) =  tmp2 * fjac(1,4,k+1)
     $                       - tmp1 * sjac(1,4,k+1)
               c(1,5,i,j,k) =  tmp2 * fjac(1,5,k+1)
     $                       - tmp1 * sjac(1,5,k+1)
c     
               c(2,1,i,j,k) =  tmp2 * fjac(2,1,k+1)
     $                       - tmp1 * sjac(2,1,k+1)
               c(2,2,i,j,k) =  tmp2 * fjac(2,2,k+1)
     $                       - tmp1 * sjac(2,2,k+1)
     $                       - tmp1 * dz2
               c(2,3,i,j,k) =  tmp2 * fjac(2,3,k+1)
     $                       - tmp1 * sjac(2,3,k+1)
               c(2,4,i,j,k) =  tmp2 * fjac(2,4,k+1)
     $                       - tmp1 * sjac(2,4,k+1)
               c(2,5,i,j,k) =  tmp2 * fjac(2,5,k+1)
     $                       - tmp1 * sjac(2,5,k+1)
c     
               c(3,1,i,j,k) =  tmp2 * fjac(3,1,k+1)
     $                       - tmp1 * sjac(3,1,k+1)
               c(3,2,i,j,k) =  tmp2 * fjac(3,2,k+1)
     $                       - tmp1 * sjac(3,2,k+1)
               c(3,3,i,j,k) =  tmp2 * fjac(3,3,k+1)
     $                       - tmp1 * sjac(3,3,k+1)
     $                       - tmp1 * dz3
               c(3,4,i,j,k) =  tmp2 * fjac(3,4,k+1)
     $                       - tmp1 * sjac(3,4,k+1)
               c(3,5,i,j,k) =  tmp2 * fjac(3,5,k+1)
     $                       - tmp1 * sjac(3,5,k+1)
c     
               c(4,1,i,j,k) =  tmp2 * fjac(4,1,k+1)
     $                       - tmp1 * sjac(4,1,k+1)
               c(4,2,i,j,k) =  tmp2 * fjac(4,2,k+1)
     $                       - tmp1 * sjac(4,2,k+1)
               c(4,3,i,j,k) =  tmp2 * fjac(4,3,k+1)
     $                       - tmp1 * sjac(4,3,k+1)
               c(4,4,i,j,k) =  tmp2 * fjac(4,4,k+1)
     $                       - tmp1 * sjac(4,4,k+1)
     $                       - tmp1 * dz4
               c(4,5,i,j,k) =  tmp2 * fjac(4,5,k+1)
     $                       - tmp1 * sjac(4,5,k+1)
c     
               c(5,1,i,j,k) =  tmp2 * fjac(5,1,k+1)
     $                       - tmp1 * sjac(5,1,k+1)
               c(5,2,i,j,k) =  tmp2 * fjac(5,2,k+1)
     $                       - tmp1 * sjac(5,2,k+1)
               c(5,3,i,j,k) =  tmp2 * fjac(5,3,k+1)
     $                       - tmp1 * sjac(5,3,k+1)
               c(5,4,i,j,k) =  tmp2 * fjac(5,4,k+1)
     $                       - tmp1 * sjac(5,4,k+1)
               c(5,5,i,j,k) =  tmp2 * fjac(5,5,k+1)
     $                       - tmp1 * sjac(5,5,k+1)
     $                       - tmp1 * dz5
c     
            end do
c
            do l = 1, 5
               do m = 1, 5
c                  
                  a(m,l,i,j,nz) = 0.0d+00
                  b(m,l,i,j,nz) = 0.0d+00
                  c(m,l,i,j,nz) = 0.0d+00
c
               end do
               b(l,l,i,j,nz) = 1.0d+00
            end do
c
         end do
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine pintgr
c
c***compute the surface integral
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension phi1(isiz1,isiz2), phi2(isiz1,isiz2)
c
      do j = ji1, ji2
c     
         do i = ii1, ii2
c
            phi1(i,j) = c2*(  u(5,i,j,ki1)
     $           - 0.50d+00 * (  u(2,i,j,ki1) ** 2
     $                         + u(3,i,j,ki1) ** 2
     $                         + u(4,i,j,ki1) ** 2 )
     $                        / u(1,i,j,ki1) )
c     
            phi2(i,j) = c2*(  u(5,i,j,ki2)
     $           - 0.50d+00 * (  u(2,i,j,ki2) ** 2
     $                         + u(3,i,j,ki2) ** 2
     $                         + u(4,i,j,ki2) ** 2 )
     $                        / u(1,i,j,ki2) )
c
         end do
c
      end do
c
      frc1 = 0.0d+00
c
      do j = ji1, ji2-1
c
         do i = ii1, ii2-1
c
            frc1 = frc1 + (  phi1(i,j)
     $                     + phi1(i+1,j)
     $                     + phi1(i,j+1)
     $                     + phi1(i+1,j+1)
     $                     + phi2(i,j)
     $                     + phi2(i+1,j)
     $                     + phi2(i,j+1)
     $                     + phi2(i+1,j+1) )
c
         end do
c
      end do
c
      frc1 = dxi * deta * frc1
c
      do k = ki1, ki2
c     
         do i = ii1, ii2
c
            phi1(i,k) = c2*(  u(5,i,ji1,k)
     $           - 0.50d+00 * (  u(2,i,ji1,k) ** 2
     $                         + u(3,i,ji1,k) ** 2
     $                         + u(4,i,ji1,k) ** 2 )
     $                        / u(1,i,ji1,k) )
c     
            phi2(i,k) = c2*(  u(5,i,ji2,k)
     $           - 0.50d+00 * (  u(2,i,ji2,k) ** 2
     $                         + u(3,i,ji2,k) ** 2
     $                         + u(4,i,ji2,k) ** 2 )
     $                        / u(1,i,ji2,k) )
c
         end do
c
      end do
c
      frc2 = 0.0d+00
c
      do k = ki1, ki2-1
c
         do i = ii1, ii2-1
c
            frc2 = frc2 + (  phi1(i,k)
     $                     + phi1(i+1,k)
     $                     + phi1(i,k+1)
     $                     + phi1(i+1,k+1)
     $                     + phi2(i,k)
     $                     + phi2(i+1,k)
     $                     + phi2(i,k+1)
     $                     + phi2(i+1,k+1) )
c
         end do
c
      end do
c
      frc2 = dxi * dzeta * frc2
c
      do k = ki1, ki2
c     
         do j = ji1, ji2
c
            phi1(j,k) = c2*(  u(5,ii1,j,k)
     $           - 0.50d+00 * (  u(2,ii1,j,k) ** 2
     $                         + u(3,ii1,j,k) ** 2
     $                         + u(4,ii1,j,k) ** 2 )
     $                        / u(1,ii1,j,k) )
c     
            phi2(j,k) = c2*(  u(5,ii2,j,k)
     $           - 0.50d+00 * (  u(2,ii2,j,k) ** 2
     $                         + u(3,ii2,j,k) ** 2
     $                         + u(4,ii2,j,k) ** 2 )
     $                        / u(1,ii2,j,k) )
c
         end do
c
      end do
c
      frc3 = 0.0d+00
c
      do k = ki1, ki2-1
c
         do j = ji1, ji2-1
c
            frc3 = frc3 + (  phi1(j,k)
     $                     + phi1(j+1,k)
     $                     + phi1(j,k+1)
     $                     + phi1(j+1,k+1)
     $                     + phi2(j,k)
     $                     + phi2(j+1,k)
     $                     + phi2(j,k+1)
     $                     + phi2(j+1,k+1) )
c
         end do
c
      end do
c
      frc3 = deta * dzeta * frc3
c
      frc = 0.25d+00 * ( frc1 + frc2 + frc3 )
c
      write (iout,1001) frc
c
      return
c
 1001 format (//5x,'surface integral = ',1pe12.5//)
c
      end
c
c
c
c
c
      subroutine rhs
c
c***compute the right hand sides
c
c Author: Sisira Weeratunga
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension flux(5,isiz1)
c
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  rsd(m,i,j,k) = - frct(m,i,j,k)
               end do
            end do
         end do
      end do
c
c***xi-direction flux differences
c
      do k = 2, nz - 1
c
         do j = 2, ny - 1
c 
            do i = 1, nx
c
               flux(1,i) = u(2,i,j,k)
c
               u21 = u(2,i,j,k) / u(1,i,j,k)
c
               q = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
     $                         + u(3,i,j,k) * u(3,i,j,k)
     $                         + u(4,i,j,k) * u(4,i,j,k) )
     $                      / u(1,i,j,k)
c
               flux(2,i) = u(2,i,j,k) * u21 + c2 * ( u(5,i,j,k) - q )
c
               flux(3,i) = u(3,i,j,k) * u21
c
               flux(4,i) = u(4,i,j,k) * u21
c
               flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
c    
            end do
c
            do i = 2, nx - 1
c
               do m = 1, 5
c
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     $                       - tx2 * ( flux(m,i+1) - flux(m,i-1) )
c
               end do
c
            end do
c
            do i = 2, nx
c
               tmp = 1.0d+00 / u(1,i,j,k)
c
               u21i = tmp * u(2,i,j,k)
               u31i = tmp * u(3,i,j,k)
               u41i = tmp * u(4,i,j,k)
               u51i = tmp * u(5,i,j,k)
c
               tmp = 1.0d+00 / u(1,i-1,j,k)
c
               u21im1 = tmp * u(2,i-1,j,k)
               u31im1 = tmp * u(3,i-1,j,k)
               u41im1 = tmp * u(4,i-1,j,k)
               u51im1 = tmp * u(5,i-1,j,k)
c
               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * ( u21i - u21im1 )
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     $                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * tx3 * ( u21i**2 - u21im1**2 )
     $              + c1 * c5 * tx3 * ( u51i - u51im1 )
c
            end do
c
            do i = 2, nx - 1
c
               rsd(1,i,j,k) = rsd(1,i,j,k)
     $              + dx1 * tx1 * (            u(1,i-1,j,k)
     $                             - 2.0d+00 * u(1,i,j,k)
     $                             +           u(1,i+1,j,k) )
c
               rsd(2,i,j,k) = rsd(2,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     $              + dx2 * tx1 * (            u(2,i-1,j,k)
     $                             - 2.0d+00 * u(2,i,j,k)
     $                             +           u(2,i+1,j,k) )
c
               rsd(3,i,j,k) = rsd(3,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     $              + dx3 * tx1 * (            u(3,i-1,j,k)
     $                             - 2.0d+00 * u(3,i,j,k)
     $                             +           u(3,i+1,j,k) )
c
               rsd(4,i,j,k) = rsd(4,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     $              + dx4 * tx1 * (            u(4,i-1,j,k)
     $                             - 2.0d+00 * u(4,i,j,k)
     $                             +           u(4,i+1,j,k) )
c
               rsd(5,i,j,k) = rsd(5,i,j,k)
     $              + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     $              + dx5 * tx1 * (            u(5,i-1,j,k)
     $                             - 2.0d+00 * u(5,i,j,k)
     $                             +           u(5,i+1,j,k) )
c
            end do
c
c***Fourth-order dissipation
c
            do m = 1, 5
c
               rsd(m,2,j,k) = rsd(m,2,j,k)
     $           - dssp * ( + 5.0d+00 * u(m,2,j,k)
     $                      - 4.0d+00 * u(m,3,j,k)
     $                      +           u(m,4,j,k) )
c
               rsd(m,3,j,k) = rsd(m,3,j,k)
     $           - dssp * ( - 4.0d+00 * u(m,2,j,k)
     $                      + 6.0d+00 * u(m,3,j,k)
     $                      - 4.0d+00 * u(m,4,j,k)
     $                      +           u(m,5,j,k) )
c
            end do
c
            do i = 4, nx - 3
c
               do m = 1, 5
c
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     $              - dssp * (            u(m,i-2,j,k)
     $                        - 4.0d+00 * u(m,i-1,j,k)
     $                        + 6.0d+00 * u(m,i,j,k)
     $                        - 4.0d+00 * u(m,i+1,j,k)
     $                        +           u(m,i+2,j,k) )
c
               end do
c
            end do
c
            do m = 1, 5
c
               rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
     $           - dssp * (             u(m,nx-4,j,k)
     $                      - 4.0d+00 * u(m,nx-3,j,k)
     $                      + 6.0d+00 * u(m,nx-2,j,k)
     $                      - 4.0d+00 * u(m,nx-1,j,k)  )
c
               rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
     $           - dssp * (             u(m,nx-3,j,k)
     $                      - 4.0d+00 * u(m,nx-2,j,k)
     $                      + 5.0d+00 * u(m,nx-1,j,k) )
c
            end do
c
         end do
c
      end do
c
c***eta-direction flux differences
c
      do k = 2, nz - 1
c
         do i = 2, nx - 1
c 
            do j = 1, ny
c
               flux(1,j) = u(3,i,j,k)
c
               u31 = u(3,i,j,k) / u(1,i,j,k)
c
               q = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
     $                         + u(3,i,j,k) * u(3,i,j,k)
     $                         + u(4,i,j,k) * u(4,i,j,k) )
     $                      / u(1,i,j,k)
c
               flux(2,j) = u(2,i,j,k) * u31 
c
               flux(3,j) = u(3,i,j,k) * u31 + c2 * ( u(5,i,j,k) - q )
c
               flux(4,j) = u(4,i,j,k) * u31
c
               flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
c    
            end do
c
            do j = 2, ny - 1
c
               do m = 1, 5
c
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     $                       - ty2 * ( flux(m,j+1) - flux(m,j-1) )
c
               end do
c
            end do
c
            do j = 2, ny
c
               tmp = 1.0d+00 / u(1,i,j,k)
c
               u21j = tmp * u(2,i,j,k)
               u31j = tmp * u(3,i,j,k)
               u41j = tmp * u(4,i,j,k)
               u51j = tmp * u(5,i,j,k)
c
               tmp = 1.0d+00 / u(1,i,j-1,k)
c
               u21jm1 = tmp * u(2,i,j-1,k)
               u31jm1 = tmp * u(3,i,j-1,k)
               u41jm1 = tmp * u(4,i,j-1,k)
               u51jm1 = tmp * u(5,i,j-1,k)
c
               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * ( u31j - u31jm1 )
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     $                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * ty3 * ( u31j**2 - u31jm1**2 )
     $              + c1 * c5 * ty3 * ( u51j - u51jm1 )
c
            end do
c
            do j = 2, ny - 1
c
               rsd(1,i,j,k) = rsd(1,i,j,k)
     $              + dy1 * ty1 * (            u(1,i,j-1,k)
     $                             - 2.0d+00 * u(1,i,j,k)
     $                             +           u(1,i,j+1,k) )
c
               rsd(2,i,j,k) = rsd(2,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     $              + dy2 * ty1 * (            u(2,i,j-1,k)
     $                             - 2.0d+00 * u(2,i,j,k)
     $                             +           u(2,i,j+1,k) )
c
               rsd(3,i,j,k) = rsd(3,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     $              + dy3 * ty1 * (            u(3,i,j-1,k)
     $                             - 2.0d+00 * u(3,i,j,k)
     $                             +           u(3,i,j+1,k) )
c
               rsd(4,i,j,k) = rsd(4,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     $              + dy4 * ty1 * (            u(4,i,j-1,k)
     $                             - 2.0d+00 * u(4,i,j,k)
     $                             +           u(4,i,j+1,k) )
c
               rsd(5,i,j,k) = rsd(5,i,j,k)
     $              + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     $              + dy5 * ty1 * (            u(5,i,j-1,k)
     $                             - 2.0d+00 * u(5,i,j,k)
     $                             +           u(5,i,j+1,k) )
c
            end do
c
c***fourth-order dissipation
c
            do m = 1, 5
c
               rsd(m,i,2,k) = rsd(m,i,2,k)
     $           - dssp * ( + 5.0d+00 * u(m,i,2,k)
     $                      - 4.0d+00 * u(m,i,3,k)
     $                      +           u(m,i,4,k) )
c
               rsd(m,i,3,k) = rsd(m,i,3,k)
     $           - dssp * ( - 4.0d+00 * u(m,i,2,k)
     $                      + 6.0d+00 * u(m,i,3,k)
     $                      - 4.0d+00 * u(m,i,4,k)
     $                      +           u(m,i,5,k) )
c
            end do
c
            do j = 4, ny - 3
c
               do m = 1, 5
c
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     $              - dssp * (            u(m,i,j-2,k)
     $                        - 4.0d+00 * u(m,i,j-1,k)
     $                        + 6.0d+00 * u(m,i,j,k)
     $                        - 4.0d+00 * u(m,i,j+1,k)
     $                        +           u(m,i,j+2,k) )
c
               end do
c
            end do
c
            do m = 1, 5
c
               rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
     $           - dssp * (             u(m,i,ny-4,k)
     $                      - 4.0d+00 * u(m,i,ny-3,k)
     $                      + 6.0d+00 * u(m,i,ny-2,k)
     $                      - 4.0d+00 * u(m,i,ny-1,k)  )
c
               rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
     $           - dssp * (             u(m,i,ny-3,k)
     $                      - 4.0d+00 * u(m,i,ny-2,k)
     $                      + 5.0d+00 * u(m,i,ny-1,k) )
c
            end do
c
         end do
c
      end do
c
c***zeta-direction flux differences
c
      do j = 2, ny - 1
c
         do i = 2, nx - 1
c
            do k = 1, nz
c
               flux(1,k) = u(4,i,j,k)
c
               u41 = u(4,i,j,k) / u(1,i,j,k)
c
               q = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
     $                         + u(3,i,j,k) * u(3,i,j,k)
     $                         + u(4,i,j,k) * u(4,i,j,k) )
     $                      / u(1,i,j,k)
c
               flux(2,k) = u(2,i,j,k) * u41 
c
               flux(3,k) = u(3,i,j,k) * u41 
c
               flux(4,k) = u(4,i,j,k) * u41 + c2 * ( u(5,i,j,k) - q )
c
               flux(5,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
c    
            end do
c
            do k = 2, nz - 1
c
               do m = 1, 5
c
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     $                       - tz2 * ( flux(m,k+1) - flux(m,k-1) )
c
               end do
c
            end do
c
            do k = 2, nz
c
               tmp = 1.0d+00 / u(1,i,j,k)
c
               u21k = tmp * u(2,i,j,k)
               u31k = tmp * u(3,i,j,k)
               u41k = tmp * u(4,i,j,k)
               u51k = tmp * u(5,i,j,k)
c
               tmp = 1.0d+00 / u(1,i,j,k-1)
c
               u21km1 = tmp * u(2,i,j,k-1)
               u31km1 = tmp * u(3,i,j,k-1)
               u41km1 = tmp * u(4,i,j,k-1)
               u51km1 = tmp * u(5,i,j,k-1)
c
               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * ( u41k - u41km1 )
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     $                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * tz3 * ( u41k**2 - u41km1**2 )
     $              + c1 * c5 * tz3 * ( u51k - u51km1 )
c
            end do
c
            do k = 2, nz - 1
c
               rsd(1,i,j,k) = rsd(1,i,j,k)
     $              + dz1 * tz1 * (            u(1,i,j,k-1)
     $                             - 2.0d+00 * u(1,i,j,k)
     $                             +           u(1,i,j,k+1) )
c
               rsd(2,i,j,k) = rsd(2,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     $              + dz2 * tz1 * (            u(2,i,j,k-1)
     $                             - 2.0d+00 * u(2,i,j,k)
     $                             +           u(2,i,j,k+1) )
c
               rsd(3,i,j,k) = rsd(3,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     $              + dz3 * tz1 * (            u(3,i,j,k-1)
     $                             - 2.0d+00 * u(3,i,j,k)
     $                             +           u(3,i,j,k+1) )
c
               rsd(4,i,j,k) = rsd(4,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     $              + dz4 * tz1 * (            u(4,i,j,k-1)
     $                             - 2.0d+00 * u(4,i,j,k)
     $                             +           u(4,i,j,k+1) )
c
               rsd(5,i,j,k) = rsd(5,i,j,k)
     $              + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     $              + dz5 * tz1 * (            u(5,i,j,k-1)
     $                             - 2.0d+00 * u(5,i,j,k)
     $                             +           u(5,i,j,k+1) )
c
            end do
c
c***fourth-order dissipation
c
            do m = 1, 5
c
               rsd(m,i,j,2) = rsd(m,i,j,2)
     $           - dssp * ( + 5.0d+00 * u(m,i,j,2)
     $                      - 4.0d+00 * u(m,i,j,3)
     $                      +           u(m,i,j,4) )
c
               rsd(m,i,j,3) = rsd(m,i,j,3)
     $           - dssp * ( - 4.0d+00 * u(m,i,j,2)
     $                      + 6.0d+00 * u(m,i,j,3)
     $                      - 4.0d+00 * u(m,i,j,4)
     $                      +           u(m,i,j,5) )
c
            end do
c
            do k = 4, nz - 3
c
               do m = 1, 5
c
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     $              - dssp * (            u(m,i,j,k-2)
     $                        - 4.0d+00 * u(m,i,j,k-1)
     $                        + 6.0d+00 * u(m,i,j,k)
     $                        - 4.0d+00 * u(m,i,j,k+1)
     $                        +           u(m,i,j,k+2) )
c
               end do
c
            end do
c
            do m = 1, 5
c
               rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
     $           - dssp * (             u(m,i,j,nz-4)
     $                      - 4.0d+00 * u(m,i,j,nz-3)
     $                      + 6.0d+00 * u(m,i,j,nz-2)
     $                      - 4.0d+00 * u(m,i,j,nz-1)  )
c
               rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
     $           - dssp * (             u(m,i,j,nz-3)
     $                      - 4.0d+00 * u(m,i,j,nz-2)
     $                      + 5.0d+00 * u(m,i,j,nz-1) )
c
            end do
c
         end do
c
      end do
c
      return
      end
c
c
c
c
c
      subroutine verify ( xcr, xce, xci )
c
c***verification routine
c
c Author: Sisira Weeratunga 
c         NASA Ames Research Center
c         (10/25/90)
c
      include 'appbt.incl'
c
      dimension xcr(5), xce(5),
     $          xrr(5), xre(5)
c
c***tolerance level
c
      epsilon = 1.0e-08
c
      if ( ( nx .eq. 12 ) .and.
     $     ( ny .eq. 12 ) .and.
     $     ( nz .eq. 12 ) ) then
c
c***Reference values of RMS-norms of residual, for the (12X12X12) grid,
c   after 60 time steps, with  DT = 1.0e-02
c
         xrr(1) = 1.7034283709541311d-01
         xrr(2) = 1.2975252070034097d-02
         xrr(3) = 3.2527926989486055d-02
         xrr(4) = 2.6436421275166801d-02
         xrr(5) = 1.9211784131744430d-01
c
c***Reference values of RMS-norms of solution error, for the (12X12X12) grid,
c   after 60 time steps, with  DT = 1.0e-02
c
         xre(1) = 4.9976913345811579d-04
         xre(2) = 4.5195666782961927d-05
         xre(3) = 7.3973765172921357d-05
         xre(4) = 7.3821238632439731d-05
         xre(5) = 8.9269630987491446d-04
c
c***Reference value of surface integral, for the (12X12X12) grid,
c   after 60 time steps, with  DT = 1.0e-02
c
         xri = 7.8414305521809488e+00
c
c***verification test for residuals
c
         do m = 1, 5
c     
            tmp = abs ( ( xcr(m) - xrr(m) ) / xrr(m) )
c     
            if ( tmp .gt. epsilon ) then
c     
               write (iout,1001) 
 1001          format(/5x,'VERIFICATION TEST FOR RESIDUALS FAILED')
c     
               go to 100
c
            end if
c     
         end do
c     
         write (iout,1002) 
 1002    format (/5x,'VERIFICATION TEST FOR RESIDUALS ',
     $        'IS SUCCESSFUL')
c     
c***verification test for solution error
c
 100     continue
c     
         do m = 1, 5
c
            tmp = abs ( ( xce(m) - xre(m) ) / xre(m) )
c     
            if ( tmp .gt. epsilon ) then
c     
               write (iout,1003) 
 1003          format(/5x,'VERIFICATION TEST FOR SOLUTION ',
     $              'ERRORS FAILED')
c     
               go to 200
c     
            end if
c     
         end do
c     
         write (iout,1004) 
 1004    format (/5x,'VERIFICATION TEST FOR SOLUTION ERRORS ',
     $        'IS SUCCESSFUL')
c
c***verification test for surface integral
c
 200     continue
c
         tmp = abs ( ( xci - xri ) / xri )
c     
         if ( tmp .gt. epsilon ) then
c     
            write (iout,1005)
 1005       format (/5x,'VERIFICATION TEST FOR SURFACE INTEGRAL FAILED')
c     
         else
c     
            write (iout,1006)
 1006       format (/5x,'VERIFICATION TEST FOR SURFACE INTEGRAL ',
     $           'IS SUCCESSFUL')
c     
         end if
c     
         write (iout,1007)
 1007    format(//10x,'CAUTION',
     $ //5x,'REFERENCE VALUES CURRENTLY IN THIS VERIFICATION ',
     $     'ROUTINE ',
     $ /5x,'ARE VALID ONLY FOR RUNS WITH THE FOLLOWING PARAMETER ',
     $     'VALUES:',
     $ //5x,'NX = 12;  NY = 12;  NZ = 12 ',
     $ //5x,'ITMAX = 60',
     $ //5x,'DT = 1.0e-02',
     $ //5x,'CHANGE IN ANY OF THE ABOVE VALUES RENDER THE REFERENCE ',
     $     'VALUES ',
     $ /5x,'INVALID AND CAUSES A FAILURE OF THE VERIFICATION TEST.')
c
      else if ( ( nx .eq. 64 ) .and.
     $          ( ny .eq. 64 ) .and.
     $          ( nz .eq. 64 ) ) then
c
c***Reference values of RMS-norms of residual, for the (64X64X64) grid,
c   after 200 time steps, with  DT = 1.0e-03
c
         xrr(1) = 1.0806346714637264d+02
         xrr(2) = 1.1319730901220813d+01
         xrr(3) = 2.5974354511582465d+01
         xrr(4) = 2.3665622544678910d+01
         xrr(5) = 2.5278963211748344d+02
c     
c***Reference values of RMS-norms of solution error, for the (64X64X64) grid,
c   after 200 time steps, with  DT = 1.0e-03
c
         xre(1) = 4.2348416040525025d+00
         xre(2) = 4.4390282496995698d-01
         xre(3) = 9.6692480136345650d-01
         xre(4) = 8.8302063039765474d-01
         xre(5) = 9.7379901770829278d+00
c
c***Reference value of surface integral, for the (64X64X64) grid,
c   after 200 time steps, with DT = 1.0e-03
c
         xri = 1.4101360012200832d+01
c
c***verification test for residuals
c
         do m = 1, 5
c     
            tmp = abs ( ( xcr(m) - xrr(m) ) / xrr(m) )
c     
            if ( tmp .gt. epsilon ) then
c     
               write (iout,1001) 
               go to 400
c
            end if
c     
         end do
c     
         write (iout,1002) 
c     
c***verification test for solution error
c
 400     continue
c     
         do m = 1, 5
c
            tmp = abs ( ( xce(m) - xre(m) ) / xre(m) )
c     
            if ( tmp .gt. epsilon ) then
c     
               write (iout,1003) 
               go to 500
c     
            end if
c     
         end do
c     
         write (iout,1004) 
c
c***verification test for surface integral
c
 500     continue
c
         tmp = abs ( ( xci - xri ) / xri )
c     
         if ( tmp .gt. epsilon ) then
c     
            write (iout,1005)
c     
         else
c     
            write (iout,1006)
c     
         end if
c     
         write (iout,1008)
 1008    format(//10x,'CAUTION',
     $ //5x,'REFERENCE VALUES CURRENTLY IN THIS VERIFICATION ',
     $     'ROUTINE ',
     $ /5x,'ARE VALID ONLY FOR RUNS WITH THE FOLLOWING PARAMETER ',
     $     'VALUES:',
     $ //5x,'NX = 64;  NY = 64;  NZ = 64 ',
     $ //5x,'ITMAX = 200',
     $ //5x,'DT = 8.0e-04',
     $ //5x,'CHANGE IN ANY OF THE ABOVE VALUES RENDER THE REFERENCE ',
     $     'VALUES ',
     $ /5x,'INVALID AND CAUSES A FAILURE OF THE VERIFICATION TEST.')
c
      else
c
         write (iout,1009) 
 1009    format (//1x,'FOR THE PROBLEM PARAMETERS IN USE ',
     $        'NO REFERENCE VALUES ARE PROVIDED'
     $        /1x,'IN THE CURRENT VERIFICATION ROUTINE - ',
     $        'NO VERIFIACTION TEST WAS PERFORMED')
c
      end if
c
      return
      end
      function timer(x)
      real x, timer
      real secs,tarray(2)
c
c   
c     sun version of timing routine
cccccccccccccccccccccccccccccccccccccc
      secs = etime(tarray)
      timer = secs
cccccccccccccccccccccccccccccccccccccc
c
c       cray version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      timer = second ()
cccccccccccccccccccccccccccccccccccccc
c
c       sgi version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      secs = secnds(0.0)
      timer = secs
cccccccccccccccccccccccccccccccccccccc
      return
      end
