1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
6 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
12 double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
16 c---------------------------------------------------------------------
17 c loop over all cells owned by this node
18 c---------------------------------------------------------------------
21 c---------------------------------------------------------------------
22 c compute the reciprocal of density, and the kinetic energy,
23 c and the speed of sound.
24 c---------------------------------------------------------------------
26 do k = -1, cell_size(3,c)
27 do j = -1, cell_size(2,c)
28 do i = -1, cell_size(1,c)
29 rho_inv = 1.0d0/u(i,j,k,1,c)
30 rho_i(i,j,k,c) = rho_inv
31 us(i,j,k,c) = u(i,j,k,2,c) * rho_inv
32 vs(i,j,k,c) = u(i,j,k,3,c) * rho_inv
33 ws(i,j,k,c) = u(i,j,k,4,c) * rho_inv
34 square(i,j,k,c) = 0.5d0* (
35 > u(i,j,k,2,c)*u(i,j,k,2,c) +
36 > u(i,j,k,3,c)*u(i,j,k,3,c) +
37 > u(i,j,k,4,c)*u(i,j,k,4,c) ) * rho_inv
38 qs(i,j,k,c) = square(i,j,k,c) * rho_inv
39 c---------------------------------------------------------------------
40 c (don't need speed and ainx until the lhs computation)
41 c---------------------------------------------------------------------
42 aux = c1c2*rho_inv* (u(i,j,k,5,c) - square(i,j,k,c))
45 ainv(i,j,k,c) = 1.0d0/aux
50 c---------------------------------------------------------------------
51 c copy the exact forcing term to the right hand side; because
52 c this forcing term is known, we can store it on the whole of every
53 c cell, including the boundary
54 c---------------------------------------------------------------------
57 do k = 0, cell_size(3,c)-1
58 do j = 0, cell_size(2,c)-1
59 do i = 0, cell_size(1,c)-1
60 rhs(i,j,k,m,c) = forcing(i,j,k,m,c)
67 c---------------------------------------------------------------------
68 c compute xi-direction fluxes
69 c---------------------------------------------------------------------
70 do k = start(3,c), cell_size(3,c)-end(3,c)-1
71 do j = start(2,c), cell_size(2,c)-end(2,c)-1
72 do i = start(1,c), cell_size(1,c)-end(1,c)-1
77 rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dx1tx1 *
78 > (u(i+1,j,k,1,c) - 2.0d0*u(i,j,k,1,c) +
80 > tx2 * (u(i+1,j,k,2,c) - u(i-1,j,k,2,c))
82 rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dx2tx1 *
83 > (u(i+1,j,k,2,c) - 2.0d0*u(i,j,k,2,c) +
85 > xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
86 > tx2 * (u(i+1,j,k,2,c)*up1 -
87 > u(i-1,j,k,2,c)*um1 +
88 > (u(i+1,j,k,5,c)- square(i+1,j,k,c)-
89 > u(i-1,j,k,5,c)+ square(i-1,j,k,c))*
92 rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dx3tx1 *
93 > (u(i+1,j,k,3,c) - 2.0d0*u(i,j,k,3,c) +
95 > xxcon2 * (vs(i+1,j,k,c) - 2.0d0*vs(i,j,k,c) +
97 > tx2 * (u(i+1,j,k,3,c)*up1 -
100 rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dx4tx1 *
101 > (u(i+1,j,k,4,c) - 2.0d0*u(i,j,k,4,c) +
103 > xxcon2 * (ws(i+1,j,k,c) - 2.0d0*ws(i,j,k,c) +
105 > tx2 * (u(i+1,j,k,4,c)*up1 -
106 > u(i-1,j,k,4,c)*um1)
108 rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dx5tx1 *
109 > (u(i+1,j,k,5,c) - 2.0d0*u(i,j,k,5,c) +
111 > xxcon3 * (qs(i+1,j,k,c) - 2.0d0*qs(i,j,k,c) +
113 > xxcon4 * (up1*up1 - 2.0d0*uijk*uijk +
115 > xxcon5 * (u(i+1,j,k,5,c)*rho_i(i+1,j,k,c) -
116 > 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
117 > u(i-1,j,k,5,c)*rho_i(i-1,j,k,c)) -
118 > tx2 * ( (c1*u(i+1,j,k,5,c) -
119 > c2*square(i+1,j,k,c))*up1 -
120 > (c1*u(i-1,j,k,5,c) -
121 > c2*square(i-1,j,k,c))*um1 )
126 c---------------------------------------------------------------------
127 c add fourth order xi-direction dissipation
128 c---------------------------------------------------------------------
129 if (start(1,c) .gt. 0) then
132 do k = start(3,c), cell_size(3,c)-end(3,c)-1
133 do j = start(2,c), cell_size(2,c)-end(2,c)-1
134 rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
135 > ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) +
143 do k = start(3,c), cell_size(3,c)-end(3,c)-1
144 do j = start(2,c), cell_size(2,c)-end(2,c)-1
145 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
146 > (-4.0d0*u(i-1,j,k,m,c) + 6.0d0*u(i,j,k,m,c) -
147 > 4.0d0*u(i+1,j,k,m,c) + u(i+2,j,k,m,c))
154 do k = start(3,c), cell_size(3,c)-end(3,c)-1
155 do j = start(2,c), cell_size(2,c)-end(2,c)-1
156 do i = 3*start(1,c),cell_size(1,c)-3*end(1,c)-1
157 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
158 > ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) +
159 > 6.0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) +
167 if (end(1,c) .gt. 0) then
170 do k = start(3,c), cell_size(3,c)-end(3,c)-1
171 do j = start(2,c), cell_size(2,c)-end(2,c)-1
172 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
173 > ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) +
174 > 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) )
181 do k = start(3,c), cell_size(3,c)-end(3,c)-1
182 do j = start(2,c), cell_size(2,c)-end(2,c)-1
183 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
184 > ( u(i-2,j,k,m,c) - 4.d0*u(i-1,j,k,m,c) +
185 > 5.d0*u(i,j,k,m,c) )
191 c---------------------------------------------------------------------
192 c compute eta-direction fluxes
193 c---------------------------------------------------------------------
194 do k = start(3,c), cell_size(3,c)-end(3,c)-1
195 do j = start(2,c), cell_size(2,c)-end(2,c)-1
196 do i = start(1,c), cell_size(1,c)-end(1,c)-1
200 rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dy1ty1 *
201 > (u(i,j+1,k,1,c) - 2.0d0*u(i,j,k,1,c) +
203 > ty2 * (u(i,j+1,k,3,c) - u(i,j-1,k,3,c))
204 rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dy2ty1 *
205 > (u(i,j+1,k,2,c) - 2.0d0*u(i,j,k,2,c) +
207 > yycon2 * (us(i,j+1,k,c) - 2.0d0*us(i,j,k,c) +
209 > ty2 * (u(i,j+1,k,2,c)*vp1 -
210 > u(i,j-1,k,2,c)*vm1)
211 rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dy3ty1 *
212 > (u(i,j+1,k,3,c) - 2.0d0*u(i,j,k,3,c) +
214 > yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
215 > ty2 * (u(i,j+1,k,3,c)*vp1 -
216 > u(i,j-1,k,3,c)*vm1 +
217 > (u(i,j+1,k,5,c) - square(i,j+1,k,c) -
218 > u(i,j-1,k,5,c) + square(i,j-1,k,c))
220 rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dy4ty1 *
221 > (u(i,j+1,k,4,c) - 2.0d0*u(i,j,k,4,c) +
223 > yycon2 * (ws(i,j+1,k,c) - 2.0d0*ws(i,j,k,c) +
225 > ty2 * (u(i,j+1,k,4,c)*vp1 -
226 > u(i,j-1,k,4,c)*vm1)
227 rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dy5ty1 *
228 > (u(i,j+1,k,5,c) - 2.0d0*u(i,j,k,5,c) +
230 > yycon3 * (qs(i,j+1,k,c) - 2.0d0*qs(i,j,k,c) +
232 > yycon4 * (vp1*vp1 - 2.0d0*vijk*vijk +
234 > yycon5 * (u(i,j+1,k,5,c)*rho_i(i,j+1,k,c) -
235 > 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
236 > u(i,j-1,k,5,c)*rho_i(i,j-1,k,c)) -
237 > ty2 * ((c1*u(i,j+1,k,5,c) -
238 > c2*square(i,j+1,k,c)) * vp1 -
239 > (c1*u(i,j-1,k,5,c) -
240 > c2*square(i,j-1,k,c)) * vm1)
245 c---------------------------------------------------------------------
246 c add fourth order eta-direction dissipation
247 c---------------------------------------------------------------------
248 if (start(2,c) .gt. 0) then
251 do k = start(3,c), cell_size(3,c)-end(3,c)-1
252 do i = start(1,c), cell_size(1,c)-end(1,c)-1
253 rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
254 > ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) +
262 do k = start(3,c), cell_size(3,c)-end(3,c)-1
263 do i = start(1,c), cell_size(1,c)-end(1,c)-1
264 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
265 > (-4.0d0*u(i,j-1,k,m,c) + 6.0d0*u(i,j,k,m,c) -
266 > 4.0d0*u(i,j+1,k,m,c) + u(i,j+2,k,m,c))
273 do k = start(3,c), cell_size(3,c)-end(3,c)-1
274 do j = 3*start(2,c), cell_size(2,c)-3*end(2,c)-1
275 do i = start(1,c),cell_size(1,c)-end(1,c)-1
276 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
277 > ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) +
278 > 6.0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) +
285 if (end(2,c) .gt. 0) then
288 do k = start(3,c), cell_size(3,c)-end(3,c)-1
289 do i = start(1,c), cell_size(1,c)-end(1,c)-1
290 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
291 > ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) +
292 > 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) )
299 do k = start(3,c), cell_size(3,c)-end(3,c)-1
300 do i = start(1,c), cell_size(1,c)-end(1,c)-1
301 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
302 > ( u(i,j-2,k,m,c) - 4.d0*u(i,j-1,k,m,c) +
303 > 5.d0*u(i,j,k,m,c) )
310 c---------------------------------------------------------------------
311 c compute zeta-direction fluxes
312 c---------------------------------------------------------------------
313 do k = start(3,c), cell_size(3,c)-end(3,c)-1
314 do j = start(2,c), cell_size(2,c)-end(2,c)-1
315 do i = start(1,c), cell_size(1,c)-end(1,c)-1
320 rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dz1tz1 *
321 > (u(i,j,k+1,1,c) - 2.0d0*u(i,j,k,1,c) +
323 > tz2 * (u(i,j,k+1,4,c) - u(i,j,k-1,4,c))
324 rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dz2tz1 *
325 > (u(i,j,k+1,2,c) - 2.0d0*u(i,j,k,2,c) +
327 > zzcon2 * (us(i,j,k+1,c) - 2.0d0*us(i,j,k,c) +
329 > tz2 * (u(i,j,k+1,2,c)*wp1 -
330 > u(i,j,k-1,2,c)*wm1)
331 rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dz3tz1 *
332 > (u(i,j,k+1,3,c) - 2.0d0*u(i,j,k,3,c) +
334 > zzcon2 * (vs(i,j,k+1,c) - 2.0d0*vs(i,j,k,c) +
336 > tz2 * (u(i,j,k+1,3,c)*wp1 -
337 > u(i,j,k-1,3,c)*wm1)
338 rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dz4tz1 *
339 > (u(i,j,k+1,4,c) - 2.0d0*u(i,j,k,4,c) +
341 > zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
342 > tz2 * (u(i,j,k+1,4,c)*wp1 -
343 > u(i,j,k-1,4,c)*wm1 +
344 > (u(i,j,k+1,5,c) - square(i,j,k+1,c) -
345 > u(i,j,k-1,5,c) + square(i,j,k-1,c))
347 rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dz5tz1 *
348 > (u(i,j,k+1,5,c) - 2.0d0*u(i,j,k,5,c) +
350 > zzcon3 * (qs(i,j,k+1,c) - 2.0d0*qs(i,j,k,c) +
352 > zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk +
354 > zzcon5 * (u(i,j,k+1,5,c)*rho_i(i,j,k+1,c) -
355 > 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
356 > u(i,j,k-1,5,c)*rho_i(i,j,k-1,c)) -
357 > tz2 * ( (c1*u(i,j,k+1,5,c) -
358 > c2*square(i,j,k+1,c))*wp1 -
359 > (c1*u(i,j,k-1,5,c) -
360 > c2*square(i,j,k-1,c))*wm1)
365 c---------------------------------------------------------------------
366 c add fourth order zeta-direction dissipation
367 c---------------------------------------------------------------------
368 if (start(3,c) .gt. 0) then
371 do j = start(2,c), cell_size(2,c)-end(2,c)-1
372 do i = start(1,c), cell_size(1,c)-end(1,c)-1
373 rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
374 > ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) +
382 do j = start(2,c), cell_size(2,c)-end(2,c)-1
383 do i = start(1,c), cell_size(1,c)-end(1,c)-1
384 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
385 > (-4.0d0*u(i,j,k-1,m,c) + 6.0d0*u(i,j,k,m,c) -
386 > 4.0d0*u(i,j,k+1,m,c) + u(i,j,k+2,m,c))
393 do k = 3*start(3,c), cell_size(3,c)-3*end(3,c)-1
394 do j = start(2,c), cell_size(2,c)-end(2,c)-1
395 do i = start(1,c),cell_size(1,c)-end(1,c)-1
396 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
397 > ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) +
398 > 6.0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) +
405 if (end(3,c) .gt. 0) then
408 do j = start(2,c), cell_size(2,c)-end(2,c)-1
409 do i = start(1,c), cell_size(1,c)-end(1,c)-1
410 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
411 > ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) +
412 > 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) )
419 do j = start(2,c), cell_size(2,c)-end(2,c)-1
420 do i = start(1,c), cell_size(1,c)-end(1,c)-1
421 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
422 > ( u(i,j,k-2,m,c) - 4.d0*u(i,j,k-1,m,c) +
423 > 5.d0*u(i,j,k,m,c) )
430 do k = start(3,c), cell_size(3,c)-end(3,c)-1
431 do j = start(2,c), cell_size(2,c)-end(2,c)-1
432 do i = start(1,c), cell_size(1,c)-end(1,c)-1
433 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) * dt