c
      subroutine adi
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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c
c
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
      dimension idmax(5), jdmax(5), kdmax(5),
     $          imax(5), jmax(5), kmax(5),
     $          delunm(5)
c
      dimension flux1(isiz1)
      dimension flux2(isiz1)
      dimension flux3(isiz1)
      dimension flux4(isiz1)
      dimension flux5(isiz1)
c
      dimension cv(isiz3), aa(isiz3), rhos(isiz3)
      dimension  rhon(isiz1)
      dimension  rhoq(isiz2)
c
      parameter ( one = 1.0d+00 )
c
      lnorm = 2
      bt = sqrt ( 0.50d+00 )
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c
		call rhs
c
c
c***compute the norms of the residuals
c
c     
c        call l2norm ( isiz1, isiz2, isiz3,
c    $                 nx, ny, nz,
c    $                 rsd, rsdnm )
c
c
         sum1 = 0.0d+00
         sum2 = 0.0d+00
         sum3 = 0.0d+00
         sum4 = 0.0d+00
         sum5 = 0.0d+00
c
c
CHPF$ INDEPENDENT
      do k = 2, nz-1
         do j = 2, ny-1
            do i = 2, nx-1
c
c
                  sum1 = sum1 + rsd(1,i,j,k) * rsd(1,i,j,k)
                  sum2 = sum2 + rsd(2,i,j,k) * rsd(2,i,j,k)
                  sum3 = sum3 + rsd(3,i,j,k) * rsd(3,i,j,k)
                  sum4 = sum4 + rsd(4,i,j,k) * rsd(4,i,j,k)
                  sum5 = sum5 + rsd(5,i,j,k) * rsd(5,i,j,k)
c
c
            end do
         end do
      end do
