36
37
38
39 USE elbufdef_mod
43
44
45
46
47
48
49
50
51
52
53
54
55
56#include "implicit_f.inc"
57
58
59
60#include "spmd_c.inc"
61#include "com01_c.inc"
62#include "com04_c.inc"
63#include "com08_c.inc"
64#include "param_c.inc"
65#include "task_c.inc"
66
67
68
69 INTEGER IXS(NIXS,*),NV46,ITRIMAT,ITASK
71 . flux(nv46,*),vfrac(*)
72 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
73
74
75
76 INTEGER K,KK,J1,J2,J
78 . vol0,av0,uav0,alphi,ualphi,aaa,ff(nv46,5),udt,phi0
79 INTEGER :: IE, MLW, IADJv, NADJv, IB, NBF, NBL, ICELL,ICELLM, MCELL, IE_M, IBM,NG,IDLOC,NADJ,IADJ
80 INTEGER :: NIN,NCELL,IBV,IFV,ICELLv, IEV
81 my_real :: volg, alph, alphv(6,5), tmpflux(nv46
83 LOGICAL :: debug_outp
84
85
86
87 IF(trimat==0)RETURN
88
89
90
91
92 IF(dt1>zero)THEN
93 udt = one/dt1
94 ELSE
95 udt = zero
96 ENDIF
97
98 nin = 1
99 nbf = 1+itask*
nb/nthread
100 nbl = (itask+1)*
nb/nthread
102
103
104
105 debug_outp = .false.
107 debug_outp = .false.
109 do ib=nbf,nbl
113 if(mlw==51)then
114 debug_outp=.true.
115 endif
116 endif
117 enddo
119 debug_outp = .true.
120 kk = 1
121 do ib=nbf,nbl
123 if(mlw/=51)then
124 kk = 0
125 endif
126 enddo
127 if (kk==0)debug_outp=.false.
128 endif
130 endif
131 if(debug_outp)then
132 print *, " |----ale51_antidiff3_int22.F-----|"
133 print *, " | THREAD INFORMATION |"
134 print *, " |--------------------------------|"
135 print *, " NCYCLE =", ncycle
136 print *, " ITRIMAT =", itrimat
137 endif
138
139
140
141 DO ib=nbf,nbl
146 icell = 0
147 DO WHILE (icell<=ncell)
148 icell = icell +1
149 IF (icell>ncell .AND. ncell/=0)icell=9
150
151 j =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
152 icellm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(2)
153 IF(j==0)THEN
154 ie_m = ie
155 ibm = ib
156 icellm = mcell
157 ELSEIF(j<=nv46)THEN
160 ELSE
161 j1 = j/10
162 j2 = mod(j,10)
163 ibv =
brick_list(nin,ib )%Adjacent_Brick(j1,4)
164 ie_m =
brick_list(nin,ibv)%Adjacent_Brick(j2,1)
165 ibm =
brick_list(nin,ibv)%Adjacent_Brick(j2,4)
166 ENDIF
169 mlw
170 IF(mlw/=51)cycle
171 alph =
brick_list(nin,ibm)%POLY(icellm)%VFRACm
172 volg = elbuf_tab(ng)%GBUF%VOL(idloc)
173 vol0 = volg*udt
174 av0 = alph * vol0
175 uav0 = vol0 - av0
176 alphi = zero
177 ualphi = zero
178 phi0 = zero
179
180
181
182
183 DO k=1,nv46
184 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
185 DO iadj=1,nadj
186 tmpflux(k,iadj) =
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_UpwFLUX(iadj)
187 IF(tmpflux(k,iadj)>zero)THEN
192 IF(icellv==0) THEN
193 alphv(k,iadj) = alph
194 ELSE
195 IF(ibv==0)THEN
196 IF(iev==0)print *, "inter22 : potential material leakage, Check domain boundaries..."
197 alphv(k,iadj) = vfrac(iev)
198 ELSE
199 alphv(k,iadj) =
brick_list(nin,ibv)%POLY(icellv)%VFRACm(itrimat)
200 ENDIF
201 ENDIF
202 ff(k,iadj)= alphv(k,iadj) * tmpflux(k,iadj)
203
204 alphi = alphi + ff(k,iadj)
205
206 phi0 = phi0 + tmpflux(k,iadj)
207 ENDIF
208 enddo
209 enddo
210
211 ualphi = phi0 - alphi
212
213
214
215 IF(alphi>av0.AND.av0>zero)THEN
216
217
218
219 aaa = av0 / alphi
220 DO k=1,nv46
221 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
222 DO iadj=1,nadj
223 IF(tmpflux(k,iadj)>zero)THEN
224 ff(k,iadj) = ff(k,iadj) * aaa
225 ENDIF
226 enddo
227 enddo
228 ELSEIF(ualphi>uav0.AND.uav0>zero)THEN
229
230
231
232 aaa = uav0/ualphi
233 DO k=1,nv46
234 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
235 DO iadj=1,nadj
236 IF(tmpflux(k,iadj)>zero)THEN
237 ff(k,iadj) = tmpflux(k,iadj) + (ff(k,iadj)-tmpflux(k
238 ENDIF
239 enddo
240 enddo
241 ENDIF
242
243
244
245 DO k=1,nv46
246 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
247 DO iadj=1,nadj
251 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_Cell(iadj
252 IF(tmpflux(k,iadj)>zero)THEN
253 ff(k,iadj) = half * ( ff(k,iadj)*(one-
ale%UPWIND%UPWSM)+alph
254
255
256
257
258 if(debug_outp)then
261
262 print *, " brique ="
263 print *, " icell ="
264 print *, " FACE =", k
265 print *, " ALPH =", alph
266 print *, " ALPHv =", alphv(k,iadj)
267 write (*,fmt=
'(A,6E26.14)')
" WAS Flux(J) =",
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_upwFLUX(iadj)
268 write (*,fmt='(A,6E26.14)')" IS Flux(J) =", ff(k,iadj)
269 print *, " ------------------------"
270 endif
271 endif
272
273
274
275 brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_UpwFLUX(iadj) = ff(k,iadj)
276
277
278 IF(icellv>0)THEN
279
280
281 IF(ibv>0)THEN
282
283 nadjv =
brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%NAdjCell
284 DO iadjv=1,nadjv
285 IF(
brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%Adjacent_Cell(iadjv)==icell)
EXIT
286
287 ENDDO
288 debug_tmp =
brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%Adjacent_UpwFLUX(iadjv)
289 brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%Adjacent_UpwFLUX(iadjv) = -ff(k,iadj)
290 ELSE
291
292 debug_tmp = flux(ifv,iev)
293 flux(ifv,iev) = -ff(k,iadj)
294 ENDIF
295
296
297 if(debug_outp)then
299 print *, " => Setting adjacent flux consequently :"
300 print *, " brique.V =", ixs(11,iev)
301 print *, " icell.V =", icellv
302 print *, " FACE.V =", ifv
303 write (*,fmt='(A,6E26.14)')
304 . " WAS Flux(J) =", debug_tmp
305 write (*,fmt='(A,6E26.14)')
306 . " IS Flux(J) =", -ff(k,iadj)
307 print *, " ---"
308 endif
309 endif
310
311 ELSE
312
313 ENDIF
314 ENDIF
315 enddo
316 enddo
317
318 enddo
319 enddo
320
321
322 RETURN
type(brick_entity), dimension(:,:), allocatable, target brick_list