2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
10 c---------------------------------------------------------------------
11 c compute the right hand side based on exact solution
12 c---------------------------------------------------------------------
16 double precision dtemp(5), xi, eta, zeta, dtpp
17 integer c, m, i, j, k, ip1, im1, jp1,
21 c---------------------------------------------------------------------
22 c loop over all cells owned by this node
23 c---------------------------------------------------------------------
26 c---------------------------------------------------------------------
28 c---------------------------------------------------------------------
29 do k= 0, cell_size(3,c)-1
30 do j = 0, cell_size(2,c)-1
31 do i = 0, cell_size(1,c)-1
33 forcing(m,i,j,k,c) = 0.0d0
39 c---------------------------------------------------------------------
40 c xi-direction flux differences
41 c---------------------------------------------------------------------
42 do k = start(3,c), cell_size(3,c)-end(3,c)-1
43 zeta = dble(k+cell_low(3,c)) * dnzm1
44 do j = start(2,c), cell_size(2,c)-end(2,c)-1
45 eta = dble(j+cell_low(2,c)) * dnym1
47 do i=-2*(1-start(1,c)), cell_size(1,c)+1-2*end(1,c)
48 xi = dble(i+cell_low(1,c)) * dnxm1
50 call exact_solution(xi, eta, zeta, dtemp)
55 dtpp = 1.0d0 / dtemp(1)
58 buf(i,m) = dtpp * dtemp(m)
61 cuf(i) = buf(i,2) * buf(i,2)
62 buf(i,1) = cuf(i) + buf(i,3) * buf(i,3) +
64 q(i) = 0.5d0*(buf(i,2)*ue(i,2) + buf(i,3)*ue(i,3) +
69 do i = start(1,c), cell_size(1,c)-end(1,c)-1
73 forcing(1,i,j,k,c) = forcing(1,i,j,k,c) -
74 > tx2*( ue(ip1,2)-ue(im1,2) )+
75 > dx1tx1*(ue(ip1,1)-2.0d0*ue(i,1)+ue(im1,1))
77 forcing(2,i,j,k,c) = forcing(2,i,j,k,c) - tx2 * (
78 > (ue(ip1,2)*buf(ip1,2)+c2*(ue(ip1,5)-q(ip1)))-
79 > (ue(im1,2)*buf(im1,2)+c2*(ue(im1,5)-q(im1))))+
80 > xxcon1*(buf(ip1,2)-2.0d0*buf(i,2)+buf(im1,2))+
81 > dx2tx1*( ue(ip1,2)-2.0d0* ue(i,2)+ue(im1,2))
83 forcing(3,i,j,k,c) = forcing(3,i,j,k,c) - tx2 * (
84 > ue(ip1,3)*buf(ip1,2)-ue(im1,3)*buf(im1,2))+
85 > xxcon2*(buf(ip1,3)-2.0d0*buf(i,3)+buf(im1,3))+
86 > dx3tx1*( ue(ip1,3)-2.0d0*ue(i,3) +ue(im1,3))
88 forcing(4,i,j,k,c) = forcing(4,i,j,k,c) - tx2*(
89 > ue(ip1,4)*buf(ip1,2)-ue(im1,4)*buf(im1,2))+
90 > xxcon2*(buf(ip1,4)-2.0d0*buf(i,4)+buf(im1,4))+
91 > dx4tx1*( ue(ip1,4)-2.0d0* ue(i,4)+ ue(im1,4))
93 forcing(5,i,j,k,c) = forcing(5,i,j,k,c) - tx2*(
94 > buf(ip1,2)*(c1*ue(ip1,5)-c2*q(ip1))-
95 > buf(im1,2)*(c1*ue(im1,5)-c2*q(im1)))+
96 > 0.5d0*xxcon3*(buf(ip1,1)-2.0d0*buf(i,1)+
98 > xxcon4*(cuf(ip1)-2.0d0*cuf(i)+cuf(im1))+
99 > xxcon5*(buf(ip1,5)-2.0d0*buf(i,5)+buf(im1,5))+
100 > dx5tx1*( ue(ip1,5)-2.0d0* ue(i,5)+ ue(im1,5))
103 c---------------------------------------------------------------------
104 c Fourth-order dissipation
105 c---------------------------------------------------------------------
106 if (start(1,c) .gt. 0) then
109 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
110 > (5.0d0*ue(i,m) - 4.0d0*ue(i+1,m) +ue(i+2,m))
112 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
113 > (-4.0d0*ue(i-1,m) + 6.0d0*ue(i,m) -
114 > 4.0d0*ue(i+1,m) + ue(i+2,m))
118 do i = start(1,c)*3, cell_size(1,c)-3*end(1,c)-1
120 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp*
121 > (ue(i-2,m) - 4.0d0*ue(i-1,m) +
122 > 6.0d0*ue(i,m) - 4.0d0*ue(i+1,m) + ue(i+2,m))
126 if (end(1,c) .gt. 0) then
129 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
130 > (ue(i-2,m) - 4.0d0*ue(i-1,m) +
131 > 6.0d0*ue(i,m) - 4.0d0*ue(i+1,m))
133 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
134 > (ue(i-2,m) - 4.0d0*ue(i-1,m) + 5.0d0*ue(i,m))
141 c---------------------------------------------------------------------
142 c eta-direction flux differences
143 c---------------------------------------------------------------------
144 do k = start(3,c), cell_size(3,c)-end(3,c)-1
145 zeta = dble(k+cell_low(3,c)) * dnzm1
146 do i=start(1,c), cell_size(1,c)-end(1,c)-1
147 xi = dble(i+cell_low(1,c)) * dnxm1
149 do j=-2*(1-start(2,c)), cell_size(2,c)+1-2*end(2,c)
150 eta = dble(j+cell_low(2,c)) * dnym1
152 call exact_solution(xi, eta, zeta, dtemp)
157 dtpp = 1.0d0/dtemp(1)
160 buf(j,m) = dtpp * dtemp(m)
163 cuf(j) = buf(j,3) * buf(j,3)
164 buf(j,1) = cuf(j) + buf(j,2) * buf(j,2) +
165 > buf(j,4) * buf(j,4)
166 q(j) = 0.5d0*(buf(j,2)*ue(j,2) + buf(j,3)*ue(j,3) +
170 do j = start(2,c), cell_size(2,c)-end(2,c)-1
174 forcing(1,i,j,k,c) = forcing(1,i,j,k,c) -
175 > ty2*( ue(jp1,3)-ue(jm1,3) )+
176 > dy1ty1*(ue(jp1,1)-2.0d0*ue(j,1)+ue(jm1,1))
178 forcing(2,i,j,k,c) = forcing(2,i,j,k,c) - ty2*(
179 > ue(jp1,2)*buf(jp1,3)-ue(jm1,2)*buf(jm1,3))+
180 > yycon2*(buf(jp1,2)-2.0d0*buf(j,2)+buf(jm1,2))+
181 > dy2ty1*( ue(jp1,2)-2.0* ue(j,2)+ ue(jm1,2))
183 forcing(3,i,j,k,c) = forcing(3,i,j,k,c) - ty2*(
184 > (ue(jp1,3)*buf(jp1,3)+c2*(ue(jp1,5)-q(jp1)))-
185 > (ue(jm1,3)*buf(jm1,3)+c2*(ue(jm1,5)-q(jm1))))+
186 > yycon1*(buf(jp1,3)-2.0d0*buf(j,3)+buf(jm1,3))+
187 > dy3ty1*( ue(jp1,3)-2.0d0*ue(j,3) +ue(jm1,3))
189 forcing(4,i,j,k,c) = forcing(4,i,j,k,c) - ty2*(
190 > ue(jp1,4)*buf(jp1,3)-ue(jm1,4)*buf(jm1,3))+
191 > yycon2*(buf(jp1,4)-2.0d0*buf(j,4)+buf(jm1,4))+
192 > dy4ty1*( ue(jp1,4)-2.0d0*ue(j,4)+ ue(jm1,4))
194 forcing(5,i,j,k,c) = forcing(5,i,j,k,c) - ty2*(
195 > buf(jp1,3)*(c1*ue(jp1,5)-c2*q(jp1))-
196 > buf(jm1,3)*(c1*ue(jm1,5)-c2*q(jm1)))+
197 > 0.5d0*yycon3*(buf(jp1,1)-2.0d0*buf(j,1)+
199 > yycon4*(cuf(jp1)-2.0d0*cuf(j)+cuf(jm1))+
200 > yycon5*(buf(jp1,5)-2.0d0*buf(j,5)+buf(jm1,5))+
201 > dy5ty1*(ue(jp1,5)-2.0d0*ue(j,5)+ue(jm1,5))
204 c---------------------------------------------------------------------
205 c Fourth-order dissipation
206 c---------------------------------------------------------------------
207 if (start(2,c) .gt. 0) then
210 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
211 > (5.0d0*ue(j,m) - 4.0d0*ue(j+1,m) +ue(j+2,m))
213 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
214 > (-4.0d0*ue(j-1,m) + 6.0d0*ue(j,m) -
215 > 4.0d0*ue(j+1,m) + ue(j+2,m))
219 do j = start(2,c)*3, cell_size(2,c)-3*end(2,c)-1
221 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp*
222 > (ue(j-2,m) - 4.0d0*ue(j-1,m) +
223 > 6.0d0*ue(j,m) - 4.0d0*ue(j+1,m) + ue(j+2,m))
227 if (end(2,c) .gt. 0) then
230 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
231 > (ue(j-2,m) - 4.0d0*ue(j-1,m) +
232 > 6.0d0*ue(j,m) - 4.0d0*ue(j+1,m))
234 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
235 > (ue(j-2,m) - 4.0d0*ue(j-1,m) + 5.0d0*ue(j,m))
243 c---------------------------------------------------------------------
244 c zeta-direction flux differences
245 c---------------------------------------------------------------------
246 do j=start(2,c), cell_size(2,c)-end(2,c)-1
247 eta = dble(j+cell_low(2,c)) * dnym1
248 do i = start(1,c), cell_size(1,c)-end(1,c)-1
249 xi = dble(i+cell_low(1,c)) * dnxm1
251 do k=-2*(1-start(3,c)), cell_size(3,c)+1-2*end(3,c)
252 zeta = dble(k+cell_low(3,c)) * dnzm1
254 call exact_solution(xi, eta, zeta, dtemp)
259 dtpp = 1.0d0/dtemp(1)
262 buf(k,m) = dtpp * dtemp(m)
265 cuf(k) = buf(k,4) * buf(k,4)
266 buf(k,1) = cuf(k) + buf(k,2) * buf(k,2) +
267 > buf(k,3) * buf(k,3)
268 q(k) = 0.5d0*(buf(k,2)*ue(k,2) + buf(k,3)*ue(k,3) +
272 do k=start(3,c), cell_size(3,c)-end(3,c)-1
276 forcing(1,i,j,k,c) = forcing(1,i,j,k,c) -
277 > tz2*( ue(kp1,4)-ue(km1,4) )+
278 > dz1tz1*(ue(kp1,1)-2.0d0*ue(k,1)+ue(km1,1))
280 forcing(2,i,j,k,c) = forcing(2,i,j,k,c) - tz2 * (
281 > ue(kp1,2)*buf(kp1,4)-ue(km1,2)*buf(km1,4))+
282 > zzcon2*(buf(kp1,2)-2.0d0*buf(k,2)+buf(km1,2))+
283 > dz2tz1*( ue(kp1,2)-2.0d0* ue(k,2)+ ue(km1,2))
285 forcing(3,i,j,k,c) = forcing(3,i,j,k,c) - tz2 * (
286 > ue(kp1,3)*buf(kp1,4)-ue(km1,3)*buf(km1,4))+
287 > zzcon2*(buf(kp1,3)-2.0d0*buf(k,3)+buf(km1,3))+
288 > dz3tz1*(ue(kp1,3)-2.0d0*ue(k,3)+ue(km1,3))
290 forcing(4,i,j,k,c) = forcing(4,i,j,k,c) - tz2 * (
291 > (ue(kp1,4)*buf(kp1,4)+c2*(ue(kp1,5)-q(kp1)))-
292 > (ue(km1,4)*buf(km1,4)+c2*(ue(km1,5)-q(km1))))+
293 > zzcon1*(buf(kp1,4)-2.0d0*buf(k,4)+buf(km1,4))+
294 > dz4tz1*( ue(kp1,4)-2.0d0*ue(k,4) +ue(km1,4))
296 forcing(5,i,j,k,c) = forcing(5,i,j,k,c) - tz2 * (
297 > buf(kp1,4)*(c1*ue(kp1,5)-c2*q(kp1))-
298 > buf(km1,4)*(c1*ue(km1,5)-c2*q(km1)))+
299 > 0.5d0*zzcon3*(buf(kp1,1)-2.0d0*buf(k,1)
301 > zzcon4*(cuf(kp1)-2.0d0*cuf(k)+cuf(km1))+
302 > zzcon5*(buf(kp1,5)-2.0d0*buf(k,5)+buf(km1,5))+
303 > dz5tz1*( ue(kp1,5)-2.0d0*ue(k,5)+ ue(km1,5))
306 c---------------------------------------------------------------------
307 c Fourth-order dissipation
308 c---------------------------------------------------------------------
309 if (start(3,c) .gt. 0) then
312 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
313 > (5.0d0*ue(k,m) - 4.0d0*ue(k+1,m) +ue(k+2,m))
315 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
316 > (-4.0d0*ue(k-1,m) + 6.0d0*ue(k,m) -
317 > 4.0d0*ue(k+1,m) + ue(k+2,m))
321 do k = start(3,c)*3, cell_size(3,c)-3*end(3,c)-1
323 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp*
324 > (ue(k-2,m) - 4.0d0*ue(k-1,m) +
325 > 6.0d0*ue(k,m) - 4.0d0*ue(k+1,m) + ue(k+2,m))
329 if (end(3,c) .gt. 0) then
332 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
333 > (ue(k-2,m) - 4.0d0*ue(k-1,m) +
334 > 6.0d0*ue(k,m) - 4.0d0*ue(k+1,m))
336 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp *
337 > (ue(k-2,m) - 4.0d0*ue(k-1,m) + 5.0d0*ue(k,m))
344 c---------------------------------------------------------------------
345 c now change the sign of the forcing function,
346 c---------------------------------------------------------------------
347 do k = start(3,c), cell_size(3,c)-end(3,c)-1
348 do j = start(2,c), cell_size(2,c)-end(2,c)-1
349 do i = start(1,c), cell_size(1,c)-end(1,c)-1
351 forcing(m,i,j,k,c) = -1.d0 * forcing(m,i,j,k,c)