c
c
         rsdnm(1) = sqrt ( sum1 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(2) = sqrt ( sum2 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(3) = sqrt ( sum3 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(4) = sqrt ( sum4 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(5) = sqrt ( sum5 / ( (nx-2) * (ny-2) * (nz-2) ) )
c
c
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
c***begin pseudo-time stepping iterations
c
      tstart = timer( )
c
      do istep = 1, itmax
c
c
         if ( ( mod ( istep, inorm ) .eq. 0 ) .and.
     $        ( ipr .eq. 1 ) ) then
c
            write ( iout, 1001 ) istep
c
         end if
c
CHPF$ INDEPENDENT
         do k = 1, nz
            do j = 1, ny
               do i = 1, nx
c
                     rsd(1,i,j,k) = dt * rsd(1,i,j,k)
                     rsd(2,i,j,k) = dt * rsd(2,i,j,k)
                     rsd(3,i,j,k) = dt * rsd(3,i,j,k)
                     rsd(4,i,j,k) = dt * rsd(4,i,j,k)
                     rsd(5,i,j,k) = dt * rsd(5,i,j,k)
c
               end do
            end do
         end do
c
c***perform 3-factor, scalar ADI iterations
c
c***perform the block diagonal inversion
c
c        call txinvr
c
CHPF$ INDEPENDENT
CAPR$ IGNORE ALL INHIBITORS
      do k = 2, nz - 1
c
         do j = 2, ny - 1
c
            do i = 2, nx - 1
c
               ru1 = 1.0d+00 / u(1,i,j,k)
c
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               ac2 = c1 * c2 * ( ru1 * u(5,i,j,k) - q )
c
               if ( ac2 .le. 0.0d+00 ) then
                  write (*,*) 'Speed of sound is zero'
                  stop
               end if
c
               ac = sqrt ( ac2 )
c
               r1 = rsd(1,i,j,k)
               r2 = rsd(2,i,j,k)
               r3 = rsd(3,i,j,k)
               r4 = rsd(4,i,j,k)
               r5 = rsd(5,i,j,k)
c
               t1 = ( c2 / ac2 ) * ( q * r1 - uu * r2
     $                                      - vv * r3
     $                                      - ww * r4
     $                                      +      r5 )
c
               rsd(1,i,j,k) = r1 - t1 
c
               rsd(2,i,j,k) = - ru1 * ( ww * r1 - r4 )
c
               rsd(3,i,j,k) =   ru1 * ( vv * r1 - r3 )
c
               t2 = bt * ru1 * ( uu * r1 - r2 )
               t3 = ( bt * ru1 * ac ) * t1
c
               rsd(4,i,j,k) = - t2 + t3
c
               rsd(5,i,j,k) =   t2 + t3
c
            end do
			end do
			end do
c
         call xisweep
c
c***perform the block diagonal inversion
c
c        call ninvr
c
CHPF$ INDEPENDENT
        do k = 2, nz -1
		  do j = 2, ny -1
            do i = 2, nx - 1
c
               r1 = rsd(1,i,j,k)
               r2 = rsd(2,i,j,k)
               r3 = rsd(3,i,j,k)
               r4 = rsd(4,i,j,k)
               r5 = rsd(5,i,j,k)
c
               rsd(1,i,j,k) = - r2
c
               rsd(2,i,j,k) =   r1
c
               rsd(3,i,j,k) = bt * ( r4 - r5 )
c
               t1 = bt * r3
c     
               t2 = 0.50d+00 * ( r4 + r5 )
c
               rsd(4,i,j,k) = - t1 + t2
c
               rsd(5,i,j,k) =   t1 + t2
c
            end do
c
         end do
c
      end do
c
       call etasweep
c
c***perform the block diagonal inversion
c
c        call pinvr
c
CHPF$ INDEPENDENT
        do k = 2, nz -1
            do j = 2, ny - 1
              do i = 2, nx -1
c
               r1 = rsd(1,i,j,k)
               r2 = rsd(2,i,j,k)
               r3 = rsd(3,i,j,k)
               r4 = rsd(4,i,j,k)
               r5 = rsd(5,i,j,k)
c
               rsd(1,i,j,k) = bt * ( r4 - r5 )
c
               rsd(2,i,j,k) = - r3
c
               rsd(3,i,j,k) = r2
c
               t1 = bt * r1
c     
               t2 = 0.50d+00 * ( r4 + r5 )
c
               rsd(4,i,j,k) = - t1 + t2
c
               rsd(5,i,j,k) =   t1 + t2
c
            end do
c
         end do
c
      end do
c
	  call ztasweep
c
c     perform the block diagonal matrix-vector multiplication
c
c        call tzetar
c
CHPF$ INDEPENDENT
      do k = 2, nz - 1
c
         do j = 2, ny - 1
c
            do i = 2, nx - 1
c
               ru1 = 1.0d+00 / u(1,i,j,k)
c
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               ac2 = c1 * c2 * ( ru1 * u(5,i,j,k) - q )
c
               ac = sqrt ( ac2 )
c
               alph = ( bt * u(1,i,j,k) ) / ac
c
               r1 = rsd(1,i,j,k)
               r2 = rsd(2,i,j,k)
               r3 = rsd(3,i,j,k)
               r4 = rsd(4,i,j,k)
               r5 = rsd(5,i,j,k)
c
               t1 = alph * ( r4 + r5 )
               t2 = r3 + t1
               t3 = bt * u(1,i,j,k) * ( r4 - r5 )
c
               rsd(1,i,j,k) = t2
c
               rsd(2,i,j,k) = - u(1,i,j,k) * r2 + uu * t2
c
               rsd(3,i,j,k) =   u(1,i,j,k) * r1 + vv * t2
c
               rsd(4,i,j,k) =   ww * t2 + t3
c
               rsd(5,i,j,k) =  u(1,i,j,k) * ( - uu * r2 + vv * r1 )
     $                       + q * t2
     $                       + ( ac2 / c2 ) * t1
     $                       + ww * t3
c
c
c***update the variables
c
c
                     u(1,  i, j, k ) =  u(1,  i, j, k )
     $                                + rsd(1,  i, j, k )
                     u(2,  i, j, k ) =  u(2,  i, j, k )
     $                                + rsd(2,  i, j, k )
                     u(3,  i, j, k ) =  u(3,  i, j, k )
     $                                + rsd(3,  i, j, k )
                     u(4,  i, j, k ) =  u(4,  i, j, k )
     $                                + rsd(4,  i, j, k )
                     u(5,  i, j, k ) =  u(5,  i, j, k )
     $                                + rsd(5,  i, j, k )
c
               end do
            end do
         end do
c
c***compute the norms of pseudo-time iteration corrections
c     
         if ( mod ( istep, inorm ) .eq. 0 ) then
c
c              call l2norm ( isiz1, isiz2, isiz3,
c    $                       nx, ny, nz,
c    $                       rsd, delunm )
c
c
         sum1 = 0.0d+00
         sum2 = 0.0d+00
         sum3 = 0.0d+00
         sum4 = 0.0d+00
         sum5 = 0.0d+00
c
c
CHPF$ INDEPENDENT
      do k = 2, nz-1
         do j = 2, ny-1
            do i = 2, nx-1
c
c
                  sum1 = sum1 + rsd(1,i,j,k) * rsd(1,i,j,k)
                  sum2 = sum2 + rsd(2,i,j,k) * rsd(2,i,j,k)
                  sum3 = sum3 + rsd(3,i,j,k) * rsd(3,i,j,k)
                  sum4 = sum4 + rsd(4,i,j,k) * rsd(4,i,j,k)
                  sum5 = sum5 + rsd(5,i,j,k) * rsd(5,i,j,k)
c
c
            end do
         end do
      end do
c
c
         delunm(1) = sqrt ( sum1 / ( (nx-2) * (ny-2) * (nz-2) ) )
         delunm(2) = sqrt ( sum2 / ( (nx-2) * (ny-2) * (nz-2) ) )
         delunm(3) = sqrt ( sum3 / ( (nx-2) * (ny-2) * (nz-2) ) )
         delunm(4) = sqrt ( sum4 / ( (nx-2) * (ny-2) * (nz-2) ) )
         delunm(5) = sqrt ( sum5 / ( (nx-2) * (ny-2) * (nz-2) ) )
c
c
               if ( ipr .eq. 1 ) then
c
                  write (iout,1006) ( delunm(m), m = 1, 5 )
c
               end if
c
c
         end if
c
c***compute the steady-state residuals
c
         call rhs
c
c
c***compute the norms of the residuals
c
         if (  istep .eq. itmax  ) then
c
c              call l2norm ( isiz1, isiz2, isiz3,
c    $                       nx, ny, nz,
c    $                       rsd, rsdnm )
c
c
         sum1 = 0.0d+00
         sum2 = 0.0d+00
         sum3 = 0.0d+00
         sum4 = 0.0d+00
         sum5 = 0.0d+00
c
c
CHPF$ INDEPENDENT
      do k = 2, nz-1
         do j = 2, ny-1
            do i = 2, nx-1
c
c
                  sum1 = sum1 + rsd(1,i,j,k) * rsd(1,i,j,k)
                  sum2 = sum2 + rsd(2,i,j,k) * rsd(2,i,j,k)
                  sum3 = sum3 + rsd(3,i,j,k) * rsd(3,i,j,k)
                  sum4 = sum4 + rsd(4,i,j,k) * rsd(4,i,j,k)
                  sum5 = sum5 + rsd(5,i,j,k) * rsd(5,i,j,k)
c
c
            end do
         end do
      end do
c
c
         rsdnm(1) = sqrt ( sum1 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(2) = sqrt ( sum2 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(3) = sqrt ( sum3 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(4) = sqrt ( sum4 / ( (nx-2) * (ny-2) * (nz-2) ) )
         rsdnm(5) = sqrt ( sum5 / ( (nx-2) * (ny-2) * (nz-2) ) )
c
c
               if ( ipr .eq. 1 )  then
c
                  write (iout,1007) ( rsdnm(m), m = 1, 5 )
c     
               end if
c
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 Scalar ADI iteration no.=',i4/)
c
 1002 format (1x/1x,'max-norm of Scalar ADI-iteration correction ',
     $ 'for first pde  = ',1pe12.5/,
     $ 59x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Scalar ADI-iteration correction ',
     $ 'for second pde = ',1pe12.5/,
     $ 59x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Scalar ADI-iteration correction ',
     $ 'for third pde  = ',1pe12.5/,
     $ 59x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Scalar ADI-iteration correction ',
     $ 'for fourth pde = ',1pe12.5/,
     $ 59x,'(',i4,',',i4,',',i4,')'/,
     $ 1x,'max-norm of Scalar ADI-iteration correction ',
     $ 'for fifth pde  = ',1pe12.5/,
     $ 59x,'(',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
 1005 format (/5x,'CPU time = ',1pe12.4,' Sec.')
c
 1006 format (1x/1x,'RMS-norm of scalar adi-iteration correction ',
     $ 'for first pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of scalar adi-iteration correction ',
     $ 'for second pde = ',1pe12.5/,
     $ 1x,'RMS-norm of scalar adi-iteration correction ',
     $ 'for third pde  = ',1pe12.5/,
     $ 1x,'RMS-norm of scalar adi-iteration correction ',
     $ 'for fourth pde = ',1pe12.5/,
     $ 1x,'RMS-norm of scalar 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
      program appsp
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 
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
c***open file for input data
c
      open (unit=5,file='inputspB.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='appsp.out',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***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 scalar approximate factorization iterations
c
      call adi
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 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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
      dimension flux(5,isiz1), ue(5,isiz1)
c
      dsspm = dssp
c
CHPF$ INDEPENDENT
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
                  frct(1, i, j, k ) = 0.0d+00
                  frct(2, i, j, k ) = 0.0d+00
                  frct(3, i, j, k ) = 0.0d+00
                  frct(4, i, j, k ) = 0.0d+00
                  frct(5, i, j, k ) = 0.0d+00
            end do
         end do
      end do
c
c***xi-direction flux differences
c
CHPF$ INDEPENDENT
      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
c
                  frct(1,i,j,k) =  frct(1,i,j,k)
     $                       - tx2 * ( flux(1,i+1) - flux(1,i-1) )
                  frct(2,i,j,k) =  frct(2,i,j,k)
     $                       - tx2 * ( flux(2,i+1) - flux(2,i-1) )
                  frct(3,i,j,k) =  frct(3,i,j,k)
     $                       - tx2 * ( flux(3,i+1) - flux(3,i-1) )
                  frct(4,i,j,k) =  frct(4,i,j,k)
     $                       - tx2 * ( flux(4,i+1) - flux(4,i-1) )
                  frct(5,i,j,k) =  frct(5,i,j,k)
     $                       - tx2 * ( flux(5,i+1) - flux(5,i-1) )
c
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
c
               frct(1,2,j,k) = frct(1,2,j,k)
     $           - dsspm * ( + 5.0d+00 * ue(1,2)
     $                       - 4.0d+00 * ue(1,3)
     $                       +           ue(1,4) )
c
               frct(1,3,j,k) = frct(1,3,j,k)
     $           - dsspm * ( - 4.0d+00 * ue(1,2)
     $                       + 6.0d+00 * ue(1,3)
     $                       - 4.0d+00 * ue(1,4)
     $                       +           ue(1,5) )
c
c
               frct(2,2,j,k) = frct(2,2,j,k)
     $           - dsspm * ( + 5.0d+00 * ue(2,2)
     $                       - 4.0d+00 * ue(2,3)
     $                       +           ue(2,4) )
c
               frct(2,3,j,k) = frct(2,3,j,k)
     $           - dsspm * ( - 4.0d+00 * ue(2,2)
     $                       + 6.0d+00 * ue(2,3)
     $                       - 4.0d+00 * ue(2,4)
     $                       +           ue(2,5) )
c
c
               frct(3,2,j,k) = frct(3,2,j,k)
     $           - dsspm * ( + 5.0d+00 * ue(3,2)
     $                       - 4.0d+00 * ue(3,3)
     $                       +           ue(3,4) )
c
               frct(3,3,j,k) = frct(3,3,j,k)
     $           - dsspm * ( - 4.0d+00 * ue(3,2)
     $                       + 6.0d+00 * ue(3,3)
     $                       - 4.0d+00 * ue(3,4)
     $                       +           ue(3,5) )
c
c
               frct(4,2,j,k) = frct(4,2,j,k)
     $           - dsspm * ( + 5.0d+00 * ue(4,2)
     $                       - 4.0d+00 * ue(4,3)
     $                       +           ue(4,4) )
c
               frct(4,3,j,k) = frct(4,3,j,k)
     $           - dsspm * ( - 4.0d+00 * ue(4,2)
     $                       + 6.0d+00 * ue(4,3)
     $                       - 4.0d+00 * ue(4,4)
     $                       +           ue(4,5) )
c
c
               frct(5,2,j,k) = frct(5,2,j,k)
     $           - dsspm * ( + 5.0d+00 * ue(5,2)
     $                       - 4.0d+00 * ue(5,3)
     $                       +           ue(5,4) )
c
               frct(5,3,j,k) = frct(5,3,j,k)
     $           - dsspm * ( - 4.0d+00 * ue(5,2)
     $                       + 6.0d+00 * ue(5,3)
     $                       - 4.0d+00 * ue(5,4)
     $                       +           ue(5,5) )
c
c
            do i = 4, nx - 3
c
c
                  frct(1,i,j,k) = frct(1,i,j,k)
     $              - dsspm * (            ue(1,i-2)
     $                         - 4.0d+00 * ue(1,i-1)
     $                         + 6.0d+00 * ue(1,i)
     $                         - 4.0d+00 * ue(1,i+1)
     $                         +           ue(1,i+2) )
c
c
                  frct(2,i,j,k) = frct(2,i,j,k)
     $              - dsspm * (            ue(2,i-2)
     $                         - 4.0d+00 * ue(2,i-1)
     $                         + 6.0d+00 * ue(2,i)
     $                         - 4.0d+00 * ue(2,i+1)
     $                         +           ue(2,i+2) )
c
c
                  frct(3,i,j,k) = frct(3,i,j,k)
     $              - dsspm * (            ue(3,i-2)
     $                         - 4.0d+00 * ue(3,i-1)
     $                         + 6.0d+00 * ue(3,i)
     $                         - 4.0d+00 * ue(3,i+1)
     $                         +           ue(3,i+2) )
c
c
                  frct(4,i,j,k) = frct(4,i,j,k)
     $              - dsspm * (            ue(4,i-2)
     $                         - 4.0d+00 * ue(4,i-1)
     $                         + 6.0d+00 * ue(4,i)
     $                         - 4.0d+00 * ue(4,i+1)
     $                         +           ue(4,i+2) )
c
c
                  frct(5,i,j,k) = frct(5,i,j,k)
     $              - dsspm * (            ue(5,i-2)
     $                         - 4.0d+00 * ue(5,i-1)
     $                         + 6.0d+00 * ue(5,i)
     $                         - 4.0d+00 * ue(5,i+1)
     $                         +           ue(5,i+2) )
c
c
            end do
c
c
               frct(1,nx-2,j,k) = frct(1,nx-2,j,k)
     $           - dsspm * (             ue(1,nx-4)
     $                       - 4.0d+00 * ue(1,nx-3)
     $                       + 6.0d+00 * ue(1,nx-2)
     $                       - 4.0d+00 * ue(1,nx-1)  )
c
               frct(1,nx-1,j,k) = frct(1,nx-1,j,k)
     $           - dsspm * (             ue(1,nx-3)
     $                       - 4.0d+00 * ue(1,nx-2)
     $                       + 5.0d+00 * ue(1,nx-1) )
c
c
               frct(2,nx-2,j,k) = frct(2,nx-2,j,k)
     $           - dsspm * (             ue(2,nx-4)
     $                       - 4.0d+00 * ue(2,nx-3)
     $                       + 6.0d+00 * ue(2,nx-2)
     $                       - 4.0d+00 * ue(2,nx-1)  )
c
               frct(2,nx-1,j,k) = frct(2,nx-1,j,k)
     $           - dsspm * (             ue(2,nx-3)
     $                       - 4.0d+00 * ue(2,nx-2)
     $                       + 5.0d+00 * ue(2,nx-1) )
c
c
               frct(3,nx-2,j,k) = frct(3,nx-2,j,k)
     $           - dsspm * (             ue(3,nx-4)
     $                       - 4.0d+00 * ue(3,nx-3)
     $                       + 6.0d+00 * ue(3,nx-2)
     $                       - 4.0d+00 * ue(3,nx-1)  )
c
               frct(3,nx-1,j,k) = frct(3,nx-1,j,k)
     $           - dsspm * (             ue(3,nx-3)
     $                       - 4.0d+00 * ue(3,nx-2)
     $                       + 5.0d+00 * ue(3,nx-1) )
c
c
               frct(4,nx-2,j,k) = frct(4,nx-2,j,k)
     $           - dsspm * (             ue(4,nx-4)
     $                       - 4.0d+00 * ue(4,nx-3)
     $                       + 6.0d+00 * ue(4,nx-2)
     $                       - 4.0d+00 * ue(4,nx-1)  )
c
               frct(4,nx-1,j,k) = frct(4,nx-1,j,k)
     $           - dsspm * (             ue(4,nx-3)
     $                       - 4.0d+00 * ue(4,nx-2)
     $                       + 5.0d+00 * ue(4,nx-1) )
c
c
               frct(5,nx-2,j,k) = frct(5,nx-2,j,k)
     $           - dsspm * (             ue(5,nx-4)
     $                       - 4.0d+00 * ue(5,nx-3)
     $                       + 6.0d+00 * ue(5,nx-2)
     $                       - 4.0d+00 * ue(5,nx-1)  )
c
               frct(5,nx-1,j,k) = frct(5,nx-1,j,k)
     $           - dsspm * (             ue(5,nx-3)
     $                       - 4.0d+00 * ue(5,nx-2)
     $                       + 5.0d+00 * ue(5,nx-1) )
c
c
         end do
c
      end do
c
c***eta-direction flux differences
c
CHPF$ INDEPENDENT
      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
c
                  frct(1,i,j,k) =  frct(1,i,j,k)
     $                       - ty2 * ( flux(1,j+1) - flux(1,j-1) )
                  frct(2,i,j,k) =  frct(2,i,j,k)
     $                       - ty2 * ( flux(2,j+1) - flux(2,j-1) )
                  frct(3,i,j,k) =  frct(3,i,j,k)
     $                       - ty2 * ( flux(3,j+1) - flux(3,j-1) )
                  frct(4,i,j,k) =  frct(4,i,j,k)
     $                       - ty2 * ( flux(4,j+1) - flux(4,j-1) )
                  frct(5,i,j,k) =  frct(5,i,j,k)
     $                       - ty2 * ( flux(5,j+1) - flux(5,j-1) )
c
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
c
               frct(1,i,2,k) = frct(1,i,2,k)
     $           - dsspm * ( + 5.0d+00 * ue(1,2)
     $                       - 4.0d+00 * ue(1,3)
     $                       +           ue(1,4) )
c
               frct(1,i,3,k) = frct(1,i,3,k)
     $           - dsspm * ( - 4.0d+00 * ue(1,2)
     $                       + 6.0d+00 * ue(1,3)
     $                       - 4.0d+00 * ue(1,4)
     $                       +           ue(1,5) )
c
c
               frct(2,i,2,k) = frct(2,i,2,k)
     $           - dsspm * ( + 5.0d+00 * ue(2,2)
     $                       - 4.0d+00 * ue(2,3)
     $                       +           ue(2,4) )
c
               frct(2,i,3,k) = frct(2,i,3,k)
     $           - dsspm * ( - 4.0d+00 * ue(2,2)
     $                       + 6.0d+00 * ue(2,3)
     $                       - 4.0d+00 * ue(2,4)
     $                       +           ue(2,5) )
c
c
               frct(3,i,2,k) = frct(3,i,2,k)
     $           - dsspm * ( + 5.0d+00 * ue(3,2)
     $                       - 4.0d+00 * ue(3,3)
     $                       +           ue(3,4) )
c
               frct(3,i,3,k) = frct(3,i,3,k)
     $           - dsspm * ( - 4.0d+00 * ue(3,2)
     $                       + 6.0d+00 * ue(3,3)
     $                       - 4.0d+00 * ue(3,4)
     $                       +           ue(3,5) )
c
c
               frct(4,i,2,k) = frct(4,i,2,k)
     $           - dsspm * ( + 5.0d+00 * ue(4,2)
     $                       - 4.0d+00 * ue(4,3)
     $                       +           ue(4,4) )
c
               frct(4,i,3,k) = frct(4,i,3,k)
     $           - dsspm * ( - 4.0d+00 * ue(4,2)
     $                       + 6.0d+00 * ue(4,3)
     $                       - 4.0d+00 * ue(4,4)
     $                       +           ue(4,5) )
c
c
               frct(5,i,2,k) = frct(5,i,2,k)
     $           - dsspm * ( + 5.0d+00 * ue(5,2)
     $                       - 4.0d+00 * ue(5,3)
     $                       +           ue(5,4) )
c
               frct(5,i,3,k) = frct(5,i,3,k)
     $           - dsspm * ( - 4.0d+00 * ue(5,2)
     $                       + 6.0d+00 * ue(5,3)
     $                       - 4.0d+00 * ue(5,4)
     $                       +           ue(5,5) )
c
c
            do j = 4, ny - 3
c
c
                  frct(1,i,j,k) = frct(1,i,j,k)
     $              - dsspm * (            ue(1,j-2)
     $                        - 4.0d+00 * ue(1,j-1)
     $                        + 6.0d+00 * ue(1,j)
     $                        - 4.0d+00 * ue(1,j+1)
     $                        +           ue(1,j+2) )
c
c
                  frct(2,i,j,k) = frct(2,i,j,k)
     $              - dsspm * (            ue(2,j-2)
     $                        - 4.0d+00 * ue(2,j-1)
     $                        + 6.0d+00 * ue(2,j)
     $                        - 4.0d+00 * ue(2,j+1)
     $                        +           ue(2,j+2) )
c
c
                  frct(3,i,j,k) = frct(3,i,j,k)
     $              - dsspm * (            ue(3,j-2)
     $                        - 4.0d+00 * ue(3,j-1)
     $                        + 6.0d+00 * ue(3,j)
     $                        - 4.0d+00 * ue(3,j+1)
     $                        +           ue(3,j+2) )
c
c
                  frct(4,i,j,k) = frct(4,i,j,k)
     $              - dsspm * (            ue(4,j-2)
     $                        - 4.0d+00 * ue(4,j-1)
     $                        + 6.0d+00 * ue(4,j)
     $                        - 4.0d+00 * ue(4,j+1)
     $                        +           ue(4,j+2) )
c
c
                  frct(5,i,j,k) = frct(5,i,j,k)
     $              - dsspm * (            ue(5,j-2)
     $                        - 4.0d+00 * ue(5,j-1)
     $                        + 6.0d+00 * ue(5,j)
     $                        - 4.0d+00 * ue(5,j+1)
     $                        +           ue(5,j+2) )
c
c
            end do
c
c
               frct(1,i,ny-2,k) = frct(1,i,ny-2,k)
     $           - dsspm * (             ue(1,ny-4)
     $                       - 4.0d+00 * ue(1,ny-3)
     $                       + 6.0d+00 * ue(1,ny-2)
     $                       - 4.0d+00 * ue(1,ny-1)  )
c
               frct(1,i,ny-1,k) = frct(1,i,ny-1,k)
     $           - dsspm * (             ue(1,ny-3)
     $                       - 4.0d+00 * ue(1,ny-2)
     $                       + 5.0d+00 * ue(1,ny-1)  )
c
c
               frct(2,i,ny-2,k) = frct(2,i,ny-2,k)
     $           - dsspm * (             ue(2,ny-4)
     $                       - 4.0d+00 * ue(2,ny-3)
     $                       + 6.0d+00 * ue(2,ny-2)
     $                       - 4.0d+00 * ue(2,ny-1)  )
c
               frct(2,i,ny-1,k) = frct(2,i,ny-1,k)
     $           - dsspm * (             ue(2,ny-3)
     $                       - 4.0d+00 * ue(2,ny-2)
     $                       + 5.0d+00 * ue(2,ny-1)  )
c
c
               frct(3,i,ny-2,k) = frct(3,i,ny-2,k)
     $           - dsspm * (             ue(3,ny-4)
     $                       - 4.0d+00 * ue(3,ny-3)
     $                       + 6.0d+00 * ue(3,ny-2)
     $                       - 4.0d+00 * ue(3,ny-1)  )
c
               frct(3,i,ny-1,k) = frct(3,i,ny-1,k)
     $           - dsspm * (             ue(3,ny-3)
     $                       - 4.0d+00 * ue(3,ny-2)
     $                       + 5.0d+00 * ue(3,ny-1)  )
c
c
               frct(4,i,ny-2,k) = frct(4,i,ny-2,k)
     $           - dsspm * (             ue(4,ny-4)
     $                       - 4.0d+00 * ue(4,ny-3)
     $                       + 6.0d+00 * ue(4,ny-2)
     $                       - 4.0d+00 * ue(4,ny-1)  )
c
               frct(4,i,ny-1,k) = frct(4,i,ny-1,k)
     $           - dsspm * (             ue(4,ny-3)
     $                       - 4.0d+00 * ue(4,ny-2)
     $                       + 5.0d+00 * ue(4,ny-1)  )
c
c
               frct(5,i,ny-2,k) = frct(5,i,ny-2,k)
     $           - dsspm * (             ue(5,ny-4)
     $                       - 4.0d+00 * ue(5,ny-3)
     $                       + 6.0d+00 * ue(5,ny-2)
     $                       - 4.0d+00 * ue(5,ny-1)  )
c
               frct(5,i,ny-1,k) = frct(5,i,ny-1,k)
     $           - dsspm * (             ue(5,ny-3)
     $                       - 4.0d+00 * ue(5,ny-2)
     $                       + 5.0d+00 * ue(5,ny-1)  )
c
c
         end do
c
      end do
c
c***zeta-direction flux differences
c
CHPF$ INDEPENDENT
      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
c
                  frct(1,i,j,k) =  frct(1,i,j,k)
     $                       - tz2 * ( flux(1,k+1) - flux(1,k-1) )
                  frct(2,i,j,k) =  frct(2,i,j,k)
     $                       - tz2 * ( flux(2,k+1) - flux(2,k-1) )
                  frct(3,i,j,k) =  frct(3,i,j,k)
     $                       - tz2 * ( flux(3,k+1) - flux(3,k-1) )
                  frct(4,i,j,k) =  frct(4,i,j,k)
     $                       - tz2 * ( flux(4,k+1) - flux(4,k-1) )
                  frct(5,i,j,k) =  frct(5,i,j,k)
     $                       - tz2 * ( flux(5,k+1) - flux(5,k-1) )
c
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
c
               frct(1,i,j,2) = frct(1,i,j,2)
     $           - dsspm * ( + 5.0d+00 * ue(1,2)
     $                       - 4.0d+00 * ue(1,3)
     $                       +           ue(1,4) )
c
               frct(1,i,j,3) = frct(1,i,j,3)
     $           - dsspm * (- 4.0d+00 * ue(1,2)
     $                      + 6.0d+00 * ue(1,3)
     $                      - 4.0d+00 * ue(1,4)
     $                      +           ue(1,5) )
c
c
               frct(2,i,j,2) = frct(2,i,j,2)
     $           - dsspm * ( + 5.0d+00 * ue(2,2)
     $                       - 4.0d+00 * ue(2,3)
     $                       +           ue(2,4) )
c
               frct(2,i,j,3) = frct(2,i,j,3)
     $           - dsspm * (- 4.0d+00 * ue(2,2)
     $                      + 6.0d+00 * ue(2,3)
     $                      - 4.0d+00 * ue(2,4)
     $                      +           ue(2,5) )
c
c
               frct(3,i,j,2) = frct(3,i,j,2)
     $           - dsspm * ( + 5.0d+00 * ue(3,2)
     $                       - 4.0d+00 * ue(3,3)
     $                       +           ue(3,4) )
c
               frct(3,i,j,3) = frct(3,i,j,3)
     $           - dsspm * (- 4.0d+00 * ue(3,2)
     $                      + 6.0d+00 * ue(3,3)
     $                      - 4.0d+00 * ue(3,4)
     $                      +           ue(3,5) )
c
c
               frct(4,i,j,2) = frct(4,i,j,2)
     $           - dsspm * ( + 5.0d+00 * ue(4,2)
     $                       - 4.0d+00 * ue(4,3)
     $                       +           ue(4,4) )
c
               frct(4,i,j,3) = frct(4,i,j,3)
     $           - dsspm * (- 4.0d+00 * ue(4,2)
     $                      + 6.0d+00 * ue(4,3)
     $                      - 4.0d+00 * ue(4,4)
     $                      +           ue(4,5) )
c
c
               frct(5,i,j,2) = frct(5,i,j,2)
     $           - dsspm * ( + 5.0d+00 * ue(5,2)
     $                       - 4.0d+00 * ue(5,3)
     $                       +           ue(5,4) )
c
               frct(5,i,j,3) = frct(5,i,j,3)
     $           - dsspm * (- 4.0d+00 * ue(5,2)
     $                      + 6.0d+00 * ue(5,3)
     $                      - 4.0d+00 * ue(5,4)
     $                      +           ue(5,5) )
c
c
            do k = 4, nz - 3
c
c
                  frct(1,i,j,k) = frct(1,i,j,k)
     $              - dsspm * (           ue(1,k-2)
     $                        - 4.0d+00 * ue(1,k-1)
     $                        + 6.0d+00 * ue(1,k)
     $                        - 4.0d+00 * ue(1,k+1)
     $                        +           ue(1,k+2) )
c
c
                  frct(2,i,j,k) = frct(2,i,j,k)
     $              - dsspm * (           ue(2,k-2)
     $                        - 4.0d+00 * ue(2,k-1)
     $                        + 6.0d+00 * ue(2,k)
     $                        - 4.0d+00 * ue(2,k+1)
     $                        +           ue(2,k+2) )
c
c
                  frct(3,i,j,k) = frct(3,i,j,k)
     $              - dsspm * (           ue(3,k-2)
     $                        - 4.0d+00 * ue(3,k-1)
     $                        + 6.0d+00 * ue(3,k)
     $                        - 4.0d+00 * ue(3,k+1)
     $                        +           ue(3,k+2) )
c
c
                  frct(4,i,j,k) = frct(4,i,j,k)
     $              - dsspm * (           ue(4,k-2)
     $                        - 4.0d+00 * ue(4,k-1)
     $                        + 6.0d+00 * ue(4,k)
     $                        - 4.0d+00 * ue(4,k+1)
     $                        +           ue(4,k+2) )
c
c
                  frct(5,i,j,k) = frct(5,i,j,k)
     $              - dsspm * (           ue(5,k-2)
     $                        - 4.0d+00 * ue(5,k-1)
     $                        + 6.0d+00 * ue(5,k)
     $                        - 4.0d+00 * ue(5,k+1)
     $                        +           ue(5,k+2) )
c
c
            end do
c
c
               frct(1,i,j,nz-2) = frct(1,i,j,nz-2)
     $           - dsspm * (            ue(1,nz-4)
     $                      - 4.0d+00 * ue(1,nz-3)
     $                      + 6.0d+00 * ue(1,nz-2)
     $                      - 4.0d+00 * ue(1,nz-1)  )
c
               frct(1,i,j,nz-1) = frct(1,i,j,nz-1)
     $           - dsspm * (             ue(1,nz-3)
     $                       - 4.0d+00 * ue(1,nz-2)
     $                       + 5.0d+00 * ue(1,nz-1)  )
c
c
               frct(2,i,j,nz-2) = frct(2,i,j,nz-2)
     $           - dsspm * (            ue(2,nz-4)
     $                      - 4.0d+00 * ue(2,nz-3)
     $                      + 6.0d+00 * ue(2,nz-2)
     $                      - 4.0d+00 * ue(2,nz-1)  )
c
               frct(2,i,j,nz-1) = frct(2,i,j,nz-1)
     $           - dsspm * (             ue(2,nz-3)
     $                       - 4.0d+00 * ue(2,nz-2)
     $                       + 5.0d+00 * ue(2,nz-1)  )
c
c
               frct(3,i,j,nz-2) = frct(3,i,j,nz-2)
     $           - dsspm * (            ue(3,nz-4)
     $                      - 4.0d+00 * ue(3,nz-3)
     $                      + 6.0d+00 * ue(3,nz-2)
     $                      - 4.0d+00 * ue(3,nz-1)  )
c
               frct(3,i,j,nz-1) = frct(3,i,j,nz-1)
     $           - dsspm * (             ue(3,nz-3)
     $                       - 4.0d+00 * ue(3,nz-2)
     $                       + 5.0d+00 * ue(3,nz-1)  )
c
c
               frct(4,i,j,nz-2) = frct(4,i,j,nz-2)
     $           - dsspm * (            ue(4,nz-4)
     $                      - 4.0d+00 * ue(4,nz-3)
     $                      + 6.0d+00 * ue(4,nz-2)
     $                      - 4.0d+00 * ue(4,nz-1)  )
c
               frct(4,i,j,nz-1) = frct(4,i,j,nz-1)
     $           - dsspm * (             ue(4,nz-3)
     $                       - 4.0d+00 * ue(4,nz-2)
     $                       + 5.0d+00 * ue(4,nz-1)  )
c
c
               frct(5,i,j,nz-2) = frct(5,i,j,nz-2)
     $           - dsspm * (            ue(5,nz-4)
     $                      - 4.0d+00 * ue(5,nz-3)
     $                      + 6.0d+00 * ue(5,nz-2)
     $                      - 4.0d+00 * ue(5,nz-1)  )
c
               frct(5,i,j,nz-1) = frct(5,i,j,nz-1)
     $           - dsspm * (             ue(5,nz-3)
     $                       - 4.0d+00 * ue(5,nz-2)
     $                       + 5.0d+00 * ue(5,nz-1)  )
c
c
         end do
c
      end do
      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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





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
c
                     tmp1 = abs ( u000ijk(m) - u(1,i,j,k) )
c
                     if ( tmp1 .gt. errmax(m) ) then
c
                        errmax(m) = tmp1
                        imax(m) = i
                        jmax(m) = j
                        kmax(m) = k
c
                     end if
c
c
                     tmp2 = abs ( u000ijk(m) - u(2,i,j,k) )
c
                     if ( tmp2 .gt. errmax(m) ) then
c
                        errmax(m) = tmp2
                        imax(m) = i
                        jmax(m) = j
                        kmax(m) = k
c
                     end if
c
c
                     tmp3 = abs ( u000ijk(m) - u(3,i,j,k) )
c
                     if ( tmp3 .gt. errmax(m) ) then
c
                        errmax(m) = tmp3
                        imax(m) = i
                        jmax(m) = j
                        kmax(m) = k
c
                     end if
c
c
                     tmp4 = abs ( u000ijk(m) - u(4,i,j,k) )
c
                     if ( tmp4 .gt. errmax(m) ) then
c
                        errmax(m) = tmp4
                        imax(m) = i
                        jmax(m) = j
                        kmax(m) = k
c
                     end if
c
c
                     tmp5 = abs ( u000ijk(m) - u(5,i,j,k) )
c
                     if ( tmp5 .gt. errmax(m) ) then
c
                        errmax(m) = tmp5
                        imax(m) = i
                        jmax(m) = j
                        kmax(m) = k
c
                     end if
c
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
c
                     tmp1 = ( u000ijk(1) - u(1,i,j,k) )
                     tmp2 = ( u000ijk(2) - u(2,i,j,k) )
                     tmp3 = ( u000ijk(3) - u(3,i,j,k) )
                     tmp4 = ( u000ijk(4) - u(4,i,j,k) )
                     tmp5 = ( u000ijk(5) - u(5,i,j,k) )
c
                     errnm(1) = errnm(1) + tmp1 ** 2
                     errnm(2) = errnm(2) + tmp2 ** 2
                     errnm(3) = errnm(3) + tmp3 ** 2
                     errnm(4) = errnm(4) + tmp4 ** 2
                     errnm(5) = errnm(5) + tmp5 ** 2
c
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
      subroutine etasweep
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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c
c
c
c
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
      dimension idmax(5), jdmax(5), kdmax(5),
     $          imax(5), jmax(5), kmax(5),
     $          delunm(5)
c
      dimension flux1(isiz1)
      dimension flux2(isiz1)
      dimension flux3(isiz1)
      dimension flux4(isiz1)
      dimension flux5(isiz1)
c
      dimension cv(isiz3), aa(isiz3), rhos(isiz3)
      dimension  rhon(isiz1)
      dimension  rhoq(isiz2)
c
      parameter ( one = 1.0d+00 )
c
c
      bt = sqrt ( 0.50d+00 )
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c
c***perform the eta-direction sweep
c
c        call jacy ( 3 ) 
c
c
CHPF$ INDEPENDENT
      do k = 2, nz - 1
c
         do i = 2, nx - 1
c
            do j = 1, ny
c
               ru1 = 1.0d+00 / u(1,i,j,k)
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(j) = vv
               aa(j) = sqrt ( c1 * c2 * ( ru1 * u(5,i,j,k) - q ) )
c
               rhoq(j) = max ( dy1,
     $                         dy2 + c34 * ru1,
     $                         dy3 + r43 * c34 * ru1,
     $                         dy4 + c34 * ru1,
     $                         dy5 + c1345 * ru1 )
c
            end do
c
            a(i,1,k) = 0.0d+00
            b(i,1,k) = 0.0d+00
            c(i,1,k) = 1.0d+00
            d(i,1,k) = 0.0d+00
            e(i,1,k) = 0.0d+00
c
            do j = 2, ny - 1
c
               a(i,j,k) = 0.0d+00
               b(i,j,k) = - dt * ty2 *  cv(j-1)
     $                    - dt * rhoq(j-1) * ty1 
               c(i,j,k) =   1.0d+00
     $                    + dt * rhoq(j) * ty1 * 2.0d+00
               d(i,j,k) =   dt * ty2 *  cv(j+1)
     $                    - dt * rhoq(j+1) * ty1
               e(i,j,k) = 0.0d+00
c
            end do
c
            a(i,ny,k) = 0.0d+00
            b(i,ny,k) = 0.0d+00
            c(i,ny,k) = 1.0d+00
            d(i,ny,k) = 0.0d+00
            e(i,ny,k) = 0.0d+00
c
c***fourth order dissipation 
c
            c(i,2,k) = c(i,2,k) + dt * dssp * ( + 5.0d+00 )
            d(i,2,k) = d(i,2,k) + dt * dssp * ( - 4.0d+00 )
            e(i,2,k) = e(i,2,k) + dt * dssp * ( + 1.0d+00 )
c
            b(i,3,k) = b(i,3,k) + dt * dssp * ( - 4.0d+00 )
            c(i,3,k) = c(i,3,k) + dt * dssp * ( + 6.0d+00 )
            d(i,3,k) = d(i,3,k) + dt * dssp * ( - 4.0d+00 )
            e(i,3,k) = e(i,3,k) + dt * dssp * ( + 1.0d+00 )
c
            do j = 4, ny - 3
c
               a(i,j,k) = a(i,j,k) + dt * dssp * ( + 1.0d+00 )
               b(i,j,k) = b(i,j,k) + dt * dssp * ( - 4.0d+00 )
               c(i,j,k) = c(i,j,k) + dt * dssp * ( + 6.0d+00 )
               d(i,j,k) = d(i,j,k) + dt * dssp * ( - 4.0d+00 )
               e(i,j,k) = e(i,j,k) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(i,ny-2,k) = a(i,ny-2,k) + dt * dssp * ( + 1.0d+00 )
            b(i,ny-2,k) = b(i,ny-2,k) + dt * dssp * ( - 4.0d+00 )
            c(i,ny-2,k) = c(i,ny-2,k) + dt * dssp * ( + 6.0d+00 )
            d(i,ny-2,k) = d(i,ny-2,k) + dt * dssp * ( - 4.0d+00 )
c
            a(i,ny-1,k) = a(i,ny-1,k) + dt * dssp * ( + 1.0d+00 )
            b(i,ny-1,k) = b(i,ny-1,k) + dt * dssp * ( - 4.0d+00 )
            c(i,ny-1,k) = c(i,ny-1,k) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentay3 ( isiz1, isiz2, isiz3,
c    $                   nx, ny, nz,
c    $                   a, b, c, d, e,
c    $                   rsd )
c
c***forward elimination
c
            do j = 2, ny-1
c
               tmp1 = b(i,j,k) / c(i,j-1,k)
               tmp2 = a(i,j+1,k) / c(i,j-1,k)
c     
               c(i,j,k) = c(i,j,k) - tmp1 * d(i,j-1,k)
               d(i,j,k) = d(i,j,k) - tmp1 * e(i,j-1,k) 
               rsd(1,i,j,k) = rsd(1,i,j,k) - tmp1 * rsd(1,i,j-1,k)
               rsd(2,i,j,k) = rsd(2,i,j,k) - tmp1 * rsd(2,i,j-1,k)
               rsd(3,i,j,k) = rsd(3,i,j,k) - tmp1 * rsd(3,i,j-1,k)
c
               b(i,j+1,k) = b(i,j+1,k) - tmp2 * d(i,j-1,k)
               c(i,j+1,k) = c(i,j+1,k) - tmp2 * e(i,j-1,k)
               rsd(1,i,j+1,k) = rsd(1,i,j+1,k) - tmp2 * rsd(1,i,j-1,k)
               rsd(2,i,j+1,k) = rsd(2,i,j+1,k) - tmp2 * rsd(2,i,j-1,k)
               rsd(3,i,j+1,k) = rsd(3,i,j+1,k) - tmp2 * rsd(3,i,j-1,k)
c     
            end do
c     
            tmp1 = b(i,ny,k) / c(i,ny-1,k)
c
            c(i,ny,k) = c(i,ny,k) - tmp1 * d(i,ny-1,k)
            d(i,ny,k) = d(i,ny,k) - tmp1 * e(i,ny-1,k)
            rsd(1,i,ny,k) = rsd(1,i,ny,k) - tmp1 * rsd(1,i,ny-1,k)
            rsd(2,i,ny,k) = rsd(2,i,ny,k) - tmp1 * rsd(2,i,ny-1,k)
            rsd(3,i,ny,k) = rsd(3,i,ny,k) - tmp1 * rsd(3,i,ny-1,k)
c
c***back-substitution phase
c
            rsd(1,i,ny,k) = rsd(1,i,ny,k) / c(i,ny,k)
      rsd(1,i,ny-1,k) = ( rsd(1,i,ny-1,k) - d(i,ny-1,k)*rsd(1,i,ny,k) )
     $           / c(i,ny-1,k)
c
            rsd(2,i,ny,k) = rsd(2,i,ny,k) / c(i,ny,k)
      rsd(2,i,ny-1,k) = ( rsd(2,i,ny-1,k) - d(i,ny-1,k)*rsd(2,i,ny,k) )
     $           / c(i,ny-1,k)
c
            rsd(3,i,ny,k) = rsd(3,i,ny,k) / c(i,ny,k)
       rsd(3,i,ny-1,k) = ( rsd(3,i,ny-1,k) - d(i,ny-1,k)*rsd(3,i,ny,k) )
     $           / c(i,ny-1,k)
c
            do j = ny-2, 1, -1
c
               rsd(1,i,j,k) = ( rsd(1,i,j,k) - d(i,j,k)*rsd(1,i,j+1,k)
     $              - e(i,j,k)*rsd(1,i,j+2,k) ) / c(i,j,k)
c
               rsd(2,i,j,k) = ( rsd(2,i,j,k) - d(i,j,k)*rsd(2,i,j+1,k)
     $              - e(i,j,k)*rsd(2,i,j+2,k) ) / c(i,j,k)
c
               rsd(3,i,j,k) = ( rsd(3,i,j,k) - d(i,j,k)*rsd(3,i,j+1,k)
     $              - e(i,j,k)*rsd(3,i,j+2,k) ) / c(i,j,k)
c     
            end do
c
c
c        call jacy ( 4 ) 
c
            do j = 1, ny
c
               ru1 = 1.0d+00 / u(1,i,j,k)
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(j) = vv
               aa(j) = sqrt ( c1 * c2 * ( ru1 * u(5,i,j,k) - q ) )
c
               rhoq(j) = max ( dy1,
     $                         dy2 + c34 * ru1,
     $                         dy3 + r43 * c34 * ru1,
     $                         dy4 + c34 * ru1,
     $                         dy5 + c1345 * ru1 )
c
            end do
c
            a(i,1,k) = 0.0d+00
            b(i,1,k) = 0.0d+00
            c(i,1,k) = 1.0d+00
            d(i,1,k) = 0.0d+00
            e(i,1,k) = 0.0d+00
c
            do j = 2, ny - 1
c
               a(i,j,k) = 0.0d+00
               b(i,j,k) = - dt * ty2 * ( cv(j-1) + aa(j-1) )
     $                    - dt * rhoq(j-1) * ty1 
               c(i,j,k) =   1.0d+00
     $                    + dt * rhoq(j) * ty1 * 2.0d+00
               d(i,j,k) =   dt * ty2 * ( cv(j+1) + aa(j+1) )
     $                    - dt * rhoq(j+1) * ty1
               e(i,j,k) = 0.0d+00
c
            end do
c
            a(i,ny,k) = 0.0d+00
            b(i,ny,k) = 0.0d+00
            c(i,ny,k) = 1.0d+00
            d(i,ny,k) = 0.0d+00
            e(i,ny,k) = 0.0d+00
c
c***fourth order dissipation 
c
            c(i,2,k) = c(i,2,k) + dt * dssp * ( + 5.0d+00 )
            d(i,2,k) = d(i,2,k) + dt * dssp * ( - 4.0d+00 )
            e(i,2,k) = e(i,2,k) + dt * dssp * ( + 1.0d+00 )
c
            b(i,3,k) = b(i,3,k) + dt * dssp * ( - 4.0d+00 )
            c(i,3,k) = c(i,3,k) + dt * dssp * ( + 6.0d+00 )
            d(i,3,k) = d(i,3,k) + dt * dssp * ( - 4.0d+00 )
            e(i,3,k) = e(i,3,k) + dt * dssp * ( + 1.0d+00 )
c
            do j = 4, ny - 3
c
               a(i,j,k) = a(i,j,k) + dt * dssp * ( + 1.0d+00 )
               b(i,j,k) = b(i,j,k) + dt * dssp * ( - 4.0d+00 )
               c(i,j,k) = c(i,j,k) + dt * dssp * ( + 6.0d+00 )
               d(i,j,k) = d(i,j,k) + dt * dssp * ( - 4.0d+00 )
               e(i,j,k) = e(i,j,k) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(i,ny-2,k) = a(i,ny-2,k) + dt * dssp * ( + 1.0d+00 )
            b(i,ny-2,k) = b(i,ny-2,k) + dt * dssp * ( - 4.0d+00 )
            c(i,ny-2,k) = c(i,ny-2,k) + dt * dssp * ( + 6.0d+00 )
            d(i,ny-2,k) = d(i,ny-2,k) + dt * dssp * ( - 4.0d+00 )
c
            a(i,ny-1,k) = a(i,ny-1,k) + dt * dssp * ( + 1.0d+00 )
            b(i,ny-1,k) = b(i,ny-1,k) + dt * dssp * ( - 4.0d+00 )
            c(i,ny-1,k) = c(i,ny-1,k) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentay ( isiz1, isiz2, isiz3,
c    $                  4, nx, ny, nz,
c    $                  a, b, c, d, e,
c    $                  rsd )
c
c***forward elimination
c
            do j = 2, ny - 1
c
               tmp1 = b(i,j,k) / c(i,j-1,k)
               tmp2 = a(i,j+1,k) / c(i,j-1,k)
c     
               c(i,j,k) = c(i,j,k) - tmp1 * d(i,j-1,k)
               d(i,j,k) = d(i,j,k) - tmp1 * e(i,j-1,k) 
               rsd(4,i,j,k) = rsd(4,i,j,k) - tmp1 * rsd(4,i,j-1,k)
c
               b(i,j+1,k) = b(i,j+1,k) - tmp2 * d(i,j-1,k)
               c(i,j+1,k) = c(i,j+1,k) - tmp2 * e(i,j-1,k)
               rsd(4,i,j+1,k) = rsd(4,i,j+1,k) - tmp2 * rsd(4,i,j-1,k)
c     
            end do
c     
            tmp1 = b(i,ny,k) / c(i,ny-1,k)
c
            c(i,ny,k) = c(i,ny,k) - tmp1 * d(i,ny-1,k)
            d(i,ny,k) = d(i,ny,k) - tmp1 * e(i,ny-1,k)
            rsd(4,i,ny,k) = rsd(4,i,ny,k) - tmp1 * rsd(4,i,ny-1,k)
c
c***back-substitution phase
c
            rsd(4,i,ny,k) = rsd(4,i,ny,k) / c(i,ny,k)
      rsd(4,i,ny-1,k) = ( rsd(4,i,ny-1,k) - d(i,ny-1,k)*rsd(4,i,ny,k) )
     $           / c(i,ny-1,k)
c
            do j = ny - 2, 1, -1
c
               rsd(4,i,j,k) = ( rsd(4,i,j,k) - d(i,j,k)*rsd(4,i,j+1,k)
     $              - e(i,j,k)*rsd(4,i,j+2,k) ) / c(i,j,k)
c     
            end do
c
c
c        call jacy ( 5 ) 
c
            do j = 1, ny
c
               ru1 = 1.0d+00 / u(1,i,j,k)
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(j) = vv
               aa(j) = sqrt ( c1 * c2 * ( ru1 * u(5,i,j,k) - q ) )
c
               rhoq(j) = max ( dy1,
     $                         dy2 + c34 * ru1,
     $                         dy3 + r43 * c34 * ru1,
     $                         dy4 + c34 * ru1,
     $                         dy5 + c1345 * ru1 )
c
            end do
c
            a(i,1,k) = 0.0d+00
            b(i,1,k) = 0.0d+00
            c(i,1,k) = 1.0d+00
            d(i,1,k) = 0.0d+00
            e(i,1,k) = 0.0d+00
c
            do j = 2, ny - 1
c
               a(i,j,k) = 0.0d+00
               b(i,j,k) = - dt * ty2 * ( cv(j-1) -  aa(j-1) )
     $                    - dt * rhoq(j-1) * ty1 
               c(i,j,k) =   1.0d+00
     $                    + dt * rhoq(j) * ty1 * 2.0d+00
               d(i,j,k) =   dt * ty2 * ( cv(j+1) - aa(j+1) )
     $                    - dt * rhoq(j+1) * ty1
               e(i,j,k) = 0.0d+00
c
            end do
c
            a(i,ny,k) = 0.0d+00
            b(i,ny,k) = 0.0d+00
            c(i,ny,k) = 1.0d+00
            d(i,ny,k) = 0.0d+00
            e(i,ny,k) = 0.0d+00
c
c***fourth order dissipation 
c
            c(i,2,k) = c(i,2,k) + dt * dssp * ( + 5.0d+00 )
            d(i,2,k) = d(i,2,k) + dt * dssp * ( - 4.0d+00 )
            e(i,2,k) = e(i,2,k) + dt * dssp * ( + 1.0d+00 )
c
            b(i,3,k) = b(i,3,k) + dt * dssp * ( - 4.0d+00 )
            c(i,3,k) = c(i,3,k) + dt * dssp * ( + 6.0d+00 )
            d(i,3,k) = d(i,3,k) + dt * dssp * ( - 4.0d+00 )
            e(i,3,k) = e(i,3,k) + dt * dssp * ( + 1.0d+00 )
c
            do j = 4, ny - 3
c
               a(i,j,k) = a(i,j,k) + dt * dssp * ( + 1.0d+00 )
               b(i,j,k) = b(i,j,k) + dt * dssp * ( - 4.0d+00 )
               c(i,j,k) = c(i,j,k) + dt * dssp * ( + 6.0d+00 )
               d(i,j,k) = d(i,j,k) + dt * dssp * ( - 4.0d+00 )
               e(i,j,k) = e(i,j,k) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(i,ny-2,k) = a(i,ny-2,k) + dt * dssp * ( + 1.0d+00 )
            b(i,ny-2,k) = b(i,ny-2,k) + dt * dssp * ( - 4.0d+00 )
            c(i,ny-2,k) = c(i,ny-2,k) + dt * dssp * ( + 6.0d+00 )
            d(i,ny-2,k) = d(i,ny-2,k) + dt * dssp * ( - 4.0d+00 )
c
            a(i,ny-1,k) = a(i,ny-1,k) + dt * dssp * ( + 1.0d+00 )
            b(i,ny-1,k) = b(i,ny-1,k) + dt * dssp * ( - 4.0d+00 )
            c(i,ny-1,k) = c(i,ny-1,k) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentay ( isiz1, isiz2, isiz3,
c    $                  5, nx, ny, nz,
c    $                  a, b, c, d, e,
c    $                  rsd )
c***forward elimination
c
            do j = 2, ny - 1
c
               tmp1 = b(i,j,k) / c(i,j-1,k)
               tmp2 = a(i,j+1,k) / c(i,j-1,k)
c     
               c(i,j,k) = c(i,j,k) - tmp1 * d(i,j-1,k)
               d(i,j,k) = d(i,j,k) - tmp1 * e(i,j-1,k) 
               rsd(5,i,j,k) = rsd(5,i,j,k) - tmp1 * rsd(5,i,j-1,k)
c
               b(i,j+1,k) = b(i,j+1,k) - tmp2 * d(i,j-1,k)
               c(i,j+1,k) = c(i,j+1,k) - tmp2 * e(i,j-1,k)
               rsd(5,i,j+1,k) = rsd(5,i,j+1,k) - tmp2 * rsd(5,i,j-1,k)
c     
            end do
c     
            tmp1 = b(i,ny,k) / c(i,ny-1,k)
c
            c(i,ny,k) = c(i,ny,k) - tmp1 * d(i,ny-1,k)
            d(i,ny,k) = d(i,ny,k) - tmp1 * e(i,ny-1,k)
            rsd(5,i,ny,k) = rsd(5,i,ny,k) - tmp1 * rsd(5,i,ny-1,k)
c
c***back-substitution phase
c
            rsd(5,i,ny,k) = rsd(5,i,ny,k) / c(i,ny,k)
      rsd(5,i,ny-1,k) = ( rsd(5,i,ny-1,k) - d(i,ny-1,k)*rsd(5,i,ny,k) )
     $           / c(i,ny-1,k)
c
            do j = ny - 2, 1, -1
c
               rsd(5,i,j,k) = ( rsd(5,i,j,k) - d(i,j,k)*rsd(5,i,j+1,k)
     $              - e(i,j,k)*rsd(5,i,j+2,k) ) / c(i,j,k)
c     
            end do
c
         end do
		 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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





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 pintgr
c
c***compute the surface integral
c
c Author: Sisira Weeratunga
c         NASA Ames Research Center
c         (10/25/90)
c
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





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
	  subroutine rhs
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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
      common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3,
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c
c
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c
 
 
 
 
 
c
      dimension idmax(5), jdmax(5), kdmax(5),
     $          imax(5), jmax(5), kmax(5),
     $          delunm(5)
c
c
      real * 8 njac
c
      dimension flux1jm1(isiz1),flux2jm1(isiz1),flux3jm1(isiz1),
     $          flux4jm1(isiz1),flux5jm1(isiz1)
      dimension flux1j(isiz1),flux2j(isiz1),flux3j(isiz1),
     $          flux4j(isiz1),flux5j(isiz1)
      dimension flux1jp1(isiz1),flux2jp1(isiz1),flux3jp1(isiz1),
     $          flux4jp1(isiz1),flux5jp1(isiz1)
      dimension flx2j(isiz1),flx3j(isiz1),
     $          flx4j(isiz1),flx5j(isiz1)
      dimension flx2jp1(isiz1),flx3jp1(isiz1),
     $          flx4jp1(isiz1),flx5jp1(isiz1)
      dimension u21jm1(isiz1),u31jm1(isiz1),u41jm1(isiz1),
     $          u51jm1(isiz1)
      dimension u21j(isiz1),u31j(isiz1),u41j(isiz1),
     $          u51j(isiz1)
      dimension u21jp1(isiz1),u31jp1(isiz1),u41jp1(isiz1),
     $          u51jp1(isiz1)
      dimension rsd1ijk(isiz1),rsd2ijk(isiz1),rsd3ijk(isiz1),
     $          rsd4ijk(isiz1),rsd5ijk(isiz1)
      dimension tm(5,5),tv(5)
      parameter ( one = 1.0d+00 )
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c
c***compute the steady-state residuals
c     
c     call rhs
c
c
      do k = 1, nz
         do j = 1, ny
               do m = 1, 5
                  rsd(m,1,j,k) = - frct(m,1,j,k)
                  rsd(m,nx,j,k) = - frct(m,nx,j,k)
               end do
         end do
      end do
c
 
c***xi-direction flux differences
 
c
	   do k = 2, nz - 1
c
c
            do i = 2, nx - 1
c
               flux1jm1(i) = u(3,i,1,k)
c
               u31 = u(3,i,1,k) / u(1,i,1,k)
c
               q = 0.50d+00 * (  u(2,i,1,k) * u(2,i,1,k)
     $                         + u(3,i,1,k) * u(3,i,1,k)
     $                         + u(4,i,1,k) * u(4,i,1,k) )
     $                      / u(1,i,1,k)
c
               flux2jm1(i) = u(2,i,1,k) * u31
c
               flux3jm1(i) = u(3,i,1,k) * u31 + c2 * ( u(5,i,1,k) - q )
c
               flux4jm1(i) = u(4,i,1,k) * u31
c
               flux5jm1(i) = ( c1 * u(5,i,1,k) - c2 * q ) * u31
c 
c
               flux1j(i) = u(3,i,2,k)
c
               u31 = u(3,i,2,k) / u(1,i,2,k)
c
               q = 0.50d+00 * (  u(2,i,2,k) * u(2,i,2,k)
     $                         + u(3,i,2,k) * u(3,i,2,k)
     $                         + u(4,i,2,k) * u(4,i,2,k) )
     $                      / u(1,i,2,k)
c
               flux2j(i) = u(2,i,2,k) * u31 
c
               flux3j(i) = u(3,i,2,k) * u31+ c2 * ( u(5,i,2,k) - q )
c
               flux4j(i) = u(4,i,2,k) * u31
c
               flux5j(i) = ( c1 * u(5,i,2,k) - c2 * q ) * u31
c 
               tmp = 1.0d+00 / u(1,i,1,k)
c
               u21jm1(i) = tmp * u(2,i,1,k)
               u31jm1(i) = tmp * u(3,i,1,k)
               u41jm1(i) = tmp * u(4,i,1,k)
               u51jm1(i) = tmp * u(5,i,1,k)
c
               tmp = 1.0d+00 / u(1,i,2,k)
c
               u21j(i) = tmp * u(2,i,2,k)
               u31j(i) = tmp * u(3,i,2,k)
               u41j(i) = tmp * u(4,i,2,k)
               u51j(i) = tmp * u(5,i,2,k)
c
      flx3j(i) = (4.0d+00/3.0d+00) * ty3 * ( u31j(i) - u31jm1(i) )
               flx2j(i) = ty3 * ( u21j(i) - u21jm1(i) )
               flx4j(i) = ty3 * ( u41j(i) - u41jm1(i) )
               flx5j(i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $  * ty3 * ( ( u21j(i)  **2 + u31j(i)  **2 + u41j(i)  **2 )
     $   - ( u21jm1(i)**2 + u31jm1(i)**2 + u41jm1(i)**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * ty3 * ( u31j(i)**2 - u31jm1(i)**2 )
     $              + c1 * c5 * ty3 * ( u51j(i) - u51jm1(i) )
c
         end do
         do j = 2, ny - 1
c
               flux1im1 = u(2,1,j,k)
c
               u21 = u(2,1,j,k) / u(1,1,j,k)
c
               q = 0.50d+00 * (  u(2,1,j,k) * u(2,1,j,k)
     $                         + u(3,1,j,k) * u(3,1,j,k)
     $                         + u(4,1,j,k) * u(4,1,j,k) )
     $                      / u(1,1,j,k)
c
               flux2im1 = u(2,1,j,k) * u21 + c2 * ( u(5,1,j,k) - q )
c
               flux3im1 = u(3,1,j,k) * u21
c
               flux4im1 = u(4,1,j,k) * u21
c
               flux5im1 = ( c1 * u(5,1,j,k) - c2 * q ) * u21
c 
c
               flux1i = u(2,2,j,k)
c
               u21 = u(2,2,j,k) / u(1,2,j,k)
c
               q = 0.50d+00 * (  u(2,2,j,k) * u(2,2,j,k)
     $                         + u(3,2,j,k) * u(3,2,j,k)
     $                         + u(4,2,j,k) * u(4,2,j,k) )
     $                      / u(1,2,j,k)
c
               flux2i = u(2,2,j,k) * u21 + c2 * ( u(5,2,j,k) - q )
c
               flux3i = u(3,2,j,k) * u21
c
               flux4i = u(4,2,j,k) * u21
c
               flux5i = ( c1 * u(5,2,j,k) - c2 * q ) * u21
c 
               tmp = 1.0d+00 / u(1,1,j,k)
c
               u21im1 = tmp * u(2,1,j,k)
               u31im1 = tmp * u(3,1,j,k)
               u41im1 = tmp * u(4,1,j,k)
               u51im1 = tmp * u(5,1,j,k)
c
               tmp = 1.0d+00 / u(1,2,j,k)
c
               u21i = tmp * u(2,2,j,k)
               u31i = tmp * u(3,2,j,k)
               u41i = tmp * u(4,2,j,k)
               u51i = tmp * u(5,2,j,k)
c
               flx2i = (4.0d+00/3.0d+00) * tx3 * ( u21i - u21im1 )
               flx3i = tx3 * ( u31i - u31im1 )
               flx4i = tx3 * ( u41i - u41im1 )
               flx5i = 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
            do i = 2, nx-1
c
               flux1ip1 = u(2,i+1,j,k)
c
               u21 = u(2,i+1,j,k) / u(1,i+1,j,k)
c
               q = 0.50d+00 * (  u(2,i+1,j,k) * u(2,i+1,j,k)
     $                         + u(3,i+1,j,k) * u(3,i+1,j,k)
     $                         + u(4,i+1,j,k) * u(4,i+1,j,k) )
     $                      / u(1,i+1,j,k)
c
               flux2ip1 = u(2,i+1,j,k) * u21 + c2 * ( u(5,i+1,j,k) - q )
c
               flux3ip1 = u(3,i+1,j,k) * u21
c
               flux4ip1 = u(4,i+1,j,k) * u21
c
               flux5ip1 = ( c1 * u(5,i+1,j,k) - c2 * q ) * u21
c
                  rsd1 =  -frct(1,i,j,k)
     $                       - tx2 * ( flux1ip1 - flux1im1 )
                  rsd2 =  -frct(2,i,j,k)
     $                       - tx2 * ( flux2ip1 - flux2im1 )
                  rsd3 =  -frct(3,i,j,k)
     $                       - tx2 * ( flux3ip1 - flux3im1 )
                  rsd4 =  -frct(4,i,j,k)
     $                       - tx2 * ( flux4ip1 - flux4im1 )
                  rsd5 =  -frct(5,i,j,k)
     $                       - tx2 * ( flux5ip1 - flux5im1 )
c
               tmp = 1.0d+00 / u(1,i+1,j,k)
c
               u21ip1 = tmp * u(2,i+1,j,k)
               u31ip1 = tmp * u(3,i+1,j,k)
               u41ip1 = tmp * u(4,i+1,j,k)
               u51ip1 = tmp * u(5,i+1,j,k)
c
               flx2ip1 = (4.0d+00/3.0d+00) * tx3 * ( u21ip1 - u21i )
               flx3ip1 = tx3 * ( u31ip1 - u31i )
               flx4ip1 = tx3 * ( u41ip1 - u41i )
               flx5ip1 = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $        * tx3 * ( ( u21ip1  **2 + u31ip1  **2 + u41ip1  **2 )
     $                      - ( u21i**2 + u31i**2 + u41i**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * tx3 * ( u21ip1**2 - u21i**2 )
     $              + c1 * c5 * tx3 * ( u51ip1 - u51i )
c
               rsd1ijk(i) = rsd1
     $              + dx1 * tx1 * (        u(1,i-1,j,k)
     $                             - 2.0d+00 * u(1,i,j,k)
     $                             +           u(1,i+1,j,k) )
c
               rsd2ijk(i) = rsd2
     $              + tx3 * c3 * c4 * ( flx2ip1 - flx2i )
     $              + dx2 * tx1 * (        u(2,i-1,j,k)
     $                             - 2.0d+00 * u(2,i,j,k)
     $                             +           u(2,i+1,j,k) )
c
               rsd3ijk(i) = rsd3
     $              + tx3 * c3 * c4 * ( flx3ip1 - flx3i )
     $              + dx3 * tx1 * (        u(3,i-1,j,k)
     $                             - 2.0d+00 * u(3,i,j,k)
     $                             +           u(3,i+1,j,k) )
c
               rsd4ijk(i) = rsd4
     $              + tx3 * c3 * c4 * ( flx4ip1 - flx4i )
     $              + dx4 * tx1 * (        u(4,i-1,j,k)
     $                             - 2.0d+00 * u(4,i,j,k)
     $                             +           u(4,i+1,j,k) )
c
               rsd5ijk(i) = rsd5
     $              + tx3 * c3 * c4 * ( flx5ip1 - flx5i )
     $              + dx5 * tx1 * (        u(5,i-1,j,k)
     $                             - 2.0d+00 * u(5,i,j,k)
     $                             +           u(5,i+1,j,k) )
c
               flux1im1 = flux1i
c
               flux2im1 = flux2i
c
               flux3im1 = flux3i
c
               flux4im1 = flux4i
c
               flux5im1 = flux5i
c
               flux1i = flux1ip1
c
               flux2i = flux2ip1
c
               flux3i = flux3ip1
c
               flux4i = flux4ip1
c
               flux5i = flux5ip1
c 
               u21im1 = u21i
               u31im1 = u31i
               u41im1 = u41i
               u51im1 = u51i
c
               u21i = u21ip1
               u31i = u31ip1
               u41i = u41ip1
               u51i = u51ip1
c
               flx2i = flx2ip1
               flx3i = flx3ip1
               flx4i = flx4ip1
               flx5i = flx5ip1
c
       end do
c
c***Fourth-order dissipation
c
          i = 2
 
          rsd1ijk(i) = rsd1ijk(i) - dssp * ( + 5.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i+1,j,k)
     $                      +           u(1,i+2,j,k) )
          rsd2ijk(i) = rsd2ijk(i) - dssp * ( + 5.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i+1,j,k)
     $                      +           u(2,i+2,j,k) )
          rsd3ijk(i) = rsd3ijk(i) - dssp * ( + 5.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i+1,j,k)
     $                      +           u(3,i+2,j,k) )
          rsd4ijk(i) = rsd4ijk(i) - dssp * ( + 5.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i+1,j,k)
     $                      +           u(4,i+2,j,k) )
          rsd5ijk(i) = rsd5ijk(i) - dssp * ( + 5.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i+1,j,k)
     $                      +           u(5,i+2,j,k) )
          i = 3
c
          rsd1ijk(i) = rsd1ijk(i) - dssp * ( - 4.0d+00 * u(1,i-1,j,k)
     $                      + 6.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i+1,j,k)
     $                      +           u(1,i+2,j,k) )
          rsd2ijk(i) = rsd2ijk(i) - dssp * ( - 4.0d+00 * u(2,i-1,j,k)
     $                      + 6.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i+1,j,k)
     $                      +           u(2,i+2,j,k) )
          rsd3ijk(i) = rsd3ijk(i) - dssp * ( - 4.0d+00 * u(3,i-1,j,k)
     $                      + 6.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i+1,j,k)
     $                      +           u(3,i+2,j,k) )
          rsd4ijk(i) = rsd4ijk(i) - dssp * ( - 4.0d+00 * u(4,i-1,j,k)
     $                      + 6.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i+1,j,k)
     $                      +           u(4,i+2,j,k) )
          rsd5ijk(i) = rsd5ijk(i) - dssp * ( - 4.0d+00 * u(5,i-1,j,k)
     $                      + 6.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i+1,j,k)
     $                      +           u(5,i+2,j,k) )
c
c
         i = nx-2
c
         rsd1ijk(i) = rsd1ijk(i) - dssp * (       u(1,i-2,j,k)
     $                      - 4.0d+00 * u(1,i-1,j,k)
     $                      + 6.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i+1,j,k)  )
         rsd2ijk(i) = rsd2ijk(i) - dssp * (       u(2,i-2,j,k)
     $                      - 4.0d+00 * u(2,i-1,j,k)
     $                      + 6.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i+1,j,k)  )
         rsd3ijk(i) = rsd3ijk(i) - dssp * (       u(3,i-2,j,k)
     $                      - 4.0d+00 * u(3,i-1,j,k)
     $                      + 6.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i+1,j,k)  )
         rsd4ijk(i) = rsd4ijk(i) - dssp * (       u(4,i-2,j,k)
     $                      - 4.0d+00 * u(4,i-1,j,k)
     $                      + 6.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i+1,j,k)  )
         rsd5ijk(i) = rsd5ijk(i) - dssp * (       u(5,i-2,j,k)
     $                      - 4.0d+00 * u(5,i-1,j,k)
     $                      + 6.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i+1,j,k)  )
