2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
10 c---------------------------------------------------------------------
11 c compute the upper triangular part of the jacobian matrix
12 c---------------------------------------------------------------------
19 c---------------------------------------------------------------------
21 c---------------------------------------------------------------------
24 c---------------------------------------------------------------------
26 c---------------------------------------------------------------------
29 double precision c1345
31 double precision tmp1, tmp2, tmp3
35 r43 = ( 4.0d+00 / 3.0d+00 )
36 c1345 = c1 * c3 * c4 * c5
42 c---------------------------------------------------------------------
43 c form the block daigonal
44 c---------------------------------------------------------------------
45 tmp1 = 1.0d+00 / u(1,i,j,k)
50 > + dt * 2.0d+00 * ( tx1 * dx1
58 d(2,1,i,j) = dt * 2.0d+00
59 > * ( tx1 * ( - r43 * c34 * tmp2 * u(2,i,j,k) )
60 > + ty1 * ( - c34 * tmp2 * u(2,i,j,k) )
61 > + tz1 * ( - c34 * tmp2 * u(2,i,j,k) ) )
64 > * ( tx1 * r43 * c34 * tmp1
66 > + tz1 * c34 * tmp1 )
67 > + dt * 2.0d+00 * ( tx1 * dx2
74 d(3,1,i,j) = dt * 2.0d+00
75 > * ( tx1 * ( - c34 * tmp2 * u(3,i,j,k) )
76 > + ty1 * ( - r43 * c34 * tmp2 * u(3,i,j,k) )
77 > + tz1 * ( - c34 * tmp2 * u(3,i,j,k) ) )
81 > * ( tx1 * c34 * tmp1
82 > + ty1 * r43 * c34 * tmp1
83 > + tz1 * c34 * tmp1 )
84 > + dt * 2.0d+00 * ( tx1 * dx3
90 d(4,1,i,j) = dt * 2.0d+00
91 > * ( tx1 * ( - c34 * tmp2 * u(4,i,j,k) )
92 > + ty1 * ( - c34 * tmp2 * u(4,i,j,k) )
93 > + tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k) ) )
98 > * ( tx1 * c34 * tmp1
100 > + tz1 * r43 * c34 * tmp1 )
101 > + dt * 2.0d+00 * ( tx1 * dx4
106 d(5,1,i,j) = dt * 2.0d+00
107 > * ( tx1 * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
108 > - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
109 > - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
110 > - ( c1345 ) * tmp2 * u(5,i,j,k) )
111 > + ty1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
112 > - ( r43*c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
113 > - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
114 > - ( c1345 ) * tmp2 * u(5,i,j,k) )
115 > + tz1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
116 > - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
117 > - ( r43*c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
118 > - ( c1345 ) * tmp2 * u(5,i,j,k) ) )
119 d(5,2,i,j) = dt * 2.0d+00
120 > * ( tx1 * ( r43*c34 - c1345 ) * tmp2 * u(2,i,j,k)
121 > + ty1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k)
122 > + tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k) )
123 d(5,3,i,j) = dt * 2.0d+00
124 > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k)
125 > + ty1 * ( r43*c34 -c1345 ) * tmp2 * u(3,i,j,k)
126 > + tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k) )
127 d(5,4,i,j) = dt * 2.0d+00
128 > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
129 > + ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
130 > + tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k) )
132 > + dt * 2.0d+00 * ( tx1 * c1345 * tmp1
133 > + ty1 * c1345 * tmp1
134 > + tz1 * c1345 * tmp1 )
135 > + dt * 2.0d+00 * ( tx1 * dx5
139 c---------------------------------------------------------------------
140 c form the first block sub-diagonal
141 c---------------------------------------------------------------------
142 tmp1 = 1.0d+00 / u(1,i+1,j,k)
146 a(1,1,i,j) = - dt * tx1 * dx1
147 a(1,2,i,j) = dt * tx2
152 a(2,1,i,j) = dt * tx2
153 > * ( - ( u(2,i+1,j,k) * tmp1 ) ** 2
154 > + c2 * 0.50d+00 * ( u(2,i+1,j,k) * u(2,i+1,j,k)
155 > + u(3,i+1,j,k) * u(3,i+1,j,k)
156 > + u(4,i+1,j,k) * u(4,i+1,j,k) ) * tmp2 )
157 > - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i+1,j,k) )
158 a(2,2,i,j) = dt * tx2
159 > * ( ( 2.0d+00 - c2 ) * ( u(2,i+1,j,k) * tmp1 ) )
160 > - dt * tx1 * ( r43 * c34 * tmp1 )
162 a(2,3,i,j) = dt * tx2
163 > * ( - c2 * ( u(3,i+1,j,k) * tmp1 ) )
164 a(2,4,i,j) = dt * tx2
165 > * ( - c2 * ( u(4,i+1,j,k) * tmp1 ) )
166 a(2,5,i,j) = dt * tx2 * c2
168 a(3,1,i,j) = dt * tx2
169 > * ( - ( u(2,i+1,j,k) * u(3,i+1,j,k) ) * tmp2 )
170 > - dt * tx1 * ( - c34 * tmp2 * u(3,i+1,j,k) )
171 a(3,2,i,j) = dt * tx2 * ( u(3,i+1,j,k) * tmp1 )
172 a(3,3,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
173 > - dt * tx1 * ( c34 * tmp1 )
178 a(4,1,i,j) = dt * tx2
179 > * ( - ( u(2,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 )
180 > - dt * tx1 * ( - c34 * tmp2 * u(4,i+1,j,k) )
181 a(4,2,i,j) = dt * tx2 * ( u(4,i+1,j,k) * tmp1 )
183 a(4,4,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
184 > - dt * tx1 * ( c34 * tmp1 )
188 a(5,1,i,j) = dt * tx2
189 > * ( ( c2 * ( u(2,i+1,j,k) * u(2,i+1,j,k)
190 > + u(3,i+1,j,k) * u(3,i+1,j,k)
191 > + u(4,i+1,j,k) * u(4,i+1,j,k) ) * tmp2
192 > - c1 * ( u(5,i+1,j,k) * tmp1 ) )
193 > * ( u(2,i+1,j,k) * tmp1 ) )
195 > * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i+1,j,k)**2 )
196 > - ( c34 - c1345 ) * tmp3 * ( u(3,i+1,j,k)**2 )
197 > - ( c34 - c1345 ) * tmp3 * ( u(4,i+1,j,k)**2 )
198 > - c1345 * tmp2 * u(5,i+1,j,k) )
199 a(5,2,i,j) = dt * tx2
200 > * ( c1 * ( u(5,i+1,j,k) * tmp1 )
202 > * ( ( 3.0d+00*u(2,i+1,j,k)*u(2,i+1,j,k)
203 > + u(3,i+1,j,k)*u(3,i+1,j,k)
204 > + u(4,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 ) )
206 > * ( r43*c34 - c1345 ) * tmp2 * u(2,i+1,j,k)
207 a(5,3,i,j) = dt * tx2
208 > * ( - c2 * ( u(3,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
210 > * ( c34 - c1345 ) * tmp2 * u(3,i+1,j,k)
211 a(5,4,i,j) = dt * tx2
212 > * ( - c2 * ( u(4,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
214 > * ( c34 - c1345 ) * tmp2 * u(4,i+1,j,k)
215 a(5,5,i,j) = dt * tx2
216 > * ( c1 * ( u(2,i+1,j,k) * tmp1 ) )
217 > - dt * tx1 * c1345 * tmp1
220 c---------------------------------------------------------------------
221 c form the second block sub-diagonal
222 c---------------------------------------------------------------------
223 tmp1 = 1.0d+00 / u(1,i,j+1,k)
227 b(1,1,i,j) = - dt * ty1 * dy1
229 b(1,3,i,j) = dt * ty2
233 b(2,1,i,j) = dt * ty2
234 > * ( - ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
235 > - dt * ty1 * ( - c34 * tmp2 * u(2,i,j+1,k) )
236 b(2,2,i,j) = dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
237 > - dt * ty1 * ( c34 * tmp1 )
239 b(2,3,i,j) = dt * ty2 * ( u(2,i,j+1,k) * tmp1 )
243 b(3,1,i,j) = dt * ty2
244 > * ( - ( u(3,i,j+1,k) * tmp1 ) ** 2
245 > + 0.50d+00 * c2 * ( ( u(2,i,j+1,k) * u(2,i,j+1,k)
246 > + u(3,i,j+1,k) * u(3,i,j+1,k)
247 > + u(4,i,j+1,k) * u(4,i,j+1,k) )
249 > - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j+1,k) )
250 b(3,2,i,j) = dt * ty2
251 > * ( - c2 * ( u(2,i,j+1,k) * tmp1 ) )
252 b(3,3,i,j) = dt * ty2 * ( ( 2.0d+00 - c2 )
253 > * ( u(3,i,j+1,k) * tmp1 ) )
254 > - dt * ty1 * ( r43 * c34 * tmp1 )
256 b(3,4,i,j) = dt * ty2
257 > * ( - c2 * ( u(4,i,j+1,k) * tmp1 ) )
258 b(3,5,i,j) = dt * ty2 * c2
260 b(4,1,i,j) = dt * ty2
261 > * ( - ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
262 > - dt * ty1 * ( - c34 * tmp2 * u(4,i,j+1,k) )
264 b(4,3,i,j) = dt * ty2 * ( u(4,i,j+1,k) * tmp1 )
265 b(4,4,i,j) = dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
266 > - dt * ty1 * ( c34 * tmp1 )
270 b(5,1,i,j) = dt * ty2
271 > * ( ( c2 * ( u(2,i,j+1,k) * u(2,i,j+1,k)
272 > + u(3,i,j+1,k) * u(3,i,j+1,k)
273 > + u(4,i,j+1,k) * u(4,i,j+1,k) ) * tmp2
274 > - c1 * ( u(5,i,j+1,k) * tmp1 ) )
275 > * ( u(3,i,j+1,k) * tmp1 ) )
277 > * ( - ( c34 - c1345 )*tmp3*(u(2,i,j+1,k)**2)
278 > - ( r43*c34 - c1345 )*tmp3*(u(3,i,j+1,k)**2)
279 > - ( c34 - c1345 )*tmp3*(u(4,i,j+1,k)**2)
280 > - c1345*tmp2*u(5,i,j+1,k) )
281 b(5,2,i,j) = dt * ty2
282 > * ( - c2 * ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
284 > * ( c34 - c1345 ) * tmp2 * u(2,i,j+1,k)
285 b(5,3,i,j) = dt * ty2
286 > * ( c1 * ( u(5,i,j+1,k) * tmp1 )
288 > * ( ( u(2,i,j+1,k)*u(2,i,j+1,k)
289 > + 3.0d+00 * u(3,i,j+1,k)*u(3,i,j+1,k)
290 > + u(4,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 ) )
292 > * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j+1,k)
293 b(5,4,i,j) = dt * ty2
294 > * ( - c2 * ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
295 > - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j+1,k)
296 b(5,5,i,j) = dt * ty2
297 > * ( c1 * ( u(3,i,j+1,k) * tmp1 ) )
298 > - dt * ty1 * c1345 * tmp1
301 c---------------------------------------------------------------------
302 c form the third block sub-diagonal
303 c---------------------------------------------------------------------
304 tmp1 = 1.0d+00 / u(1,i,j,k+1)
308 c(1,1,i,j) = - dt * tz1 * dz1
311 c(1,4,i,j) = dt * tz2
314 c(2,1,i,j) = dt * tz2
315 > * ( - ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
316 > - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k+1) )
317 c(2,2,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
318 > - dt * tz1 * c34 * tmp1
321 c(2,4,i,j) = dt * tz2 * ( u(2,i,j,k+1) * tmp1 )
324 c(3,1,i,j) = dt * tz2
325 > * ( - ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
326 > - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k+1) )
328 c(3,3,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
329 > - dt * tz1 * ( c34 * tmp1 )
331 c(3,4,i,j) = dt * tz2 * ( u(3,i,j,k+1) * tmp1 )
334 c(4,1,i,j) = dt * tz2
335 > * ( - ( u(4,i,j,k+1) * tmp1 ) ** 2
337 > * ( ( u(2,i,j,k+1) * u(2,i,j,k+1)
338 > + u(3,i,j,k+1) * u(3,i,j,k+1)
339 > + u(4,i,j,k+1) * u(4,i,j,k+1) ) * tmp2 ) )
340 > - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k+1) )
341 c(4,2,i,j) = dt * tz2
342 > * ( - c2 * ( u(2,i,j,k+1) * tmp1 ) )
343 c(4,3,i,j) = dt * tz2
344 > * ( - c2 * ( u(3,i,j,k+1) * tmp1 ) )
345 c(4,4,i,j) = dt * tz2 * ( 2.0d+00 - c2 )
346 > * ( u(4,i,j,k+1) * tmp1 )
347 > - dt * tz1 * ( r43 * c34 * tmp1 )
349 c(4,5,i,j) = dt * tz2 * c2
351 c(5,1,i,j) = dt * tz2
352 > * ( ( c2 * ( u(2,i,j,k+1) * u(2,i,j,k+1)
353 > + u(3,i,j,k+1) * u(3,i,j,k+1)
354 > + u(4,i,j,k+1) * u(4,i,j,k+1) ) * tmp2
355 > - c1 * ( u(5,i,j,k+1) * tmp1 ) )
356 > * ( u(4,i,j,k+1) * tmp1 ) )
358 > * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k+1)**2)
359 > - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k+1)**2)
360 > - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k+1)**2)
361 > - c1345 * tmp2 * u(5,i,j,k+1) )
362 c(5,2,i,j) = dt * tz2
363 > * ( - c2 * ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
364 > - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k+1)
365 c(5,3,i,j) = dt * tz2
366 > * ( - c2 * ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
367 > - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k+1)
368 c(5,4,i,j) = dt * tz2
369 > * ( c1 * ( u(5,i,j,k+1) * tmp1 )
371 > * ( ( u(2,i,j,k+1)*u(2,i,j,k+1)
372 > + u(3,i,j,k+1)*u(3,i,j,k+1)
373 > + 3.0d+00*u(4,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 ) )
374 > - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k+1)
375 c(5,5,i,j) = dt * tz2
376 > * ( c1 * ( u(4,i,j,k+1) * tmp1 ) )
377 > - dt * tz1 * c1345 * tmp1