1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
6 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
12 double precision 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---------------------------------------------------------------------
25 do k = -1, cell_size(3,c)
26 do j = -1, cell_size(2,c)
27 do i = -1, cell_size(1,c)
28 rho_inv = 1.0d0/u(1,i,j,k,c)
29 rho_i(i,j,k,c) = rho_inv
30 us(i,j,k,c) = u(2,i,j,k,c) * rho_inv
31 vs(i,j,k,c) = u(3,i,j,k,c) * rho_inv
32 ws(i,j,k,c) = u(4,i,j,k,c) * rho_inv
33 square(i,j,k,c) = 0.5d0* (
34 > u(2,i,j,k,c)*u(2,i,j,k,c) +
35 > u(3,i,j,k,c)*u(3,i,j,k,c) +
36 > u(4,i,j,k,c)*u(4,i,j,k,c) ) * rho_inv
37 qs(i,j,k,c) = square(i,j,k,c) * rho_inv
42 c---------------------------------------------------------------------
43 c copy the exact forcing term to the right hand side; because
44 c this forcing term is known, we can store it on the whole of every
45 c cell, including the boundary
46 c---------------------------------------------------------------------
48 do k = 0, cell_size(3,c)-1
49 do j = 0, cell_size(2,c)-1
50 do i = 0, cell_size(1,c)-1
52 rhs(m,i,j,k,c) = forcing(m,i,j,k,c)
59 c---------------------------------------------------------------------
60 c compute xi-direction fluxes
61 c---------------------------------------------------------------------
62 do k = start(3,c), cell_size(3,c)-end(3,c)-1
63 do j = start(2,c), cell_size(2,c)-end(2,c)-1
64 do i = start(1,c), cell_size(1,c)-end(1,c)-1
69 rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dx1tx1 *
70 > (u(1,i+1,j,k,c) - 2.0d0*u(1,i,j,k,c) +
72 > tx2 * (u(2,i+1,j,k,c) - u(2,i-1,j,k,c))
74 rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dx2tx1 *
75 > (u(2,i+1,j,k,c) - 2.0d0*u(2,i,j,k,c) +
77 > xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
78 > tx2 * (u(2,i+1,j,k,c)*up1 -
79 > u(2,i-1,j,k,c)*um1 +
80 > (u(5,i+1,j,k,c)- square(i+1,j,k,c)-
81 > u(5,i-1,j,k,c)+ square(i-1,j,k,c))*
84 rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dx3tx1 *
85 > (u(3,i+1,j,k,c) - 2.0d0*u(3,i,j,k,c) +
87 > xxcon2 * (vs(i+1,j,k,c) - 2.0d0*vs(i,j,k,c) +
89 > tx2 * (u(3,i+1,j,k,c)*up1 -
92 rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dx4tx1 *
93 > (u(4,i+1,j,k,c) - 2.0d0*u(4,i,j,k,c) +
95 > xxcon2 * (ws(i+1,j,k,c) - 2.0d0*ws(i,j,k,c) +
97 > tx2 * (u(4,i+1,j,k,c)*up1 -
100 rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dx5tx1 *
101 > (u(5,i+1,j,k,c) - 2.0d0*u(5,i,j,k,c) +
103 > xxcon3 * (qs(i+1,j,k,c) - 2.0d0*qs(i,j,k,c) +
105 > xxcon4 * (up1*up1 - 2.0d0*uijk*uijk +
107 > xxcon5 * (u(5,i+1,j,k,c)*rho_i(i+1,j,k,c) -
108 > 2.0d0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
109 > u(5,i-1,j,k,c)*rho_i(i-1,j,k,c)) -
110 > tx2 * ( (c1*u(5,i+1,j,k,c) -
111 > c2*square(i+1,j,k,c))*up1 -
112 > (c1*u(5,i-1,j,k,c) -
113 > c2*square(i-1,j,k,c))*um1 )
118 c---------------------------------------------------------------------
119 c add fourth order xi-direction dissipation
120 c---------------------------------------------------------------------
121 if (start(1,c) .gt. 0) then
122 do k = start(3,c), cell_size(3,c)-end(3,c)-1
123 do j = start(2,c), cell_size(2,c)-end(2,c)-1
126 rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp *
127 > ( 5.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i+1,j,k,c) +
133 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
134 > (-4.0d0*u(m,i-1,j,k,c) + 6.0d0*u(m,i,j,k,c) -
135 > 4.0d0*u(m,i+1,j,k,c) + u(m,i+2,j,k,c))
141 do k = start(3,c), cell_size(3,c)-end(3,c)-1
142 do j = start(2,c), cell_size(2,c)-end(2,c)-1
143 do i = 3*start(1,c),cell_size(1,c)-3*end(1,c)-1
145 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
146 > ( u(m,i-2,j,k,c) - 4.0d0*u(m,i-1,j,k,c) +
147 > 6.0*u(m,i,j,k,c) - 4.0d0*u(m,i+1,j,k,c) +
155 if (end(1,c) .gt. 0) then
156 do k = start(3,c), cell_size(3,c)-end(3,c)-1
157 do j = start(2,c), cell_size(2,c)-end(2,c)-1
160 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
161 > ( u(m,i-2,j,k,c) - 4.0d0*u(m,i-1,j,k,c) +
162 > 6.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i+1,j,k,c) )
167 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
168 > ( u(m,i-2,j,k,c) - 4.d0*u(m,i-1,j,k,c) +
169 > 5.d0*u(m,i,j,k,c) )
175 c---------------------------------------------------------------------
176 c compute eta-direction fluxes
177 c---------------------------------------------------------------------
178 do k = start(3,c), cell_size(3,c)-end(3,c)-1
179 do j = start(2,c), cell_size(2,c)-end(2,c)-1
180 do i = start(1,c), cell_size(1,c)-end(1,c)-1
184 rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dy1ty1 *
185 > (u(1,i,j+1,k,c) - 2.0d0*u(1,i,j,k,c) +
187 > ty2 * (u(3,i,j+1,k,c) - u(3,i,j-1,k,c))
188 rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dy2ty1 *
189 > (u(2,i,j+1,k,c) - 2.0d0*u(2,i,j,k,c) +
191 > yycon2 * (us(i,j+1,k,c) - 2.0d0*us(i,j,k,c) +
193 > ty2 * (u(2,i,j+1,k,c)*vp1 -
194 > u(2,i,j-1,k,c)*vm1)
195 rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dy3ty1 *
196 > (u(3,i,j+1,k,c) - 2.0d0*u(3,i,j,k,c) +
198 > yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
199 > ty2 * (u(3,i,j+1,k,c)*vp1 -
200 > u(3,i,j-1,k,c)*vm1 +
201 > (u(5,i,j+1,k,c) - square(i,j+1,k,c) -
202 > u(5,i,j-1,k,c) + square(i,j-1,k,c))
204 rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dy4ty1 *
205 > (u(4,i,j+1,k,c) - 2.0d0*u(4,i,j,k,c) +
207 > yycon2 * (ws(i,j+1,k,c) - 2.0d0*ws(i,j,k,c) +
209 > ty2 * (u(4,i,j+1,k,c)*vp1 -
210 > u(4,i,j-1,k,c)*vm1)
211 rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dy5ty1 *
212 > (u(5,i,j+1,k,c) - 2.0d0*u(5,i,j,k,c) +
214 > yycon3 * (qs(i,j+1,k,c) - 2.0d0*qs(i,j,k,c) +
216 > yycon4 * (vp1*vp1 - 2.0d0*vijk*vijk +
218 > yycon5 * (u(5,i,j+1,k,c)*rho_i(i,j+1,k,c) -
219 > 2.0d0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
220 > u(5,i,j-1,k,c)*rho_i(i,j-1,k,c)) -
221 > ty2 * ((c1*u(5,i,j+1,k,c) -
222 > c2*square(i,j+1,k,c)) * vp1 -
223 > (c1*u(5,i,j-1,k,c) -
224 > c2*square(i,j-1,k,c)) * vm1)
229 c---------------------------------------------------------------------
230 c add fourth order eta-direction dissipation
231 c---------------------------------------------------------------------
232 if (start(2,c) .gt. 0) then
233 do k = start(3,c), cell_size(3,c)-end(3,c)-1
235 do i = start(1,c), cell_size(1,c)-end(1,c)-1
237 rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp *
238 > ( 5.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j+1,k,c) +
244 do i = start(1,c), cell_size(1,c)-end(1,c)-1
246 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
247 > (-4.0d0*u(m,i,j-1,k,c) + 6.0d0*u(m,i,j,k,c) -
248 > 4.0d0*u(m,i,j+1,k,c) + u(m,i,j+2,k,c))
254 do k = start(3,c), cell_size(3,c)-end(3,c)-1
255 do j = 3*start(2,c), cell_size(2,c)-3*end(2,c)-1
256 do i = start(1,c),cell_size(1,c)-end(1,c)-1
258 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
259 > ( u(m,i,j-2,k,c) - 4.0d0*u(m,i,j-1,k,c) +
260 > 6.0*u(m,i,j,k,c) - 4.0d0*u(m,i,j+1,k,c) +
267 if (end(2,c) .gt. 0) then
268 do k = start(3,c), cell_size(3,c)-end(3,c)-1
270 do i = start(1,c), cell_size(1,c)-end(1,c)-1
272 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
273 > ( u(m,i,j-2,k,c) - 4.0d0*u(m,i,j-1,k,c) +
274 > 6.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j+1,k,c) )
279 do i = start(1,c), cell_size(1,c)-end(1,c)-1
281 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
282 > ( u(m,i,j-2,k,c) - 4.d0*u(m,i,j-1,k,c) +
283 > 5.d0*u(m,i,j,k,c) )
289 c---------------------------------------------------------------------
290 c compute zeta-direction fluxes
291 c---------------------------------------------------------------------
292 do k = start(3,c), cell_size(3,c)-end(3,c)-1
293 do j = start(2,c), cell_size(2,c)-end(2,c)-1
294 do i = start(1,c), cell_size(1,c)-end(1,c)-1
299 rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dz1tz1 *
300 > (u(1,i,j,k+1,c) - 2.0d0*u(1,i,j,k,c) +
302 > tz2 * (u(4,i,j,k+1,c) - u(4,i,j,k-1,c))
303 rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dz2tz1 *
304 > (u(2,i,j,k+1,c) - 2.0d0*u(2,i,j,k,c) +
306 > zzcon2 * (us(i,j,k+1,c) - 2.0d0*us(i,j,k,c) +
308 > tz2 * (u(2,i,j,k+1,c)*wp1 -
309 > u(2,i,j,k-1,c)*wm1)
310 rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dz3tz1 *
311 > (u(3,i,j,k+1,c) - 2.0d0*u(3,i,j,k,c) +
313 > zzcon2 * (vs(i,j,k+1,c) - 2.0d0*vs(i,j,k,c) +
315 > tz2 * (u(3,i,j,k+1,c)*wp1 -
316 > u(3,i,j,k-1,c)*wm1)
317 rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dz4tz1 *
318 > (u(4,i,j,k+1,c) - 2.0d0*u(4,i,j,k,c) +
320 > zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
321 > tz2 * (u(4,i,j,k+1,c)*wp1 -
322 > u(4,i,j,k-1,c)*wm1 +
323 > (u(5,i,j,k+1,c) - square(i,j,k+1,c) -
324 > u(5,i,j,k-1,c) + square(i,j,k-1,c))
326 rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dz5tz1 *
327 > (u(5,i,j,k+1,c) - 2.0d0*u(5,i,j,k,c) +
329 > zzcon3 * (qs(i,j,k+1,c) - 2.0d0*qs(i,j,k,c) +
331 > zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk +
333 > zzcon5 * (u(5,i,j,k+1,c)*rho_i(i,j,k+1,c) -
334 > 2.0d0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
335 > u(5,i,j,k-1,c)*rho_i(i,j,k-1,c)) -
336 > tz2 * ( (c1*u(5,i,j,k+1,c) -
337 > c2*square(i,j,k+1,c))*wp1 -
338 > (c1*u(5,i,j,k-1,c) -
339 > c2*square(i,j,k-1,c))*wm1)
344 c---------------------------------------------------------------------
345 c add fourth order zeta-direction dissipation
346 c---------------------------------------------------------------------
347 if (start(3,c) .gt. 0) then
349 do j = start(2,c), cell_size(2,c)-end(2,c)-1
350 do i = start(1,c), cell_size(1,c)-end(1,c)-1
352 rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp *
353 > ( 5.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j,k+1,c) +
360 do j = start(2,c), cell_size(2,c)-end(2,c)-1
361 do i = start(1,c), cell_size(1,c)-end(1,c)-1
363 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
364 > (-4.0d0*u(m,i,j,k-1,c) + 6.0d0*u(m,i,j,k,c) -
365 > 4.0d0*u(m,i,j,k+1,c) + u(m,i,j,k+2,c))
371 do k = 3*start(3,c), cell_size(3,c)-3*end(3,c)-1
372 do j = start(2,c), cell_size(2,c)-end(2,c)-1
373 do i = start(1,c),cell_size(1,c)-end(1,c)-1
375 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
376 > ( u(m,i,j,k-2,c) - 4.0d0*u(m,i,j,k-1,c) +
377 > 6.0*u(m,i,j,k,c) - 4.0d0*u(m,i,j,k+1,c) +
384 if (end(3,c) .gt. 0) then
386 do j = start(2,c), cell_size(2,c)-end(2,c)-1
387 do i = start(1,c), cell_size(1,c)-end(1,c)-1
389 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
390 > ( u(m,i,j,k-2,c) - 4.0d0*u(m,i,j,k-1,c) +
391 > 6.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j,k+1,c) )
397 do j = start(2,c), cell_size(2,c)-end(2,c)-1
398 do i = start(1,c), cell_size(1,c)-end(1,c)-1
400 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
401 > ( u(m,i,j,k-2,c) - 4.d0*u(m,i,j,k-1,c) +
402 > 5.d0*u(m,i,j,k,c) )
408 do k = start(3,c), cell_size(3,c)-end(3,c)-1
409 do j = start(2,c), cell_size(2,c)-end(2,c)-1
410 do i = start(1,c), cell_size(1,c)-end(1,c)-1
412 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) * dt