c
         i = nx -1
c
         rsd1ijk(i) = rsd1ijk(i) - dssp * (         u(1,i-2,j,k)
     $                      - 4.0d+00 * u(1,i-1,j,k)
     $                      + 5.0d+00 * u(1,i,j,k) )
         rsd2ijk(i) = rsd2ijk(i) - dssp * (         u(2,i-2,j,k)
     $                      - 4.0d+00 * u(2,i-1,j,k)
     $                      + 5.0d+00 * u(2,i,j,k) )
         rsd3ijk(i) = rsd3ijk(i) - dssp * (         u(3,i-2,j,k)
     $                      - 4.0d+00 * u(3,i-1,j,k)
     $                      + 5.0d+00 * u(3,i,j,k) )
         rsd4ijk(i) = rsd4ijk(i) - dssp * (         u(4,i-2,j,k)
     $                      - 4.0d+00 * u(4,i-1,j,k)
     $                      + 5.0d+00 * u(4,i,j,k) )
         rsd5ijk(i) = rsd5ijk(i) - dssp * (         u(5,i-2,j,k)
     $                      - 4.0d+00 * u(5,i-1,j,k)
     $                      + 5.0d+00 * u(5,i,j,k) )
c
        do i = 4, nx -3
c
         rsd1ijk(i) = rsd1ijk(i) - dssp * (        u(1,i-2,j,k)
     $                        - 4.0d+00 * u(1,i-1,j,k)
     $                        + 6.0d+00 * u(1,i,j,k)
     $                        - 4.0d+00 * u(1,i+1,j,k)
     $                        +           u(1,i+2,j,k) )
