--- /dev/null
+
+c---------------------------------------------------------------------
+c---------------------------------------------------------------------
+
+ subroutine jacu(k)
+
+c---------------------------------------------------------------------
+c---------------------------------------------------------------------
+
+c---------------------------------------------------------------------
+c compute the upper triangular part of the jacobian matrix
+c---------------------------------------------------------------------
+
+
+ implicit none
+
+ include 'applu.incl'
+
+c---------------------------------------------------------------------
+c input parameters
+c---------------------------------------------------------------------
+ integer k
+
+c---------------------------------------------------------------------
+c local variables
+c---------------------------------------------------------------------
+ integer i, j
+ double precision r43
+ double precision c1345
+ double precision c34
+ double precision tmp1, tmp2, tmp3
+
+
+
+ r43 = ( 4.0d+00 / 3.0d+00 )
+ c1345 = c1 * c3 * c4 * c5
+ c34 = c3 * c4
+
+ do j = jst, jend
+ do i = ist, iend
+
+c---------------------------------------------------------------------
+c form the block daigonal
+c---------------------------------------------------------------------
+ tmp1 = 1.0d+00 / u(1,i,j,k)
+ tmp2 = tmp1 * tmp1
+ tmp3 = tmp1 * tmp2
+
+ d(1,1,i,j) = 1.0d+00
+ > + dt * 2.0d+00 * ( tx1 * dx1
+ > + ty1 * dy1
+ > + tz1 * dz1 )
+ d(1,2,i,j) = 0.0d+00
+ d(1,3,i,j) = 0.0d+00
+ d(1,4,i,j) = 0.0d+00
+ d(1,5,i,j) = 0.0d+00
+
+ d(2,1,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( - r43 * c34 * tmp2 * u(2,i,j,k) )
+ > + ty1 * ( - c34 * tmp2 * u(2,i,j,k) )
+ > + tz1 * ( - c34 * tmp2 * u(2,i,j,k) ) )
+ d(2,2,i,j) = 1.0d+00
+ > + dt * 2.0d+00
+ > * ( tx1 * r43 * c34 * tmp1
+ > + ty1 * c34 * tmp1
+ > + tz1 * c34 * tmp1 )
+ > + dt * 2.0d+00 * ( tx1 * dx2
+ > + ty1 * dy2
+ > + tz1 * dz2 )
+ d(2,3,i,j) = 0.0d+00
+ d(2,4,i,j) = 0.0d+00
+ d(2,5,i,j) = 0.0d+00
+
+ d(3,1,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( - c34 * tmp2 * u(3,i,j,k) )
+ > + ty1 * ( - r43 * c34 * tmp2 * u(3,i,j,k) )
+ > + tz1 * ( - c34 * tmp2 * u(3,i,j,k) ) )
+ d(3,2,i,j) = 0.0d+00
+ d(3,3,i,j) = 1.0d+00
+ > + dt * 2.0d+00
+ > * ( tx1 * c34 * tmp1
+ > + ty1 * r43 * c34 * tmp1
+ > + tz1 * c34 * tmp1 )
+ > + dt * 2.0d+00 * ( tx1 * dx3
+ > + ty1 * dy3
+ > + tz1 * dz3 )
+ d(3,4,i,j) = 0.0d+00
+ d(3,5,i,j) = 0.0d+00
+
+ d(4,1,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( - c34 * tmp2 * u(4,i,j,k) )
+ > + ty1 * ( - c34 * tmp2 * u(4,i,j,k) )
+ > + tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k) ) )
+ d(4,2,i,j) = 0.0d+00
+ d(4,3,i,j) = 0.0d+00
+ d(4,4,i,j) = 1.0d+00
+ > + dt * 2.0d+00
+ > * ( tx1 * c34 * tmp1
+ > + ty1 * c34 * tmp1
+ > + tz1 * r43 * c34 * tmp1 )
+ > + dt * 2.0d+00 * ( tx1 * dx4
+ > + ty1 * dy4
+ > + tz1 * dz4 )
+ d(4,5,i,j) = 0.0d+00
+
+ d(5,1,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
+ > - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
+ > - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
+ > - ( c1345 ) * tmp2 * u(5,i,j,k) )
+ > + ty1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
+ > - ( r43*c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
+ > - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
+ > - ( c1345 ) * tmp2 * u(5,i,j,k) )
+ > + tz1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
+ > - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
+ > - ( r43*c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
+ > - ( c1345 ) * tmp2 * u(5,i,j,k) ) )
+ d(5,2,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( r43*c34 - c1345 ) * tmp2 * u(2,i,j,k)
+ > + ty1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k)
+ > + tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k) )
+ d(5,3,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k)
+ > + ty1 * ( r43*c34 -c1345 ) * tmp2 * u(3,i,j,k)
+ > + tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k) )
+ d(5,4,i,j) = dt * 2.0d+00
+ > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
+ > + ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
+ > + tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k) )
+ d(5,5,i,j) = 1.0d+00
+ > + dt * 2.0d+00 * ( tx1 * c1345 * tmp1
+ > + ty1 * c1345 * tmp1
+ > + tz1 * c1345 * tmp1 )
+ > + dt * 2.0d+00 * ( tx1 * dx5
+ > + ty1 * dy5
+ > + tz1 * dz5 )
+
+c---------------------------------------------------------------------
+c form the first block sub-diagonal
+c---------------------------------------------------------------------
+ tmp1 = 1.0d+00 / u(1,i+1,j,k)
+ tmp2 = tmp1 * tmp1
+ tmp3 = tmp1 * tmp2
+
+ a(1,1,i,j) = - dt * tx1 * dx1
+ a(1,2,i,j) = dt * tx2
+ a(1,3,i,j) = 0.0d+00
+ a(1,4,i,j) = 0.0d+00
+ a(1,5,i,j) = 0.0d+00
+
+ a(2,1,i,j) = dt * tx2
+ > * ( - ( u(2,i+1,j,k) * tmp1 ) ** 2
+ > + c2 * 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) ) * tmp2 )
+ > - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i+1,j,k) )
+ a(2,2,i,j) = dt * tx2
+ > * ( ( 2.0d+00 - c2 ) * ( u(2,i+1,j,k) * tmp1 ) )
+ > - dt * tx1 * ( r43 * c34 * tmp1 )
+ > - dt * tx1 * dx2
+ a(2,3,i,j) = dt * tx2
+ > * ( - c2 * ( u(3,i+1,j,k) * tmp1 ) )
+ a(2,4,i,j) = dt * tx2
+ > * ( - c2 * ( u(4,i+1,j,k) * tmp1 ) )
+ a(2,5,i,j) = dt * tx2 * c2
+
+ a(3,1,i,j) = dt * tx2
+ > * ( - ( u(2,i+1,j,k) * u(3,i+1,j,k) ) * tmp2 )
+ > - dt * tx1 * ( - c34 * tmp2 * u(3,i+1,j,k) )
+ a(3,2,i,j) = dt * tx2 * ( u(3,i+1,j,k) * tmp1 )
+ a(3,3,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
+ > - dt * tx1 * ( c34 * tmp1 )
+ > - dt * tx1 * dx3
+ a(3,4,i,j) = 0.0d+00
+ a(3,5,i,j) = 0.0d+00
+
+ a(4,1,i,j) = dt * tx2
+ > * ( - ( u(2,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 )
+ > - dt * tx1 * ( - c34 * tmp2 * u(4,i+1,j,k) )
+ a(4,2,i,j) = dt * tx2 * ( u(4,i+1,j,k) * tmp1 )
+ a(4,3,i,j) = 0.0d+00
+ a(4,4,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
+ > - dt * tx1 * ( c34 * tmp1 )
+ > - dt * tx1 * dx4
+ a(4,5,i,j) = 0.0d+00
+
+ a(5,1,i,j) = dt * tx2
+ > * ( ( c2 * ( 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) ) * tmp2
+ > - c1 * ( u(5,i+1,j,k) * tmp1 ) )
+ > * ( u(2,i+1,j,k) * tmp1 ) )
+ > - dt * tx1
+ > * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i+1,j,k)**2 )
+ > - ( c34 - c1345 ) * tmp3 * ( u(3,i+1,j,k)**2 )
+ > - ( c34 - c1345 ) * tmp3 * ( u(4,i+1,j,k)**2 )
+ > - c1345 * tmp2 * u(5,i+1,j,k) )
+ a(5,2,i,j) = dt * tx2
+ > * ( c1 * ( u(5,i+1,j,k) * tmp1 )
+ > - 0.50d+00 * c2
+ > * ( ( 3.0d+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) ) * tmp2 ) )
+ > - dt * tx1
+ > * ( r43*c34 - c1345 ) * tmp2 * u(2,i+1,j,k)
+ a(5,3,i,j) = dt * tx2
+ > * ( - c2 * ( u(3,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
+ > - dt * tx1
+ > * ( c34 - c1345 ) * tmp2 * u(3,i+1,j,k)
+ a(5,4,i,j) = dt * tx2
+ > * ( - c2 * ( u(4,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
+ > - dt * tx1
+ > * ( c34 - c1345 ) * tmp2 * u(4,i+1,j,k)
+ a(5,5,i,j) = dt * tx2
+ > * ( c1 * ( u(2,i+1,j,k) * tmp1 ) )
+ > - dt * tx1 * c1345 * tmp1
+ > - dt * tx1 * dx5
+
+c---------------------------------------------------------------------
+c form the second block sub-diagonal
+c---------------------------------------------------------------------
+ tmp1 = 1.0d+00 / u(1,i,j+1,k)
+ tmp2 = tmp1 * tmp1
+ tmp3 = tmp1 * tmp2
+
+ b(1,1,i,j) = - dt * ty1 * dy1
+ b(1,2,i,j) = 0.0d+00
+ b(1,3,i,j) = dt * ty2
+ b(1,4,i,j) = 0.0d+00
+ b(1,5,i,j) = 0.0d+00
+
+ b(2,1,i,j) = dt * ty2
+ > * ( - ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
+ > - dt * ty1 * ( - c34 * tmp2 * u(2,i,j+1,k) )
+ b(2,2,i,j) = dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
+ > - dt * ty1 * ( c34 * tmp1 )
+ > - dt * ty1 * dy2
+ b(2,3,i,j) = dt * ty2 * ( u(2,i,j+1,k) * tmp1 )
+ b(2,4,i,j) = 0.0d+00
+ b(2,5,i,j) = 0.0d+00
+
+ b(3,1,i,j) = dt * ty2
+ > * ( - ( u(3,i,j+1,k) * tmp1 ) ** 2
+ > + 0.50d+00 * c2 * ( ( 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) )
+ > * tmp2 ) )
+ > - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j+1,k) )
+ b(3,2,i,j) = dt * ty2
+ > * ( - c2 * ( u(2,i,j+1,k) * tmp1 ) )
+ b(3,3,i,j) = dt * ty2 * ( ( 2.0d+00 - c2 )
+ > * ( u(3,i,j+1,k) * tmp1 ) )
+ > - dt * ty1 * ( r43 * c34 * tmp1 )
+ > - dt * ty1 * dy3
+ b(3,4,i,j) = dt * ty2
+ > * ( - c2 * ( u(4,i,j+1,k) * tmp1 ) )
+ b(3,5,i,j) = dt * ty2 * c2
+
+ b(4,1,i,j) = dt * ty2
+ > * ( - ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
+ > - dt * ty1 * ( - c34 * tmp2 * u(4,i,j+1,k) )
+ b(4,2,i,j) = 0.0d+00
+ b(4,3,i,j) = dt * ty2 * ( u(4,i,j+1,k) * tmp1 )
+ b(4,4,i,j) = dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
+ > - dt * ty1 * ( c34 * tmp1 )
+ > - dt * ty1 * dy4
+ b(4,5,i,j) = 0.0d+00
+
+ b(5,1,i,j) = dt * ty2
+ > * ( ( c2 * ( 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) ) * tmp2
+ > - c1 * ( u(5,i,j+1,k) * tmp1 ) )
+ > * ( u(3,i,j+1,k) * tmp1 ) )
+ > - dt * ty1
+ > * ( - ( c34 - c1345 )*tmp3*(u(2,i,j+1,k)**2)
+ > - ( r43*c34 - c1345 )*tmp3*(u(3,i,j+1,k)**2)
+ > - ( c34 - c1345 )*tmp3*(u(4,i,j+1,k)**2)
+ > - c1345*tmp2*u(5,i,j+1,k) )
+ b(5,2,i,j) = dt * ty2
+ > * ( - c2 * ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
+ > - dt * ty1
+ > * ( c34 - c1345 ) * tmp2 * u(2,i,j+1,k)
+ b(5,3,i,j) = dt * ty2
+ > * ( c1 * ( u(5,i,j+1,k) * tmp1 )
+ > - 0.50d+00 * c2
+ > * ( ( u(2,i,j+1,k)*u(2,i,j+1,k)
+ > + 3.0d+00 * u(3,i,j+1,k)*u(3,i,j+1,k)
+ > + u(4,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 ) )
+ > - dt * ty1
+ > * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j+1,k)
+ b(5,4,i,j) = dt * ty2
+ > * ( - c2 * ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
+ > - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j+1,k)
+ b(5,5,i,j) = dt * ty2
+ > * ( c1 * ( u(3,i,j+1,k) * tmp1 ) )
+ > - dt * ty1 * c1345 * tmp1
+ > - dt * ty1 * dy5
+
+c---------------------------------------------------------------------
+c form the third block sub-diagonal
+c---------------------------------------------------------------------
+ tmp1 = 1.0d+00 / u(1,i,j,k+1)
+ tmp2 = tmp1 * tmp1
+ tmp3 = tmp1 * tmp2
+
+ c(1,1,i,j) = - dt * tz1 * dz1
+ c(1,2,i,j) = 0.0d+00
+ c(1,3,i,j) = 0.0d+00
+ c(1,4,i,j) = dt * tz2
+ c(1,5,i,j) = 0.0d+00
+
+ c(2,1,i,j) = dt * tz2
+ > * ( - ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
+ > - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k+1) )
+ c(2,2,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
+ > - dt * tz1 * c34 * tmp1
+ > - dt * tz1 * dz2
+ c(2,3,i,j) = 0.0d+00
+ c(2,4,i,j) = dt * tz2 * ( u(2,i,j,k+1) * tmp1 )
+ c(2,5,i,j) = 0.0d+00
+
+ c(3,1,i,j) = dt * tz2
+ > * ( - ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
+ > - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k+1) )
+ c(3,2,i,j) = 0.0d+00
+ c(3,3,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
+ > - dt * tz1 * ( c34 * tmp1 )
+ > - dt * tz1 * dz3
+ c(3,4,i,j) = dt * tz2 * ( u(3,i,j,k+1) * tmp1 )
+ c(3,5,i,j) = 0.0d+00
+
+ c(4,1,i,j) = dt * tz2
+ > * ( - ( u(4,i,j,k+1) * tmp1 ) ** 2
+ > + 0.50d+00 * c2
+ > * ( ( 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) ) * tmp2 ) )
+ > - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k+1) )
+ c(4,2,i,j) = dt * tz2
+ > * ( - c2 * ( u(2,i,j,k+1) * tmp1 ) )
+ c(4,3,i,j) = dt * tz2
+ > * ( - c2 * ( u(3,i,j,k+1) * tmp1 ) )
+ c(4,4,i,j) = dt * tz2 * ( 2.0d+00 - c2 )
+ > * ( u(4,i,j,k+1) * tmp1 )
+ > - dt * tz1 * ( r43 * c34 * tmp1 )
+ > - dt * tz1 * dz4
+ c(4,5,i,j) = dt * tz2 * c2
+
+ c(5,1,i,j) = dt * tz2
+ > * ( ( c2 * ( 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) ) * tmp2
+ > - c1 * ( u(5,i,j,k+1) * tmp1 ) )
+ > * ( u(4,i,j,k+1) * tmp1 ) )
+ > - dt * tz1
+ > * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k+1)**2)
+ > - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k+1)**2)
+ > - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k+1)**2)
+ > - c1345 * tmp2 * u(5,i,j,k+1) )
+ c(5,2,i,j) = dt * tz2
+ > * ( - c2 * ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
+ > - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k+1)
+ c(5,3,i,j) = dt * tz2
+ > * ( - c2 * ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
+ > - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k+1)
+ c(5,4,i,j) = dt * tz2
+ > * ( c1 * ( u(5,i,j,k+1) * tmp1 )
+ > - 0.50d+00 * c2
+ > * ( ( u(2,i,j,k+1)*u(2,i,j,k+1)
+ > + u(3,i,j,k+1)*u(3,i,j,k+1)
+ > + 3.0d+00*u(4,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 ) )
+ > - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k+1)
+ c(5,5,i,j) = dt * tz2
+ > * ( c1 * ( u(4,i,j,k+1) * tmp1 ) )
+ > - dt * tz1 * c1345 * tmp1
+ > - dt * tz1 * dz5
+
+ end do
+ end do
+
+ return
+ end
