1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
6 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
9 c---------------------------------------------------------------------
11 c compute the right hand side based on exact solution
13 c---------------------------------------------------------------------
19 c---------------------------------------------------------------------
21 c---------------------------------------------------------------------
28 double precision dsspm
29 double precision xi, eta, zeta
31 double precision u21, u31, u41
33 double precision u21i, u31i, u41i, u51i
34 double precision u21j, u31j, u41j, u51j
35 double precision u21k, u31k, u41k, u51k
36 double precision u21im1, u31im1, u41im1, u51im1
37 double precision u21jm1, u31jm1, u41jm1, u51jm1
38 double precision u21km1, u31km1, u41km1, u51km1
47 frct( m, i, j, k ) = 0.0d+00
54 zeta = ( dble(k-1) ) / ( nz - 1 )
57 eta = ( dble(jglob-1) ) / ( ny0 - 1 )
60 xi = ( dble(iglob-1) ) / ( nx0 - 1 )
62 rsd(m,i,j,k) = ce(m,1)
67 > + ce(m,6) * eta * eta
68 > + ce(m,7) * zeta * zeta
69 > + ce(m,8) * xi * xi * xi
70 > + ce(m,9) * eta * eta * eta
71 > + ce(m,10) * zeta * zeta * zeta
72 > + ce(m,11) * xi * xi * xi * xi
73 > + ce(m,12) * eta * eta * eta * eta
74 > + ce(m,13) * zeta * zeta * zeta * zeta
80 c---------------------------------------------------------------------
81 c xi-direction flux differences
82 c---------------------------------------------------------------------
84 c iex = flag : iex = 0 north/south communication
85 c : iex = 1 east/west communication
87 c---------------------------------------------------------------------
90 c---------------------------------------------------------------------
91 c communicate and receive/send two rows of data
92 c---------------------------------------------------------------------
93 call exchange_3 (rsd,iex)
96 if (north.eq.-1) L1 = 1
98 if (south.eq.-1) L2 = nx
102 if (north.eq.-1) ist1 = 4
103 if (south.eq.-1) iend1 = nx - 3
108 flux(1,i,j,k) = rsd(2,i,j,k)
109 u21 = rsd(2,i,j,k) / rsd(1,i,j,k)
110 q = 0.50d+00 * ( rsd(2,i,j,k) * rsd(2,i,j,k)
111 > + rsd(3,i,j,k) * rsd(3,i,j,k)
112 > + rsd(4,i,j,k) * rsd(4,i,j,k) )
114 flux(2,i,j,k) = rsd(2,i,j,k) * u21 + c2 *
115 > ( rsd(5,i,j,k) - q )
116 flux(3,i,j,k) = rsd(3,i,j,k) * u21
117 flux(4,i,j,k) = rsd(4,i,j,k) * u21
118 flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u21
127 frct(m,i,j,k) = frct(m,i,j,k)
128 > - tx2 * ( flux(m,i+1,j,k) - flux(m,i-1,j,k) )
132 tmp = 1.0d+00 / rsd(1,i,j,k)
134 u21i = tmp * rsd(2,i,j,k)
135 u31i = tmp * rsd(3,i,j,k)
136 u41i = tmp * rsd(4,i,j,k)
137 u51i = tmp * rsd(5,i,j,k)
139 tmp = 1.0d+00 / rsd(1,i-1,j,k)
141 u21im1 = tmp * rsd(2,i-1,j,k)
142 u31im1 = tmp * rsd(3,i-1,j,k)
143 u41im1 = tmp * rsd(4,i-1,j,k)
144 u51im1 = tmp * rsd(5,i-1,j,k)
146 flux(2,i,j,k) = (4.0d+00/3.0d+00) * tx3 *
148 flux(3,i,j,k) = tx3 * ( u31i - u31im1 )
149 flux(4,i,j,k) = tx3 * ( u41i - u41im1 )
150 flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
151 > * tx3 * ( ( u21i **2 + u31i **2 + u41i **2 )
152 > - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
153 > + (1.0d+00/6.0d+00)
154 > * tx3 * ( u21i**2 - u21im1**2 )
155 > + c1 * c5 * tx3 * ( u51i - u51im1 )
159 frct(1,i,j,k) = frct(1,i,j,k)
160 > + dx1 * tx1 * ( rsd(1,i-1,j,k)
161 > - 2.0d+00 * rsd(1,i,j,k)
163 frct(2,i,j,k) = frct(2,i,j,k)
164 > + tx3 * c3 * c4 * ( flux(2,i+1,j,k) - flux(2,i,j,k) )
165 > + dx2 * tx1 * ( rsd(2,i-1,j,k)
166 > - 2.0d+00 * rsd(2,i,j,k)
168 frct(3,i,j,k) = frct(3,i,j,k)
169 > + tx3 * c3 * c4 * ( flux(3,i+1,j,k) - flux(3,i,j,k) )
170 > + dx3 * tx1 * ( rsd(3,i-1,j,k)
171 > - 2.0d+00 * rsd(3,i,j,k)
173 frct(4,i,j,k) = frct(4,i,j,k)
174 > + tx3 * c3 * c4 * ( flux(4,i+1,j,k) - flux(4,i,j,k) )
175 > + dx4 * tx1 * ( rsd(4,i-1,j,k)
176 > - 2.0d+00 * rsd(4,i,j,k)
178 frct(5,i,j,k) = frct(5,i,j,k)
179 > + tx3 * c3 * c4 * ( flux(5,i+1,j,k) - flux(5,i,j,k) )
180 > + dx5 * tx1 * ( rsd(5,i-1,j,k)
181 > - 2.0d+00 * rsd(5,i,j,k)
185 c---------------------------------------------------------------------
186 c Fourth-order dissipation
187 c---------------------------------------------------------------------
188 IF (north.eq.-1) then
190 frct(m,2,j,k) = frct(m,2,j,k)
191 > - dsspm * ( + 5.0d+00 * rsd(m,2,j,k)
192 > - 4.0d+00 * rsd(m,3,j,k)
194 frct(m,3,j,k) = frct(m,3,j,k)
195 > - dsspm * ( - 4.0d+00 * rsd(m,2,j,k)
196 > + 6.0d+00 * rsd(m,3,j,k)
197 > - 4.0d+00 * rsd(m,4,j,k)
204 frct(m,i,j,k) = frct(m,i,j,k)
205 > - dsspm * ( rsd(m,i-2,j,k)
206 > - 4.0d+00 * rsd(m,i-1,j,k)
207 > + 6.0d+00 * rsd(m,i,j,k)
208 > - 4.0d+00 * rsd(m,i+1,j,k)
213 IF (south.eq.-1) then
215 frct(m,nx-2,j,k) = frct(m,nx-2,j,k)
216 > - dsspm * ( rsd(m,nx-4,j,k)
217 > - 4.0d+00 * rsd(m,nx-3,j,k)
218 > + 6.0d+00 * rsd(m,nx-2,j,k)
219 > - 4.0d+00 * rsd(m,nx-1,j,k) )
220 frct(m,nx-1,j,k) = frct(m,nx-1,j,k)
221 > - dsspm * ( rsd(m,nx-3,j,k)
222 > - 4.0d+00 * rsd(m,nx-2,j,k)
223 > + 5.0d+00 * rsd(m,nx-1,j,k) )
230 c---------------------------------------------------------------------
231 c eta-direction flux differences
232 c---------------------------------------------------------------------
234 c iex = flag : iex = 0 north/south communication
235 c : iex = 1 east/west communication
237 c---------------------------------------------------------------------
240 c---------------------------------------------------------------------
241 c communicate and receive/send two rows of data
242 c---------------------------------------------------------------------
243 call exchange_3 (rsd,iex)
246 if (west.eq.-1) L1 = 1
248 if (east.eq.-1) L2 = ny
252 if (west.eq.-1) jst1 = 4
253 if (east.eq.-1) jend1 = ny - 3
258 flux(1,i,j,k) = rsd(3,i,j,k)
259 u31 = rsd(3,i,j,k) / rsd(1,i,j,k)
260 q = 0.50d+00 * ( rsd(2,i,j,k) * rsd(2,i,j,k)
261 > + rsd(3,i,j,k) * rsd(3,i,j,k)
262 > + rsd(4,i,j,k) * rsd(4,i,j,k) )
264 flux(2,i,j,k) = rsd(2,i,j,k) * u31
265 flux(3,i,j,k) = rsd(3,i,j,k) * u31 + c2 *
266 > ( rsd(5,i,j,k) - q )
267 flux(4,i,j,k) = rsd(4,i,j,k) * u31
268 flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u31
277 frct(m,i,j,k) = frct(m,i,j,k)
278 > - ty2 * ( flux(m,i,j+1,k) - flux(m,i,j-1,k) )
285 tmp = 1.0d+00 / rsd(1,i,j,k)
287 u21j = tmp * rsd(2,i,j,k)
288 u31j = tmp * rsd(3,i,j,k)
289 u41j = tmp * rsd(4,i,j,k)
290 u51j = tmp * rsd(5,i,j,k)
292 tmp = 1.0d+00 / rsd(1,i,j-1,k)
294 u21jm1 = tmp * rsd(2,i,j-1,k)
295 u31jm1 = tmp * rsd(3,i,j-1,k)
296 u41jm1 = tmp * rsd(4,i,j-1,k)
297 u51jm1 = tmp * rsd(5,i,j-1,k)
299 flux(2,i,j,k) = ty3 * ( u21j - u21jm1 )
300 flux(3,i,j,k) = (4.0d+00/3.0d+00) * ty3 *
302 flux(4,i,j,k) = ty3 * ( u41j - u41jm1 )
303 flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
304 > * ty3 * ( ( u21j **2 + u31j **2 + u41j **2 )
305 > - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
306 > + (1.0d+00/6.0d+00)
307 > * ty3 * ( u31j**2 - u31jm1**2 )
308 > + c1 * c5 * ty3 * ( u51j - u51jm1 )
314 frct(1,i,j,k) = frct(1,i,j,k)
315 > + dy1 * ty1 * ( rsd(1,i,j-1,k)
316 > - 2.0d+00 * rsd(1,i,j,k)
318 frct(2,i,j,k) = frct(2,i,j,k)
319 > + ty3 * c3 * c4 * ( flux(2,i,j+1,k) - flux(2,i,j,k) )
320 > + dy2 * ty1 * ( rsd(2,i,j-1,k)
321 > - 2.0d+00 * rsd(2,i,j,k)
323 frct(3,i,j,k) = frct(3,i,j,k)
324 > + ty3 * c3 * c4 * ( flux(3,i,j+1,k) - flux(3,i,j,k) )
325 > + dy3 * ty1 * ( rsd(3,i,j-1,k)
326 > - 2.0d+00 * rsd(3,i,j,k)
328 frct(4,i,j,k) = frct(4,i,j,k)
329 > + ty3 * c3 * c4 * ( flux(4,i,j+1,k) - flux(4,i,j,k) )
330 > + dy4 * ty1 * ( rsd(4,i,j-1,k)
331 > - 2.0d+00 * rsd(4,i,j,k)
333 frct(5,i,j,k) = frct(5,i,j,k)
334 > + ty3 * c3 * c4 * ( flux(5,i,j+1,k) - flux(5,i,j,k) )
335 > + dy5 * ty1 * ( rsd(5,i,j-1,k)
336 > - 2.0d+00 * rsd(5,i,j,k)
341 c---------------------------------------------------------------------
342 c fourth-order dissipation
343 c---------------------------------------------------------------------
347 frct(m,i,2,k) = frct(m,i,2,k)
348 > - dsspm * ( + 5.0d+00 * rsd(m,i,2,k)
349 > - 4.0d+00 * rsd(m,i,3,k)
351 frct(m,i,3,k) = frct(m,i,3,k)
352 > - dsspm * ( - 4.0d+00 * rsd(m,i,2,k)
353 > + 6.0d+00 * rsd(m,i,3,k)
354 > - 4.0d+00 * rsd(m,i,4,k)
363 frct(m,i,j,k) = frct(m,i,j,k)
364 > - dsspm * ( rsd(m,i,j-2,k)
365 > - 4.0d+00 * rsd(m,i,j-1,k)
366 > + 6.0d+00 * rsd(m,i,j,k)
367 > - 4.0d+00 * rsd(m,i,j+1,k)
376 frct(m,i,ny-2,k) = frct(m,i,ny-2,k)
377 > - dsspm * ( rsd(m,i,ny-4,k)
378 > - 4.0d+00 * rsd(m,i,ny-3,k)
379 > + 6.0d+00 * rsd(m,i,ny-2,k)
380 > - 4.0d+00 * rsd(m,i,ny-1,k) )
381 frct(m,i,ny-1,k) = frct(m,i,ny-1,k)
382 > - dsspm * ( rsd(m,i,ny-3,k)
383 > - 4.0d+00 * rsd(m,i,ny-2,k)
384 > + 5.0d+00 * rsd(m,i,ny-1,k) )
391 c---------------------------------------------------------------------
392 c zeta-direction flux differences
393 c---------------------------------------------------------------------
397 flux(1,i,j,k) = rsd(4,i,j,k)
398 u41 = rsd(4,i,j,k) / rsd(1,i,j,k)
399 q = 0.50d+00 * ( rsd(2,i,j,k) * rsd(2,i,j,k)
400 > + rsd(3,i,j,k) * rsd(3,i,j,k)
401 > + rsd(4,i,j,k) * rsd(4,i,j,k) )
403 flux(2,i,j,k) = rsd(2,i,j,k) * u41
404 flux(3,i,j,k) = rsd(3,i,j,k) * u41
405 flux(4,i,j,k) = rsd(4,i,j,k) * u41 + c2 *
406 > ( rsd(5,i,j,k) - q )
407 flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u41
416 frct(m,i,j,k) = frct(m,i,j,k)
417 > - tz2 * ( flux(m,i,j,k+1) - flux(m,i,j,k-1) )
426 tmp = 1.0d+00 / rsd(1,i,j,k)
428 u21k = tmp * rsd(2,i,j,k)
429 u31k = tmp * rsd(3,i,j,k)
430 u41k = tmp * rsd(4,i,j,k)
431 u51k = tmp * rsd(5,i,j,k)
433 tmp = 1.0d+00 / rsd(1,i,j,k-1)
435 u21km1 = tmp * rsd(2,i,j,k-1)
436 u31km1 = tmp * rsd(3,i,j,k-1)
437 u41km1 = tmp * rsd(4,i,j,k-1)
438 u51km1 = tmp * rsd(5,i,j,k-1)
440 flux(2,i,j,k) = tz3 * ( u21k - u21km1 )
441 flux(3,i,j,k) = tz3 * ( u31k - u31km1 )
442 flux(4,i,j,k) = (4.0d+00/3.0d+00) * tz3 * ( u41k
444 flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
445 > * tz3 * ( ( u21k **2 + u31k **2 + u41k **2 )
446 > - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
447 > + (1.0d+00/6.0d+00)
448 > * tz3 * ( u41k**2 - u41km1**2 )
449 > + c1 * c5 * tz3 * ( u51k - u51km1 )
457 frct(1,i,j,k) = frct(1,i,j,k)
458 > + dz1 * tz1 * ( rsd(1,i,j,k+1)
459 > - 2.0d+00 * rsd(1,i,j,k)
461 frct(2,i,j,k) = frct(2,i,j,k)
462 > + tz3 * c3 * c4 * ( flux(2,i,j,k+1) - flux(2,i,j,k) )
463 > + dz2 * tz1 * ( rsd(2,i,j,k+1)
464 > - 2.0d+00 * rsd(2,i,j,k)
466 frct(3,i,j,k) = frct(3,i,j,k)
467 > + tz3 * c3 * c4 * ( flux(3,i,j,k+1) - flux(3,i,j,k) )
468 > + dz3 * tz1 * ( rsd(3,i,j,k+1)
469 > - 2.0d+00 * rsd(3,i,j,k)
471 frct(4,i,j,k) = frct(4,i,j,k)
472 > + tz3 * c3 * c4 * ( flux(4,i,j,k+1) - flux(4,i,j,k) )
473 > + dz4 * tz1 * ( rsd(4,i,j,k+1)
474 > - 2.0d+00 * rsd(4,i,j,k)
476 frct(5,i,j,k) = frct(5,i,j,k)
477 > + tz3 * c3 * c4 * ( flux(5,i,j,k+1) - flux(5,i,j,k) )
478 > + dz5 * tz1 * ( rsd(5,i,j,k+1)
479 > - 2.0d+00 * rsd(5,i,j,k)
485 c---------------------------------------------------------------------
486 c fourth-order dissipation
487 c---------------------------------------------------------------------
491 frct(m,i,j,2) = frct(m,i,j,2)
492 > - dsspm * ( + 5.0d+00 * rsd(m,i,j,2)
493 > - 4.0d+00 * rsd(m,i,j,3)
495 frct(m,i,j,3) = frct(m,i,j,3)
496 > - dsspm * (- 4.0d+00 * rsd(m,i,j,2)
497 > + 6.0d+00 * rsd(m,i,j,3)
498 > - 4.0d+00 * rsd(m,i,j,4)
508 frct(m,i,j,k) = frct(m,i,j,k)
509 > - dsspm * ( rsd(m,i,j,k-2)
510 > - 4.0d+00 * rsd(m,i,j,k-1)
511 > + 6.0d+00 * rsd(m,i,j,k)
512 > - 4.0d+00 * rsd(m,i,j,k+1)
522 frct(m,i,j,nz-2) = frct(m,i,j,nz-2)
523 > - dsspm * ( rsd(m,i,j,nz-4)
524 > - 4.0d+00 * rsd(m,i,j,nz-3)
525 > + 6.0d+00 * rsd(m,i,j,nz-2)
526 > - 4.0d+00 * rsd(m,i,j,nz-1) )
527 frct(m,i,j,nz-1) = frct(m,i,j,nz-1)
528 > - dsspm * ( rsd(m,i,j,nz-3)
529 > - 4.0d+00 * rsd(m,i,j,nz-2)
530 > + 5.0d+00 * rsd(m,i,j,nz-1) )