5 &(flow,phiel,nseg,nelem,nelmax,eltseg,oriseg,ikle,iniflo,iopt,fn,
6 & yaflulim,flulim,yaflulimebe,flulimebe)
74 INTEGER,
INTENT(IN) :: NSEG,IOPT,NELEM,NELMAX
75 INTEGER,
INTENT(IN) :: ELTSEG(nelmax,3)
76 INTEGER,
INTENT(IN) :: ORISEG(nelmax,3)
77 INTEGER,
INTENT(IN) :: IKLE(nelmax,3)
78 DOUBLE PRECISION,
INTENT(INOUT) :: FLOW(nseg)
79 DOUBLE PRECISION,
INTENT(IN) :: PHIEL(nelmax,3)
80 LOGICAL,
INTENT(IN) :: INIFLO
81 TYPE(bief_obj),
INTENT(IN),
OPTIONAL :: FN
82 DOUBLE PRECISION,
INTENT(IN),
OPTIONAL :: FLULIMEBE(*)
83 DOUBLE PRECISION,
INTENT(IN),
OPTIONAL :: FLULIM(*)
84 LOGICAL,
INTENT(IN),
OPTIONAL :: YAFLULIM,YAFLULIMEBE
88 INTEGER IELEM,ISEG,I1,I2,I3
89 DOUBLE PRECISION A1,A2,A3,K1,K2,K3,THIRD,F12,F23,F31,F1,F2,F3
90 DOUBLE PRECISION FN1,FN2,FN3,F21,F32,F13
91 DOUBLE PRECISION BETA1FI,BETA2FI,BETA3FI,FI
92 LOGICAL TFLULIM,TFLULIMEBE
100 IF(
PRESENT(yaflulim))
THEN 102 IF(tflulim.AND..NOT.
PRESENT(flulim))
THEN 103 WRITE(
lu,*)
'FLUX_EF_VF:' 104 WRITE(
lu,*)
'VARIABLE FLULIM MISSING' 111 IF(
PRESENT(yaflulimebe))
THEN 112 tflulimebe=yaflulimebe
113 IF(tflulimebe.AND..NOT.
PRESENT(flulimebe))
THEN 114 WRITE(
lu,*)
'FLUX_EF_VF:' 115 WRITE(
lu,*)
'VARIABLE FLULIMEBE MISSING' 148 iseg = eltseg(ielem,1)
149 IF(oriseg(ielem,1).EQ.1)
THEN 150 flow(iseg) = flow(iseg) +
151 & phiel(ielem,1)*flulimebe((1-1)*nelmax+ielem)
153 flow(iseg) = flow(iseg) -
154 & phiel(ielem,1)*flulimebe((1-1)*nelmax+ielem)
157 iseg = eltseg(ielem,2)
158 IF(oriseg(ielem,2).EQ.1)
THEN 159 flow(iseg) = flow(iseg) +
160 & phiel(ielem,2)*flulimebe((2-1)*nelmax+ielem)
162 flow(iseg) = flow(iseg) -
163 & phiel(ielem,2)*flulimebe((2-1)*nelmax+ielem)
166 iseg = eltseg(ielem,3)
167 IF(oriseg(ielem,3).EQ.1)
THEN 168 flow(iseg) = flow(iseg) +
169 & phiel(ielem,3)*flulimebe((3-1)*nelmax+ielem)
171 flow(iseg) = flow(iseg) -
172 & phiel(ielem,3)*flulimebe((3-1)*nelmax+ielem)
179 iseg = eltseg(ielem,1)
180 IF(oriseg(ielem,1).EQ.1)
THEN 181 flow(iseg) = flow(iseg) + phiel(ielem,1)
183 flow(iseg) = flow(iseg) - phiel(ielem,1)
186 iseg = eltseg(ielem,2)
187 IF(oriseg(ielem,2).EQ.1)
THEN 188 flow(iseg) = flow(iseg) + phiel(ielem,2)
190 flow(iseg) = flow(iseg) - phiel(ielem,2)
193 iseg = eltseg(ielem,3)
194 IF(oriseg(ielem,3).EQ.1)
THEN 195 flow(iseg) = flow(iseg) + phiel(ielem,3)
197 flow(iseg) = flow(iseg) - phiel(ielem,3)
204 ELSEIF(iopt.EQ.0)
THEN 217 iseg = eltseg(ielem,1)
218 IF(oriseg(ielem,1).EQ.1)
THEN 219 flow(iseg) = flow(iseg) +
220 & third*(f1-f2)*flulimebe((1-1)*nelmax + ielem)
222 flow(iseg) = flow(iseg) -
223 & third*(f1-f2)*flulimebe((1-1)*nelmax + ielem)
226 iseg = eltseg(ielem,2)
227 IF(oriseg(ielem,2).EQ.1)
THEN 228 flow(iseg) = flow(iseg) +
229 & third*(f2-f3)*flulimebe((2-1)*nelmax + ielem)
231 flow(iseg) = flow(iseg) -
232 & third*(f2-f3)*flulimebe((2-1)*nelmax + ielem)
235 iseg = eltseg(ielem,3)
236 IF(oriseg(ielem,3).EQ.1)
THEN 237 flow(iseg) = flow(iseg) +
238 & third*(f3-f1)*flulimebe((3-1)*nelmax + ielem)
240 flow(iseg) = flow(iseg) -
241 & third*(f3-f1)*flulimebe((3-1)*nelmax + ielem)
251 iseg = eltseg(ielem,1)
252 IF(oriseg(ielem,1).EQ.1)
THEN 253 flow(iseg) = flow(iseg) + third*(f1-f2)
255 flow(iseg) = flow(iseg) - third*(f1-f2)
258 iseg = eltseg(ielem,2)
259 IF(oriseg(ielem,2).EQ.1)
THEN 260 flow(iseg) = flow(iseg) + third*(f2-f3)
262 flow(iseg) = flow(iseg) - third*(f2-f3)
265 iseg = eltseg(ielem,3)
266 IF(oriseg(ielem,3).EQ.1)
THEN 267 flow(iseg) = flow(iseg) + third*(f3-f1)
269 flow(iseg) = flow(iseg) - third*(f3-f1)
276 ELSEIF(iopt.EQ.2)
THEN 290 IF(a1.GE.a2.AND.a1.GE.a3)
THEN 292 iseg = eltseg(ielem,1)
293 IF(oriseg(ielem,1).EQ.1)
THEN 294 flow(iseg) = flow(iseg) -
295 & f2*flulimebe((1-1)*nelmax + ielem)
297 flow(iseg) = flow(iseg) +
298 & f2*flulimebe((1-1)*nelmax + ielem)
300 iseg = eltseg(ielem,3)
301 IF(oriseg(ielem,3).EQ.2)
THEN 302 flow(iseg) = flow(iseg) -
303 & f3*flulimebe((3-1)*nelmax + ielem)
305 flow(iseg) = flow(iseg) +
306 & f3*flulimebe((3-1)*nelmax + ielem)
308 ELSEIF(a2.GE.a1.AND.a2.GE.a3)
THEN 310 iseg = eltseg(ielem,1)
311 IF(oriseg(ielem,1).EQ.2)
THEN 312 flow(iseg) = flow(iseg) -
313 & f1*flulimebe((1-1)*nelmax + ielem)
315 flow(iseg) = flow(iseg) +
316 & f1*flulimebe((1-1)*nelmax + ielem)
318 iseg = eltseg(ielem,2)
319 IF(oriseg(ielem,2).EQ.1)
THEN 320 flow(iseg) = flow(iseg) -
321 & f3*flulimebe((2-1)*nelmax + ielem)
323 flow(iseg) = flow(iseg) +
324 & f3*flulimebe((2-1)*nelmax + ielem)
328 iseg = eltseg(ielem,2)
329 IF(oriseg(ielem,2).EQ.2)
THEN 330 flow(iseg) = flow(iseg) -
331 & f2*flulimebe((2-1)*nelmax + ielem)
333 flow(iseg) = flow(iseg) +
334 & f2*flulimebe((2-1)*nelmax + ielem)
336 iseg = eltseg(ielem,3)
337 IF(oriseg(ielem,3).EQ.1)
THEN 338 flow(iseg) = flow(iseg) -
339 & f1*flulimebe((3-1)*nelmax + ielem)
341 flow(iseg) = flow(iseg) +
342 & f1*flulimebe((3-1)*nelmax + ielem)
354 IF(a1.GE.a2.AND.a1.GE.a3)
THEN 356 iseg = eltseg(ielem,1)
357 IF(oriseg(ielem,1).EQ.1)
THEN 358 flow(iseg) = flow(iseg) - f2
360 flow(iseg) = flow(iseg) + f2
362 iseg = eltseg(ielem,3)
363 IF(oriseg(ielem,3).EQ.2)
THEN 364 flow(iseg) = flow(iseg) - f3
366 flow(iseg) = flow(iseg) + f3
368 ELSEIF(a2.GE.a1.AND.a2.GE.a3)
THEN 370 iseg = eltseg(ielem,1)
371 IF(oriseg(ielem,1).EQ.2)
THEN 372 flow(iseg) = flow(iseg) - f1
374 flow(iseg) = flow(iseg) + f1
376 iseg = eltseg(ielem,2)
377 IF(oriseg(ielem,2).EQ.1)
THEN 378 flow(iseg) = flow(iseg) - f3
380 flow(iseg) = flow(iseg) + f3
384 iseg = eltseg(ielem,2)
385 IF(oriseg(ielem,2).EQ.2)
THEN 386 flow(iseg) = flow(iseg) - f2
388 flow(iseg) = flow(iseg) + f2
390 iseg = eltseg(ielem,3)
391 IF(oriseg(ielem,3).EQ.1)
THEN 392 flow(iseg) = flow(iseg) - f1
394 flow(iseg) = flow(iseg) + f1
402 ELSEIF(iopt.EQ.3.AND.
PRESENT(fn))
THEN 418 f12=max(min(k1,-k2),0.d0)
419 f23=max(min(k2,-k3),0.d0)
420 f31=max(min(k3,-k1),0.d0)
421 f21=max(min(k2,-k1),0.d0)
422 f32=max(min(k3,-k2),0.d0)
423 f13=max(min(k1,-k3),0.d0)
432 beta1fi=f12*(fn1-fn2)+f13*(fn1-fn3)
433 beta2fi=f21*(fn2-fn1)+f23*(fn2-fn3)
434 beta3fi=f31*(fn3-fn1)+f32*(fn3-fn2)
436 fi=beta1fi+beta2fi+beta3fi
444 IF(beta1fi.GT.fi)
THEN 447 ELSEIF(beta1fi.LT.0.d0)
THEN 451 IF(beta2fi.GT.fi)
THEN 454 ELSEIF(beta2fi.LT.0.d0)
THEN 458 IF(beta3fi.GT.fi)
THEN 461 ELSEIF(beta3fi.LT.0.d0)
THEN 465 ELSEIF(fi.LT.0.d0)
THEN 466 IF(beta1fi.LT.fi)
THEN 469 ELSEIF(beta1fi.GT.0.d0)
THEN 473 IF(beta2fi.LT.fi)
THEN 476 ELSEIF(beta2fi.GT.0.d0)
THEN 480 IF(beta3fi.LT.fi)
THEN 483 ELSEIF(beta3fi.GT.0.d0)
THEN 502 iseg = eltseg(ielem,1)
503 IF(oriseg(ielem,1).EQ.1)
THEN 504 flow(iseg) = flow(iseg) +
505 & (+ f21 - f12)*flulimebe((1-1)*nelmax + ielem)
507 flow(iseg) = flow(iseg) +
508 & (- f21 + f12)*flulimebe((1-1)*nelmax + ielem)
511 iseg = eltseg(ielem,2)
512 IF(oriseg(ielem,2).EQ.1)
THEN 513 flow(iseg) = flow(iseg) +
514 & (+ f32 - f23)*flulimebe((2-1)*nelmax + ielem)
516 flow(iseg) = flow(iseg) +
517 & (- f32 + f23)*flulimebe((2-1)*nelmax + ielem)
520 iseg = eltseg(ielem,3)
521 IF(oriseg(ielem,3).EQ.1)
THEN 522 flow(iseg) = flow(iseg) +
523 & (+ f13 - f31)*flulimebe((3-1)*nelmax + ielem)
525 flow(iseg) = flow(iseg) +
526 & (- f13 + f31)*flulimebe((3-1)*nelmax + ielem)
539 f12=max(min(k1,-k2),0.d0)
540 f23=max(min(k2,-k3),0.d0)
541 f31=max(min(k3,-k1),0.d0)
542 f21=max(min(k2,-k1),0.d0)
543 f32=max(min(k3,-k2),0.d0)
544 f13=max(min(k1,-k3),0.d0)
553 beta1fi=f12*(fn1-fn2)+f13*(fn1-fn3)
554 beta2fi=f21*(fn2-fn1)+f23*(fn2-fn3)
555 beta3fi=f31*(fn3-fn1)+f32*(fn3-fn2)
557 fi=beta1fi+beta2fi+beta3fi
565 IF(beta1fi.GT.fi)
THEN 568 ELSEIF(beta1fi.LT.0.d0)
THEN 572 IF(beta2fi.GT.fi)
THEN 575 ELSEIF(beta2fi.LT.0.d0)
THEN 579 IF(beta3fi.GT.fi)
THEN 582 ELSEIF(beta3fi.LT.0.d0)
THEN 586 ELSEIF(fi.LT.0.d0)
THEN 587 IF(beta1fi.LT.fi)
THEN 590 ELSEIF(beta1fi.GT.0.d0)
THEN 594 IF(beta2fi.LT.fi)
THEN 597 ELSEIF(beta2fi.GT.0.d0)
THEN 601 IF(beta3fi.LT.fi)
THEN 604 ELSEIF(beta3fi.GT.0.d0)
THEN 623 iseg = eltseg(ielem,1)
624 IF(oriseg(ielem,1).EQ.1)
THEN 625 flow(iseg) = flow(iseg) + f21 - f12
627 flow(iseg) = flow(iseg) - f21 + f12
630 iseg = eltseg(ielem,2)
631 IF(oriseg(ielem,2).EQ.1)
THEN 632 flow(iseg) = flow(iseg) + f32 - f23
634 flow(iseg) = flow(iseg) - f32 + f23
637 iseg = eltseg(ielem,3)
638 IF(oriseg(ielem,3).EQ.1)
THEN 639 flow(iseg) = flow(iseg) + f13 - f31
641 flow(iseg) = flow(iseg) - f13 + f31
650 WRITE(
lu,*)
'FLUX_EF_VF:' 651 WRITE(
lu,*)
'UNKNOWN OPTION: ',iopt
654 WRITE(
lu,*)
'OPTION 3: ADVECTED FUNCTION REQUIRED' 665 IF(tflulim.AND..NOT.tflulimebe)
THEN 667 flow(iseg) = flow(iseg)*flulim(iseg)
subroutine flux_ef_vf(FLOW, PHIEL, NSEG, NELEM, NELMAX, ELTSEG, ORISEG, IKLE, INIFLO, IOPT, FN, YAFLULIM, FLULIM, YAFLULIMEBE, FLULIMEBE)