c
c
         rsd2ijk(i) = rsd2ijk(i) - dssp * (        u(2,i-2,j,k)
     $                        - 4.0d+00 * u(2,i-1,j,k)
     $                        + 6.0d+00 * u(2,i,j,k)
     $                        - 4.0d+00 * u(2,i+1,j,k)
     $                        +           u(2,i+2,j,k) )
c
c
         rsd3ijk(i) = rsd3ijk(i) - dssp * (        u(3,i-2,j,k)
     $                        - 4.0d+00 * u(3,i-1,j,k)
     $                        + 6.0d+00 * u(3,i,j,k)
     $                        - 4.0d+00 * u(3,i+1,j,k)
     $                        +           u(3,i+2,j,k) )
c
c
         rsd4ijk(i) = rsd4ijk(i) - dssp * (        u(4,i-2,j,k)
     $                        - 4.0d+00 * u(4,i-1,j,k)
     $                        + 6.0d+00 * u(4,i,j,k)
     $                        - 4.0d+00 * u(4,i+1,j,k)
     $                        +           u(4,i+2,j,k) )
c
c
         rsd5ijk(i) = rsd5ijk(i) - dssp * (        u(5,i-2,j,k)
     $                        - 4.0d+00 * u(5,i-1,j,k)
     $                        + 6.0d+00 * u(5,i,j,k)
     $                        - 4.0d+00 * u(5,i+1,j,k)
     $                        +           u(5,i+2,j,k) )
