2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
10 c---------------------------------------------------------------------
11 c compute the right hand sides
12 c---------------------------------------------------------------------
18 c---------------------------------------------------------------------
20 c---------------------------------------------------------------------
27 double precision u21, u31, u41
29 double precision u21i, u31i, u41i, u51i
30 double precision u21j, u31j, u41j, u51j
31 double precision u21k, u31k, u41k, u51k
32 double precision u21im1, u31im1, u41im1, u51im1
33 double precision u21jm1, u31jm1, u41jm1, u51jm1
34 double precision u21km1, u31km1, u41km1, u51km1
42 rsd(m,i,j,k) = - frct(m,i,j,k)
48 c---------------------------------------------------------------------
49 c xi-direction flux differences
50 c---------------------------------------------------------------------
52 c---------------------------------------------------------------------
53 c iex = flag : iex = 0 north/south communication
54 c : iex = 1 east/west communication
55 c---------------------------------------------------------------------
58 c---------------------------------------------------------------------
59 c communicate and receive/send two rows of data
60 c---------------------------------------------------------------------
61 call exchange_3(u,iex)
64 if (north.eq.-1) L1 = 1
66 if (south.eq.-1) L2 = nx
70 if (north.eq.-1) ist1 = 4
71 if (south.eq.-1) iend1 = nx - 3
76 flux(1,i,j,k) = u(2,i,j,k)
77 u21 = u(2,i,j,k) / u(1,i,j,k)
79 q = 0.50d+00 * ( u(2,i,j,k) * u(2,i,j,k)
80 > + u(3,i,j,k) * u(3,i,j,k)
81 > + u(4,i,j,k) * u(4,i,j,k) )
84 flux(2,i,j,k) = u(2,i,j,k) * u21 + c2 *
86 flux(3,i,j,k) = u(3,i,j,k) * u21
87 flux(4,i,j,k) = u(4,i,j,k) * u21
88 flux(5,i,j,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
93 rsd(m,i,j,k) = rsd(m,i,j,k)
94 > - tx2 * ( flux(m,i+1,j,k) - flux(m,i-1,j,k) )
99 tmp = 1.0d+00 / u(1,i,j,k)
101 u21i = tmp * u(2,i,j,k)
102 u31i = tmp * u(3,i,j,k)
103 u41i = tmp * u(4,i,j,k)
104 u51i = tmp * u(5,i,j,k)
106 tmp = 1.0d+00 / u(1,i-1,j,k)
108 u21im1 = tmp * u(2,i-1,j,k)
109 u31im1 = tmp * u(3,i-1,j,k)
110 u41im1 = tmp * u(4,i-1,j,k)
111 u51im1 = tmp * u(5,i-1,j,k)
113 flux(2,i,j,k) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
114 flux(3,i,j,k) = tx3 * ( u31i - u31im1 )
115 flux(4,i,j,k) = tx3 * ( u41i - u41im1 )
116 flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
117 > * tx3 * ( ( u21i **2 + u31i **2 + u41i **2 )
118 > - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
119 > + (1.0d+00/6.0d+00)
120 > * tx3 * ( u21i**2 - u21im1**2 )
121 > + c1 * c5 * tx3 * ( u51i - u51im1 )
125 rsd(1,i,j,k) = rsd(1,i,j,k)
126 > + dx1 * tx1 * ( u(1,i-1,j,k)
127 > - 2.0d+00 * u(1,i,j,k)
129 rsd(2,i,j,k) = rsd(2,i,j,k)
130 > + tx3 * c3 * c4 * ( flux(2,i+1,j,k) - flux(2,i,j,k) )
131 > + dx2 * tx1 * ( u(2,i-1,j,k)
132 > - 2.0d+00 * u(2,i,j,k)
134 rsd(3,i,j,k) = rsd(3,i,j,k)
135 > + tx3 * c3 * c4 * ( flux(3,i+1,j,k) - flux(3,i,j,k) )
136 > + dx3 * tx1 * ( u(3,i-1,j,k)
137 > - 2.0d+00 * u(3,i,j,k)
139 rsd(4,i,j,k) = rsd(4,i,j,k)
140 > + tx3 * c3 * c4 * ( flux(4,i+1,j,k) - flux(4,i,j,k) )
141 > + dx4 * tx1 * ( u(4,i-1,j,k)
142 > - 2.0d+00 * u(4,i,j,k)
144 rsd(5,i,j,k) = rsd(5,i,j,k)
145 > + tx3 * c3 * c4 * ( flux(5,i+1,j,k) - flux(5,i,j,k) )
146 > + dx5 * tx1 * ( u(5,i-1,j,k)
147 > - 2.0d+00 * u(5,i,j,k)
151 c---------------------------------------------------------------------
152 c Fourth-order dissipation
153 c---------------------------------------------------------------------
154 IF (north.eq.-1) then
156 rsd(m,2,j,k) = rsd(m,2,j,k)
157 > - dssp * ( + 5.0d+00 * u(m,2,j,k)
158 > - 4.0d+00 * u(m,3,j,k)
160 rsd(m,3,j,k) = rsd(m,3,j,k)
161 > - dssp * ( - 4.0d+00 * u(m,2,j,k)
162 > + 6.0d+00 * u(m,3,j,k)
163 > - 4.0d+00 * u(m,4,j,k)
170 rsd(m,i,j,k) = rsd(m,i,j,k)
171 > - dssp * ( u(m,i-2,j,k)
172 > - 4.0d+00 * u(m,i-1,j,k)
173 > + 6.0d+00 * u(m,i,j,k)
174 > - 4.0d+00 * u(m,i+1,j,k)
179 IF (south.eq.-1) then
181 rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
182 > - dssp * ( u(m,nx-4,j,k)
183 > - 4.0d+00 * u(m,nx-3,j,k)
184 > + 6.0d+00 * u(m,nx-2,j,k)
185 > - 4.0d+00 * u(m,nx-1,j,k) )
186 rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
187 > - dssp * ( u(m,nx-3,j,k)
188 > - 4.0d+00 * u(m,nx-2,j,k)
189 > + 5.0d+00 * u(m,nx-1,j,k) )
196 c---------------------------------------------------------------------
197 c eta-direction flux differences
198 c---------------------------------------------------------------------
200 c---------------------------------------------------------------------
201 c iex = flag : iex = 0 north/south communication
202 c---------------------------------------------------------------------
205 c---------------------------------------------------------------------
206 c communicate and receive/send two rows of data
207 c---------------------------------------------------------------------
208 call exchange_3(u,iex)
211 if (west.eq.-1) L1 = 1
213 if (east.eq.-1) L2 = ny
217 if (west.eq.-1) jst1 = 4
218 if (east.eq.-1) jend1 = ny - 3
223 flux(1,i,j,k) = u(3,i,j,k)
224 u31 = u(3,i,j,k) / u(1,i,j,k)
226 q = 0.50d+00 * ( u(2,i,j,k) * u(2,i,j,k)
227 > + u(3,i,j,k) * u(3,i,j,k)
228 > + u(4,i,j,k) * u(4,i,j,k) )
231 flux(2,i,j,k) = u(2,i,j,k) * u31
232 flux(3,i,j,k) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
233 flux(4,i,j,k) = u(4,i,j,k) * u31
234 flux(5,i,j,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
241 rsd(m,i,j,k) = rsd(m,i,j,k)
242 > - ty2 * ( flux(m,i,j+1,k) - flux(m,i,j-1,k) )
249 tmp = 1.0d+00 / u(1,i,j,k)
251 u21j = tmp * u(2,i,j,k)
252 u31j = tmp * u(3,i,j,k)
253 u41j = tmp * u(4,i,j,k)
254 u51j = tmp * u(5,i,j,k)
256 tmp = 1.0d+00 / u(1,i,j-1,k)
257 u21jm1 = tmp * u(2,i,j-1,k)
258 u31jm1 = tmp * u(3,i,j-1,k)
259 u41jm1 = tmp * u(4,i,j-1,k)
260 u51jm1 = tmp * u(5,i,j-1,k)
262 flux(2,i,j,k) = ty3 * ( u21j - u21jm1 )
263 flux(3,i,j,k) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
264 flux(4,i,j,k) = ty3 * ( u41j - u41jm1 )
265 flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
266 > * ty3 * ( ( u21j **2 + u31j **2 + u41j **2 )
267 > - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
268 > + (1.0d+00/6.0d+00)
269 > * ty3 * ( u31j**2 - u31jm1**2 )
270 > + c1 * c5 * ty3 * ( u51j - u51jm1 )
277 rsd(1,i,j,k) = rsd(1,i,j,k)
278 > + dy1 * ty1 * ( u(1,i,j-1,k)
279 > - 2.0d+00 * u(1,i,j,k)
282 rsd(2,i,j,k) = rsd(2,i,j,k)
283 > + ty3 * c3 * c4 * ( flux(2,i,j+1,k) - flux(2,i,j,k) )
284 > + dy2 * ty1 * ( u(2,i,j-1,k)
285 > - 2.0d+00 * u(2,i,j,k)
288 rsd(3,i,j,k) = rsd(3,i,j,k)
289 > + ty3 * c3 * c4 * ( flux(3,i,j+1,k) - flux(3,i,j,k) )
290 > + dy3 * ty1 * ( u(3,i,j-1,k)
291 > - 2.0d+00 * u(3,i,j,k)
294 rsd(4,i,j,k) = rsd(4,i,j,k)
295 > + ty3 * c3 * c4 * ( flux(4,i,j+1,k) - flux(4,i,j,k) )
296 > + dy4 * ty1 * ( u(4,i,j-1,k)
297 > - 2.0d+00 * u(4,i,j,k)
300 rsd(5,i,j,k) = rsd(5,i,j,k)
301 > + ty3 * c3 * c4 * ( flux(5,i,j+1,k) - flux(5,i,j,k) )
302 > + dy5 * ty1 * ( u(5,i,j-1,k)
303 > - 2.0d+00 * u(5,i,j,k)
309 c---------------------------------------------------------------------
310 c fourth-order dissipation
311 c---------------------------------------------------------------------
315 rsd(m,i,2,k) = rsd(m,i,2,k)
316 > - dssp * ( + 5.0d+00 * u(m,i,2,k)
317 > - 4.0d+00 * u(m,i,3,k)
319 rsd(m,i,3,k) = rsd(m,i,3,k)
320 > - dssp * ( - 4.0d+00 * u(m,i,2,k)
321 > + 6.0d+00 * u(m,i,3,k)
322 > - 4.0d+00 * u(m,i,4,k)
331 rsd(m,i,j,k) = rsd(m,i,j,k)
332 > - dssp * ( u(m,i,j-2,k)
333 > - 4.0d+00 * u(m,i,j-1,k)
334 > + 6.0d+00 * u(m,i,j,k)
335 > - 4.0d+00 * u(m,i,j+1,k)
344 rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
345 > - dssp * ( u(m,i,ny-4,k)
346 > - 4.0d+00 * u(m,i,ny-3,k)
347 > + 6.0d+00 * u(m,i,ny-2,k)
348 > - 4.0d+00 * u(m,i,ny-1,k) )
349 rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
350 > - dssp * ( u(m,i,ny-3,k)
351 > - 4.0d+00 * u(m,i,ny-2,k)
352 > + 5.0d+00 * u(m,i,ny-1,k) )
359 c---------------------------------------------------------------------
360 c zeta-direction flux differences
361 c---------------------------------------------------------------------
365 flux(1,i,j,k) = u(4,i,j,k)
366 u41 = u(4,i,j,k) / u(1,i,j,k)
368 q = 0.50d+00 * ( u(2,i,j,k) * u(2,i,j,k)
369 > + u(3,i,j,k) * u(3,i,j,k)
370 > + u(4,i,j,k) * u(4,i,j,k) )
373 flux(2,i,j,k) = u(2,i,j,k) * u41
374 flux(3,i,j,k) = u(3,i,j,k) * u41
375 flux(4,i,j,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q)
376 flux(5,i,j,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
385 rsd(m,i,j,k) = rsd(m,i,j,k)
386 > - tz2 * ( flux(m,i,j,k+1) - flux(m,i,j,k-1) )
395 tmp = 1.0d+00 / u(1,i,j,k)
397 u21k = tmp * u(2,i,j,k)
398 u31k = tmp * u(3,i,j,k)
399 u41k = tmp * u(4,i,j,k)
400 u51k = tmp * u(5,i,j,k)
402 tmp = 1.0d+00 / u(1,i,j,k-1)
404 u21km1 = tmp * u(2,i,j,k-1)
405 u31km1 = tmp * u(3,i,j,k-1)
406 u41km1 = tmp * u(4,i,j,k-1)
407 u51km1 = tmp * u(5,i,j,k-1)
409 flux(2,i,j,k) = tz3 * ( u21k - u21km1 )
410 flux(3,i,j,k) = tz3 * ( u31k - u31km1 )
411 flux(4,i,j,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
412 flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
413 > * tz3 * ( ( u21k **2 + u31k **2 + u41k **2 )
414 > - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
415 > + (1.0d+00/6.0d+00)
416 > * tz3 * ( u41k**2 - u41km1**2 )
417 > + c1 * c5 * tz3 * ( u51k - u51km1 )
425 rsd(1,i,j,k) = rsd(1,i,j,k)
426 > + dz1 * tz1 * ( u(1,i,j,k-1)
427 > - 2.0d+00 * u(1,i,j,k)
429 rsd(2,i,j,k) = rsd(2,i,j,k)
430 > + tz3 * c3 * c4 * ( flux(2,i,j,k+1) - flux(2,i,j,k) )
431 > + dz2 * tz1 * ( u(2,i,j,k-1)
432 > - 2.0d+00 * u(2,i,j,k)
434 rsd(3,i,j,k) = rsd(3,i,j,k)
435 > + tz3 * c3 * c4 * ( flux(3,i,j,k+1) - flux(3,i,j,k) )
436 > + dz3 * tz1 * ( u(3,i,j,k-1)
437 > - 2.0d+00 * u(3,i,j,k)
439 rsd(4,i,j,k) = rsd(4,i,j,k)
440 > + tz3 * c3 * c4 * ( flux(4,i,j,k+1) - flux(4,i,j,k) )
441 > + dz4 * tz1 * ( u(4,i,j,k-1)
442 > - 2.0d+00 * u(4,i,j,k)
444 rsd(5,i,j,k) = rsd(5,i,j,k)
445 > + tz3 * c3 * c4 * ( flux(5,i,j,k+1) - flux(5,i,j,k) )
446 > + dz5 * tz1 * ( u(5,i,j,k-1)
447 > - 2.0d+00 * u(5,i,j,k)
453 c---------------------------------------------------------------------
454 c fourth-order dissipation
455 c---------------------------------------------------------------------
459 rsd(m,i,j,2) = rsd(m,i,j,2)
460 > - dssp * ( + 5.0d+00 * u(m,i,j,2)
461 > - 4.0d+00 * u(m,i,j,3)
463 rsd(m,i,j,3) = rsd(m,i,j,3)
464 > - dssp * ( - 4.0d+00 * u(m,i,j,2)
465 > + 6.0d+00 * u(m,i,j,3)
466 > - 4.0d+00 * u(m,i,j,4)
476 rsd(m,i,j,k) = rsd(m,i,j,k)
477 > - dssp * ( u(m,i,j,k-2)
478 > - 4.0d+00 * u(m,i,j,k-1)
479 > + 6.0d+00 * u(m,i,j,k)
480 > - 4.0d+00 * u(m,i,j,k+1)
490 rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
491 > - dssp * ( u(m,i,j,nz-4)
492 > - 4.0d+00 * u(m,i,j,nz-3)
493 > + 6.0d+00 * u(m,i,j,nz-2)
494 > - 4.0d+00 * u(m,i,j,nz-1) )
495 rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
496 > - dssp * ( u(m,i,j,nz-3)
497 > - 4.0d+00 * u(m,i,j,nz-2)
498 > + 5.0d+00 * u(m,i,j,nz-1) )