49
50
51
53
54
55
56#include "implicit_f.inc"
57
58
59
60#include "mvsiz_p.inc"
61
62
63
64#include "inter22.inc"
65
66
67
68 INTEGER, INTENT(IN) :: NEL
70 . a0(mvsiz,3),a1(mvsiz,3),a2(mvsiz,3),a3(mvsiz,3),
71 . ajv(*),r(mvsiz,3),s(mvsiz,3),t(mvsiz,3),
72 . vx1(*),vx2(*),vx3(*),vx0(*),vdx0(*),
73 . vy1(*),vy2(*),vy3(*),vy0(*),vdy0(*),
74 . vz1(*),vz2(*),vz3(*),vz0(*),vdz0(*),
75 . f11(*),f21(*),f31(*),f12(*),f22(*),f32(*),
76 . f13(*),f23(*),f33(*),f14(*),f24(*),f34(*),
77 . f15(*),f25(*),f35(*),f16(*),f26(*),f36(*),
78 . f17(*),f27(*),f37(*),f18(*),f28(*),f38(*),
79 . rho(*),deltax(*),vis(*),vol(*),iad22(*)
80 INTEGER :: NC1(*), NC2(*), NC3(*), NC4(*),
81 . NC5(*), NC6(*), NC7(*), NC8(*), NALE(*)
82
83
84
86 . gam(mvsiz),
87 . dxxp(mvsiz),dxyp(mvsiz),dxzp(mvsiz),
88 . dyxp(mvsiz),dyyp(mvsiz),dyzp(mvsiz),
89 . dzxp(mvsiz),dzyp(mvsiz),dzzp(mvsiz),
90 . up0(mvsiz),up1(mvsiz),up2(mvsiz),up3(mvsiz),
91 . f1p(mvsiz),f2p(mvsiz),f3p(mvsiz),
92 . f1(mvsiz),f2(mvsiz),f3(mvsiz),fac,
93 . off1(mvsiz), off2(mvsiz), off3(mvsiz), off4(mvsiz),
94 . off5(mvsiz), off6(mvsiz), off7(mvsiz), off8(mvsiz)
95 INTEGER I
96
98 1 rho, vis, vdx0, vdy0,
99 2 vdz0, r, s, t,
100 3 deltax, gam, nel)
101 off1(1:mvsiz) = one
102 off2(1:mvsiz) = one
103 off3(1:mvsiz) = one
104 off4(1:mvsiz) = one
105 off5(1:mvsiz) = one
106 off6(1:mvsiz) = one
107 off7(1:mvsiz) = one
108 off8(1:mvsiz) = one
109 DO i=1,nel
110 IF (nale(nc1(i)) == 0) off1(i) = zero
111 IF (nale(nc2(i)) == 0) off2(i) = zero
112 IF (nale(nc3(i)) == 0) off3(i) = zero
113 IF (nale(nc4(i)) == 0) off4(i) = zero
114 IF (nale(nc5(i)) == 0) off5(i) = zero
115 IF (nale(nc6(i)) == 0) off6(i) = zero
116 IF (nale(nc7(i)) == 0) off7(i) = zero
117 IF (nale(nc8(i)) == 0) off8(i) = zero
118 ENDDO
119 IF (
ale%UPWIND%UPWM==2)
THEN
120
121
122
123
124
125 DO i=1,nel
126 gam(i)=gam(i)/ajv(i)
127
128 dxxp(i)=a0(i,1)*vx0(i)+a1(i,1)*vx1(i)+a2(i,1)*vx2(i)+a3(i,1)*vx3(i)
129 dxyp(i)=a0(i,2)*vx0(i)+a1(i,2)*vx1(i)+a2(i,2)*vx2(i)+a3(i,2)*vx3(i)
130 dxzp(i)=a0(i,3)*vx0(i)+a1(i,3)*vx1(i)+a2(i,3)*vx2(i)+a3(i,3)*vx3(i)
131 dyxp(i)=a0(i,1)*vy0(i)+a1(i,1)*vy1(i)+a2(i,1)*vy2(i)+a3(i,1)*vy3(i)
132 dyyp(i)=a0(i,2)*vy0(i)+a1(i,2)*vy1(i)+a2(i,2)*vy2(i)+a3(i,2)*vy3(i)
133 dyzp(i)=a0(i,3)*vy0(i)+a1(i,3)*vy1(i)+a2(i,3)*vy2(i)+a3(i,3)*vy3(i)
134 dzxp(i)=a0(i,1)*vz0(i)+a1(i,1)*vz1(i)+a2(i,1)*vz2(i)+a3(i,1)*vz3(i)
135 dzyp(i)=a0(i,2)*vz0(i)+a1(i,2)*vz1(i)+a2(i,2)*vz2(i)+a3(i,2)*vz3(i)
136 dzzp(i)=a0(i,3)*vz0(i)+a1(i,3)*vz1(i)+a2(i,3)*vz2(i)+a3(i,3)*vz3(i)
137
138 fac=eight*gam(i)
139 up0(i)=one +fac*(a0(i,1)*vdx0(i)+a0(i,2)*vdy0(i)+a0(i,3)*vdz0(i))
140 up1(i)=one +fac*(a1(i,1)*vdx0(i)+a1(i,2)*vdy0(i)+a1(i,3)*vdz0(i))
141 up2(i)=one +fac*(a2(i,1)*vdx0(i)+a2(i,2)*vdy0(i)+a2(i,3)*vdz0(i))
142 up3(i)=one +fac*(a3(i,1)*vdx0(i)+a3(i,2)*vdy0(i)+a3(i,3)*vdz0(i))
143
144 fac=rho(i)*vol(i)/ajv(i)/sixty4
145 f1p(i)= fac*(vdx0(i)*dxxp(i)+vdy0(i)*dxyp(i)+vdz0(i)*dxzp(i))
146 f2p(i)= fac*(vdx0(i)*dyxp(i)+vdy0(i)*dyyp(i)+vdz0(i)*dyzp(i))
147 f3p(i)= fac*(vdx0(i)*dzxp(i)+vdy0(i)*dzyp(i)+vdz0(i)*dzzp(i))
148
149 f11(i) = f11(i) - up0(i)*f1p(i)*off1(i)
150 f21(i) = f21(i) - up0(i)*f2p(i)*off1(i)
151 f31(i) = f31(i) - up0(i)*f3p(i)*off1(i)
152 f12(i) = f12(i) - up1(i)*f1p(i)*off2(i)
153 f22(i) = f22(i) - up1(i)*f2p(i)*off2(i)
154 f32(i) = f32(i) - up1(i)*f3p(i)*off2(i)
155 f13(i) = f13(i) - up2(i)*f1p(i)*off3(i)
156 f23(i) = f23(i) - up2(i)*f2p(i)*off3(i)
157 f33(i) = f33(i) - up2(i)*f3p(i)*off3(i)
158 f14(i) = f14(i) - up3(i)*f1p(i)*off4(i)
159 f24(i) = f24(i) - up3(i)*f2p(i)*off4(i)
160 f34(i) = f34(i) - up3(i)*f3p(i)*off4(i)
161 f15(i) = f15(i) - f1p(i)*off5(i)
162 f25(i) = f25(i) - f2p(i)*off5(i)
163 f35(i) = f35(i) - f3p(i)*off5(i)
164 f16(i) = f16(i) - f1p(i)*off6(i)
165 f26(i) = f26(i) - f2p(i)*off6(i)
166 f36(i) = f36(i) - f3p(i)*off6(i)
167 f17(i) = f17(i) - f1p(i)*off7(i)
168 f27(i) = f27(i) - f2p(i)*off7(i)
169 f37(i) = f37(i) - f3p(i)*off7(i)
170 f18(i) = f18(i) - f1p(i)*off8(i)
171 f28(i) = f28(i) - f2p(i)*off8(i)
172 f38(i) = f38(i) - f3p(i)*off8(i)
173 ENDDO
174
175
176 IF(int22>0)THEN
177 DO i=1,nel
178 IF(nint(iad22(i))==0)THEN
179 !currently unplugged only
for cut cell buffer
180 f1p(i) = zero
181 f2p(i) = zero
182 f3p(i) = zero
183 ENDIF
184 ENDDO
185 ENDIF
186
187 DO i=1,nel
188 f1(i)=-gam(i)*(-f1p(i)*(dyyp(i)+dzzp(i))+f2p(i)*dxyp(i)+f3p(i)*dxzp(i))
189 f2(i)=-gam(i)*( f1p(i)*dyxp(i)-f2p(i)*(dxxp(i)+dzzp(i))+f3p(i)*dyzp(i))
190 f3(i)=-gam(i)*( f1p(i)*dzxp(i)+f2p(i)*dzyp(i)-f3p(i)*(dyyp(i)+dxxp(i)))
191 f11(i) = f11(i) - f1(i)*off1(i)
192 f21(i) = f21(i) - f2(i)*off1(i)
193 f31(i) = f31(i) - f3(i)*off1(i)
194 f12(i) = f12(i) - f1(i)*off2(i)
195 f22(i) = f22(i) - f2(i)*off2(i)
196 f32(i) = f32(i) - f3(i)*off2(i)
197 f13(i) = f13(i) - f1(i)*off3(i)
198 f23(i) = f23(i) - f2(i)*off3(i)
199 f33(i) = f33(i) - f3(i)*off3(i)
200 f14(i) = f14(i) - f1(i)*off4(i)
201 f24(i) = f24(i) - f2(i)*off4(i)
202 f34(i) = f34(i) - f3(i)*off4(i)
203 f15(i) = f15(i) - f1(i)*off5(i)
204 f25(i) = f25(i) - f2(i)*off5(i)
205 f35(i) = f35(i) - f3(i)*off5(i)
206 f16(i) = f16(i) - f1(i)*off6(i)
207 f26(i) = f26(i) - f2(i)*off6(i)
208 f36(i) = f36(i) - f3(i)*off6(i)
209 f17(i) = f17(i) - f1(i)*off7(i)
210 f27(i) = f27(i) - f2(i)*off7(i)
211 f37(i) = f37(i) - f3(i)*off7(i)
212 f18(i) = f18(i) - f1(i)*off8(i)
213 f28(i) = f28(i) - f2(i)*off8(i)
214 f38(i) = f38(i) - f3(i)*off8(i)
215 ENDDO
216
217 ELSE
218
219
220
221
222
223 DO i=1,nel
224 gam(i)=gam(i)/ajv(i)
225 dxxp(i)=a0(i,1)*vx0(i)+a1(i,1)*vx1(i)+a2(i,1)*vx2(i)+a3(i,1)*vx3(i)
226 dxyp(i)=a0(i,2)*vx0(i)+a1(i,2)*vx1(i)+a2(i,2)*vx2(i)+a3(i,2)*vx3(i)
227 dxzp(i)=a0(i,3)*vx0(i)+a1(i,3)*vx1(i)+a2(i,3)*vx2(i)+a3(i,3)*vx3(i)
228 dyxp(i)=a0(i,1)*vy0(i)+a1(i,1)*vy1(i)+a2(i,1)*vy2(i)+a3(i,1)*vy3(i)
229 dyyp(i)=a0(i,2)*vy0(i)+a1(i,2)*vy1(i)+a2(i,2)*vy2(i)+a3(i,2)*vy3(i)
230 dyzp(i)=a0(i,3)*vy0(i)+a1(i,3)*vy1(i)+a2(i,3)*vy2(i)+a3(i,3)*vy3(i)
231 dzxp(i)=a0(i,1)*vz0(i)+a1(i,1)*vz1(i)+a2(i,1)*vz2(i)+a3(i,1)*vz3(i)
232 dzyp(i)=a0(i,2)*vz0(i)+a1(i,2)*vz1(i)+a2(i,2)*vz2(i)+a3(i,2)*vz3(i)
233 dzzp(i)=a0(i,3)*vz0(i)+a1(i,3)*vz1(i)+a2(i,3)*vz2(i)+a3(i,3)*vz3(i)
234
235 fac=eight*gam(i)
236 up0(i)=one+fac*(a0(i,1)*vdx0(i)+a0(i,2)*vdy0(i)+a0(i,3)*vdz0(i))
237 up1(i)=one+fac*(a1(i,1)*vdx0(i)+a1(i,2)*vdy0(i)+a1(i,3)*vdz0(i))
238 up2(i)=one+fac*(a2(i,1)*vdx0(i)+a2(i,2)*vdy0(i)+a2(i,3)*vdz0(i))
239 up3(i)=one+fac*(a3(i,1)*vdx0(i)+a3(i,2)*vdy0(i)+a3(i,3)*vdz0(i))
240
241 fac=rho(i)*vol(i)/ajv(i)/sixty4
242 f1p(i)= fac*(vdx0(i)*dxxp(i)+vdy0(i)*dxyp(i)+vdz0(i)*dxzp(i))
243 f2p(i)= fac*(vdx0(i)*dyxp(i)+vdy0(i)*dyyp(i)+vdz0(i)*dyzp(i))
244 f3p(i)= fac*(vdx0(i)*dzxp(i)+vdy0(i)*dzyp(i)+vdz0(i)*dzzp(i))
245 ENDDO
246
247
248 IF(int22>0)THEN
249 DO i=1,nel
250 IF(nint(iad22(i))==0)THEN
251
252 f1p(i) = zero
253 f2p(i) = zero
254 f3p(i) = zero
255 ENDIF
256 ENDDO
257 ENDIF
258
259 DO i=1,nel
260 f11(i) = f11(i) - up0(i)*f1p(i)*off1(i)
261 f21(i) = f21(i) - up0(i)*f2p(i)*off1(i)
262 f31(i) = f31(i) - up0(i)*f3p(i)*off1(i)
263 f12(i) = f12(i) - up1(i)*f1p(i)*off2(i)
264 f22(i) = f22(i) - up1(i)*f2p(i)*off2(i)
265 f32(i) = f32(i) - up1(i)*f3p(i)*off2(i)
266 f13(i) = f13(i) - up2(i)*f1p(i)*off3(i)
267 f23(i) = f23(i) - up2(i)*f2p(i)*off3(i)
268 f33(i) = f33(i) - up2(i)*f3p(i)*off3(i)
269 f14(i) = f14(i) - up3(i)*f1p(i)*off4(i)
270 f24(i) = f24(i) - up3(i)*f2p(i)*off4(i)
271 f34(i) = f34(i) - up3(i)*f3p(i)*off4(i)
272 f15(i) = f15(i) - f1p(i)*off5(i)
273 f25(i) = f25(i) - f2p(i)*off5(i)
274 f35(i) = f35(i) - f3p(i)*off5(i)
275 f16(i) = f16(i) - f1p(i)*off6(i)
276 f26(i) = f26(i) - f2p(i)*off6(i)
277 f36(i) = f36(i) - f3p(i)*off6(i)
278 f17(i) = f17(i) - f1p(i)*off7(i)
279 f27(i) = f27(i) - f2p(i)*off7(i)
280 f37(i) = f37(i) - f3p(i)*off7(i)
281 f18(i) = f18(i) - f1p(i)*off8(i)
282 f28(i) = f28(i) - f2p(i)*off8(i)
283 f38(i) = f38(i) - f3p(i)*off8(i)
284 ENDDO
285
286 ENDIF
287
288
289 RETURN
for(i8=*sizetab-1;i8 >=0;i8--)
subroutine upwind_v(rho, vis, vdx, vdy, vdz, r, s, t, deltax, gam, nel)
1.15.0