c
        end do
c
c
c
            do i = 2, nx - 1
c
c
c***eta-direction flux differences
c
               flux1jp1(i) = u(3,i,j+1,k)
c
               u31 = u(3,i,j+1,k) / u(1,i,j+1,k)
c
               q = 0.50d+00 * (  u(2,i,j+1,k) * u(2,i,j+1,k)
     $                         + u(3,i,j+1,k) * u(3,i,j+1,k)
     $                         + u(4,i,j+1,k) * u(4,i,j+1,k) )
     $                      / u(1,i,j+1,k)
c
               flux2jp1(i) = u(2,i,j+1,k) * u31 
c
       flux3jp1(i) = u(3,i,j+1,k) * u31 + c2 * ( u(5,i,j+1,k) - q )
c
               flux4jp1(i) = u(4,i,j+1,k) * u31
c
               flux5jp1(i) = ( c1 * u(5,i,j+1,k) - c2 * q ) * u31
               tmp = 1.0d+00 / u(1,i,j+1,k)
c
               u21jp1(i) = tmp * u(2,i,j+1,k)
               u31jp1(i) = tmp * u(3,i,j+1,k)
               u41jp1(i) = tmp * u(4,i,j+1,k)
               u51jp1(i) = tmp * u(5,i,j+1,k)
c
               flx2jp1(i) = ty3 * ( u21jp1(i) - u21j(i) )
       flx3jp1(i) = (4.0d+00/3.0d+00) * ty3 * ( u31jp1(i) - u31j(i) )
               flx4jp1(i) = ty3 * ( u41jp1(i) - u41j(i) )
               flx5jp1(i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $   * ty3 * ( ( u21jp1(i)  **2 + u31jp1(i)  **2 + u41jp1(i)  **2 )
     $      - ( u21j(i)**2 + u31j(i)**2 + u41j(i)**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * ty3 * ( u31jp1(i)**2 - u31j(i)**2 )
     $              + c1 * c5 * ty3 * ( u51jp1(i) - u51j(i) )
c
                  rsd1 =  
     $                      - ty2 * ( flux1jp1(i) - flux1jm1(i) )
                  rsd2 = 
     $                      - ty2 * ( flux2jp1(i) - flux2jm1(i) )
                  rsd3 = 
     $                      - ty2 * ( flux3jp1(i) - flux3jm1(i) )
                  rsd4 = 
     $                      - ty2 * ( flux4jp1(i) - flux4jm1(i) )
                  rsd5 = 
     $                      - ty2 * ( flux5jp1(i) - flux5jm1(i) )
c
c
                rsd1ijk(i)= rsd1ijk(i) + rsd1
     $              + dy1 * ty1 * (        u(1,i,j-1,k)
     $                             - 2.0d+00 * u(1,i,j,k)
     $                             +           u(1,i,j+1,k) )
c
               rsd2ijk(i) = rsd2ijk(i) + rsd2
     $              + ty3 * c3 * c4 * ( flx2jp1(i) - flx2j(i) )
     $              + dy2 * ty1 * (        u(2,i,j-1,k)
     $                             - 2.0d+00 * u(2,i,j,k)
     $                             +           u(2,i,j+1,k) )
c
               rsd3ijk(i) = rsd3ijk(i) + rsd3
     $              + ty3 * c3 * c4 * ( flx3jp1(i) - flx3j(i) )
     $              + dy3 * ty1 * (        u(3,i,j-1,k)
     $                             - 2.0d+00 * u(3,i,j,k)
     $                             +           u(3,i,j+1,k) )
c
               rsd4ijk(i) = rsd4ijk(i) + rsd4
     $              + ty3 * c3 * c4 * ( flx4jp1(i) - flx4j(i) )
     $              + dy4 * ty1 * (        u(4,i,j-1,k)
     $                             - 2.0d+00 * u(4,i,j,k)
     $                             +           u(4,i,j+1,k) )
c
               rsd5ijk(i) = rsd5ijk(i) + rsd5
     $              + ty3 * c3 * c4 * ( flx5jp1(i) - flx5j(i) )
     $              + dy5 * ty1 * (        u(5,i,j-1,k)
     $                             - 2.0d+00 * u(5,i,j,k)
     $                             +           u(5,i,j+1,k) )
      end do
c
c***fourth-order dissipation
c
      if ( j .eq. 2) then
      do i = 2,nx -1
c
      rsd1ijk(i)=rsd1ijk(i)      - dssp * ( + 5.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i,j+1,k)
     $                      +           u(1,i,j+2,k) )
 
      rsd2ijk(i)=rsd2ijk(i)      - dssp * ( + 5.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i,j+1,k)
     $                      +           u(2,i,j+2,k) )
 
      rsd3ijk(i)=rsd3ijk(i)      - dssp * ( + 5.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i,j+1,k)
     $                      +           u(3,i,j+2,k) )
 
      rsd4ijk(i)=rsd4ijk(i)      - dssp * ( + 5.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i,j+1,k)
     $                      +           u(4,i,j+2,k) )
 
      rsd5ijk(i)=rsd5ijk(i)      - dssp * ( + 5.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i,j+1,k)
     $                      +           u(5,i,j+2,k) )
c
      end do
      else if ( j .eq. 3) then
      do i = 2,nx -1
c
      rsd1ijk(i)=rsd1ijk(i)      - dssp * ( - 4.0d+00 * u(1,i,j-1,k)
     $                      + 6.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i,j+1,k)
     $                      +           u(1,i,j+2,k) )
 
      rsd2ijk(i)=rsd2ijk(i)      - dssp * ( - 4.0d+00 * u(2,i,j-1,k)
     $                      + 6.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i,j+1,k)
     $                      +           u(2,i,j+2,k) )
 
      rsd3ijk(i)=rsd3ijk(i)      - dssp * ( - 4.0d+00 * u(3,i,j-1,k)
     $                      + 6.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i,j+1,k)
     $                      +           u(3,i,j+2,k) )
 
      rsd4ijk(i)=rsd4ijk(i)      - dssp * ( - 4.0d+00 * u(4,i,j-1,k)
     $                      + 6.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i,j+1,k)
     $                      +           u(4,i,j+2,k) )
 
      rsd5ijk(i)=rsd5ijk(i)      - dssp * ( - 4.0d+00 * u(5,i,j-1,k)
     $                      + 6.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i,j+1,k)
     $                      +           u(5,i,j+2,k) )
c
      end do
      else if ( j.eq. ny-2) then
      do i = 2,nx -1
c
      rsd1ijk(i)=rsd1ijk(i)      - dssp * (         u(1,i,j-2,k)
     $                      - 4.0d+00 * u(1,i,j-1,k)
     $                      + 6.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i,j+1,k)  )
 
      rsd2ijk(i)=rsd2ijk(i)      - dssp * (         u(2,i,j-2,k)
     $                      - 4.0d+00 * u(2,i,j-1,k)
     $                      + 6.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i,j+1,k)  )
 
      rsd3ijk(i)=rsd3ijk(i)      - dssp * (         u(3,i,j-2,k)
     $                      - 4.0d+00 * u(3,i,j-1,k)
     $                      + 6.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i,j+1,k)  )
 
      rsd4ijk(i)=rsd4ijk(i)      - dssp * (         u(4,i,j-2,k)
     $                      - 4.0d+00 * u(4,i,j-1,k)
     $                      + 6.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i,j+1,k)  )
 
      rsd5ijk(i)=rsd5ijk(i)      - dssp * (         u(5,i,j-2,k)
     $                      - 4.0d+00 * u(5,i,j-1,k)
     $                      + 6.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i,j+1,k)  )
c
      end do
      else if ( j.eq. ny-1) then
      do i = 2,nx -1
c
      rsd1ijk(i)=rsd1ijk(i)      - dssp * (         u(1,i,j-2,k)
     $                      - 4.0d+00 * u(1,i,j-1,k)
     $                      + 5.0d+00 * u(1,i,j,k) )
 
      rsd2ijk(i)=rsd2ijk(i)      - dssp * (         u(2,i,j-2,k)
     $                      - 4.0d+00 * u(2,i,j-1,k)
     $                      + 5.0d+00 * u(2,i,j,k) )
 
      rsd3ijk(i)=rsd3ijk(i)      - dssp * (         u(3,i,j-2,k)
     $                      - 4.0d+00 * u(3,i,j-1,k)
     $                      + 5.0d+00 * u(3,i,j,k) )
 
      rsd4ijk(i)=rsd4ijk(i)      - dssp * (         u(4,i,j-2,k)
     $                      - 4.0d+00 * u(4,i,j-1,k)
     $                      + 5.0d+00 * u(4,i,j,k) )
 
      rsd5ijk(i)=rsd5ijk(i)      - dssp * (         u(5,i,j-2,k)
     $                      - 4.0d+00 * u(5,i,j-1,k)
     $                      + 5.0d+00 * u(5,i,j,k) )
c
      end do
       else
      do i = 2,nx -1
c
 
      rsd1ijk(i)=rsd1ijk(i)         - dssp * (        u(1,i,j-2,k)
     $                        - 4.0d+00 * u(1,i,j-1,k)
     $                        + 6.0d+00 * u(1,i,j,k)
     $                        - 4.0d+00 * u(1,i,j+1,k)
     $                        +           u(1,i,j+2,k) )
 
      rsd2ijk(i)=rsd2ijk(i)         - dssp * (        u(2,i,j-2,k)
     $                        - 4.0d+00 * u(2,i,j-1,k)
     $                        + 6.0d+00 * u(2,i,j,k)
     $                        - 4.0d+00 * u(2,i,j+1,k)
     $                        +           u(2,i,j+2,k) ) 
      rsd3ijk(i)=rsd3ijk(i)         - dssp * (        u(3,i,j-2,k)
     $                        - 4.0d+00 * u(3,i,j-1,k)
     $                        + 6.0d+00 * u(3,i,j,k)
     $                        - 4.0d+00 * u(3,i,j+1,k)
     $                        +           u(3,i,j+2,k) ) 
      rsd4ijk(i)=rsd4ijk(i)         - dssp * (        u(4,i,j-2,k)
     $                        - 4.0d+00 * u(4,i,j-1,k)
     $                        + 6.0d+00 * u(4,i,j,k)
     $                        - 4.0d+00 * u(4,i,j+1,k)
     $                        +           u(4,i,j+2,k) ) 
      rsd5ijk(i)=rsd5ijk(i)         - dssp * (        u(5,i,j-2,k)
     $                        - 4.0d+00 * u(5,i,j-1,k)
     $                        + 6.0d+00 * u(5,i,j,k)
     $                        - 4.0d+00 * u(5,i,j+1,k)
     $                        +           u(5,i,j+2,k) )
c
         end do
       end if
c
            do i = 2, nx - 1
c
               flux1km1 = u(4,i,j,k-1)
c
               u41 = u(4,i,j,k-1) / u(1,i,j,k-1)
c
               q = 0.50d+00 * (  u(2,i,j,k-1) * u(2,i,j,k-1)
     $                         + u(3,i,j,k-1) * u(3,i,j,k-1)
     $                         + u(4,i,j,k-1) * u(4,i,j,k-1) )
     $                      / u(1,i,j,k-1)
c
               flux2km1 = u(2,i,j,k-1) * u41 
c
               flux3km1 = u(3,i,j,k-1) * u41
c
               flux4km1 = u(4,i,j,k-1) * u41+ c2 * ( u(5,i,j,k-1) - q )
c
               flux5km1 = ( c1 * u(5,i,j,k-1) - c2 * q ) * u41
c 
c
               flux1k = 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
               flux2k = u(2,i,j,k) * u41 
c
               flux3k = u(3,i,j,k) * u41
c
               flux4k = u(4,i,j,k) * u41+ c2 * ( u(5,i,j,k) - q )
c
               flux5k = ( c1 * u(5,i,j,k) - c2 * q ) * u41
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
               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
               flx4k = (4.0d+00/3.0d+00)*tz3*
     $      (u41k-u41km1)
               flx3k = tz3 * ( u31k - u31km1 )
               flx2k = tz3 * ( u21k - u21km1 )
               flx5k = 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
c***zeta-direction flux differences
c
               flux1kp1 = u(4,i,j,k+1)
c
               u41 = u(4,i,j,k+1) / u(1,i,j,k+1)
c
               q = 0.50d+00 * (  u(2,i,j,k+1) * u(2,i,j,k+1)
     $                         + u(3,i,j,k+1) * u(3,i,j,k+1)
     $                         + u(4,i,j,k+1) * u(4,i,j,k+1) )
     $                      / u(1,i,j,k+1)
c
               flux2kp1 = u(2,i,j,k+1) * u41 
c
               flux3kp1 = u(3,i,j,k+1) * u41 
c
       flux4kp1 = u(4,i,j,k+1) * u41 + c2 * ( u(5,i,j,k+1) - q )
c
               flux5kp1 = ( c1 * u(5,i,j,k+1) - c2 * q ) * u41
               tmp = 1.0d+00 / u(1,i,j,k+1)
c
               u21kp1 = tmp * u(2,i,j,k+1)
               u31kp1 = tmp * u(3,i,j,k+1)
               u41kp1 = tmp * u(4,i,j,k+1)
               u51kp1 = tmp * u(5,i,j,k+1)
c
               flx2kp1 = tz3 * ( u21kp1 - u21k )
               flx3kp1 = tz3 * ( u31kp1 - u31k )
      flx4kp1=(4.0d+00/3.0d+00)*tz3*(u41kp1-u41k )
      flx5kp1 = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     $  * tz3*((u21kp1**2 + u31kp1**2+u41kp1  **2 )
     $  - ( u21k**2 + u31k**2 + u41k**2 ) )
     $              + (1.0d+00/6.0d+00)
     $              * tz3 * ( u41kp1**2 - u41k**2 )
     $              + c1 * c5 * tz3 * ( u51kp1 - u51k )
c    
                  rsd1 = 
     $                       - tz2 * ( flux1kp1 - flux1km1 )
                  rsd2 =
     $                       - tz2 * ( flux2kp1 - flux2km1 )
                  rsd3 =
     $                       - tz2 * ( flux3kp1 - flux3km1 )
                  rsd4 = 
     $                       - tz2 * ( flux4kp1 - flux4km1 )
                  rsd5 =
     $                       - tz2 * ( flux5kp1 - flux5km1 )
c
c
               rsd1ijk(i) = rsd1ijk(i) + rsd1
     $              + dz1 * tz1 * (        u(1,i,j,k-1)
     $                             - 2.0d+00 * u(1,i,j,k)
     $                             +           u(1,i,j,k+1) )
c
               rsd2ijk(i) = rsd2ijk(i) + rsd2
     $            + tz3 * c3 * c4 * ( flx2kp1 - flx2k )
     $              + dz2 * tz1 * (        u(2,i,j,k-1)
     $                             - 2.0d+00 * u(2,i,j,k)
     $                             +           u(2,i,j,k+1) )
c
               rsd3ijk(i) = rsd3ijk(i) + rsd3
     $              + tz3 * c3 * c4 * ( flx3kp1 - flx3k )
     $              + dz3 * tz1 * (        u(3,i,j,k-1)
     $                             - 2.0d+00 * u(3,i,j,k)
     $                             +           u(3,i,j,k+1) )
c
               rsd4ijk(i) = rsd4ijk(i) + rsd4
     $              + tz3 * c3 * c4 * ( flx4kp1 - flx4k )
     $              + dz4 * tz1 * (        u(4,i,j,k-1)
     $                             - 2.0d+00 * u(4,i,j,k)
     $                             +           u(4,i,j,k+1) )
c
               rsd5ijk(i) = rsd5ijk(i) + rsd5
     $              + tz3 * c3 * c4 * ( flx5kp1 - flx5k )
     $              + 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***fourth-order dissipation
c
      if ( k.eq.2 ) then
      do i = 2,nx-1
      rsd1ijk(i)=rsd1ijk(i)      - dssp * ( + 5.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i,j,k+1)
     $                      +           u(1,i,j,k+2) )
      rsd2ijk(i)=rsd2ijk(i)      - dssp * ( + 5.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i,j,k+1)
     $                      +           u(2,i,j,k+2) )
      rsd3ijk(i)=rsd3ijk(i)      - dssp * ( + 5.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i,j,k+1)
     $                      +           u(3,i,j,k+2) )
      rsd4ijk(i)=rsd4ijk(i)      - dssp * ( + 5.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i,j,k+1)
     $                      +           u(4,i,j,k+2) )
      rsd5ijk(i)=rsd5ijk(i)      - dssp * ( + 5.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i,j,k+1)
     $                      +           u(5,i,j,k+2) )
c
      end do
      else if ( k.eq.3 ) then
      do i = 2,nx-1
       rsd1ijk(i)=rsd1ijk(i)      - dssp * ( - 4.0d+00 * u(1,i,j,k-1)
     $                      + 6.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i,j,k+1)
     $                      +           u(1,i,j,k+2) )
       rsd2ijk(i)=rsd2ijk(i)      - dssp * ( - 4.0d+00 * u(2,i,j,k-1)
     $                      + 6.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i,j,k+1)
     $                      +           u(2,i,j,k+2) )
       rsd3ijk(i)=rsd3ijk(i)      - dssp * ( - 4.0d+00 * u(3,i,j,k-1)
     $                      + 6.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i,j,k+1)
     $                      +           u(3,i,j,k+2) )
       rsd4ijk(i)=rsd4ijk(i)      - dssp * ( - 4.0d+00 * u(4,i,j,k-1)
     $                      + 6.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i,j,k+1)
     $                      +           u(4,i,j,k+2) )
       rsd5ijk(i)=rsd5ijk(i)      - dssp * ( - 4.0d+00 * u(5,i,j,k-1)
     $                      + 6.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i,j,k+1)
     $                      +           u(5,i,j,k+2) )
