2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
11 c---------------------------------------------------------------------
12 c compute the lower triangular part of the jacobian matrix
13 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,j,k-1)
146 a(1,1,i,j) = - dt * tz1 * dz1
149 a(1,4,i,j) = - dt * tz2
152 a(2,1,i,j) = - dt * tz2
153 > * ( - ( u(2,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
154 > - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k-1) )
155 a(2,2,i,j) = - dt * tz2 * ( u(4,i,j,k-1) * tmp1 )
156 > - dt * tz1 * c34 * tmp1
159 a(2,4,i,j) = - dt * tz2 * ( u(2,i,j,k-1) * tmp1 )
162 a(3,1,i,j) = - dt * tz2
163 > * ( - ( u(3,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
164 > - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k-1) )
166 a(3,3,i,j) = - dt * tz2 * ( u(4,i,j,k-1) * tmp1 )
167 > - dt * tz1 * ( c34 * tmp1 )
169 a(3,4,i,j) = - dt * tz2 * ( u(3,i,j,k-1) * tmp1 )
172 a(4,1,i,j) = - dt * tz2
173 > * ( - ( u(4,i,j,k-1) * tmp1 ) ** 2
175 > * ( ( u(2,i,j,k-1) * u(2,i,j,k-1)
176 > + u(3,i,j,k-1) * u(3,i,j,k-1)
177 > + u(4,i,j,k-1) * u(4,i,j,k-1) ) * tmp2 ) )
178 > - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k-1) )
179 a(4,2,i,j) = - dt * tz2
180 > * ( - c2 * ( u(2,i,j,k-1) * tmp1 ) )
181 a(4,3,i,j) = - dt * tz2
182 > * ( - c2 * ( u(3,i,j,k-1) * tmp1 ) )
183 a(4,4,i,j) = - dt * tz2 * ( 2.0d+00 - c2 )
184 > * ( u(4,i,j,k-1) * tmp1 )
185 > - dt * tz1 * ( r43 * c34 * tmp1 )
187 a(4,5,i,j) = - dt * tz2 * c2
189 a(5,1,i,j) = - dt * tz2
190 > * ( ( c2 * ( u(2,i,j,k-1) * u(2,i,j,k-1)
191 > + u(3,i,j,k-1) * u(3,i,j,k-1)
192 > + u(4,i,j,k-1) * u(4,i,j,k-1) ) * tmp2
193 > - c1 * ( u(5,i,j,k-1) * tmp1 ) )
194 > * ( u(4,i,j,k-1) * tmp1 ) )
196 > * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k-1)**2)
197 > - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k-1)**2)
198 > - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k-1)**2)
199 > - c1345 * tmp2 * u(5,i,j,k-1) )
200 a(5,2,i,j) = - dt * tz2
201 > * ( - c2 * ( u(2,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
202 > - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k-1)
203 a(5,3,i,j) = - dt * tz2
204 > * ( - c2 * ( u(3,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
205 > - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k-1)
206 a(5,4,i,j) = - dt * tz2
207 > * ( c1 * ( u(5,i,j,k-1) * tmp1 )
209 > * ( ( u(2,i,j,k-1)*u(2,i,j,k-1)
210 > + u(3,i,j,k-1)*u(3,i,j,k-1)
211 > + 3.0d+00*u(4,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 ) )
212 > - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k-1)
213 a(5,5,i,j) = - dt * tz2
214 > * ( c1 * ( u(4,i,j,k-1) * tmp1 ) )
215 > - dt * tz1 * c1345 * tmp1
218 c---------------------------------------------------------------------
219 c form the second block sub-diagonal
220 c---------------------------------------------------------------------
221 tmp1 = 1.0d+00 / u(1,i,j-1,k)
225 b(1,1,i,j) = - dt * ty1 * dy1
227 b(1,3,i,j) = - dt * ty2
231 b(2,1,i,j) = - dt * ty2
232 > * ( - ( u(2,i,j-1,k)*u(3,i,j-1,k) ) * tmp2 )
233 > - dt * ty1 * ( - c34 * tmp2 * u(2,i,j-1,k) )
234 b(2,2,i,j) = - dt * ty2 * ( u(3,i,j-1,k) * tmp1 )
235 > - dt * ty1 * ( c34 * tmp1 )
237 b(2,3,i,j) = - dt * ty2 * ( u(2,i,j-1,k) * tmp1 )
241 b(3,1,i,j) = - dt * ty2
242 > * ( - ( u(3,i,j-1,k) * tmp1 ) ** 2
243 > + 0.50d+00 * c2 * ( ( u(2,i,j-1,k) * u(2,i,j-1,k)
244 > + u(3,i,j-1,k) * u(3,i,j-1,k)
245 > + u(4,i,j-1,k) * u(4,i,j-1,k) )
247 > - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j-1,k) )
248 b(3,2,i,j) = - dt * ty2
249 > * ( - c2 * ( u(2,i,j-1,k) * tmp1 ) )
250 b(3,3,i,j) = - dt * ty2 * ( ( 2.0d+00 - c2 )
251 > * ( u(3,i,j-1,k) * tmp1 ) )
252 > - dt * ty1 * ( r43 * c34 * tmp1 )
254 b(3,4,i,j) = - dt * ty2
255 > * ( - c2 * ( u(4,i,j-1,k) * tmp1 ) )
256 b(3,5,i,j) = - dt * ty2 * c2
258 b(4,1,i,j) = - dt * ty2
259 > * ( - ( u(3,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 )
260 > - dt * ty1 * ( - c34 * tmp2 * u(4,i,j-1,k) )
262 b(4,3,i,j) = - dt * ty2 * ( u(4,i,j-1,k) * tmp1 )
263 b(4,4,i,j) = - dt * ty2 * ( u(3,i,j-1,k) * tmp1 )
264 > - dt * ty1 * ( c34 * tmp1 )
268 b(5,1,i,j) = - dt * ty2
269 > * ( ( c2 * ( u(2,i,j-1,k) * u(2,i,j-1,k)
270 > + u(3,i,j-1,k) * u(3,i,j-1,k)
271 > + u(4,i,j-1,k) * u(4,i,j-1,k) ) * tmp2
272 > - c1 * ( u(5,i,j-1,k) * tmp1 ) )
273 > * ( u(3,i,j-1,k) * tmp1 ) )
275 > * ( - ( c34 - c1345 )*tmp3*(u(2,i,j-1,k)**2)
276 > - ( r43*c34 - c1345 )*tmp3*(u(3,i,j-1,k)**2)
277 > - ( c34 - c1345 )*tmp3*(u(4,i,j-1,k)**2)
278 > - c1345*tmp2*u(5,i,j-1,k) )
279 b(5,2,i,j) = - dt * ty2
280 > * ( - c2 * ( u(2,i,j-1,k)*u(3,i,j-1,k) ) * tmp2 )
282 > * ( c34 - c1345 ) * tmp2 * u(2,i,j-1,k)
283 b(5,3,i,j) = - dt * ty2
284 > * ( c1 * ( u(5,i,j-1,k) * tmp1 )
286 > * ( ( u(2,i,j-1,k)*u(2,i,j-1,k)
287 > + 3.0d+00 * u(3,i,j-1,k)*u(3,i,j-1,k)
288 > + u(4,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 ) )
290 > * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j-1,k)
291 b(5,4,i,j) = - dt * ty2
292 > * ( - c2 * ( u(3,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 )
293 > - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j-1,k)
294 b(5,5,i,j) = - dt * ty2
295 > * ( c1 * ( u(3,i,j-1,k) * tmp1 ) )
296 > - dt * ty1 * c1345 * tmp1
299 c---------------------------------------------------------------------
300 c form the third block sub-diagonal
301 c---------------------------------------------------------------------
302 tmp1 = 1.0d+00 / u(1,i-1,j,k)
306 c(1,1,i,j) = - dt * tx1 * dx1
307 c(1,2,i,j) = - dt * tx2
312 c(2,1,i,j) = - dt * tx2
313 > * ( - ( u(2,i-1,j,k) * tmp1 ) ** 2
314 > + c2 * 0.50d+00 * ( u(2,i-1,j,k) * u(2,i-1,j,k)
315 > + u(3,i-1,j,k) * u(3,i-1,j,k)
316 > + u(4,i-1,j,k) * u(4,i-1,j,k) ) * tmp2 )
317 > - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i-1,j,k) )
318 c(2,2,i,j) = - dt * tx2
319 > * ( ( 2.0d+00 - c2 ) * ( u(2,i-1,j,k) * tmp1 ) )
320 > - dt * tx1 * ( r43 * c34 * tmp1 )
322 c(2,3,i,j) = - dt * tx2
323 > * ( - c2 * ( u(3,i-1,j,k) * tmp1 ) )
324 c(2,4,i,j) = - dt * tx2
325 > * ( - c2 * ( u(4,i-1,j,k) * tmp1 ) )
326 c(2,5,i,j) = - dt * tx2 * c2
328 c(3,1,i,j) = - dt * tx2
329 > * ( - ( u(2,i-1,j,k) * u(3,i-1,j,k) ) * tmp2 )
330 > - dt * tx1 * ( - c34 * tmp2 * u(3,i-1,j,k) )
331 c(3,2,i,j) = - dt * tx2 * ( u(3,i-1,j,k) * tmp1 )
332 c(3,3,i,j) = - dt * tx2 * ( u(2,i-1,j,k) * tmp1 )
333 > - dt * tx1 * ( c34 * tmp1 )
338 c(4,1,i,j) = - dt * tx2
339 > * ( - ( u(2,i-1,j,k)*u(4,i-1,j,k) ) * tmp2 )
340 > - dt * tx1 * ( - c34 * tmp2 * u(4,i-1,j,k) )
341 c(4,2,i,j) = - dt * tx2 * ( u(4,i-1,j,k) * tmp1 )
343 c(4,4,i,j) = - dt * tx2 * ( u(2,i-1,j,k) * tmp1 )
344 > - dt * tx1 * ( c34 * tmp1 )
348 c(5,1,i,j) = - dt * tx2
349 > * ( ( c2 * ( u(2,i-1,j,k) * u(2,i-1,j,k)
350 > + u(3,i-1,j,k) * u(3,i-1,j,k)
351 > + u(4,i-1,j,k) * u(4,i-1,j,k) ) * tmp2
352 > - c1 * ( u(5,i-1,j,k) * tmp1 ) )
353 > * ( u(2,i-1,j,k) * tmp1 ) )
355 > * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i-1,j,k)**2 )
356 > - ( c34 - c1345 ) * tmp3 * ( u(3,i-1,j,k)**2 )
357 > - ( c34 - c1345 ) * tmp3 * ( u(4,i-1,j,k)**2 )
358 > - c1345 * tmp2 * u(5,i-1,j,k) )
359 c(5,2,i,j) = - dt * tx2
360 > * ( c1 * ( u(5,i-1,j,k) * tmp1 )
362 > * ( ( 3.0d+00*u(2,i-1,j,k)*u(2,i-1,j,k)
363 > + u(3,i-1,j,k)*u(3,i-1,j,k)
364 > + u(4,i-1,j,k)*u(4,i-1,j,k) ) * tmp2 ) )
366 > * ( r43*c34 - c1345 ) * tmp2 * u(2,i-1,j,k)
367 c(5,3,i,j) = - dt * tx2
368 > * ( - c2 * ( u(3,i-1,j,k)*u(2,i-1,j,k) ) * tmp2 )
370 > * ( c34 - c1345 ) * tmp2 * u(3,i-1,j,k)
371 c(5,4,i,j) = - dt * tx2
372 > * ( - c2 * ( u(4,i-1,j,k)*u(2,i-1,j,k) ) * tmp2 )
374 > * ( c34 - c1345 ) * tmp2 * u(4,i-1,j,k)
375 c(5,5,i,j) = - dt * tx2
376 > * ( c1 * ( u(2,i-1,j,k) * tmp1 ) )
377 > - dt * tx1 * c1345 * tmp1