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