c
      end do
       else if( k.eq.nz-2) then
      do i = 2,nx-1
       rsd1ijk(i)=rsd1ijk(i)      - dssp * (         u(1,i,j,k-2)
     $                      - 4.0d+00 * u(1,i,j,k-1)
     $                      + 6.0d+00 * u(1,i,j,k)
     $                      - 4.0d+00 * u(1,i,j,k+1)  )
       rsd2ijk(i)=rsd2ijk(i)      - dssp * (         u(2,i,j,k-2)
     $                      - 4.0d+00 * u(2,i,j,k-1)
     $                      + 6.0d+00 * u(2,i,j,k)
     $                      - 4.0d+00 * u(2,i,j,k+1)  )
       rsd3ijk(i)=rsd3ijk(i)      - dssp * (         u(3,i,j,k-2)
     $                      - 4.0d+00 * u(3,i,j,k-1)
     $                      + 6.0d+00 * u(3,i,j,k)
     $                      - 4.0d+00 * u(3,i,j,k+1)  )
       rsd4ijk(i)=rsd4ijk(i)      - dssp * (         u(4,i,j,k-2)
     $                      - 4.0d+00 * u(4,i,j,k-1)
     $                      + 6.0d+00 * u(4,i,j,k)
     $                      - 4.0d+00 * u(4,i,j,k+1)  )
       rsd5ijk(i)=rsd5ijk(i)      - dssp * (         u(5,i,j,k-2)
     $                      - 4.0d+00 * u(5,i,j,k-1)
     $                      + 6.0d+00 * u(5,i,j,k)
     $                      - 4.0d+00 * u(5,i,j,k+1)  )
c
       end do
       else if( k.eq.nz-1 ) then
       do i = 2,nx-1
       rsd1ijk(i)=rsd1ijk(i)      - dssp * (         u(1,i,j,k-2)
     $                      - 4.0d+00 * u(1,i,j,k-1)
     $                      + 5.0d+00 * u(1,i,j,k) )
       rsd2ijk(i)=rsd2ijk(i)      - dssp * (         u(2,i,j,k-2)
     $                      - 4.0d+00 * u(2,i,j,k-1)
     $                      + 5.0d+00 * u(2,i,j,k) )
       rsd3ijk(i)=rsd3ijk(i)      - dssp * (         u(3,i,j,k-2)
     $                      - 4.0d+00 * u(3,i,j,k-1)
     $                      + 5.0d+00 * u(3,i,j,k) )
       rsd4ijk(i)=rsd4ijk(i)      - dssp * (         u(4,i,j,k-2)
     $                      - 4.0d+00 * u(4,i,j,k-1)
     $                      + 5.0d+00 * u(4,i,j,k) )
       rsd5ijk(i)=rsd5ijk(i)      - dssp * (         u(5,i,j,k-2)
     $                      - 4.0d+00 * u(5,i,j,k-1)
     $                      + 5.0d+00 * u(5,i,j,k) )
c
      end do
       else
      do i = 2,nx-1
       rsd1ijk(i)=rsd1ijk(i)         - dssp * (        u(1,i,j,k-2)
     $                        - 4.0d+00 * u(1,i,j,k-1)
     $                        + 6.0d+00 * u(1,i,j,k)
     $                        - 4.0d+00 * u(1,i,j,k+1)
     $                        +           u(1,i,j,k+2) )
       rsd2ijk(i)=rsd2ijk(i)         - dssp * (        u(2,i,j,k-2)
     $                        - 4.0d+00 * u(2,i,j,k-1)
     $                        + 6.0d+00 * u(2,i,j,k)
     $                        - 4.0d+00 * u(2,i,j,k+1)
     $                        +           u(2,i,j,k+2) )
       rsd3ijk(i)=rsd3ijk(i)         - dssp * (        u(3,i,j,k-2)
     $                        - 4.0d+00 * u(3,i,j,k-1)
     $                        + 6.0d+00 * u(3,i,j,k)
     $                        - 4.0d+00 * u(3,i,j,k+1)
     $                        +           u(3,i,j,k+2) )
       rsd4ijk(i)=rsd4ijk(i)         - dssp * (        u(4,i,j,k-2)
     $                        - 4.0d+00 * u(4,i,j,k-1)
     $                        + 6.0d+00 * u(4,i,j,k)
     $                        - 4.0d+00 * u(4,i,j,k+1)
     $                        +           u(4,i,j,k+2) )
       rsd5ijk(i)=rsd5ijk(i)         - dssp * (        u(5,i,j,k-2)
     $                        - 4.0d+00 * u(5,i,j,k-1)
     $                        + 6.0d+00 * u(5,i,j,k)
     $                        - 4.0d+00 * u(5,i,j,k+1)
     $                        +           u(5,i,j,k+2) )
c
       end do
      end if
         do i = 2,nx-1
               rsd(1,i,j,k) =  rsd1ijk(i)
               rsd(2,i,j,k) =  rsd2ijk(i)
               rsd(3,i,j,k) =  rsd3ijk(i)
               rsd(4,i,j,k) =  rsd4ijk(i)
               rsd(5,i,j,k) =  rsd5ijk(i)
c
c
c   end of major i loop
c
         end do
            do i = 2, nx - 1
c
               flux1jm1(i) = flux1j(i)
c
               flux2jm1(i) = flux2j(i)
c
               flux3jm1(i) = flux3j(i)
c
               flux4jm1(i) = flux4j(i)
c
               flux5jm1(i) = flux5j(i)
c
               flux1j(i) = flux1jp1(i)
c
               flux2j(i) = flux2jp1(i)
c
               flux3j(i) = flux3jp1(i)
c
               flux4j(i) = flux4jp1(i)
c
               flux5j(i) = flux5jp1(i)
c
               u21jm1(i) = u21j(i)
               u31jm1(i) = u31j(i)
               u41jm1(i) = u41j(i)
               u51jm1(i) = u51j(i)
c
               u21j(i) = u21jp1(i)
               u31j(i) = u31jp1(i)
               u41j(i) = u41jp1(i)
               u51j(i) = u51jp1(i)
c
               flx2j(i) = flx2jp1(i)
               flx3j(i) = flx3jp1(i)
               flx4j(i) = flx4jp1(i)
               flx5j(i) = flx5jp1(i)
c
            end do
c
c   end of major j loop
         end do
c   end of major k loop
c
      end do
c
      return
      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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c

      dimension t1(5),t2(5)




c
c***set the dependent variable values along the top and bottom faces
c
CHPF$ INDEPENDENT
      do j = 1, ny
         do i = 1, nx
c
            call exact ( i, j, 1, t1)
            call exact ( i, j, nz, t2)
         u(1, i, j, 1 ) = t1(1)  
         u(2, i, j, 1 ) = t1(2)  
         u(3, i, j, 1 ) = t1(3)  
         u(4, i, j, 1 ) = t1(4)  
         u(5, i, j, 1 ) = t1(5)  
         u(1, i, j, nz ) = t2(1)
         u(2, i, j, nz ) = t2(2)
         u(3, i, j, nz ) = t2(3)
         u(4, i, j, nz ) = t2(4)
         u(5, i, j, nz ) = t2(5)
c
         end do
      end do
c
c***set the dependent variable values along north and south faces
c
CHPF$ INDEPENDENT
      do k = 1, nz
         do i = 1, nx
c
            call exact ( i, 1, k, t1)
            call exact ( i, ny, k, t2)
         u(1, i, 1,k ) = t1(1)  
         u(2, i, 1,k ) = t1(2)  
         u(3, i, 1,k ) = t1(3)  
         u(4, i, 1,k ) = t1(4)  
         u(5, i, 1,k ) = t1(5)  
         u(1, i, ny,k ) = t2(1)
         u(2, i, ny,k ) = t2(2)
         u(3, i, ny,k ) = t2(3)
         u(4, i, ny,k ) = t2(4)
         u(5, i, ny,k ) = t2(5)
c
         end do
      end do
c
c***set the dependent variable values along east and west faces
c
CHPF$ INDEPENDENT
      do k = 1, nz
         do j = 1, ny
c
            call exact ( 1, j, k, t1)
            call exact ( nx, j, k, t2)
         u(1, 1,j,k ) = t1(1)  
         u(2, 1,j,k ) = t1(2)  
         u(3, 1,j,k ) = t1(3)  
         u(4, 1,j,k ) = t1(4)  
         u(5, 1,j,k ) = t1(5)  
         u(1, nx,j,k ) = t2(1)
         u(2, nx,j,k ) = t2(2)
         u(3, nx,j,k ) = t2(3)
         u(4, nx,j,k ) = t2(4)
         u(5, nx,j,k ) = t2(5)
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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
c Author: Sisira Weeratunga
c         NASA Ames Research Center
c         (10/25/90)
c
CHPF$ INDEPENDENT
      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
c
                  pxi =   ( 1.0d+00 - xi ) * u(1,1,j,k)
     $                              + xi   * u(1,nx,j,k)
c
                  peta =  ( 1.0d+00 - eta ) * u(1,i,1,k)
     $                              + eta   * u(1,i,ny,k)
c
                  pzeta = ( 1.0d+00 - zeta ) * u(1,i,j,1)
     $                              + zeta   * u(1,i,j,nz)
c
                  u(1, i, j, k ) = pxi + peta + pzeta
     $                 - pxi * peta - peta * pzeta - pzeta * pxi
     $                 + pxi * peta * pzeta
c
c
                  pxi =   ( 1.0d+00 - xi ) * u(2,1,j,k)
     $                              + xi   * u(2,nx,j,k)
c
                  peta =  ( 1.0d+00 - eta ) * u(2,i,1,k)
     $                              + eta   * u(2,i,ny,k)
c
                  pzeta = ( 1.0d+00 - zeta ) * u(2,i,j,1)
     $                              + zeta   * u(2,i,j,nz)
c
                  u(2, i, j, k ) = pxi + peta + pzeta
     $                 - pxi * peta - peta * pzeta - pzeta * pxi
     $                 + pxi * peta * pzeta
c
c
                  pxi =   ( 1.0d+00 - xi ) * u(3,1,j,k)
     $                              + xi   * u(3,nx,j,k)
c
                  peta =  ( 1.0d+00 - eta ) * u(3,i,1,k)
     $                              + eta   * u(3,i,ny,k)
c
                  pzeta = ( 1.0d+00 - zeta ) * u(3,i,j,1)
     $                              + zeta   * u(3,i,j,nz)
c
                  u(3, i, j, k ) = pxi + peta + pzeta
     $                 - pxi * peta - peta * pzeta - pzeta * pxi
     $                 + pxi * peta * pzeta
c
c
                  pxi =   ( 1.0d+00 - xi ) * u(4,1,j,k)
     $                              + xi   * u(4,nx,j,k)
c
                  peta =  ( 1.0d+00 - eta ) * u(4,i,1,k)
     $                              + eta   * u(4,i,ny,k)
c
                  pzeta = ( 1.0d+00 - zeta ) * u(4,i,j,1)
     $                              + zeta   * u(4,i,j,nz)
c
                  u(4, i, j, k ) = pxi + peta + pzeta
     $                 - pxi * peta - peta * pzeta - pzeta * pxi
     $                 + pxi * peta * pzeta
c
c
                  pxi =   ( 1.0d+00 - xi ) * u(5,1,j,k)
     $                              + xi   * u(5,nx,j,k)
c
                  peta =  ( 1.0d+00 - eta ) * u(5,i,1,k)
     $                              + eta   * u(5,i,ny,k)
c
                  pzeta = ( 1.0d+00 - zeta ) * u(5,i,j,1)
     $                              + zeta   * u(5,i,j,nz)
c
                  u(5, i, j, k ) = pxi + peta + pzeta
     $                 - pxi * peta - peta * pzeta - pzeta * pxi
     $                 + pxi * peta * pzeta
c
c
            end do
c
         end do
c
      end do
c
      return
      end
      function timer(x)
      real *8 second, timer
c     real secs,tarray(2)
c
c   
c     sun version of timing routine
cccccccccccccccccccccccccccccccccccccc
c     secs = etime(tarray)
c     timer = secs
cccccccccccccccccccccccccccccccccccccc
c
c       cray version of timing routine
cccccccccccccccccccccccccccccccccccccc
       timer = second ()
cccccccccccccccccccccccccccccccccccccc
c
c       sgi version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      secs = secnds(0.0)
c     timer = secs
cccccccccccccccccccccccccccccccccccccc
      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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
      dimension xcr(5), xce(5),
     $          xrr(5), xre(5)
c
c***tolerance level
c
      epsilon = 1.0d-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 100 time steps, with  DT = 1.50e-02
c
         xrr(1) = 2.7470315451339479d-02
         xrr(2) = 1.0360746705285417d-02
         xrr(3) = 1.6235745065095532d-02
         xrr(4) = 1.5840557224455615d-02
         xrr(5) = 3.4849040609362460d-02
c
c***Reference values of RMS-norms of solution error, for the (12X12X12) grid,
c   after 100 time steps, with  DT = 1.50e-02
c
         xre(1) = 2.7289258557377227d-05
         xre(2) = 1.0364446640837285d-05
         xre(3) = 1.6154798287166471d-05
         xre(4) = 1.5750704994480102d-05
         xre(5) = 3.4177666183390531d-05
c
c***Reference value of surface integral, for the (12X12X12) grid,
c   after 100 time steps, with  DT = 1.50e-02
c
         xri = 7.8406293309126962e+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 = 100',
     $ //5x,'DT = 1.5e-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 400 time steps, with  DT = 1.5e-03
c
         xrr(1) = 2.4799822399300195d+00
         xrr(2) = 1.1276337964368832d+00
         xrr(3) = 1.5028977888770491d+00
         xrr(4) = 1.4217816211695179d+00
         xrr(5) = 2.1292113035138280d+00
c     
c***Reference values of RMS-norms of solution error, for the (64X64X64) grid,
c   after 400 time steps, with  DT = 1.5e-03
c
         xre(1) = 1.0900140297820550d-04
         xre(2) = 3.7343951769282091d-05
         xre(3) = 5.0092785406541633d-05
         xre(4) = 4.7671093939528255d-05
         xre(5) = 1.3621613399213001d-04
c
c***Reference value of surface integral, for the (64X64X64) grid,
c   after 400 time steps, with DT = 1.5e-03
c
         xri = 1.2080439685038630d+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 = 400',
     $ //5x,'DT = 1.5e-03',
     $ //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
c
c
      subroutine xisweep
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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c
c
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
      dimension idmax(5), jdmax(5), kdmax(5),
     $          imax(5), jmax(5), kmax(5),
     $          delunm(5)
c
      dimension flux1(isiz1)
      dimension flux2(isiz1)
      dimension flux3(isiz1)
      dimension flux4(isiz1)
      dimension flux5(isiz1)
c
      dimension cv(isiz3), aa(isiz3), rhos(isiz3)
      dimension  rhon(isiz1)
      dimension  rhoq(isiz2)
c
      parameter ( one = 1.0d+00 )
c
      lnorm = 2
      bt = sqrt ( 0.50d+00 )
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c***perform the xsi-direction sweep
c
c        call jacx ( 3 )
c
CHPF$ INDEPENDENT
         do k = 2, nz -1
		  do j = 2, ny -1
            do i = 1, nx
c
               ru1 = 1.0d+00 / u(1,i,j,k)
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(i) = uu
               aa(i) = sqrt ( c1 * c2 * ( ru1 * u(5,i,j,k) - q ) )
c
               rhon(i) = max ( dx1,
     $                         dx2 + r43 * c34 * ru1,
     $                         dx3 + c34 * ru1,
     $                         dx4 + c34 * ru1,
     $                         dx5 + c1345 * ru1 )
c
            end do
c
            a(1,j,k) = 0.0d+00
            b(1,j,k) = 0.0d+00
            c(1,j,k) = 1.0d+00
            d(1,j,k) = 0.0d+00
            e(1,j,k) = 0.0d+00
c
            do i = 2, nx - 1
c
               a(i,j,k) =   0.0d+00
               b(i,j,k) = - dt * tx2 *  cv(i-1) 
     $                    - dt * rhon(i-1) * tx1 
               c(i,j,k) =   1.0d+00
     $                    + dt * rhon(i) * tx1 * 2.0d+00
               d(i,j,k) =   dt * tx2 *  cv(i+1)
     $                    - dt * rhon(i+1) * tx1
               e(i,j,k) =   0.0d+00
c
            end do
c
            a(nx,j,k) = 0.0d+00
            b(nx,j,k) = 0.0d+00
            c(nx,j,k) = 1.0d+00
            d(nx,j,k) = 0.0d+00
            e(nx,j,k) = 0.0d+00
c
c***fourth order dissipation 
c
            c(2,j,k) = c(2,j,k) + dt * dssp * ( + 5.0d+00 )
            d(2,j,k) = d(2,j,k) + dt * dssp * ( - 4.0d+00 )
            e(2,j,k) = e(2,j,k) + dt * dssp * ( + 1.0d+00 )
c
            b(3,j,k) = b(3,j,k) + dt * dssp * ( - 4.0d+00 )
            c(3,j,k) = c(3,j,k) + dt * dssp * ( + 6.0d+00 )
            d(3,j,k) = d(3,j,k) + dt * dssp * ( - 4.0d+00 )
            e(3,j,k) = e(3,j,k) + dt * dssp * ( + 1.0d+00 )
c
            do i = 4, nx - 3
c
               a(i,j,k) = a(i,j,k) + dt * dssp * ( + 1.0d+00 )
               b(i,j,k) = b(i,j,k) + dt * dssp * ( - 4.0d+00 )
               c(i,j,k) = c(i,j,k) + dt * dssp * ( + 6.0d+00 )
               d(i,j,k) = d(i,j,k) + dt * dssp * ( - 4.0d+00 )
               e(i,j,k) = e(i,j,k) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(nx-2,j,k) = a(nx-2,j,k) + dt * dssp * ( + 1.0d+00 )
            b(nx-2,j,k) = b(nx-2,j,k) + dt * dssp * ( - 4.0d+00 )
            c(nx-2,j,k) = c(nx-2,j,k) + dt * dssp * ( + 6.0d+00 )
            d(nx-2,j,k) = d(nx-2,j,k) + dt * dssp * ( - 4.0d+00 )
c
            a(nx-1,j,k) = a(nx-1,j,k) + dt * dssp * ( + 1.0d+00 )
            b(nx-1,j,k) = b(nx-1,j,k) + dt * dssp * ( - 4.0d+00 )
            c(nx-1,j,k) = c(nx-1,j,k) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentax3 ( isiz1, isiz2, isiz3,
c    $                   nx, ny, nz,
c    $                   a, b, c, d, e,
c    $                   rsd )
c
c***forward elimination
c
            do i = 2, nx-1
