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