52
53
54
56
57
58
59#include "implicit_f.inc"
60#include "comlock.inc"
61
62
63
64#include "mvsiz_p.inc"
65
66
67
68#include "param_c.inc"
69#include "com04_c.inc"
70#include "com06_c.inc"
71#include "com08_c.inc"
72#include "scr14_c.inc"
73#include "scr16_c.inc"
74#include "parit_c.inc"
75#include "inter22.inc"
76
77
78
79 INTEGER, INTENT(IN) :: NEL
80 INTEGER, INTENT(IN) :: NFT
81 INTEGER, INTENT(IN) :: MTN
82 INTEGER, INTENT(IN) :: ISMSTR
83 INTEGER, INTENT(IN) :: JLAG
84 INTEGER, INTENT(IN) :: JHBE
85 INTEGER MAT(*),(*),IPARTS(*), IPARG1(*)
87 . pm(npropm,nummat),geo(npropg,numgeo), rho(*),off(*),
88 . vx1(*),vx2(*),vx3(*),vx4(*),vx5(*),vx6(*),vx7(*),vx8(*),
89 . vy1(*),vy2(*),vy3(*),vy4(*),vy5(*),vy6(*),vy7(*),vy8(*),
90 . vz1(*),vz2(*),vz3(*),vz4(*),vz5(*),vz6(*),vz7(*),vz8(*),
91 . f11(*),f21(*),f31(*),f12(*),f22(*),f32(*),
92 . f13(*),f23(*),f33(*),f14(*),f24(*),f34(*),
93 . f15(*),f25(*),f35(*),f16(*),f26(*),f36(*),
94 . f17(*),f27(*),f37(*),f18(*),f28(*),f38(*),
95 . px1h1(*), px1h2(*), px1h3(*),
96 . px2h1(*), px2h2(*), px2h3(*),
97 . px3h1(*), px3h2(*), px3h3(*),
98 . px4h1(*), px4h2(*), px4h3(*),eani(*),partsav(npsav,*),
99 . vol(*),cxx(*),vis(*),vd2(*),deltax(*) ,offg(*),vol0(*)
100
101
102
103 INTEGER I, FLUID,MX, J, II, IC, JST(MVSIZ+1), MT,
105 . caq(mvsiz), fcl(mvsiz), fcq(mvsiz),ehou(mvsiz),
106 . hx1(mvsiz), hx2(mvsiz), hx3(mvsiz), hx4(mvsiz),
107 . hy1(mvsiz), hy2(mvsiz), hy3(mvsiz), hy4(mvsiz),
108 . hz1(mvsiz), hz2(mvsiz), hz3(mvsiz), hz4(mvsiz),
109 . g11(mvsiz),g21(mvsiz),g31(mvsiz),g41(mvsiz),
110 . g51(mvsiz),g61(mvsiz),g71(mvsiz),g81(mvsiz),
111 . g12(mvsiz),g22(mvsiz),g32(mvsiz),g42(mvsiz),
112 . g52(mvsiz),g62(mvsiz),g72(mvsiz),g82(mvsiz),
113 . g13(mvsiz),g23(mvsiz),g33(mvsiz),g43(mvsiz),
114 . g53(mvsiz),g63(mvsiz),g73(mvsiz),g83(mvsiz),
115 . len(mvsiz),rho0(mvsiz),
116 . ehourt,aa
117
118 REAL(kind=8) ::
119 . hgx1(mvsiz), hgx2(mvsiz), hgx3(mvsiz), hgx4(mvsiz),
120 . hgy1(mvsiz), hgy2(mvsiz), hgy3(mvsiz), hgy4(mvsiz),
121 . hgz1(mvsiz), hgz2(mvsiz), hgz3(mvsiz), hgz4(mvsiz)
122 REAL(kind=8) ::
123 . vx3478, vx2358, vx1467, vx1256,
124 . vy3478, vy2358, vy1467, vy1256,
125 . vz3478, vz2358, vz1467, vz1256
126
127
128 IF(iparg1(64)==1 .OR. (mtn==17 .AND.
ale%UPWIND%UPWM<2) .OR. int22>0 .OR. ifvm_skip==1)
THEN
129
130 DO i=1,nel
131
132 f11(i) =zero
133 f12(i) =zero
134 f13(i) =zero
135 f14(i) =zero
136 f15(i) =zero
137 f16(i) =zero
138 f17(i) =zero
139 f18(i) =zero
140
141 f21(i) =zero
142 f22(i) =zero
143 f23(i) =zero
144 f24(i) =zero
145 f25(i) =zero
146 f26(i) =zero
147 f27(i) =zero
148 f28(i) =zero
149
150 f31(i) =zero
151 f32(i) =zero
152 f33(i) =zero
153 f34(i) =zero
154 f35(i) =zero
155 f36(i) =zero
156 f37(i) =zero
157 f38(i) =zero
158 ENDDO
159
160 RETURN
161 ENDIF
162
163 IF(invstr>=35)THEN
164 mt = pid(1)
165 DO i=1,nel
166 caq(i)=fourth*off(i)*geo(13,mt)
167 ENDDO
168 ELSE
169 mx = mat(1)
170 DO i=1,nel
171 caq(i)=fourth*off(i)*pm(4,mx)
172 ENDDO
173 ENDIF
174
175 fluid=iparg1(63)
176
177 IF(fluid==1)THEN
178 IF(
ale%UPWIND%UPWM>1 .OR. jlag==1)
THEN
179 DO i=1,nel
180 IF(vis(i)>zero)THEN
181 fcq(i)=zero
182 fcl(i)=ten*caq(i)*vis(i)*deltax(i)
183 ELSE
184 fcq(i)=zero
185 fcl(i)=caq(i)*rho(i)*cxx(i)*deltax(i)**2
186 ENDIF
187 ENDDO
188 ELSEIF(
ale%UPWIND%UPWM==0)
THEN
189 DO i=1,nel
190 fcl(i)=caq(i)*rho(i)*cxx(i)*vol(i)**two_third
191 fcq(i)=zero
192 ENDDO
193 ELSEIF(
ale%UPWIND%UPWM==1)
THEN
194 DO i=1,nel
195 fcl(i)=caq(i)*rho(i)*deltax(i)**2
196 fcl(i)=
min(fcl(i)*cxx(i),
max(ten*caq(i)*vis(i)*deltax(i),fcl(i)*sqrt(vd2(i))))
197 fcq(i)=zero
198 ENDDO
199 ENDIF
200 ELSE
201
202 IF(ismstr == 1) THEN
203 mx = mat(1)
204 DO i=1,nel
205 rho0(i) = pm(1,mx)
206 END DO
207 DO i=1,nel
208 fcl(i)=caq(i)*rho0(i)*vol(i)**two_third
209 fcq(i)=fcl(i)*caq(i)*hundred
210 fcl(i)=fcl(i)*cxx(i)
211 ENDDO
212 ELSEIF(ismstr == 2) THEN
213 mx = mat(1)
214 DO i=1,nel
215 IF(offg(i) > one) THEN
216 rho0(i) = pm(1,mx)
217 END IF
218 END DO
219 DO i=1,nel
220 IF(offg(i) > one) THEN
221 aa = rho0(i)*vol0(i)/
max(em20,vol(i))
222 fcl(i)=caq(i)*aa*vol(i)**two_third
223 fcq(i)=fcl(i)*caq(i)*hundred
224 fcl(i)=fcl(i)*cxx(i)
225 ELSE
226 fcl(i)=caq(i)*rho(i)*vol(i)**two_third
227 fcq(i)=fcl(i)*caq(i)*hundred
228 fcl(i)=fcl(i)*cxx(i)
229 ENDIF
230 ENDDO
231 ELSE
232 DO i=1,nel
233 fcl(i)=caq(i)*rho(i)*vol(i)**two_third
234 fcq(i)=fcl(i)*caq(i)*hundred
235 fcl(i)=fcl(i)*cxx(i)
236 ENDDO
237 ENDIF
238 END IF
239 IF(jhbe==0)THEN
240 DO i=1,nel
241 vx3478=vx3(i)-vx4(i)-vx7(i)+vx8(i)
242 vx2358=vx2(i)-vx3(i)-vx5(i)+vx8(i)
243 vx1467=vx1(i)-vx4(i)-vx6(i)+vx7(i)
244 vx1256=vx1(i)-vx2(i)-vx5(i)+vx6(i)
245
246 vy3478=vy3(i)-vy4(i)-vy7(i)+vy8(i)
247 vy2358=vy2(i)-vy3(i)-vy5(i)+vy8(i)
248 vy1467=vy1(i)-vy4(i)-vy6(i)+vy7(i)
249 vy1256=vy1(i)-vy2(i)-vy5(i)+vy6(i)
250
251 vz3478=vz3(i)-vz4(i)-vz7(i)+vz8(i)
252 vz2358=vz2(i)-vz3(i)-vz5(i)+vz8(i)
253 vz1467=vz1(i)-vz4(i)-vz6(i)+vz7(i)
254 vz1256=vz1(i)-vz2(i)-vz5(i)+vz6(i)
255
256 hgx1(i)=vx1467-vx2358
257 hgx2(i)=vx1467+vx2358
258 hgx3(i)=vx1256-vx3478
259 hgx4(i)=vx1256+vx3478
260
261 hgy1(i)=vy1467-vy2358
262 hgy2(i)=vy1467+vy2358
263 hgy3(i)=vy1256-vy3478
264 hgy4(i)=vy1256+vy3478
265
266 hgz1(i)=vz1467-vz2358
267 hgz2(i)=vz1467+vz2358
268 hgz3(i)=vz1256-vz3478
269 hgz4(i)=vz1256+vz3478
270 ENDDO
271
272 DO i=1,nel
273 hx1(i)=hgx1(i)*(fcl(i)+abs(hgx1(i))*fcq(i))
274 hx2(i)=hgx2(i)*(fcl(i)+abs(hgx2(i))*fcq(i))
275 hx3(i)=hgx3(i)*(fcl(i)+abs(hgx3(i))*fcq(i))
276 hx4(i)=hgx4(i)*(fcl(i)+abs(hgx4(i))*fcq(i))
277
278 hy1(i)=hgy1(i)*(fcl(i)+abs(hgy1(i))*fcq(i))
279 hy2(i)=hgy2(i)*(fcl(i)+abs(hgy2(i))*fcq(i))
280 hy3(i)=hgy3(i)*(fcl(i)+abs(hgy3(i))*fcq(i))
281 hy4(i)=hgy4(i)*(fcl(i)+abs(hgy4(i))*fcq(i))
282
283 hz1(i)=hgz1(i)*(fcl(i)+abs(hgz1(i))*fcq(i))
284 hz2(i)=hgz2(i)*(fcl(i)+abs(hgz2(i))*fcq(i))
285 hz3(i)=hgz3(i)*(fcl(i)+abs(hgz3(i))*fcq(i))
286 hz4(i)=hgz4(i)*(fcl(i)+abs(hgz4(i))*fcq(i))
287 ENDDO
288
289 DO i=1,nel
290 f11(i) =-hx1(i)-hx2(i)-hx3(i)-hx4(i)
291 f12(i) = hx1(i)-hx2(i)+hx3(i)+hx4(i)
292 f13(i) =-hx1(i)+hx2(i)+hx3(i)-hx4(i)
293 f14(i) = hx1(i)+hx2(i)-hx3(i)+hx4(i)
294 f15(i) =-hx1(i)+hx2(i)+hx3(i)+hx4(i)
295 f16(i) = hx1(i)+hx2(i)-hx3(i)-hx4(i)
296 f17(i) =-hx1(i)-hx2(i)-hx3(i)+hx4(i)
297 f18(i) = hx1(i)-hx2(i)+hx3(i)-hx4(i)
298
299 f21(i) =-hy1(i)-hy2(i)-hy3(i)-hy4(i)
300 f22(i) = hy1(i)-hy2(i)+hy3(i)+hy4(i)
301 f23(i) =-hy1(i)+hy2(i)+hy3(i)-hy4(i)
302 f24(i) = hy1(i)+hy2(i)-hy3(i)+hy4(i)
303 f25(i) =-hy1(i)+hy2(i)+hy3(i)+hy4(i)
304 f26(i) = hy1(i)+hy2(i)-hy3(i)-hy4(i)
305 f27(i) =-hy1(i)-hy2(i)-hy3(i)+hy4(i)
306 f28(i) = hy1(i)-hy2(i)+hy3(i)-hy4(i)
307
308 f31(i) =-hz1(i)-hz2(i)-hz3(i)-hz4(i)
309 f32(i) = hz1(i)-hz2(i)+hz3(i)+hz4(i)
310 f33(i) =-hz1(i)+hz2(i)+hz3(i)-hz4(i)
311 f34(i) = hz1(i)+hz2(i)-hz3(i)+hz4(i)
312 f35(i) =-hz1(i)+hz2(i)+hz3(i)+hz4(i)
313 f36(i) = hz1(i)+hz2(i)-hz3(i)-hz4(i)
314 f37(i) =-hz1(i)-hz2(i)-hz3(i)+hz4(i)
315 f38(i) = hz1(i)-hz2(i)+hz3(i)-hz4(i)
316 ENDDO
317 ELSEIF(jhbe>=1) THEN
318 DO i=1,nel
319 g11(i)= one-px1h1(i)
320 g21(i)=-one-px2h1(i)
321 g31(i)= one-px3h1(i)
322 g41(i)=-one-px4h1(i)
323 g51(i)= one+px3h1(i)
324 g61(i)=-one+px4h1(i)
325 g71(i)= one+px1h1(i)
326 g81(i)=-one+px2h1(i)
327 hgx1(i) = g11(i)*vx1(i)+g21(i)*vx2(i)+g31(i)*vx3(i)+g41(i)*vx4(i)+g51(i)*vx5(i)+g61(i)*vx6(i)+g71(i)*vx7(i)+g81(i)*vx8(i)
328 hgy1(i) = g11(i)*vy1(i)+g21(i)*vy2(i)+g31(i)*vy3(i)+g41(i)*vy4(i)+g51(i)*vy5(i)+g61(i)*vy6(i)+g71(i)*vy7(i)+g81(i)*vy8(i)
329 hgz1(i) = g11(i)*vz1(i)+g21(i)*vz2(i)+g31(i)*vz3(i)+g41(i)*vz4(i)+g51(i)*vz5(i)+g61(i)*vz6(i)+g71(i)*vz7(i)+g81(i)*vz8(i)
330 ENDDO
331
332 DO i=1,nel
333 g12(i)= one-px1h2(i)
334 g22(i)= one-px2h2(i)
335 g32(i)=-one-px3h2(i)
336 g42(i)=-one-px4h2(i)
337 g52(i)=-one+px3h2(i)
338 g62(i)=-one+px4h2(i)
339 g72(i)= one+px1h2(i)
340 g82(i)= one+px2h2(i)
341 hgx2(i)=g12(i)*vx1(i)+g22(i)*vx2(i)+g32(i)*vx3(i)+g42(i)*vx4(i)+g52(i)*vx5(i)+g62(i)*vx6(i)+g72(i)*vx7(i)+g82(i)*vx8(i)
342 hgy2(i)=g12(i)*vy1(i)+g22(i)*vy2(i)+g32(i)*vy3(i)+g42(i)*vy4(i)+g52(i)*vy5(i)+g62(i)*vy6(i)+g72(i)*vy7(i)+g82(i)*vy8(i)
343 hgz2(i)=g12(i)*vz1(i)+g22(i)*vz2(i)+g32(i)*vz3(i)+g42(i)*vz4(i)+g52(i)*vz5(i)+g62(i)*vz6(i)+g72(i)*vz7(i)+g82(i)*vz8(i)
344 ENDDO
345
346 DO i=1,nel
347 g13(i)= one-px1h3(i)
348 g23(i)=-one-px2h3(i)
349 g33(i)=-one-px3h3(i)
350 g43(i)= one-px4h3(i)
351 g53(i)=-one+px3h3(i)
352 g63(i)= one+px4h3(i)
353 g73(i)= one+px1h3(i)
354 g83(i)=-one+px2h3(i)
355 hgx3(i)=g13(i)*vx1(i)+g23(i)*vx2(i)+g33(i)*vx3(i)+g43(i)*vx4(i
356 hgy3(i)=g13(i)*vy1(i)+g23(i)*vy2(i)+g33(i)*vy3(i)+g43(i)*vy4(i)+g53(i)*vy5(i)+g63(i)*vy6(i)+g73(i)*vy7(i)+g83(i)*vy8(i)
357 hgz3(i)=g13(i)*vz1(i)+g23(i)*vz2(i)+g33(i)*vz3(i)+g43(i)*vz4(i)+g53(i)*vz5(i)+g63(i)*vz6(i)+g73(i)*vz7(i)+g83(i)*vz8(i)
358 ENDDO
359
360 DO i=1,nel
361
362 hgx4(i)=vx1(i)-vx2(i)+vx3(i)-vx4(i)-vx5(i)+vx6(i)-vx7(i)+vx8(i)
363 hgy4(i)=vy1(i)-vy2(i)+vy3(i)-vy4(i)-vy5(i)+vy6(i)-vy7(i)+vy8(i)
364 hgz4(i)=vz1(i)-vz2(i)+vz3(i)-vz4(i)-vz5(i)+vz6(i)-vz7(i)+vz8(i)
365 ENDDO
366
367 DO i=1,nel
368 hx1(i)=hgx1(i)*(fcl(i)+abs(hgx1(i))*fcq(i))
369 hx2(i)=hgx2(i)*(fcl(i)+abs(hgx2(i))*fcq(i))
370 hx3(i)=hgx3(i)*(fcl(i)+abs(hgx3(i))*fcq(i))
371 hx4(i)=hgx4(i)*(fcl(i)+abs(hgx4(i))*fcq(i))
372
373 hy1(i)=hgy1(i)*(fcl(i)+abs(hgy1(i))*fcq(i))
374 hy2(i)=hgy2(i)*(fcl(i)+abs(hgy2(i))*fcq(i))
375 hy3(i)=hgy3(i)*(fcl(i)+abs(hgy3(i))*fcq(i))
376 hy4(i)=hgy4(i)*(fcl(i)+abs(hgy4(i))*fcq(i))
377
378 hz1(i)=hgz1(i)*(fcl(i)+abs(hgz1(i))*fcq(i
379 hz2(i)=hgz2(i)*(fcl(i)+abs(hgz2(i))*fcq(i))
380 hz3(i)=hgz3(i)*(fcl(i)+abs(hgz3
381 hz4(i)=hgz4(i)*(fcl(i)+abs(hgz4(i))*fcq(i))
382 ENDDO
383
384 DO i=1,nel
385 f11(i) =-g11(i)*hx1(i)-g12(i)*hx2(i)-g13(i)*hx3(i)-hx4(i
386 f12(i) =-g21(i)*hx1(i)-g22(i)*hx2(i)-g23(i)*hx3(i)+hx4(i)
387 f13(i) =-g31(i)*hx1(i)-g32(i)*hx2(i)-g33(i)*hx3(i)-hx4(i)
388 f14(i) =-g41(i)*hx1(i)-g42(i)*hx2(i)-g43(i)*hx3(i)+hx4(i)
389 f15(i) =-g51(i)*hx1(i)-g52(i)*hx2(i)-g53(i)*hx3(i)+hx4(i)
390 f16(i) =-g61(i)*hx1(i)-g62(i)*hx2(i)-g63(i)*hx3(i)-hx4(i)
391 f17(i) =-g71(i)*hx1(i)-g72(i)*hx2(i)-g73(i)*hx3(i)+hx4(i)
392 f18(i) =-g81(i)*hx1(i)-g82(i)*hx2(i)-g83(i)*hx3(i)-hx4(i)
393
394 f21(i) =-g11(i)*hy1(i)-g12(i)*hy2(i)-g13(i)*hy3(i)-hy4(i)
395 f22(i) =-g21(i)*hy1(i)-g22(i)*hy2(i)-g23(i)*hy3(i)+hy4(i)
396 f23(i) =-g31(i)*hy1(i)-g32(i)*hy2(i)-g33(i)*hy3(i)-hy4(i)
397 f24(i) =-g41(i)*hy1(i)-g42(i)*hy2(i)-g43(i)*hy3(i)+hy4(i)
398 f25(i) =-g51(i)*hy1(i)-g52(i)*hy2(i)-g53(i)*hy3(i)+hy4(i)
399 f26(i) =-g61(i)*hy1(i)-g62(i)*hy2(i)-g63(i)*hy3(i)-hy4(i)
400 f27(i) =-g71(i)*hy1(i)-g72(i)*hy2(i)-g73(i)*hy3(i)+hy4(i)
401 f28(i) =-g81(i)*hy1(i)-g82(i)*hy2(i)-g83(i)*hy3(i)-hy4(i)
402
403 f31(i) =-g11(i)*hz1(i)-g12(i)*hz2(i)-g13(i)*hz3(i)-hz4(i)
404 f32(i) =-g21(i)*hz1(i)-g22(i)*hz2(i)-g23(i)*hz3(i)+hz4(i)
405 f33(i) =-g31(i)*hz1(i)-g32(i)*hz2(i)-g33(i)*hz3(i)-hz4(i)
406 f34(i) =-g41(i)*hz1(i)-g42(i)*hz2(i)-g43(i)*hz3(i)+hz4(i)
407 f35(i) =-g51(i)*hz1(i)-g52(i)*hz2(i)-g53(i)*hz3(i)+hz4(i)
408 f36(i) =-g61(i)*hz1(i)-g62(i)*hz2(i)-g63(i)*hz3(i)-hz4(i)
409 f37(i) =-g71(i)*hz1(i)-g72(i)*hz2(i)-g73(i)*hz3(i)+hz4(i)
410 f38(i) =-g81(i)*hz1(i)-g82(i)*hz2(i)-g83(i)*hz3(i)-hz4(i)
411 ENDDO
412 ENDIF
413
414 DO i=1,nel
415 ehou(i)= dt1*(
416 & hz1(i)*hgz1(i) + hz2(i)*hgz2(i) +
417 & hz3(i)*hgz3(i) + hz4(i)*hgz4(i) +
418 & hx1(i)*hgx1(i) + hx2(i)*hgx2(i) +
419 & hx3(i)*hgx3(i) + hx4(i)*hgx4(i) +
420 & hy1(i)*hgy1(i) + hy2(i)*hgy2(i) +
421 & hy3(i)*hgy3(i) + hy4(i)*hgy4(i) )
422 END DO
423
424 ehourt = zero
425 DO i=1,nel
426 ehourt= ehourt+ehou(i)
427 ENDDO
428
429
430
431 IF(jlag==1)THEN
432
433 IF (ivector==0) THEN
434 DO i=1,nel
435 mx = iparts(i)
436 partsav(8,mx)=partsav(8,mx) + ehou(i)
437 ENDDO
438 ELSE
439 ic=1
440 jst(ic)=1
441 DO j=1+1,nel
442 IF (iparts(j)/=iparts(j-1)) THEN
443 ic=ic+1
444 jst(ic)=j
445 ENDIF
446 ENDDO
447 jst(ic+1)=nel+1
448 DO ii=1,ic
449 mx=iparts(jst(ii))
450 IF (jst(ii+1)-jst(ii)>15) THEN
451#include "vectorize.inc"
452 DO j=jst(ii),jst(ii+1)-1
453 partsav(8,mx)=partsav(8,mx) + ehou(j)
454 ENDDO
455 ELSE
456 DO j=jst(ii),jst(ii+1)-1
457 partsav(8,mx)=partsav(8,mx) + ehou(j)
458 ENDDO
459 ENDIF
460 ENDDO
461 ENDIF
462
463 ehour = ehour + ehourt
464 ENDIF
465
466#include "vectorize.inc"
467 DO i=1,nel
468 eani(nft+i) = eani(nft+i)+ehou(i)/
max(em30,rho(i)*vol(i))
469 ENDDO
470
471 RETURN