c
               tmp1 = b(i,j,k) / c(i-1,j,k)
               tmp2 = a(i+1,j,k) / c(i-1,j,k)
c     
               c(i,j,k) = c(i,j,k) - tmp1 * d(i-1,j,k)
               d(i,j,k) = d(i,j,k) - tmp1 * e(i-1,j,k) 
               rsd(1,i,j,k) = rsd(1,i,j,k) - tmp1 * rsd(1,i-1,j,k)
               rsd(2,i,j,k) = rsd(2,i,j,k) - tmp1 * rsd(2,i-1,j,k)
               rsd(3,i,j,k) = rsd(3,i,j,k) - tmp1 * rsd(3,i-1,j,k)
c
               b(i+1,j,k) = b(i+1,j,k) - tmp2 * d(i-1,j,k)
               c(i+1,j,k) = c(i+1,j,k) - tmp2 * e(i-1,j,k)
               rsd(1,i+1,j,k) = rsd(1,i+1,j,k) - tmp2 * rsd(1,i-1,j,k)
               rsd(2,i+1,j,k) = rsd(2,i+1,j,k) - tmp2 * rsd(2,i-1,j,k)
               rsd(3,i+1,j,k) = rsd(3,i+1,j,k) - tmp2 * rsd(3,i-1,j,k)
c     
            end do
c     
            tmp1 = b(nx,j,k) / c(nx-1,j,k)
c
            c(nx,j,k) = c(nx,j,k) - tmp1 * d(nx-1,j,k)
            d(nx,j,k) = d(nx,j,k) - tmp1 * e(nx-1,j,k)
            rsd(1,nx,j,k) = rsd(1,nx,j,k) - tmp1 * rsd(1,nx-1,j,k)
            rsd(2,nx,j,k) = rsd(2,nx,j,k) - tmp1 * rsd(2,nx-1,j,k)
            rsd(3,nx,j,k) = rsd(3,nx,j,k) - tmp1 * rsd(3,nx-1,j,k)
c
c***back-substitution phase
c
            rsd(1,nx,j,k) = rsd(1,nx,j,k) / c(nx,j,k)
      rsd(1,nx-1,j,k) = ( rsd(1,nx-1,j,k) - d(nx-1,j,k)*rsd(1,nx,j,k) )
     $           / c(nx-1,j,k)
c
            rsd(2,nx,j,k) = rsd(2,nx,j,k) / c(nx,j,k)
      rsd(2,nx-1,j,k) = ( rsd(2,nx-1,j,k) - d(nx-1,j,k)*rsd(2,nx,j,k) )
     $           / c(nx-1,j,k)
c
            rsd(3,nx,j,k) = rsd(3,nx,j,k) / c(nx,j,k)
      rsd(3,nx-1,j,k) = ( rsd(3,nx-1,j,k) - d(nx-1,j,k)*rsd(3,nx,j,k) )
     $           / c(nx-1,j,k)
c
            do i = nx-2, 1, -1
c
               rsd(1,i,j,k) = ( rsd(1,i,j,k) - d(i,j,k)*rsd(1,i+1,j,k)
     $              - e(i,j,k)*rsd(1,i+2,j,k) ) / c(i,j,k)
c
               rsd(2,i,j,k) = ( rsd(2,i,j,k) - d(i,j,k)*rsd(2,i+1,j,k)
     $              - e(i,j,k)*rsd(2,i+2,j,k) ) / c(i,j,k)
c
               rsd(3,i,j,k) = ( rsd(3,i,j,k) - d(i,j,k)*rsd(3,i+1,j,k)
     $              - e(i,j,k)*rsd(3,i+2,j,k) ) / c(i,j,k)
c     
            end do
c
c
c        call jacx ( 4 )
c
            do i = 1, nx
c
               ru1 = 1.0d+00 / u(1,i,j,k)
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(i) = uu
               aa(i) = sqrt ( c1 * c2 * ( ru1 * u(5,i,j,k) - q ) )
c
               rhon(i) = max ( dx1,
     $                         dx2 + r43 * c34 * ru1,
     $                         dx3 + c34 * ru1,
     $                         dx4 + c34 * ru1,
     $                         dx5 + c1345 * ru1 )
c
            end do
c
            a(1,j,k) = 0.0d+00
            b(1,j,k) = 0.0d+00
            c(1,j,k) = 1.0d+00
            d(1,j,k) = 0.0d+00
            e(1,j,k) = 0.0d+00
c
            do i = 2, nx - 1
c
               a(i,j,k) =   0.0d+00
               b(i,j,k) = - dt * tx2 * ( cv(i-1) +  aa(i-1) )
     $                    - dt * rhon(i-1) * tx1 
               c(i,j,k) =   1.0d+00
     $                    + dt * rhon(i) * tx1 * 2.0d+00
               d(i,j,k) =   dt * tx2 * ( cv(i+1) +  aa(i+1) )
     $                    - dt * rhon(i+1) * tx1
               e(i,j,k) =   0.0d+00
c
            end do
c
            a(nx,j,k) = 0.0d+00
            b(nx,j,k) = 0.0d+00
            c(nx,j,k) = 1.0d+00
            d(nx,j,k) = 0.0d+00
            e(nx,j,k) = 0.0d+00
c
c***fourth order dissipation 
c
            c(2,j,k) = c(2,j,k) + dt * dssp * ( + 5.0d+00 )
            d(2,j,k) = d(2,j,k) + dt * dssp * ( - 4.0d+00 )
            e(2,j,k) = e(2,j,k) + dt * dssp * ( + 1.0d+00 )
c
            b(3,j,k) = b(3,j,k) + dt * dssp * ( - 4.0d+00 )
            c(3,j,k) = c(3,j,k) + dt * dssp * ( + 6.0d+00 )
            d(3,j,k) = d(3,j,k) + dt * dssp * ( - 4.0d+00 )
            e(3,j,k) = e(3,j,k) + dt * dssp * ( + 1.0d+00 )
c
            do i = 4, nx - 3
c
               a(i,j,k) = a(i,j,k) + dt * dssp * ( + 1.0d+00 )
               b(i,j,k) = b(i,j,k) + dt * dssp * ( - 4.0d+00 )
               c(i,j,k) = c(i,j,k) + dt * dssp * ( + 6.0d+00 )
               d(i,j,k) = d(i,j,k) + dt * dssp * ( - 4.0d+00 )
               e(i,j,k) = e(i,j,k) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(nx-2,j,k) = a(nx-2,j,k) + dt * dssp * ( + 1.0d+00 )
            b(nx-2,j,k) = b(nx-2,j,k) + dt * dssp * ( - 4.0d+00 )
            c(nx-2,j,k) = c(nx-2,j,k) + dt * dssp * ( + 6.0d+00 )
            d(nx-2,j,k) = d(nx-2,j,k) + dt * dssp * ( - 4.0d+00 )
c
            a(nx-1,j,k) = a(nx-1,j,k) + dt * dssp * ( + 1.0d+00 )
            b(nx-1,j,k) = b(nx-1,j,k) + dt * dssp * ( - 4.0d+00 )
            c(nx-1,j,k) = c(nx-1,j,k) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentax ( isiz1, isiz2, isiz3,
c    $                  4, nx, ny, nz,
c    $                  a, b, c, d, e,
c    $                  rsd )
c
c***forward elimination
c
            do i = 2, nx-1
c
               tmp1 = b(i,j,k) / c(i-1,j,k)
               tmp2 = a(i+1,j,k) / c(i-1,j,k)
c     
               c(i,j,k) = c(i,j,k) - tmp1 * d(i-1,j,k)
               d(i,j,k) = d(i,j,k) - tmp1 * e(i-1,j,k) 
               rsd(4,i,j,k) = rsd(4,i,j,k) - tmp1 * rsd(4,i-1,j,k)
c
               b(i+1,j,k) = b(i+1,j,k) - tmp2 * d(i-1,j,k)
               c(i+1,j,k) = c(i+1,j,k) - tmp2 * e(i-1,j,k)
               rsd(4,i+1,j,k) = rsd(4,i+1,j,k) - tmp2 * rsd(4,i-1,j,k)
c     
            end do
c     
            tmp1 = b(nx,j,k) / c(nx-1,j,k)
c
            c(nx,j,k) = c(nx,j,k) - tmp1 * d(nx-1,j,k)
            d(nx,j,k) = d(nx,j,k) - tmp1 * e(nx-1,j,k)
            rsd(4,nx,j,k) = rsd(4,nx,j,k) - tmp1 * rsd(4,nx-1,j,k)
c
c***back-substitution phase
c
            rsd(4,nx,j,k) = rsd(4,nx,j,k) / c(nx,j,k)
      rsd(4,nx-1,j,k) = ( rsd(4,nx-1,j,k) - d(nx-1,j,k)*rsd(4,nx,j,k) )
     $           / c(nx-1,j,k)
c
            do i = nx-2, 1, -1
c
      rsd(4,i,j,k) = ( rsd(4,i,j,k) - d(i,j,k)*rsd(4,i+1,j,k)
     $              - e(i,j,k)*rsd(4,i+2,j,k) ) / c(i,j,k)
c     
            end do
c
c
c        call jacx ( 5 )
c
            do i = 1, nx
c
               ru1 = 1.0d+00 / u(1,i,j,k)
               uu = ru1 * u(2,i,j,k)
               vv = ru1 * u(3,i,j,k)
               ww = ru1 * u(4,i,j,k)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(i) = uu
               aa(i) = sqrt ( c1 * c2 * ( ru1 * u(5,i,j,k) - q ) )
c
               rhon(i) = max ( dx1,
     $                         dx2 + r43 * c34 * ru1,
     $                         dx3 + c34 * ru1,
     $                         dx4 + c34 * ru1,
     $                         dx5 + c1345 * ru1 )
c
            end do
c
            a(1,j,k) = 0.0d+00
            b(1,j,k) = 0.0d+00
            c(1,j,k) = 1.0d+00
            d(1,j,k) = 0.0d+00
            e(1,j,k) = 0.0d+00
c
            do i = 2, nx - 1
c
               a(i,j,k) =   0.0d+00
               b(i,j,k) = - dt * tx2 * ( cv(i-1) -  aa(i-1) )
     $                    - dt * rhon(i-1) * tx1 
               c(i,j,k) =   1.0d+00
     $                    + dt * rhon(i) * tx1 * 2.0d+00
               d(i,j,k) =   dt * tx2 * ( cv(i+1) - aa(i+1) )
     $                    - dt * rhon(i+1) * tx1
               e(i,j,k) =   0.0d+00
c
            end do
c
            a(nx,j,k) = 0.0d+00
            b(nx,j,k) = 0.0d+00
            c(nx,j,k) = 1.0d+00
            d(nx,j,k) = 0.0d+00
            e(nx,j,k) = 0.0d+00
c
c***fourth order dissipation 
c
            c(2,j,k) = c(2,j,k) + dt * dssp * ( + 5.0d+00 )
            d(2,j,k) = d(2,j,k) + dt * dssp * ( - 4.0d+00 )
            e(2,j,k) = e(2,j,k) + dt * dssp * ( + 1.0d+00 )
c
            b(3,j,k) = b(3,j,k) + dt * dssp * ( - 4.0d+00 )
            c(3,j,k) = c(3,j,k) + dt * dssp * ( + 6.0d+00 )
            d(3,j,k) = d(3,j,k) + dt * dssp * ( - 4.0d+00 )
            e(3,j,k) = e(3,j,k) + dt * dssp * ( + 1.0d+00 )
c
            do i = 4, nx - 3
c
               a(i,j,k) = a(i,j,k) + dt * dssp * ( + 1.0d+00 )
               b(i,j,k) = b(i,j,k) + dt * dssp * ( - 4.0d+00 )
               c(i,j,k) = c(i,j,k) + dt * dssp * ( + 6.0d+00 )
               d(i,j,k) = d(i,j,k) + dt * dssp * ( - 4.0d+00 )
               e(i,j,k) = e(i,j,k) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(nx-2,j,k) = a(nx-2,j,k) + dt * dssp * ( + 1.0d+00 )
            b(nx-2,j,k) = b(nx-2,j,k) + dt * dssp * ( - 4.0d+00 )
            c(nx-2,j,k) = c(nx-2,j,k) + dt * dssp * ( + 6.0d+00 )
            d(nx-2,j,k) = d(nx-2,j,k) + dt * dssp * ( - 4.0d+00 )
c
            a(nx-1,j,k) = a(nx-1,j,k) + dt * dssp * ( + 1.0d+00 )
            b(nx-1,j,k) = b(nx-1,j,k) + dt * dssp * ( - 4.0d+00 )
            c(nx-1,j,k) = c(nx-1,j,k) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentax ( isiz1, isiz2, isiz3,
c    $                  5, nx, ny, nz,
c    $                  a, b, c, d, e,
c    $                  rsd )
c
c***forward elimination
c
            do i = 2, nx-1
c
               tmp1 = b(i,j,k) / c(i-1,j,k)
               tmp2 = a(i+1,j,k) / c(i-1,j,k)
c     
               c(i,j,k) = c(i,j,k) - tmp1 * d(i-1,j,k)
               d(i,j,k) = d(i,j,k) - tmp1 * e(i-1,j,k) 
               rsd(5,i,j,k) = rsd(5,i,j,k) - tmp1 * rsd(5,i-1,j,k)
c
               b(i+1,j,k) = b(i+1,j,k) - tmp2 * d(i-1,j,k)
               c(i+1,j,k) = c(i+1,j,k) - tmp2 * e(i-1,j,k)
               rsd(5,i+1,j,k) = rsd(5,i+1,j,k) - tmp2 * rsd(5,i-1,j,k)
c     
            end do
c     
            tmp1 = b(nx,j,k) / c(nx-1,j,k)
c
            c(nx,j,k) = c(nx,j,k) - tmp1 * d(nx-1,j,k)
            d(nx,j,k) = d(nx,j,k) - tmp1 * e(nx-1,j,k)
            rsd(5,nx,j,k) = rsd(5,nx,j,k) - tmp1 * rsd(5,nx-1,j,k)
c
c***back-substitution phase
c
            rsd(5,nx,j,k) = rsd(5,nx,j,k) / c(nx,j,k)
      rsd(5,nx-1,j,k) = ( rsd(5,nx-1,j,k) - d(nx-1,j,k)*rsd(5,nx,j,k) )
     $           / c(nx-1,j,k)
c
            do i = nx-2, 1, -1
c
           rsd(5,i,j,k) = ( rsd(5,i,j,k) - d(i,j,k)*rsd(5,i+1,j,k)
     $              - e(i,j,k)*rsd(5,i+2,j,k) ) / c(i,j,k)
c     
            end do
c
         end do
		 end do
		 return
		 end
c
c
      subroutine ztasweep
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
      implicit real*8 (a-h,o-z)
      include 'appspB.incl'
c
      parameter ( c1 = 1.40d+00, c2 = 0.40d+00,
     $            c3 = 1.00d-01, c4 = 1.00d+00,
     $            c5 = 1.40d+00 )
      real *8 timer, tstart, tend
c
c***grid
c
      common/cgcon/ nx, ny, nz,
     $              ii1, ii2, ji1, ji2, ki1, ki2
	  common/cgcon2/
     $              dxi, deta, dzeta,
     $              tx1, tx2, tx3, 
     $              ty1, ty2, ty3, 
     $              tz1, tz2, tz3
c
c***dissipation
c
      common/disp/ dx1,dx2,dx3,dx4,dx5,
     $             dy1,dy2,dy3,dy4,dy5,
     $             dz1,dz2,dz3,dz4,dz5,
     $             dssp
c
c***field variables and residuals
c
c
c***output control parameters
c
      common/cprcon/ ipr, iout, inorm
c
c***newton-raphson iteration control parameters
c
      common/ctscon/ itmax, invert,
     $               dt, omega, tolrsd(5),
     $               rsdnm(5), errnm(5), frc, ttotal
c
c***global jacobian matrix
c***coefficients of the exact solution
c
      common/cexact/ ce(5,13)
c





c
      dimension idmax(5), jdmax(5), kdmax(5),
     $          imax(5), jmax(5), kmax(5),
     $          delunm(5)
c
      dimension flux1(isiz1)
      dimension flux2(isiz1)
      dimension flux3(isiz1)
      dimension flux4(isiz1)
      dimension flux5(isiz1)
c
      dimension cv(isiz3), aa(isiz3), rhos(isiz3)
      dimension  rhon(isiz1)
      dimension  rhoq(isiz2)
c
      parameter ( one = 1.0d+00 )
c
      lnorm = 2
      bt = sqrt ( 0.50d+00 )
      r43 = 4.0d+00 / 3.0d+00
      c34 = c3 * c4
      c1345 = c1 * c3 * c4 * c5
c***perform the zeta-direction sweep
c
c        call jacz ( 3 )
c
c
c***perform the zeta-direction sweep
c
c        call jacz ( 3 )
c
CHPF$ INDEPENDENT
      do k = 1, nz
c
         do j = 1, ny
c
            do i = 1, nx
               rsdt(1,k,i,j) = rsd(1,i,j,k)
               rsdt(2,k,i,j) = rsd(2,i,j,k)
               rsdt(3,k,i,j) = rsd(3,i,j,k)
               rsdt(4,k,i,j) = rsd(4,i,j,k)
               rsdt(5,k,i,j) = rsd(5,i,j,k)
               ut(1,k,i,j) = u(1,i,j,k)
               ut(2,k,i,j) = u(2,i,j,k)
               ut(3,k,i,j) = u(3,i,j,k)
               ut(4,k,i,j) = u(4,i,j,k)
               ut(5,k,i,j) = u(5,i,j,k)
             end do
c
         end do
c
      end do
c
CHPF$ INDEPENDENT
      do j = 2, ny - 1
c
         do i = 2, nx - 1
c
            do k = 1, nz
c
               ru1 = 1.0d+00 / ut(1,k,i,j)
               uu = ru1 * ut(2,k,i,j)
               vv = ru1 * ut(3,k,i,j)
               ww = ru1 * ut(4,k,i,j)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(k) = ww
               aa(k) = sqrt ( c1 * c2 * ( ru1 * ut(5,k,i,j) - q ) )
c
               rhos(k) = max ( dz1,
     $                         dz2 + c34 * ru1 ,
     $                         dz3 + c34 * ru1 ,
     $                         dz4 + r43 * c34 * ru1,
     $                         dz5 + c1345 * ru1 )
c
            end do
c
            a(1,i,j) = 0.0d+00
            b(1,i,j) = 0.0d+00
            c(1,i,j) = 1.0d+00
            d(1,i,j) = 0.0d+00
            e(1,i,j) = 0.0d+00
c
            do k = 2, nz - 1
