41
42
43
44
45
46
47
48
49
50
52
53
54
55#include "implicit_f.inc"
56
57
58
59#include "mvsiz_p.inc"
60
61
62
63#include "com04_c.inc"
64#include "vect01_c.inc"
65#include "param_c.inc"
66
67
68
70 . pm(npropm,nummat), w(3,numnod), rho(*),
71 . y1(*), y2(*), y3(*), y4(*), z1(*), z2(*), z3(*), z4(*),
72 . t11(*), t12(*), t13(*), t14(*), t21(*), t22(*), t23(*),t24(*),
73 . py1(*), py2(*), pz1(*), pz2(*), airem(*),
74 . vy1(*), vy2(*), vy3(*), vy4(*),
75 . vz1(*), vz2(*), vz3(*), vz4(*),
76 . dyy(*), dzz(*), dyz(*), dzy(*), off(*),qmv(8,*),bufmat(*), deltax(*), vis(*)
77 INTEGER MAT(*), NC1(*), NC2(*), NC3(*), NC4(*)
78 INTEGER,INTENT(IN) :: IPM(NPROPMI,NUMMAT)
79
80
81
82 INTEGER :: I,IADBUF,IFLG
83 my_real :: f1(mvsiz), f2(mvsiz),a1(mvsiz), a2(mvsiz), xms(mvsiz), gamma(mvsiz),vdy(mvsiz), vdz(mvsiz)
85 my_real :: fvy1(mvsiz), fvz1(mvsiz), fvy2(mvsiz), fvz2(mvsiz), fvy3(mvsiz), fvz3(mvsiz), fvy4(mvsiz), fvz4(mvsiz)
88 . y12(mvsiz), z12(mvsiz), y23(mvsiz), z23(mvsiz),
89 . y34(mvsiz), z34(mvsiz), y41(mvsiz), z41(mvsiz),
90 . y24(mvsiz), z24(mvsiz), y13(mvsiz), z13(mvsiz),
91 . dyyp(mvsiz), dzzp(mvsiz), dyzp(mvsiz), dzyp(mvsiz),
92 . f1p(mvsiz), f2p(mvsiz), a3(mvsiz),
93 . s(3,mvsiz), t(3,mvsiz), tr(mvsiz)
94
95
96
97
98
99
100
101 DO i=lft,llt
102
103 fvy1(i) = vy1(i)
104 vy1(i) = vy1(i) - w(2,nc1(i))
105 ! Save fluid velocity
106 fvz1(i) = vz1(i)
107 vz1(i) = vz1(i) - w(3,nc1(i))
108
109 fvy2(i) = vy2(i)
110 vy2(i) = vy2(i) - w(2,nc2(i))
111
112 fvz2(i) = vz2(i)
113 vz2(i) = vz2(i) - w(3,nc2(i))
114
115 fvy3(i) = vy3(i)
116 vy3(i) = vy3(i) - w(2,nc3(i))
117
118 fvz3(i) = vz3(i)
119 vz3(i) = vz3(i) - w(3,nc3(i))
120
121 fvy4(i) = vy4(i)
122 vy4(i) = vy4(i) - w(2,nc4(i))
123
124 fvz4(i) = vz4(i)
125 vz4(i) = vz4(i) - w(3,nc4(i))
126 ENDDO
127
128
129
130
131 DO i=lft,llt
132 xms(i) = fourth*rho(i)*airem(i)
133 gamma(i)= pm(15,mat(i))
134 ENDDO
135 DO i=lft,llt
136 vdy(i) = fourth*(vy1(i)+vy2(i)+vy3(i)+vy4(i))
137 vdz(i) = fourth*(vz1(i)+vz2(i)+vz3(i)+vz4(i))
138 ENDDO
139
140
141
142
143
144 IF(mtn == 51 .AND.
ale%UPWIND%UPWM /= 3)
THEN
145 iadbuf = ipm(27,mat(1))
146 iflg = nint(bufmat(31+iadbuf-1))
147 IF(iflg > 1)RETURN
148 IF(
ale%UPWIND%UPWM == 0.)
THEN
149 DO i=lft,llt
150 gamma(i)= pm(15,mat(i))
151 ENDDO
152 ELSE
153 DO i=lft,llt
154 gamma(i)=
ale%UPWIND%CUPWM
155 ENDDO
156 ENDIF
157 DO i=lft,llt
158 aaa = qmv(1,i)+qmv(2,i)+qmv(3,i)+qmv(4,i)
159 aaa = fourth * aaa
160 qmv(1,i) = fourth * (qmv(1,i) - aaa)
161 qmv(2,i) = fourth * (qmv(2,i) - aaa)
162 qmv(3,i) = fourth * (qmv(3,i) - aaa)
163 qmv(4,i) = fourth * (qmv(4,i) - aaa)
164 dmy = zero
165 dmz = zero
166 dm = qmv(4,i)+qmv(1,i)
167 dmy = dmy + vy1(i)*dm
168 dmz = dmz + vz1(i)*dm
169 dm = qmv(1,i)+qmv(2,i)
170 dmy = dmy + vy2(i)*dm
171 dmz = dmz + vz2(i)*dm
172 dm = qmv(2,i)+qmv(3,i)
173 dmy = dmy + vy3(i)*dm
174 dmz = dmz + vz3(i)*dm
175 dm = qmv(3,i)+qmv(4,i)
176 dmy = dmy + vy4(i)*dm
177 dmz = dmz + vz4(i)*dm
178 dmy = -fourth * dmy
179 dmz = -fourth * dmz
180 a1(i) = py1(i)*vdy(i)+pz1(i)*vdz(i)
181 a2(i) = py2(i)*vdy(i)+pz2(i)*vdz(i)
182 a1(i) = sign(gamma(i),a1(i))
183 a2(i) = sign(gamma(i),a2(i))
184 t11(i) = t11(i)+(one+a1(i))*dmy
185 t12(i) = t12(i)+(one+a2(i))*dmy
186 t13(i) = t13(i)+(one-a1(i))*dmy
187 t14(i) = t14(i)+(one-a2(i))*dmy
188 t21(i) = t21(i)+(one+a1(i))*dmz
189 t22(i) = t22(i)+(one+a2(i))*dmz
190 t23(i) = t23(i)+(one-a1(i))*dmz
191 t24(i) = t24(i)+(one-a2(i))*dmz
192 ENDDO
193 ELSEIF (mtn /= 51 .AND.
ale%UPWIND%UPWM /= 3)
THEN
194
195
196
197 DO i=lft,llt
198 f1(i) = (vdy(i)*dyy(i)+vdz(i)*dyz(i))*xms(i)*off(i)
199 f2(i) = (vdy(i)*dzy(i)+vdz(i)*dzz(i))*xms(i)*off(i)
200 ENDDO
201 DO i=lft,llt
202 a1(i) = py1(i)*vdy(i)+pz1(i)*vdz(i)
203 a2(i) = py2(i)*vdy(i)+pz2(i)*vdz(i)
204 ENDDO
205 DO i=lft,llt
206 a1(i) = sign(gamma(i),a1(i))
207 a2(i) = sign(gamma(i),a2(i))
208 ENDDO
209 DO i=lft,llt
210 t11(i) = t11(i)+(one+a1(i))*f1(i)
211 t12(i) = t12(i)+(one+a2(i))*f1(i)
212 t13(i) = t13(i)+(one-a1(i))*f1(i)
213 t14(i) = t14(i)+(one-a2(i))*f1(i)
214
215 t21(i) = t21(i)+(one+a1(i))*f2(i)
216 t22(i) = t22(i)+(one+a2(i))*f2(i)
217 t23(i) = t23(i)+(one-a1(i))*f2(i)
218 t24(i) = t24(i)+(one-a2(i))*f2(i)
219 ENDDO
220
221 ELSEIF (
ale%UPWIND%UPWM == 3)
THEN
222 pc=fourth
223 pd=one_over_16
224 ps=one_over_16
225 pp=one_over_16
226 DO i=lft,llt
227 y13(i)=y3(i)-y1(i)
228 z13(i)=z3(i)-z1(i)
229 y24(i)=y4(i)-y2(i)
230 z24(i)=z4(i)-z2(i)
231 y12(i)=y2(i)-y1(i)
232 z12(i)=z2(i)-z1(i)
233 y23(i)=y3(i)-y2(i)
234 z23(i)=z3(i)-z2(i)
235 y34(i)=y4(i)-y3(i)
236 z34(i)=z4(i)-z3(i)
237 y41(i)=y1(i)-y4(i)
238 z41(i)=z1(i)-z4(i)
239 ENDDO
240 DO i=lft,llt
241 s(2,i)=half*(y12(i)-y34(i))
242 s(3,i)=half*(z12(i)-z34(i))
243 t(2,i)=half*(-y41(i)+y23(i))
244 t(3,i)=half*(-z41(i)+z23(i))
245 ENDDO
246
247
248
249
251 1 rho, vis, vy1, vy1,
252 2 vz1, s, s, t,
253 3 deltax, gamma, llt)
254 DO i=lft,llt
255 tr(i)=(y41(i)*z12(i)-y12(i)*z41(i))
256 dyyp(i)=(-z24(i)*fvy1(i)-z41(i)*fvy2(i)-z12(i)*fvy4(i))
257 dyzp(i)=( y24(i)*fvy1(i)+y41(i)*fvy2(i)+y12(i)*fvy4(i))
258 dzzp(i)=( y24(i)*fvz1(i)+y41(i)*fvz2(i)+y12(i)*fvz4(i))
259 dzyp(i)=(-z24(i)*fvz1(i)-z41(i)*fvz2(i)-z12(i)*fvz4(i))
260 f1p(i)=rho(i)*(vy1(i)*dyyp(i)+vz1(i)*dyzp(i))
261 f2p(i)=rho(i)*(vy1(i)*dzyp(i)+vz1(i)*dzzp(i))
262 fac=zero
263 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
264 a1(i)=fac*(-z24(i)*vy1(i)+y24(i)*vz1(i))
265 a2(i)=fac*(-z41(i)*vy1(i)+y41(i)*vz1(i))
266 a3(i)=fac*(-z12(i)*vy1(i)+y12(i)*vz1(i))
267
268 t11(i) = t11(i)+(pp+a1(i))*f1p(i)
269 t12(i) = t12(i)+(ps+a2(i))*f1p(i)
270 t13(i) = t13(i)+pd*f1p(i)
271 t14(i) = t14(i)+(ps+a3(i))*f1p(i)
272 t21(i) = t21(i)+(pp+a1(i))*f2p(i)
273 t22(i) = t22(i)+(ps+a2(i))*f2p(i)
274 t23(i) = t23(i)+pd*f2p(i)
275 t24(i) = t24(i)+(ps+a3(i))*f2p(i)
276 ENDDO
277
278
279
280
282 1 rho, vis, vy2, vy2,
283 2 vz2, s, s, t,
284 3 deltax, gamma, llt)
285 DO i=lft,llt
286 tr(i)=(y12(i)*z23(i)-y23(i)*z12(i))
287 dyyp(i)=(-z23(i)*fvy1(i)+z13(i)*fvy2(i)-z12(i)*fvy3(i))
288 dyzp(i)=( y23(i)*fvy1(i)-y13(i)*fvy2(i)+y12(i)*fvy3(i))
289 dzzp(i)=( y23(i)*fvz1(i)-y13(i)*fvz2(i)+y12(i)*fvz3(i))
290 dzyp(i)=(-z23(i)*fvz1(i)+z13(i)*fvz2(i)-z12(i)*fvz3(i))
291 f1p(i)=rho(i)*(vy2(i)*dyyp(i)+vz2(i)*dyzp(i))
292 f2p(i)=rho(i)*(vy2(i)*dzyp(i)+vz2(i)*dzzp(i))
293 fac=zero
294 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
295 a1(i)=fac*(-z23(i)*vy2(i)+y23(i)*vz2(i))
296 a2(i)=fac*( z13(i)*vy2(i)-y13(i)*vz2(i))
297 a3(i)=fac*(-z12(i)*vy2(i)+y12(i)*vz2(i))
298 t11(i) = t11(i)+(ps+a1(i))*f1p(i)
299 t12(i) = t12(i)+(pp+a2(i))*f1p(i)
300 t13(i) = t13(i)+(ps+a3(i))*f1p(i)
301 t14(i) = t14(i)+pd*f1p(i)
302 t21(i) = t21(i)+(ps+a1(i))*f2p(i)
303 t22(i) = t22(i)+(pp+a2(i))*f2p(i)
304 t23(i) = t23(i)+(ps+a3(i))*f2p(i)
305 t24(i) = t24(i)+pd*f2p(i)
306 ENDDO
307
308
309
310
312 1 rho, vis, vy3, vy3,
313 2 vz3, s, s, t,
314 3 deltax, gamma, llt)
315 DO i=lft,llt
316 tr(i)=(y23(i)*z34(i)-y34(i)*z23(i))
317 dyyp(i)=(-z34(i)*fvy2(i)+z24(i)*fvy3(i)-z23(i)*fvy4(i))
318 dyzp(i)=( y34(i)*fvy2(i)-y24(i)*fvy3(i)+y23(i)*fvy4(i))
319 dzzp(i)=( y34(i)*fvz2(i)-y24(i)*fvz3(i)+y23(i)*fvz4(i))
320 dzyp(i)=(-z34(i)*fvz2(i)+z24(i)*fvz3(i)-z23(i)*fvz4(i))
321 f1p(i)=rho(i)*(vy3(i)*dyyp(i)+vz3(i)*dyzp(i))
322 f2p(i)=rho(i)*(vy3(i)*dzyp(i)+vz3(i)*dzzp(i))
323 fac=zero
324 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
325 a1(i)=fac*(-z34(i)*vy3(i)+y34(i)*vz3(i))
326 a2(i)=fac*( z24(i)*vy3(i)-y24(i)*vz3(i))
327 a3(i)=fac*(-z23(i)*vy3(i)+y23(i)*vz3(i))
328 t11(i) = t11(i)+pd*f1p(i)
329 t12(i) = t12(i)+(ps+a1(i))*f1p(i)
330 t13(i) = t13(i)+(pp+a2(i))*f1p(i)
331 t14(i) = t14(i)+(ps+a3(i))*f1p(i)
332 t21(i) = t21(i)+pd*f2p(i)
333 t22(i) = t22(i)+(ps+a1(i))*f2p(i)
334 t23(i) = t23(i)+(pp+a2(i))*f2p(i)
335 t24(i) = t24(i)+(ps+a3(i))*f2p(i)
336 ENDDO
337
338
339
340
342 1 rho, vis, vy4, vy4,
343 2 vz4, s, s, t,
344 3 deltax, gamma, llt)
345 DO i=lft,llt
346 tr(i)=(y34(i)*z41(i)-y41(i)*z34(i))
347 dyyp(i)=(-z34(i)*fvy1(i)-z41(i)*fvy3(i)-z13(i)*fvy4(i))
348 dyzp(i)=( y34(i)*fvy1(i)+y41(i)*fvy3(i)+y13(i)*fvy4(i))
349 dzzp(i)=( y34(i)*fvz1(i)+y41(i)*fvz3(i)+y13(i)*fvz4(i))
350 dzyp(i)=(-z34(i)*fvz1(i)-z41(i)*fvz3(i)-z13(i)*fvz4(i))
351 f1p(i)=rho(i)*(vy4(i)*dyyp(i)+vz4(i)*dyzp(i))
352 f2p(i)=rho(i)*(vy4(i)*dzyp(i)+vz4(i)*dzzp(i))
353 fac=zero
354 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
355 a1(i)=fac*(-z34(i)*vy4(i)+y34(i)*vz4(i))
356 a2(i)=fac*(-z41(i)*vy4(i)+y41(i)*vz4(i))
357 a3(i)=fac*(-z13(i)*vy4(i)+y13(i)*vz4(i))
358 t11(i) = t11(i)+(ps+a1(i))*f1p(i)
359 t12(i) = t12(i)+pd*f1p(i)
360 t13(i) = t13(i)+(ps+a2(i))*f1p(i)
361 t14(i) = t14(i)+(pp+a3(i))*f1p(i)
362 t21(i) = t21(i)+(ps+a1(i))*f2p(i)
363 t22(i) = t22(i)+pd*f2p(i)
364 t23(i) = t23(i)+(ps+a2(i))*f2p(i)
365 t24(i) = t24(i)+(pp+a3(i))*f2p(i)
366 ENDDO
367 ENDIF
368
369 RETURN
subroutine upwind(rho, vis, vdx, vdy, vdz, r, s, t, deltax, gam, nel)