File: /nfs/home/0/users/jenkins/mfix.git/model/xsi_mod.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 MODULE XSI
21
22 CONTAINS
23
24 SUBROUTINE CALC_XSI(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T,incr)
25
26
27
28
29
30
31
32
33 USE chischeme
34 USE compar
35 USE discretization
36 USE functions
37 USE geometry
38 USE indices
39 USE param
40 USE param1
41 USE run
42 USE sendrecv
43 USE vshear
44
45 IMPLICIT NONE
46
47
48
49
50
51
52
53
54 INTEGER DISCR
55
56
57 DOUBLE PRECISION PHI(DIMENSION_3)
58
59
60 DOUBLE PRECISION U(DIMENSION_3), V(DIMENSION_3), W(DIMENSION_3)
61
62
63 DOUBLE PRECISION XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),&
64 XSI_t(DIMENSION_3)
65
66
67 INTEGER IJK, IJKC, IJKD, IJKU, I, J, K
68
69
70 CHARACTER(LEN=80) LINE(1)
71
72 DOUBLE PRECISION PHI_C
73
74
75 DOUBLE PRECISION dwf
76
77
78 DOUBLE PRECISION cf
79
80
81 DOUBLE PRECISION oDXc, oDXuc, oDYc, oDYuc, oDZc, oDZuc
82
83 INTEGER incr
84
85 IF (SHEAR) THEN
86
87
88 call CXS(incr,DISCR,U,V,W,PHI,XSI_E,XSI_N,XSI_T)
89
90
91 ELSE
92 SELECT CASE (DISCR)
93 CASE (:1)
94
95
96 DO IJK = ijkstart3, ijkend3
97 XSI_E(IJK) = XSI_func(U(IJK),ZERO)
98 XSI_N(IJK) = XSI_func(V(IJK),ZERO)
99 IF (DO_K) XSI_T(IJK) = XSI_func(W(IJK),ZERO)
100 END DO
101 CASE (2)
102
103
104 DO IJK = ijkstart3, ijkend3
105 IF (U(IJK) >= ZERO) THEN
106 IJKC = IJK
107 IJKD = EAST_OF(IJK)
108 IJKU = WEST_OF(IJKC)
109 ELSE
110 IJKC = EAST_OF(IJK)
111 IJKD = IJK
112 IJKU = EAST_OF(IJKC)
113 ENDIF
114
115 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
116 DWF = SUPERBEE(PHI_C)
117 XSI_E(IJK) = XSI_func(U(IJK),DWF)
118
119 IF (V(IJK) >= ZERO) THEN
120 IJKC = IJK
121 IJKD = NORTH_OF(IJK)
122 IJKU = SOUTH_OF(IJKC)
123 ELSE
124 IJKC = NORTH_OF(IJK)
125 IJKD = IJK
126 IJKU = NORTH_OF(IJKC)
127 ENDIF
128
129 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
130 DWF = SUPERBEE(PHI_C)
131 XSI_N(IJK) = XSI_func(V(IJK),DWF)
132
133 IF (DO_K) THEN
134 IF (W(IJK) >= ZERO) THEN
135 IJKC = IJK
136 IJKD = TOP_OF(IJK)
137 IJKU = BOTTOM_OF(IJKC)
138 ELSE
139 IJKC = TOP_OF(IJK)
140 IJKD = IJK
141 IJKU = TOP_OF(IJKC)
142 ENDIF
143
144 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
145
146 DWF = SUPERBEE(PHI_C)
147
148 XSI_T(IJK) = XSI_func(W(IJK),DWF)
149 ENDIF
150 END DO
151 CASE (3)
152
153
154 DO IJK = ijkstart3, ijkend3
155 IF (U(IJK) >= ZERO) THEN
156 IJKC = IJK
157 IJKD = EAST_OF(IJK)
158 IJKU = WEST_OF(IJKC)
159 ELSE
160 IJKC = EAST_OF(IJK)
161 IJKD = IJK
162 IJKU = EAST_OF(IJKC)
163 ENDIF
164
165 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
166
167 if(Chi_flag) then
168 DWF = Chi_SMART(PHI_C, CHI_e(IJK))
169 else
170 DWF = SMART(PHI_C)
171 endif
172
173 XSI_E(IJK) = XSI_func(U(IJK),DWF)
174
175 IF (V(IJK) >= ZERO) THEN
176 IJKC = IJK
177 IJKD = NORTH_OF(IJK)
178 IJKU = SOUTH_OF(IJKC)
179 ELSE
180 IJKC = NORTH_OF(IJK)
181 IJKD = IJK
182 IJKU = NORTH_OF(IJKC)
183 ENDIF
184
185 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
186
187 if(Chi_flag) then
188 DWF = Chi_SMART(PHI_C, CHI_n(IJK))
189 else
190 DWF = SMART(PHI_C)
191 endif
192
193 XSI_N(IJK) = XSI_func(V(IJK),DWF)
194
195 IF (DO_K) THEN
196 IF (W(IJK) >= ZERO) THEN
197 IJKC = IJK
198 IJKD = TOP_OF(IJK)
199 IJKU = BOTTOM_OF(IJKC)
200 ELSE
201 IJKC = TOP_OF(IJK)
202 IJKD = IJK
203 IJKU = TOP_OF(IJKC)
204 ENDIF
205
206 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
207
208 if(Chi_flag) then
209 DWF = Chi_SMART(PHI_C, CHI_t(IJK))
210 else
211 DWF = SMART(PHI_C)
212 endif
213
214 XSI_T(IJK) = XSI_func(W(IJK),DWF)
215 ENDIF
216 END DO
217 CASE (4)
218
219
220 DO IJK = ijkstart3, ijkend3
221 I = I_OF(IJK)
222 IF (U(IJK) >= ZERO) THEN
223 IJKC = IJK
224 IJKD = EAST_OF(IJK)
225 IJKU = WEST_OF(IJKC)
226 ELSE
227 IJKC = EAST_OF(IJK)
228 IJKD = IJK
229 IJKU = EAST_OF(IJKC)
230 ENDIF
231
232 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
233 CF = ABS(U(IJK))*DT*ODX_E(I)
234 DWF = ULTRA_QUICK(PHI_C,CF)
235 XSI_E(IJK) = XSI_func(U(IJK),DWF)
236
237 J = J_OF(IJK)
238 IF (V(IJK) >= ZERO) THEN
239 IJKC = IJK
240 IJKD = NORTH_OF(IJK)
241 IJKU = SOUTH_OF(IJKC)
242 ELSE
243 IJKC = NORTH_OF(IJK)
244 IJKD = IJK
245 IJKU = NORTH_OF(IJKC)
246 ENDIF
247
248 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
249 CF = ABS(V(IJK))*DT*ODY_N(J)
250 DWF = ULTRA_QUICK(PHI_C,CF)
251 XSI_N(IJK) = XSI_func(V(IJK),DWF)
252
253 IF (DO_K) THEN
254 K = K_OF(IJK)
255 IF (W(IJK) >= ZERO) THEN
256 IJKC = IJK
257 IJKD = TOP_OF(IJK)
258 IJKU = BOTTOM_OF(IJKC)
259 ELSE
260 IJKC = TOP_OF(IJK)
261 IJKD = IJK
262 IJKU = TOP_OF(IJKC)
263 ENDIF
264
265 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
266 CF = ABS(W(IJK))*DT*OX(I)*ODZ_T(K)
267 DWF = ULTRA_QUICK(PHI_C,CF)
268
269 XSI_T(IJK) = XSI_func(W(IJK),DWF)
270 ENDIF
271 END DO
272 CASE (5)
273
274
275
276
277
278 DO IJK = ijkstart3, ijkend3
279 I = I_OF(IJK)
280 IF (U(IJK) >= ZERO) THEN
281 IJKC = IJK
282 IJKD = EAST_OF(IJK)
283 IJKU = WEST_OF(IJKC)
284 ODXC = ODX(I)
285 ODXUC = ODX_E(IM1(I))
286 ELSE
287 IJKC = EAST_OF(IJK)
288 IJKD = IJK
289 IJKU = EAST_OF(IJKC)
290 ODXC = ODX(IP1(I))
291 ODXUC = ODX_E(IP1(I))
292 ENDIF
293
294 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
295 CF = ABS(U(IJK))*DT*ODX_E(I)
296 DWF = QUICKEST(PHI_C,CF,ODXC,ODXUC,ODX_E(I))
297 XSI_E(IJK) = XSI_func(U(IJK),DWF)
298
299 J = J_OF(IJK)
300 IF (V(IJK) >= ZERO) THEN
301 IJKC = IJK
302 IJKD = NORTH_OF(IJK)
303 IJKU = SOUTH_OF(IJKC)
304 ODYC = ODY(J)
305 ODYUC = ODY_N(JM1(J))
306 ELSE
307 IJKC = NORTH_OF(IJK)
308 IJKD = IJK
309 IJKU = NORTH_OF(IJKC)
310 ODYC = ODY(JP1(J))
311 ODYUC = ODY_N(JP1(J))
312 ENDIF
313
314 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
315 CF = ABS(V(IJK))*DT*ODY_N(J)
316 DWF = QUICKEST(PHI_C,CF,ODYC,ODYUC,ODY_N(J))
317 XSI_N(IJK) = XSI_func(V(IJK),DWF)
318
319 IF (DO_K) THEN
320 K = K_OF(IJK)
321 IF (W(IJK) >= ZERO) THEN
322 IJKC = IJK
323 IJKD = TOP_OF(IJK)
324 IJKU = BOTTOM_OF(IJKC)
325 ODZC = ODZ(K)
326 ODZUC = ODZ_T(KM1(K))
327 ELSE
328 IJKC = TOP_OF(IJK)
329 IJKD = IJK
330 IJKU = TOP_OF(IJKC)
331 ODZC = ODZ(KP1(K))
332 ODZUC = ODZ_T(KP1(K))
333 ENDIF
334
335 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
336 CF = ABS(W(IJK))*DT*OX(I)*ODZ_T(K)
337 DWF = QUICKEST(PHI_C,CF,ODZC,ODZUC,ODZ_T(K))
338
339 XSI_T(IJK) = XSI_func(W(IJK),DWF)
340 ENDIF
341 END DO
342 CASE (6)
343
344
345 DO IJK = ijkstart3, ijkend3
346 IF (U(IJK) >= ZERO) THEN
347 IJKC = IJK
348 IJKD = EAST_OF(IJK)
349 IJKU = WEST_OF(IJKC)
350 ELSE
351 IJKC = EAST_OF(IJK)
352 IJKD = IJK
353 IJKU = EAST_OF(IJKC)
354 ENDIF
355
356 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
357
358 if(Chi_flag) then
359 DWF = Chi_MUSCL(PHI_C, CHI_e(IJK))
360 else
361 DWF = MUSCL(PHI_C)
362 endif
363
364 XSI_E(IJK) = XSI_func(U(IJK),DWF)
365
366 IF (V(IJK) >= ZERO) THEN
367 IJKC = IJK
368 IJKD = NORTH_OF(IJK)
369 IJKU = SOUTH_OF(IJKC)
370 ELSE
371 IJKC = NORTH_OF(IJK)
372 IJKD = IJK
373 IJKU = NORTH_OF(IJKC)
374 ENDIF
375
376 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
377
378 if(Chi_flag) then
379 DWF = Chi_MUSCL(PHI_C, CHI_n(IJK))
380 else
381 DWF = MUSCL(PHI_C)
382 endif
383
384 XSI_N(IJK) = XSI_func(V(IJK),DWF)
385
386 IF (DO_K) THEN
387 IF (W(IJK) >= ZERO) THEN
388 IJKC = IJK
389 IJKD = TOP_OF(IJK)
390 IJKU = BOTTOM_OF(IJKC)
391 ELSE
392 IJKC = TOP_OF(IJK)
393 IJKD = IJK
394 IJKU = TOP_OF(IJKC)
395 ENDIF
396
397 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
398
399 if(Chi_flag) then
400 DWF = Chi_MUSCL(PHI_C, CHI_t(IJK))
401 else
402 DWF = MUSCL(PHI_C)
403 endif
404
405 XSI_T(IJK) = XSI_func(W(IJK),DWF)
406 ENDIF
407 END DO
408 CASE (7)
409
410
411
412 DO IJK = ijkstart3, ijkend3
413 IF (U(IJK) >= ZERO) THEN
414 IJKC = IJK
415 IJKD = EAST_OF(IJK)
416 IJKU = WEST_OF(IJKC)
417 ELSE
418 IJKC = EAST_OF(IJK)
419 IJKD = IJK
420 IJKU = EAST_OF(IJKC)
421 ENDIF
422
423 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
424 DWF = VANLEER(PHI_C)
425 XSI_E(IJK) = XSI_func(U(IJK),DWF)
426
427 IF (V(IJK) >= ZERO) THEN
428 IJKC = IJK
429 IJKD = NORTH_OF(IJK)
430 IJKU = SOUTH_OF(IJKC)
431 ELSE
432 IJKC = NORTH_OF(IJK)
433 IJKD = IJK
434 IJKU = NORTH_OF(IJKC)
435 ENDIF
436
437 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
438 DWF = VANLEER(PHI_C)
439 XSI_N(IJK) = XSI_func(V(IJK),DWF)
440
441 IF (DO_K) THEN
442 IF (W(IJK) >= ZERO) THEN
443 IJKC = IJK
444 IJKD = TOP_OF(IJK)
445 IJKU = BOTTOM_OF(IJKC)
446 ELSE
447 IJKC = TOP_OF(IJK)
448 IJKD = IJK
449 IJKU = TOP_OF(IJKC)
450 ENDIF
451
452 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
453 DWF = VANLEER(PHI_C)
454 XSI_T(IJK) = XSI_func(W(IJK),DWF)
455 ENDIF
456 END DO
457 CASE (8)
458
459
460 DO IJK = ijkstart3, ijkend3
461 IF (U(IJK) >= ZERO) THEN
462 IJKC = IJK
463 IJKD = EAST_OF(IJK)
464 IJKU = WEST_OF(IJKC)
465 ELSE
466 IJKC = EAST_OF(IJK)
467 IJKD = IJK
468 IJKU = EAST_OF(IJKC)
469 ENDIF
470
471 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
472 DWF = MINMOD(PHI_C)
473 XSI_E(IJK) = XSI_func(U(IJK),DWF)
474
475 IF (V(IJK) >= ZERO) THEN
476 IJKC = IJK
477 IJKD = NORTH_OF(IJK)
478 IJKU = SOUTH_OF(IJKC)
479 ELSE
480 IJKC = NORTH_OF(IJK)
481 IJKD = IJK
482 IJKU = NORTH_OF(IJKC)
483 ENDIF
484
485 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
486 DWF = MINMOD(PHI_C)
487 XSI_N(IJK) = XSI_func(V(IJK),DWF)
488
489 IF (DO_K) THEN
490 IF (W(IJK) >= ZERO) THEN
491 IJKC = IJK
492 IJKD = TOP_OF(IJK)
493 IJKU = BOTTOM_OF(IJKC)
494 ELSE
495 IJKC = TOP_OF(IJK)
496 IJKD = IJK
497 IJKU = TOP_OF(IJKC)
498 ENDIF
499
500 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
501 DWF = MINMOD(PHI_C)
502 XSI_T(IJK) = XSI_func(W(IJK),DWF)
503 ENDIF
504 END DO
505
506
507 CASE (9)
508
509
510 DO IJK = ijkstart3, ijkend3
511 IF (U(IJK) >= ZERO) THEN
512 IJKC = IJK
513 IJKD = EAST_OF(IJK)
514 IJKU = WEST_OF(IJKC)
515 ELSE
516 IJKC = EAST_OF(IJK)
517 IJKD = IJK
518 IJKU = EAST_OF(IJKC)
519 ENDIF
520
521 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
522 DWF = CENTRAL_SCHEME(PHI_C)
523 XSI_E(IJK) = XSI_func(U(IJK),DWF)
524
525
526 IF (V(IJK) >= ZERO) THEN
527 IJKC = IJK
528 IJKD = NORTH_OF(IJK)
529 IJKU = SOUTH_OF(IJKC)
530 ELSE
531 IJKC = NORTH_OF(IJK)
532 IJKD = IJK
533 IJKU = NORTH_OF(IJKC)
534 ENDIF
535
536 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
537 DWF = CENTRAL_SCHEME(PHI_C)
538 XSI_N(IJK) = XSI_func(V(IJK),DWF)
539
540
541 IF (DO_K) THEN
542 IF (W(IJK) >= ZERO) THEN
543 IJKC = IJK
544 IJKD = TOP_OF(IJK)
545 IJKU = BOTTOM_OF(IJKC)
546 ELSE
547 IJKC = TOP_OF(IJK)
548 IJKD = IJK
549 IJKU = TOP_OF(IJKC)
550 ENDIF
551
552 PHI_C = PHI_C_OF(PHI(IJKU),PHI(IJKC),PHI(IJKD))
553 DWF = CENTRAL_SCHEME(PHI_C)
554 XSI_T(IJK) = XSI_func(W(IJK),DWF)
555 ENDIF
556 END DO
557
558
559 CASE DEFAULT
560 WRITE (LINE, '(A,I2,A)') 'DISCRETIZE = ', DISCR, ' not supported.'
561 CALL WRITE_ERROR ('CALC_XSI', LINE, 1)
562 CALL MFIX_EXIT(myPE)
563 END SELECT
564
565 ENDIF
566
567 call send_recv(XSI_E,2)
568 call send_recv(XSI_N,2)
569 call send_recv(XSI_T,2)
570
571 RETURN
572 END SUBROUTINE CALC_XSI
573
574
575
576
577
578
579
580 SUBROUTINE CXS(INCR,DISCR,U,V,W,PHI,XSI_E,XSI_N,XSI_T)
581
582
583
584 USE param
585 USE param1
586 USE run
587 USE geometry
588 USE indices
589 USE vshear
590 USE compar
591 USE sendrecv
592 USE functions
593 IMPLICIT NONE
594
595 INTEGER incr,IJK,IJKC,IJKD,IJKU,I,DISCR
596 DOUBLE PRECISION V(DIMENSION_3), U(DIMENSION_3), W(DIMENSION_3)
597 DOUBLE PRECISION PHI(DIMENSION_3),DWFE,DWFN,DWFT
598 DOUBLE PRECISION PHICU,PHIDU,PHIUU,PHICV,PHIDV,PHIUV
599 DOUBLE PRECISION PHICW,PHIDW,PHIUW
600 DOUBLE PRECISION XSI_E(DIMENSION_3), XSI_N(DIMENSION_3)
601 DOUBLE PRECISION XSI_T(DIMENSION_3),SRT
602
603 IF (INCR .eq. 2) THEN
604 =(2d0*V_sh/XLENGTH)
605
606
607
608
609 DO IJK = ijkstart3, ijkend3
610
611 I=I_OF(IJK)
612 V(IJK)=V(IJK)+VSH(IJK)
613
614 IF (U(IJK) >= ZERO) THEN
615 IJKC = IJK
616 IJKD = EAST_OF(IJK)
617 IJKU = WEST_OF(IJKC)
618 PHICU=PHI(IJKC)
619 PHIDU=PHI(IJKD)+SRT*1d0/oDX_E(I)
620 PHIUU=PHI(IJKU)-SRT*1d0/oDX_E(IM1(I))
621
622
623
624 ELSE
625 IJKC = EAST_OF(IJK)
626 IJKD = IJK
627 IJKU = EAST_OF(IJKC)
628 PHICU=PHI(IJKC)+SRT*1d0/oDX_E(I)
629 PHIDU=PHI(IJKD)
630 PHIUU=PHI(IJKU)+SRT*1d0/oDX_E(I)&
631 +SRT*1d0/oDX_E(IP1(I))
632 ENDIF
633
634 IF (V(IJK) >= ZERO) THEN
635 IJKC = IJK
636 IJKD = NORTH_OF(IJK)
637 IJKU= SOUTH_OF(IJKC)
638 PHICV=PHI(IJKC)
639 PHIDV=PHI(IJKD)
640 PHIUV=PHI(IJKU)
641 ELSE
642 IJKC= NORTH_OF(IJK)
643 IJKD= IJK
644 IJKU= NORTH_OF(IJKC)
645 PHICV=PHI(IJKC)
646 PHIDV=PHI(IJKD)
647 PHIUV=PHI(IJKU)
648 ENDIF
649
650 IF (DO_K) THEN
651 IF (W(IJK) >= ZERO) THEN
652 IJKC = IJK
653 IJKD = TOP_OF(IJK)
654 IJKU = BOTTOM_OF(IJKC)
655 ELSE
656 IJKC = TOP_OF(IJK)
657 IJKD = IJK
658 IJKU = TOP_OF(IJKC)
659 ENDIF
660 PHICW=PHI(IJKC)
661 PHIDW=PHI(IJKD)
662 PHIUW=PHI(IJKU)
663
664 ELSE
665 PHICW=0d0
666 PHIDW=0d0
667 PHIUW=0d0
668 W(IJK)=0d0
669 DWFT=0d0
670 END IF
671
672 CALL DW(U(IJK),V(IJK),W(IJK),IJK,PHICU,PHIDU,PHIUU,PHICV,&
673 PHIDV,PHIUV,PHICW,PHIDW,PHIUW,DWFE,DWFN,DWFT,DISCR)
674 XSI_E(IJK) = XSI_func(U(IJK),DWFE)
675 XSI_N(IJK) = XSI_func(V(IJK),DWFN)
676 XSI_T(IJK) = XSI_func(W(IJK),DWFT)
677
678 V(IJK)=V(IJK)-VSH(IJK)
679 END DO
680
681
682 ELSE IF (INCR .eq. 1) THEN
683
684
685
686
687 DO IJK = ijkstart3, ijkend3
688
689 V(IJK)=V(IJK)+VSHE(IJK)
690 IF (U(IJK) >= ZERO) THEN
691 IJKC = IJK
692 IJKD = EAST_OF(IJK)
693 IJKU = WEST_OF(IJKC)
694 PHICU=PHI(IJKC)
695 PHIDU=PHI(IJKD)
696 PHIUU=PHI(IJKU)
697 ELSE
698 IJKC = EAST_OF(IJK)
699 IJKD = IJK
700 IJKU = EAST_OF(IJKC)
701 PHICU=PHI(IJKC)
702 PHIDU=PHI(IJKD)
703 PHIUU=PHI(IJKU)
704 ENDIF
705
706 IF (V(IJK) >= ZERO) THEN
707 IJKC = IJK
708 IJKD = NORTH_OF(IJK)
709 IJKU= SOUTH_OF(IJKC)
710 PHICV=PHI(IJKC)
711 PHIDV=PHI(IJKD)
712 PHIUV=PHI(IJKU)
713 ELSE
714 IJKC= NORTH_OF(IJK)
715 IJKD= IJK
716 IJKU= NORTH_OF(IJKC)
717 PHICV=PHI(IJKC)
718 PHIDV=PHI(IJKD)
719 PHIUV=PHI(IJKU)
720 ENDIF
721
722 IF (DO_K) THEN
723 IF (W(IJK) >= ZERO) THEN
724 IJKC = IJK
725 IJKD = TOP_OF(IJK)
726 IJKU = BOTTOM_OF(IJKC)
727 ELSE
728 IJKC = TOP_OF(IJK)
729 IJKD = IJK
730 IJKU = TOP_OF(IJKC)
731 ENDIF
732 PHICW=PHI(IJKC)
733 PHIDW=PHI(IJKD)
734 PHIUW=PHI(IJKU)
735 ELSE
736 PHICW=0d0
737 PHIDW=0d0
738 PHIUW=0d0
739 W(IJK)=0d0
740 DWFT=0d0
741 END IF
742
743 CALL DW(U(IJK),V(IJK),W(IJK),IJK,PHICU,PHIDU,PHIUU,PHICV,&
744 PHIDV,PHIUV,PHICW,PHIDW,PHIUW,DWFE,DWFN,DWFT,DISCR)
745
746 XSI_E(IJK) = XSI_func(U(IJK),DWFE)
747 XSI_N(IJK) = XSI_func(V(IJK),DWFN)
748 XSI_T(IJK) = XSI_func(W(IJK),DWFT)
749
750 V(IJK)=V(IJK)-VSHE(IJK)
751 END DO
752
753
754 ELSE IF (INCR .eq. 0) THEN
755
756
757
758
759 DO IJK = ijkstart3, ijkend3
760
761 V(IJK)=V(IJK)+VSH(IJK)
762 IF (U(IJK) >= ZERO) THEN
763 IJKC = IJK
764 IJKD = EAST_OF(IJK)
765 IJKU = WEST_OF(IJKC)
766 PHICU=PHI(IJKC)
767 PHIDU=PHI(IJKD)
768 PHIUU=PHI(IJKU)
769
770 ELSE
771 IJKC = EAST_OF(IJK)
772 IJKD = IJK
773 IJKU = EAST_OF(IJKC)
774 PHICU=PHI(IJKC)
775 PHIDU=PHI(IJKD)
776 PHIUU=PHI(IJKU)
777 ENDIF
778
779 IF (V(IJK) >= ZERO) THEN
780 IJKC = IJK
781 IJKD = NORTH_OF(IJK)
782 IJKU= SOUTH_OF(IJKC)
783 PHICV=PHI(IJKC)
784 PHIDV=PHI(IJKD)
785 PHIUV=PHI(IJKU)
786 ELSE
787 IJKC= NORTH_OF(IJK)
788 IJKD= IJK
789 IJKU= NORTH_OF(IJKC)
790 PHICV=PHI(IJKC)
791 PHIDV=PHI(IJKD)
792 PHIUV=PHI(IJKU)
793 ENDIF
794
795 IF (DO_K) THEN
796 IF (W(IJK) >= ZERO) THEN
797 IJKC = IJK
798 IJKD = TOP_OF(IJK)
799 IJKU = BOTTOM_OF(IJKC)
800 ELSE
801 IJKC = TOP_OF(IJK)
802 IJKD = IJK
803 IJKU = TOP_OF(IJKC)
804 ENDIF
805 PHICW=PHI(IJKC)
806 PHIDW=PHI(IJKD)
807 PHIUW=PHI(IJKU)
808 ELSE
809 PHICW=0d0
810 PHIDW=0d0
811 PHIUW=0d0
812 W(IJK)=0d0
813 DWFT=0d0
814 END IF
815
816
817
818 CALL DW(U(IJK),V(IJK),W(IJK),IJK,PHICU,PHIDU,PHIUU,PHICV,&
819 PHIDV,PHIUV,PHICW,PHIDW,PHIUW,DWFE,DWFN,DWFT,DISCR)
820
821
822 XSI_E(IJK) = XSI_func(U(IJK),DWFE)
823 XSI_N(IJK) = XSI_func(V(IJK),DWFN)
824 XSI_T(IJK) = XSI_func(W(IJK),DWFT)
825
826 V(IJK)=V(IJK)-VSH(IJK)
827 END DO
828
829 ELSE
830 write(*,*) 'INCR ERROR'
831 END IF
832
833 call send_recv(XSI_E,2)
834 call send_recv(XSI_N,2)
835 call send_recv(XSI_T,2)
836
837
838
839
840
841 RETURN
842 END SUBROUTINE CXS
843
844
845
846
847
848
849
850
851 SUBROUTINE DW(UU,VV,WW,IJK,PHICU,PHIDU,PHIUU,PHICV,&
852 PHIDV,PHIUV,PHICW,PHIDW,PHIUW,DWFE,DWFN,DWFT,DISCR)
853
854
855
856
857
858 USE compar
859 USE discretization
860 USE functions
861 USE geometry
862 USE indices
863 USE param
864 USE param1
865 USE run
866 USE sendrecv
867 USE vshear
868
869 IMPLICIT NONE
870
871
872 INTEGER DISCR
873
874
875 INTEGER IJK, I, J, K
876
877 DOUBLE PRECISION PHI_C
878
879
880 DOUBLE PRECISION PHICU, PHIDU, PHIUU, DWFE, UU
881 DOUBLE PRECISION PHICV, PHIDV, PHIUV, DWFN, VV
882 DOUBLE PRECISION PHICW, PHIDW, PHIUW, DWFT, WW
883
884
885 DOUBLE PRECISION oDXc, oDXuc, oDYc, oDYuc, oDZc, oDZuc
886 DOUBLE PRECISION CF,DZUC
887
888
889 SELECT CASE (DISCR)
890 CASE (:1)
891
892 DWFE=0d0
893 DWFN=0d0
894 DWFT=0d0
895 CASE (2)
896 = PHI_C_OF(PHIUU,PHICU,PHIDU)
897 DWFE = SUPERBEE(PHI_C)
898
899 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
900 DWFN = SUPERBEE(PHI_C)
901
902 IF (DO_K) THEN
903 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
904 DWFT = SUPERBEE(PHI_C)
905 END IF
906
907 CASE (3)
908
909 = PHI_C_OF(PHIUU,PHICU,PHIDU)
910 DWFE = SMART(PHI_C)
911
912 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
913 DWFN = SMART(PHI_C)
914
915 IF (DO_K) THEN
916 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
917 DWFT = SMART(PHI_C)
918 END IF
919 CASE (4)
920
921
922 =I_OF(IJK)
923 J=J_OF(IJK)
924
925
926 PHI_C = PHI_C_OF(PHIUU,PHICU,PHIDU)
927 CF = ABS(UU)*DT*ODX_E(I)
928 DWFE = ULTRA_QUICK(PHI_C,CF)
929
930 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
931 CF = ABS(VV)*DT*ODY_N(J)
932 DWFN = ULTRA_QUICK(PHI_C,CF)
933
934
935 IF (DO_K) THEN
936 K=K_OF(IJK)
937 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
938 CF = ABS(WW)*DT*OX(I)*ODZ_T(K)
939 DWFT = ULTRA_QUICK(PHI_C,CF)
940 END IF
941 CASE(5)
942
943 =I_OF(IJK)
944 J=J_OF(IJK)
945
946 IF (UU >= ZERO) THEN
947 ODXC = ODX(I)
948 ODXUC = ODX_E(IM1(I))
949 ELSE
950 ODXC = ODX(IP1(I))
951 ODXUC = ODX_E(IP1(I))
952 ENDIF
953
954 PHI_C = PHI_C_OF(PHIUU,PHICU,PHIDU)
955 CF = ABS(UU)*DT*ODX_E(I)
956 DWFE = QUICKEST(PHI_C,CF,ODXC,ODXUC,ODX_E(I))
957
958 IF (VV >= ZERO) THEN
959
960 ODYC = ODY(J)
961 ODYUC = ODY_N(JM1(J))
962 ELSE
963 ODYC = ODY(JP1(J))
964 ODYUC = ODY_N(JP1(J))
965 ENDIF
966
967 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
968 CF = ABS(VV)*DT*ODY_N(J)
969 DWFN = QUICKEST(PHI_C,CF,ODYC,ODYUC,ODY_N(J))
970
971 IF (DO_K) THEN
972 K=K_OF(IJK)
973 IF (WW >= ZERO) THEN
974 ODZC = ODZ(K)
975 ODZUC = ODZ_T(KM1(K))
976 ELSE
977 ODZC = ODZ(KP1(K))
978 DZUC = ODZ_T(KP1(K))
979 ENDIF
980
981 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
982 CF = ABS(WW)*DT*OX(I)*ODZ_T(K)
983 DWFT = QUICKEST(PHI_C,CF,ODZC,ODZUC,ODZ_T(K))
984 END IF
985
986 CASE (6)
987 = PHI_C_OF(PHIUU,PHICU,PHIDU)
988 DWFE = MUSCL(PHI_C)
989
990 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
991 DWFN = MUSCL(PHI_C)
992
993 IF (DO_K) THEN
994 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
995 DWFT = MUSCL(PHI_C)
996 END IF
997
998 CASE (7)
999
1000 = PHI_C_OF(PHIUU,PHICU,PHIDU)
1001 DWFE = VANLEER(PHI_C)
1002
1003 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
1004 DWFN = VANLEER(PHI_C)
1005
1006 IF (DO_K) THEN
1007 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
1008 DWFT = VANLEER(PHI_C)
1009 END IF
1010
1011 CASE (8)
1012 = PHI_C_OF(PHIUU,PHICU,PHIDU)
1013 DWFE = MINMOD(PHI_C)
1014
1015 PHI_C = PHI_C_OF(PHIUV,PHICV,PHIDV)
1016 DWFN = MINMOD(PHI_C)
1017
1018 IF (DO_K) THEN
1019 PHI_C = PHI_C_OF(PHIUW,PHICW,PHIDW)
1020 DWFT = MINMOD(PHI_C)
1021 END IF
1022 CASE DEFAULT
1023 WRITE (*,*) 'DISCRETIZE = ', DISCR, ' not supported.'
1024
1025 CALL MFIX_EXIT(myPE)
1026
1027 END SELECT
1028 RETURN
1029 END SUBROUTINE DW
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040 DOUBLE PRECISION FUNCTION XSI_func(XXXv,XXXdwf)
1041 IMPLICIT NONE
1042 DOUBLE PRECISION XXXv,XXXdwf
1043 XSI_func = (sign(1d0,(-XXXv))+1d0)/(2d0) + &
1044 sign(1d0,XXXv)*XXXdwf
1045 END FUNCTION XSI_func
1046
1047 END MODULE XSI
1048