c
               a(k,i,j) = 0.0d+00
               b(k,i,j) = - dt * tz2 *  cv(k-1)
     $                    - dt * rhos(k-1) * tz1 
               c(k,i,j) =   1.0d+00
     $                    + dt * rhos(k) * tz1 * 2.0d+00
               d(k,i,j) =   dt * tz2 *  cv(k+1)
     $                    - dt * rhos(k+1) * tz1
               e(k,i,j) = 0.0d+00
c
            end do
c
            a(nz,i,j) = 0.0d+00
            b(nz,i,j) = 0.0d+00
            c(nz,i,j) = 1.0d+00
            d(nz,i,j) = 0.0d+00
            e(nz,i,j) = 0.0d+00
c
c***fourth order dissipation 
c
            c(2,i,j) = c(2,i,j) + dt * dssp * ( + 5.0d+00 )
            d(2,i,j) = d(2,i,j) + dt * dssp * ( - 4.0d+00 )
            e(2,i,j) = e(2,i,j) + dt * dssp * ( + 1.0d+00 )
c
            b(3,i,j) = b(3,i,j) + dt * dssp * ( - 4.0d+00 )
            c(3,i,j) = c(3,i,j) + dt * dssp * ( + 6.0d+00 )
            d(3,i,j) = d(3,i,j) + dt * dssp * ( - 4.0d+00 )
            e(3,i,j) = e(3,i,j) + dt * dssp * ( + 1.0d+00 )
c
            do k = 4, nz - 3
c
               a(k,i,j) = a(k,i,j) + dt * dssp * ( + 1.0d+00 )
               b(k,i,j) = b(k,i,j) + dt * dssp * ( - 4.0d+00 )
               c(k,i,j) = c(k,i,j) + dt * dssp * ( + 6.0d+00 )
               d(k,i,j) = d(k,i,j) + dt * dssp * ( - 4.0d+00 )
               e(k,i,j) = e(k,i,j) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(nz-2,i,j) = a(nz-2,i,j) + dt * dssp * ( + 1.0d+00 )
            b(nz-2,i,j) = b(nz-2,i,j) + dt * dssp * ( - 4.0d+00 )
            c(nz-2,i,j) = c(nz-2,i,j) + dt * dssp * ( + 6.0d+00 )
            d(nz-2,i,j) = d(nz-2,i,j) + dt * dssp * ( - 4.0d+00 )
c
            a(nz-1,i,j) = a(nz-1,i,j) + dt * dssp * ( + 1.0d+00 )
            b(nz-1,i,j) = b(nz-1,i,j) + dt * dssp * ( - 4.0d+00 )
            c(nz-1,i,j) = c(nz-1,i,j) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentaz3 ( isiz1, isiz2, isiz3,
c    $                   nx, ny, nz,
c    $                   a, b, c, d, e,
c    $                   rsd )
c
c***forward elimination
c
            do k = 2, nz - 1
c
               tmp1 = b(k,i,j) / c(k-1,i,j)
               tmp2 = a(k+1,i,j) / c(k-1,i,j)
c     
               c(k,i,j) = c(k,i,j) - tmp1 * d(k-1,i,j)
               d(k,i,j) = d(k,i,j) - tmp1 * e(k-1,i,j) 
               rsdt(1,k,i,j) = rsdt(1,k,i,j) - tmp1 * rsdt(1,k-1,i,j)
               rsdt(2,k,i,j) = rsdt(2,k,i,j) - tmp1 * rsdt(2,k-1,i,j)
               rsdt(3,k,i,j) = rsdt(3,k,i,j) - tmp1 * rsdt(3,k-1,i,j)
c
               b(k+1,i,j) = b(k+1,i,j) - tmp2 * d(k-1,i,j)
               c(k+1,i,j) = c(k+1,i,j) - tmp2 * e(k-1,i,j)
       rsdt(1,k+1,i,j) = rsdt(1,k+1,i,j) - tmp2 * rsdt(1,k-1,i,j)
      rsdt(2,k+1,i,j) = rsdt(2,k+1,i,j) - tmp2 * rsdt(2,k-1,i,j)
       rsdt(3,k+1,i,j) = rsdt(3,k+1,i,j) - tmp2 * rsdt(3,k-1,i,j)
c     
            end do
c     
            tmp1 = b(nz,i,j) / c(nz-1,i,j)
c
            c(nz,i,j) = c(nz,i,j) - tmp1 * d(nz-1,i,j)
            d(nz,i,j) = d(nz,i,j) - tmp1 * e(nz-1,i,j)
            rsdt(1,nz,i,j) = rsdt(1,nz,i,j) - tmp1 * rsdt(1,nz-1,i,j)
            rsdt(2,nz,i,j) = rsdt(2,nz,i,j) - tmp1 * rsdt(2,nz-1,i,j)
            rsdt(3,nz,i,j) = rsdt(3,nz,i,j) - tmp1 * rsdt(3,nz-1,i,j)
c
c***back-substitution phase
c
            rsdt(1,nz,i,j) = rsdt(1,nz,i,j) / c(nz,i,j)
       rsdt(1,nz-1,i,j)=(rsdt(1,nz-1,i,j)-d(nz-1,i,j)*rsdt(1,nz,i,j) )
     $           / c(nz-1,i,j)
c
            rsdt(2,nz,i,j) = rsdt(2,nz,i,j) / c(nz,i,j)
       rsdt(2,nz-1,i,j)=(rsdt(2,nz-1,i,j)-d(nz-1,i,j)*rsdt(2,nz,i,j))
     $           / c(nz-1,i,j)
c
            rsdt(3,nz,i,j) = rsdt(3,nz,i,j) / c(nz,i,j)
       rsdt(3,nz-1,i,j)=(rsdt(3,nz-1,i,j)-d(nz-1,i,j)*rsdt(3,nz,i,j) )
     $           / c(nz-1,i,j)
c
            do k = nz - 2, 1, -1
c
        rsdt(1,k,i,j) = ( rsdt(1,k,i,j) - d(k,i,j)*rsdt(1,k+1,i,j)
     $              - e(k,i,j)*rsdt(1,k+2,i,j) ) / c(k,i,j)
c
      rsdt(2,k,i,j) = ( rsdt(2,k,i,j) - d(k,i,j)*rsdt(2,k+1,i,j)
     $              - e(k,i,j)*rsdt(2,k+2,i,j) ) / c(k,i,j)
c
       rsdt(3,k,i,j) = ( rsdt(3,k,i,j) - d(k,i,j)*rsdt(3,k+1,i,j)
     $              - e(k,i,j)*rsdt(3,k+2,i,j) ) / c(k,i,j)
c     
            end do
c
c
c        call jacz ( 4 )
c
            do k = 1, nz
c
               ru1 = 1.0d+00 / ut(1,k,i,j)
               uu = ru1 * ut(2,k,i,j)
               vv = ru1 * ut(3,k,i,j)
               ww = ru1 * ut(4,k,i,j)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(k) = ww
               aa(k) = sqrt ( c1 * c2 * ( ru1 * ut(5,k,i,j) - q ) )
c
               rhos(k) = max ( dz1,
     $                         dz2 + c34 * ru1 ,
     $                         dz3 + c34 * ru1 ,
     $                         dz4 + r43 * c34 * ru1,
     $                         dz5 + c1345 * ru1 )
c
            end do
c
            a(1,i,j) = 0.0d+00
            b(1,i,j) = 0.0d+00
            c(1,i,j) = 1.0d+00
            d(1,i,j) = 0.0d+00
            e(1,i,j) = 0.0d+00
c
            do k = 2, nz - 1
c
               a(k,i,j) = 0.0d+00
               b(k,i,j) = - dt * tz2 * ( cv(k-1) + aa(k-1) )
     $                    - dt * rhos(k-1) * tz1 
               c(k,i,j) =   1.0d+00
     $                    + dt * rhos(k) * tz1 * 2.0d+00
               d(k,i,j) =   dt * tz2 * ( cv(k+1) + aa(k+1) )
     $                    - dt * rhos(k+1) * tz1
               e(k,i,j) = 0.0d+00
c
            end do
c
            a(nz,i,j) = 0.0d+00
            b(nz,i,j) = 0.0d+00
            c(nz,i,j) = 1.0d+00
            d(nz,i,j) = 0.0d+00
            e(nz,i,j) = 0.0d+00
c
c***fourth order dissipation 
c
            c(2,i,j) = c(2,i,j) + dt * dssp * ( + 5.0d+00 )
            d(2,i,j) = d(2,i,j) + dt * dssp * ( - 4.0d+00 )
            e(2,i,j) = e(2,i,j) + dt * dssp * ( + 1.0d+00 )
c
            b(3,i,j) = b(3,i,j) + dt * dssp * ( - 4.0d+00 )
            c(3,i,j) = c(3,i,j) + dt * dssp * ( + 6.0d+00 )
            d(3,i,j) = d(3,i,j) + dt * dssp * ( - 4.0d+00 )
            e(3,i,j) = e(3,i,j) + dt * dssp * ( + 1.0d+00 )
c
            do k = 4, nz - 3
c
               a(k,i,j) = a(k,i,j) + dt * dssp * ( + 1.0d+00 )
               b(k,i,j) = b(k,i,j) + dt * dssp * ( - 4.0d+00 )
               c(k,i,j) = c(k,i,j) + dt * dssp * ( + 6.0d+00 )
               d(k,i,j) = d(k,i,j) + dt * dssp * ( - 4.0d+00 )
               e(k,i,j) = e(k,i,j) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(nz-2,i,j) = a(nz-2,i,j) + dt * dssp * ( + 1.0d+00 )
            b(nz-2,i,j) = b(nz-2,i,j) + dt * dssp * ( - 4.0d+00 )
            c(nz-2,i,j) = c(nz-2,i,j) + dt * dssp * ( + 6.0d+00 )
            d(nz-2,i,j) = d(nz-2,i,j) + dt * dssp * ( - 4.0d+00 )
c
            a(nz-1,i,j) = a(nz-1,i,j) + dt * dssp * ( + 1.0d+00 )
            b(nz-1,i,j) = b(nz-1,i,j) + dt * dssp * ( - 4.0d+00 )
            c(nz-1,i,j) = c(nz-1,i,j) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentaz ( isiz1, isiz2, isiz3,
c    $                  4, nx, ny, nz,
c    $                  a, b, c, d, e,
c    $                  rsd )
c
c***forward elimination
c
            do k = 2, nz - 1
c
               tmp1 = b(k,i,j) / c(k-1,i,j)
               tmp2 = a(k+1,i,j) / c(k-1,i,j)
c     
               c(k,i,j) = c(k,i,j) - tmp1 * d(k-1,i,j)
               d(k,i,j) = d(k,i,j) - tmp1 * e(k-1,i,j) 
               rsdt(4,k,i,j) = rsdt(4,k,i,j) - tmp1 * rsdt(4,k-1,i,j)
c
               b(k+1,i,j) = b(k+1,i,j) - tmp2 * d(k-1,i,j)
               c(k+1,i,j) = c(k+1,i,j) - tmp2 * e(k-1,i,j)
      rsdt(4,k+1,i,j) = rsdt(4,k+1,i,j) - tmp2 * rsdt(4,k-1,i,j)
c     
            end do
c     
            tmp1 = b(nz,i,j) / c(nz-1,i,j)
c
            c(nz,i,j) = c(nz,i,j) - tmp1 * d(nz-1,i,j)
            d(nz,i,j) = d(nz,i,j) - tmp1 * e(nz-1,i,j)
            rsdt(4,nz,i,j) = rsdt(4,nz,i,j) - tmp1 * rsdt(4,nz-1,i,j)
c
c***back-substitution phase
c
            rsdt(4,nz,i,j) = rsdt(4,nz,i,j) / c(nz,i,j)
      rsdt(4,nz-1,i,j)=(rsdt(4,nz-1,i,j)-d(nz-1,i,j)*rsdt(4,nz,i,j) )
     $           / c(nz-1,i,j)
c
            do k = nz - 2, 1, -1
c
      rsdt(4,k,i,j) = ( rsdt(4,k,i,j) - d(k,i,j)*rsdt(4,k+1,i,j)
     $              - e(k,i,j)*rsdt(4,k+2,i,j) ) / c(k,i,j)
c     
            end do
c
c
c        call jacz ( 5 )
c
            do k = 1, nz
c
               ru1 = 1.0d+00 / ut(1,k,i,j)
               uu = ru1 * ut(2,k,i,j)
               vv = ru1 * ut(3,k,i,j)
               ww = ru1 * ut(4,k,i,j)
c
               q = 0.50d+00 * (  uu ** 2
     $                         + vv ** 2
     $                         + ww ** 2 )
c
               cv(k) = ww
               aa(k) = sqrt ( c1 * c2 * ( ru1 * ut(5,k,i,j) - q ) )
c
               rhos(k) = max ( dz1,
     $                         dz2 + c34 * ru1 ,
     $                         dz3 + c34 * ru1 ,
     $                         dz4 + r43 * c34 * ru1,
     $                         dz5 + c1345 * ru1 )
c
            end do
c
            a(1,i,j) = 0.0d+00
            b(1,i,j) = 0.0d+00
            c(1,i,j) = 1.0d+00
            d(1,i,j) = 0.0d+00
            e(1,i,j) = 0.0d+00
c
            do k = 2, nz - 1
c
               a(k,i,j) = 0.0d+00
               b(k,i,j) = - dt * tz2 * ( cv(k-1) - aa(k-1) )
     $                    - dt * rhos(k-1) * tz1 
               c(k,i,j) =   1.0d+00
     $                    + dt * rhos(k) * tz1 * 2.0d+00
               d(k,i,j) =   dt * tz2 * ( cv(k+1) - aa(k+1) )
     $                    - dt * rhos(k+1) * tz1
               e(k,i,j) = 0.0d+00
c
            end do
c
            a(nz,i,j) = 0.0d+00
            b(nz,i,j) = 0.0d+00
            c(nz,i,j) = 1.0d+00
            d(nz,i,j) = 0.0d+00
            e(nz,i,j) = 0.0d+00
c
c***fourth order dissipation 
c
            c(2,i,j) = c(2,i,j) + dt * dssp * ( + 5.0d+00 )
            d(2,i,j) = d(2,i,j) + dt * dssp * ( - 4.0d+00 )
            e(2,i,j) = e(2,i,j) + dt * dssp * ( + 1.0d+00 )
c
            b(3,i,j) = b(3,i,j) + dt * dssp * ( - 4.0d+00 )
            c(3,i,j) = c(3,i,j) + dt * dssp * ( + 6.0d+00 )
            d(3,i,j) = d(3,i,j) + dt * dssp * ( - 4.0d+00 )
            e(3,i,j) = e(3,i,j) + dt * dssp * ( + 1.0d+00 )
c
            do k = 4, nz - 3
c
               a(k,i,j) = a(k,i,j) + dt * dssp * ( + 1.0d+00 )
               b(k,i,j) = b(k,i,j) + dt * dssp * ( - 4.0d+00 )
               c(k,i,j) = c(k,i,j) + dt * dssp * ( + 6.0d+00 )
               d(k,i,j) = d(k,i,j) + dt * dssp * ( - 4.0d+00 )
               e(k,i,j) = e(k,i,j) + dt * dssp * ( + 1.0d+00 )
c
            end do
c
            a(nz-2,i,j) = a(nz-2,i,j) + dt * dssp * ( + 1.0d+00 )
            b(nz-2,i,j) = b(nz-2,i,j) + dt * dssp * ( - 4.0d+00 )
            c(nz-2,i,j) = c(nz-2,i,j) + dt * dssp * ( + 6.0d+00 )
            d(nz-2,i,j) = d(nz-2,i,j) + dt * dssp * ( - 4.0d+00 )
c
            a(nz-1,i,j) = a(nz-1,i,j) + dt * dssp * ( + 1.0d+00 )
            b(nz-1,i,j) = b(nz-1,i,j) + dt * dssp * ( - 4.0d+00 )
            c(nz-1,i,j) = c(nz-1,i,j) + dt * dssp * ( + 5.0d+00 )
c
c
c        call spentaz ( isiz1, isiz2, isiz3,
c    $                  5, nx, ny, nz,
c    $                  a, b, c, d, e,
c    $                  rsd )
c
c***forward elimination
c
            do k = 2, nz - 1
c
               tmp1 = b(k,i,j) / c(k-1,i,j)
               tmp2 = a(k+1,i,j) / c(k-1,i,j)
c     
               c(k,i,j) = c(k,i,j) - tmp1 * d(k-1,i,j)
               d(k,i,j) = d(k,i,j) - tmp1 * e(k-1,i,j) 
               rsdt(5,k,i,j) = rsdt(5,k,i,j) - tmp1 * rsdt(5,k-1,i,j)
c
               b(k+1,i,j) = b(k+1,i,j) - tmp2 * d(k-1,i,j)
               c(k+1,i,j) = c(k+1,i,j) - tmp2 * e(k-1,i,j)
       rsdt(5,k+1,i,j) = rsdt(5,k+1,i,j) - tmp2 * rsdt(5,k-1,i,j)
c     
            end do
c     
            tmp1 = b(nz,i,j) / c(nz-1,i,j)
c
            c(nz,i,j) = c(nz,i,j) - tmp1 * d(nz-1,i,j)
            d(nz,i,j) = d(nz,i,j) - tmp1 * e(nz-1,i,j)
            rsdt(5,nz,i,j) = rsdt(5,nz,i,j) - tmp1 * rsdt(5,nz-1,i,j)
c
c***back-substitution phase
c
            rsdt(5,nz,i,j) = rsdt(5,nz,i,j) / c(nz,i,j)
      rsdt(5,nz-1,i,j)=(rsdt(5,nz-1,i,j)-d(nz-1,i,j)*rsdt(5,nz,i,j) )
     $           / c(nz-1,i,j)
c
            do k = nz - 2, 1, -1
c
        rsdt(5,k,i,j) = ( rsdt(5,k,i,j) - d(k,i,j)*rsdt(5,k+1,i,j)
     $              - e(k,i,j)*rsdt(5,k+2,i,j) ) / c(k,i,j)
c     
            end do
c
         end do
c
      end do
c
c
CHPF$ INDEPENDENT
      do k = 1, nz
c
         do j = 1, ny
c
            do i = 1, nx
               rsd(1,i,j,k) = rsdt(1,k,i,j)
               rsd(2,i,j,k) = rsdt(2,k,i,j)
               rsd(3,i,j,k) = rsdt(3,k,i,j)
               rsd(4,i,j,k) = rsdt(4,k,i,j)
               rsd(5,i,j,k) = rsdt(5,k,i,j)
c              u(1,i,j,k) = ut(1,k,i,j)
c              u(2,i,j,k) = ut(2,k,i,j)
c              u(3,i,j,k) = ut(3,k,i,j)
c              u(4,i,j,k) = ut(4,k,i,j)
c              u(5,i,j,k) = ut(5,k,i,j)
             end do
c
         end do
c
      end do
	   return
	   end
