File: /nfs/home/0/users/jenkins/mfix.git/model/conv_dif_w_g.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 SUBROUTINE CONV_DIF_W_G(A_M, B_M, IER)
23
24
25
26
27
28
29
30
31 USE param
32 USE param1
33 USE parallel
34 USE matrix
35 USE geometry
36 USE indices
37 USE run
38 USE visc_g
39 USE compar
40
41 IMPLICIT NONE
42
43
44
45
46
47
48
49
50
51
52 INTEGER IER
53
54
55 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
56
57
58 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
59
60
61 IF (.NOT.MOMENTUM_Z_EQ(0)) RETURN
62
63
64
65 IF (DEF_COR) THEN
66 CALL STORE_A_W_G0 (A_M(1,-3,0), IER)
67 IF (DISCRETIZE(5) > 1)CALL STORE_A_W_GDC (A_M(1,-3,0), B_M(1,0))
68 ELSE
69
70 IF (DISCRETIZE(5) == 0) THEN
71 CALL STORE_A_W_G0 (A_M(1,-3,0), IER)
72 ELSE
73 CALL STORE_A_W_G1 (A_M(1,-3,0))
74 ENDIF
75 ENDIF
76
77 CALL DIF_W_IS (MU_GT, A_M, B_M, 0, IER)
78
79
80 RETURN
81 END SUBROUTINE CONV_DIF_W_G
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107 SUBROUTINE STORE_A_W_G0(A_W_G, IER)
108
109
110
111
112
113
114
115
116 USE param
117 USE param1
118 USE parallel
119 USE matrix
120 USE geometry
121 USE indices
122 USE run
123 USE visc_g
124 USE toleranc
125 USE physprop
126 USE fldvar
127 USE output
128 USE compar
129 USE mflux
130 USE fun_avg
131 USE functions
132
133
134
135 USE cutcell
136
137
138
139 IMPLICIT NONE
140
141
142
143
144
145
146
147
148
149
150 INTEGER IER
151
152
153 INTEGER I, J, K, IPJK, IJPK, IJKN, IJKC, KP, IJKE,&
154 IJKTE, IJKP, IJKT, IJKTN, IJK
155 INTEGER IMJK, IM, IJKW, IJKWT, IMJKP
156 INTEGER IJMK, JM, IJMKP, IJKS, IJKST
157 INTEGER IJKM, KM, IJKB
158
159
160 INTEGER M
161
162
163 DOUBLE PRECISION Flux
164
165
166 DOUBLE PRECISION D_f
167
168
169 DOUBLE PRECISION A_W_g(DIMENSION_3, -3:3)
170
171
172
173 DOUBLE PRECISION :: AW,HW,VELW
174
175
176
177
178
179
180
181 = 0
182
183
184
185
186
187
188
189
190
191
192
193
194 DO IJK = ijkstart3, ijkend3
195
196 IF (FLOW_AT_T(IJK)) THEN
197 I = I_OF(IJK)
198 J = J_OF(IJK)
199 K = K_OF(IJK)
200 IPJK = IP_OF(IJK)
201 IJPK = JP_OF(IJK)
202 IJKN = NORTH_OF(IJK )
203 IJKT = TOP_OF(IJK)
204 IF (WALL_AT(IJK)) THEN
205 IJKC = IJKT
206 ELSE
207 IJKC = IJK
208 ENDIF
209 KP = KP1(K)
210 IJKE = EAST_OF(IJK)
211 IJKP = KP_OF(IJK)
212 IJKTN = NORTH_OF(IJKT)
213 IJKTE = EAST_OF(IJKT)
214
215
216
217
218
219 IF(CUT_W_TREATMENT_AT(IJK)) THEN
220 Flux = (Theta_W_be(IJK) * Flux_gE(IJK) + Theta_W_te(IJK) * Flux_gE(IJKP))
221 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
222 Flux = Flux * AW
223 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
224 ),MU_GT(IJKTE),I),K)*ONEoDX_E_W(IJK)*AYZ_W(IJK)
225 ELSE
226 = HALF * (Flux_gE(IJK) + Flux_gE(IJKP))
227 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
228 ),MU_GT(IJKTE),I),K)*ODX_E(I)*AYZ_W(IJK)
229 ENDIF
230
231
232
233 IF (Flux >= ZERO) THEN
234 A_W_G(IJK,E) = D_F
235 A_W_G(IPJK,W) = D_F + Flux
236 ELSE
237 A_W_G(IJK,E) = D_F - Flux
238 A_W_G(IPJK,W) = D_F
239 ENDIF
240
241
242
243
244
245 IF(CUT_W_TREATMENT_AT(IJK)) THEN
246 Flux = (Theta_W_bn(IJK) * Flux_gN(IJK) + Theta_W_tn(IJK) * Flux_gN(IJKP))
247 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
248 Flux = Flux * AW
249 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
250 ),MU_GT(IJKTN),J),K)*ONEoDY_N_W(IJK)*AXZ_W(IJK)
251 ELSE
252 = HALF * (Flux_gN(IJK) + Flux_gN(IJKP))
253 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
254 ),MU_GT(IJKTN),J),K)*ODY_N(J)*AXZ_W(IJK)
255 ENDIF
256
257
258
259 IF (Flux >= ZERO) THEN
260 A_W_G(IJK,N) = D_F
261 A_W_G(IJPK,S) = D_F + Flux
262 ELSE
263 A_W_G(IJK,N) = D_F - Flux
264 A_W_G(IJPK,S) = D_F
265 ENDIF
266
267
268
269
270
271 IF(CUT_W_TREATMENT_AT(IJK)) THEN
272 Flux = (Theta_Wt_bar(IJK) * Flux_gT(IJK) + Theta_Wt(IJK) * Flux_gT(IJKP))
273 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
274 Flux = Flux * AW
275 D_F = MU_GT(IJKT)*ONEoDZ_T_W(IJK)*AXY_W(IJK)
276 ELSE
277 = HALF * (Flux_gT(IJK) + Flux_gT(IJKP))
278 D_F = MU_GT(IJKT)*OX(I)*ODZ(KP)*AXY_W(IJK)
279 ENDIF
280
281
282
283 IF (Flux >= ZERO) THEN
284 A_W_G(IJK,T) = D_F
285 A_W_G(IJKP,B) = D_F + Flux
286 ELSE
287 A_W_G(IJK,T) = D_F - Flux
288 A_W_G(IJKP,B) = D_F
289 ENDIF
290
291
292 = IM_OF(IJK)
293 IF (.NOT.FLOW_AT_T(IMJK)) THEN
294 IM = IM1(I)
295 IJKW = WEST_OF(IJK)
296 IJKWT = TOP_OF(IJKW)
297 IMJKP = KP_OF(IMJK)
298
299
300
301 IF(CUT_W_TREATMENT_AT(IJK)) THEN
302 Flux = (Theta_W_be(IMJK) * Flux_gE(IMJK) + Theta_W_te(IMJK) * Flux_gE(IMJKP))
303 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IMJK),AW,HW,VELW)
304 Flux = Flux * AW
305 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
306 IJKWT),MU_GT(IJKT),IM),K)*ONEoDX_E_W(IMJK)*AYZ_W(IMJK)
307 ELSE
308 = HALF * (Flux_gE(IMJK) + Flux_gE(IMJKP))
309 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
310 IJKWT),MU_GT(IJKT),IM),K)*ODX_E(IM)*AYZ_W(IMJK)
311 ENDIF
312
313
314
315 IF (Flux >= ZERO) THEN
316 A_W_G(IJK,W) = D_F + Flux
317 ELSE
318 A_W_G(IJK,W) = D_F
319 ENDIF
320 ENDIF
321
322
323 = JM_OF(IJK)
324 IF (.NOT.FLOW_AT_T(IJMK)) THEN
325 JM = JM1(J)
326 IJMKP = KP_OF(IJMK)
327 IJKS = SOUTH_OF(IJK)
328 IJKST = TOP_OF(IJKS)
329
330
331
332 IF(CUT_W_TREATMENT_AT(IJK)) THEN
333 Flux = (Theta_W_bn(IJMK) * Flux_gN(IJMK) + Theta_W_tn(IJMK) * Flux_gN(IJMKP))
334 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJMK),AW,HW,VELW)
335 Flux = Flux * AW
336 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
337 IJKST),MU_GT(IJKT),JM),K)*ONEoDY_N_W(IJMK)*AXZ_W(IJMK)
338 ELSE
339 = HALF * (Flux_gN(IJMK) + Flux_gN(IJMKP))
340 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
341 IJKST),MU_GT(IJKT),JM),K)*ODY_N(JM)*AXZ_W(IJMK)
342 ENDIF
343
344
345
346 IF (Flux >= ZERO) THEN
347 A_W_G(IJK,S) = D_F + Flux
348 ELSE
349 A_W_G(IJK,S) = D_F
350 ENDIF
351 ENDIF
352
353
354 = KM_OF(IJK)
355 IF (.NOT.FLOW_AT_T(IJKM)) THEN
356 KM = KM1(K)
357 IJKB = BOTTOM_OF(IJK)
358
359
360
361 IF(CUT_W_TREATMENT_AT(IJK)) THEN
362 Flux = (Theta_Wt_bar(IJKM) * Flux_gT(IJKM) + Theta_Wt(IJKM) * Flux_gT(IJK))
363 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJKM),AW,HW,VELW)
364 Flux = Flux * AW
365 D_F = MU_GT(IJK)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
366 ELSE
367 = HALF * (Flux_gT(IJKM) + Flux_gT(IJK))
368 D_F = MU_GT(IJK)*OX(I)*ODZ(K)*AXY_W(IJKM)
369 ENDIF
370
371
372
373 IF (Flux >= ZERO) THEN
374 A_W_G(IJK,B) = D_F + Flux
375 ELSE
376 A_W_G(IJK,B) = D_F
377 ENDIF
378 ENDIF
379 ENDIF
380 END DO
381
382
383 RETURN
384 END SUBROUTINE STORE_A_W_G0
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411 SUBROUTINE STORE_A_W_GDC(A_W_G, B_M)
412
413
414
415
416
417
418
419
420 USE compar
421 USE discretization, ONLY: fpfoi_of
422 USE fldvar
423 USE fun_avg
424 USE function3
425 USE functions
426 USE geometry
427 USE indices
428 USE matrix
429 USE mflux
430 USE output
431 USE parallel
432 USE param
433 USE param1
434 USE physprop
435 USE run
436 USE sendrecv
437 USE sendrecv3
438 USE toleranc
439 USE visc_g
440 USE xsi
441 USE xsi_array
442 Use tmp_array, U => Array1, V => Array2, WW => Array3
443
444
445
446 USE cutcell
447
448
449
450 IMPLICIT NONE
451
452
453
454
455
456
457
458
459 INTEGER I, J, K, IPJK, IJPK, IJKN, IJKC, KP, IJKE,&
460 IJKTE, IJKP, IJKT, IJKTN, IJK
461 INTEGER IMJK, IM, IJKW, IJKWT, IMJKP
462 INTEGER IJMK, JM, IJMKP, IJKS, IJKST
463 INTEGER IJKM, KM, IJKB
464 INTEGER IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
465 INTEGER IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
466
467
468 INTEGER incr
469
470
471
472 DOUBLE PRECISION A_W_g(DIMENSION_3, -3:3)
473
474
475 DOUBLE PRECISION B_m(DIMENSION_3)
476
477
478 DOUBLE PRECISION MOM_HO
479
480
481 DOUBLE PRECISION MOM_LO
482
483
484 DOUBLE PRECISION Flux
485
486
487 DOUBLE PRECISION EAST_DC
488 DOUBLE PRECISION WEST_DC
489 DOUBLE PRECISION NORTH_DC
490 DOUBLE PRECISION SOUTH_DC
491 DOUBLE PRECISION TOP_DC
492 DOUBLE PRECISION BOTTOM_DC
493
494
495
496 DOUBLE PRECISION :: AW,HW,VELW
497
498
499
500
501 call lock_tmp4_array
502 call lock_tmp_array
503 call lock_xsi_array
504
505
506
507
508
509 IF ( FPFOI ) THEN
510 Do IJK = ijkstart3, ijkend3
511 I = I_OF(IJK)
512 J = J_OF(IJK)
513 K = K_OF(IJK)
514 IJK4 = funijk3(I,J,K)
515 TMP4(IJK4) = W_G(IJK)
516 ENDDO
517 CALL send_recv3(tmp4)
518 ENDIF
519
520
521
522 DO IJK = ijkstart3, ijkend3
523 K = K_OF(IJK)
524 IJKT = TOP_OF(IJK)
525 IJKP = KP_OF(IJK)
526
527
528
529
530
531 IF(CUT_W_TREATMENT_AT(IJK)) THEN
532 U(IJK) = (Theta_W_be(IJK) * U_G(IJK) + Theta_W_te(IJK) * U_G(IJKP))
533 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
534 U(IJK) = U(IJK) * AW
535 ELSE
536 (IJK) = AVG_Z(U_G(IJK),U_G(IJKP),K)
537 ENDIF
538
539
540
541 IF(CUT_W_TREATMENT_AT(IJK)) THEN
542 V(IJK) = (Theta_W_bn(IJK) * V_G(IJK) + Theta_W_tn(IJK) * V_G(IJKP))
543 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
544 V(IJK) = V(IJK) * AW
545 ELSE
546 (IJK) = AVG_Z(V_G(IJK),V_G(IJKP),K)
547 ENDIF
548
549
550
551 IF(CUT_W_TREATMENT_AT(IJK)) THEN
552 WW(IJK) = (Theta_Wt_bar(IJK) * W_G(IJK) + Theta_Wt(IJK) * W_G(IJKP))
553 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
554 WW(IJK) = WW(IJK) * AW
555 ELSE
556 (IJK) = AVG_Z_T(W_G(IJK),W_G(IJKP))
557 ENDIF
558
559
560
561 END DO
562
563
564 =0
565
566
567 CALL CALC_XSI (DISCRETIZE(5), W_G, U, V, WW, XSI_E, XSI_N,&
568 XSI_T,incr)
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584 DO IJK = ijkstart3, ijkend3
585
586 IF (FLOW_AT_T(IJK)) THEN
587 I = I_OF(IJK)
588 J = J_OF(IJK)
589 K = K_OF(IJK)
590 IPJK = IP_OF(IJK)
591 IMJK = IM_OF(IJK)
592 IJPK = JP_OF(IJK)
593 IJMK = JM_OF(IJK)
594 IJKP = KP_OF(IJK)
595 IJKM = KM_OF(IJK)
596 IJKN = NORTH_OF(IJK)
597 IJKT = TOP_OF(IJK)
598 IF (WALL_AT(IJK)) THEN
599 IJKC = IJKT
600 ELSE
601 IJKC = IJK
602 ENDIF
603 KP = KP1(K)
604 IJKE = EAST_OF(IJK)
605 IJKP = KP_OF(IJK)
606 IJKTN = NORTH_OF(IJKT)
607 IJKTE = EAST_OF(IJKT)
608
609
610 = IP_OF(IP_OF(IPJK))
611 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
612 IMMM = IM_OF(IM_OF(IMJK))
613 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
614 JPPP = JP_OF(JP_OF(IJPK))
615 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
616 JMMM = JM_OF(JM_OF(IJMK))
617 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
618 KPPP = KP_OF(KP_OF(IJKP))
619 KPPP4 = funijk3(K_OF(IPPP), J_OF(KPPP), K_OF(KPPP))
620 KMMM = KM_OF(KM_OF(IJKM))
621 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
622
623
624
625
626
627
628 IF(U(IJK) >= ZERO)THEN
629 MOM_LO = W_G(IJK)
630 IF ( FPFOI ) &
631 MOM_HO = FPFOI_OF(W_G(IPJK), W_G(IJK), &
632 W_G(IMJK), W_G(IM_OF(IMJK)))
633 ELSE
634 MOM_LO = W_G(IPJK)
635 IF ( FPFOI ) &
636 MOM_HO = FPFOI_OF(W_G(IJK), W_G(IPJK), &
637 W_G(IP_OF(IPJK)), TMP4(IPPP4))
638 ENDIF
639 IF (.NOT. FPFOI ) &
640 MOM_HO = XSI_E(IJK)*W_G(IPJK)+ &
641 (1.0-XSI_E(IJK))*W_G(IJK)
642
643
644
645 IF(CUT_W_TREATMENT_AT(IJK)) THEN
646 Flux = (Theta_W_be(IJK) * Flux_gE(IJK) + Theta_W_te(IJK) * Flux_gE(IJKP))
647 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
648 Flux = Flux * AW
649 ELSE
650 = HALF * (Flux_gE(IJK) + Flux_gE(IJKP))
651 ENDIF
652
653
654
655 = Flux*(MOM_LO-MOM_HO)
656
657
658
659 IF(V(IJK) >= ZERO)THEN
660 MOM_LO = W_G(IJK)
661 IF ( FPFOI ) &
662 MOM_HO = FPFOI_OF(W_G(IJPK), W_G(IJK), &
663 W_G(IJMK), W_G(JM_OF(IJMK)))
664 ELSE
665 MOM_LO = W_G(IJPK)
666 IF ( FPFOI ) &
667 MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJPK), &
668 W_G(JP_OF(IJPK)), TMP4(JPPP4))
669 ENDIF
670 IF (.NOT. FPFOI ) &
671 MOM_HO = XSI_N(IJK)*W_G(IJPK)+ &
672 (1.0-XSI_N(IJK))*W_G(IJK)
673
674
675
676 IF(CUT_W_TREATMENT_AT(IJK)) THEN
677 Flux = (Theta_W_bn(IJK) * Flux_gN(IJK) + Theta_W_tn(IJK) * Flux_gN(IJKP))
678 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
679 Flux = Flux * AW
680 ELSE
681 = HALF * (Flux_gN(IJK) + Flux_gN(IJKP))
682 ENDIF
683
684
685
686 = Flux*(MOM_LO-MOM_HO)
687
688
689
690 IF(WW(IJK) >= ZERO)THEN
691 MOM_LO = W_G(IJK)
692 IF ( FPFOI ) &
693 MOM_HO = FPFOI_OF(W_G(IJKP), W_G(IJK), &
694 W_G(IJKM), W_G(KM_OF(IJKM)))
695 ELSE
696 MOM_LO = W_G(IJKP)
697 IF ( FPFOI ) &
698 MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJKP), &
699 W_G(KP_OF(IJKP)), TMP4(KPPP4))
700 ENDIF
701 IF (.NOT. FPFOI ) &
702 MOM_HO = XSI_T(IJK)*W_G(IJKP)+ &
703 (1.0-XSI_T(IJK))*W_G(IJK)
704
705
706
707 IF(CUT_W_TREATMENT_AT(IJK)) THEN
708 Flux = (Theta_Wt_bar(IJK) * Flux_gT(IJK) + Theta_Wt(IJK) * Flux_gT(IJKP))
709 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
710 Flux = Flux * AW
711 ELSE
712 = HALF * (Flux_gT(IJK) + Flux_gT(IJKP))
713 ENDIF
714
715
716
717 = Flux*(MOM_LO-MOM_HO)
718
719
720
721 = IM_OF(IJK)
722 IM = IM1(I)
723 IJKW = WEST_OF(IJK)
724 IJKWT = TOP_OF(IJKW)
725 IMJKP = KP_OF(IMJK)
726 IF(U(IMJK) >= ZERO)THEN
727 MOM_LO = W_G(IMJK)
728 IF ( FPFOI ) &
729 MOM_HO = FPFOI_OF(W_G(IJK), W_G(IMJK), &
730 W_G(IM_OF(IMJK)), TMP4(IMMM4))
731 ELSE
732 MOM_LO = W_G(IJK)
733 IF ( FPFOI ) &
734 MOM_HO = FPFOI_OF(W_G(IMJK), W_G(IJK), &
735 W_G(IPJK), W_G(IP_OF(IPJK)))
736 ENDIF
737 IF (.NOT. FPFOI ) &
738 MOM_HO = XSI_E(IMJK)*W_G(IJK)+ &
739 (1.0-XSI_E(IMJK))*W_G(IMJK)
740
741
742
743 IF(CUT_W_TREATMENT_AT(IJK)) THEN
744 Flux = (Theta_W_be(IMJK) * Flux_gE(IMJK) + Theta_W_te(IMJK) * Flux_gE(IMJKP))
745 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IMJK),AW,HW,VELW)
746 Flux = Flux * AW
747 ELSE
748 = HALF * (Flux_gE(IMJK) + Flux_gE(IMJKP))
749 ENDIF
750
751
752
753 = Flux*(MOM_LO-MOM_HO)
754
755
756
757 = JM_OF(IJK)
758 JM = JM1(J)
759 IJMKP = KP_OF(IJMK)
760 IJKS = SOUTH_OF(IJK)
761 IJKST = TOP_OF(IJKS)
762 IF(V(IJMK) >= ZERO)THEN
763 MOM_LO = W_G(IJMK)
764 IF ( FPFOI ) &
765 MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJMK), &
766 W_G(JM_OF(IJMK)), TMP4(JMMM4))
767 ELSE
768 MOM_LO = W_G(IJK)
769 IF ( FPFOI ) &
770 MOM_HO = FPFOI_OF(W_G(IJMK), W_G(IJK), &
771 W_G(IJPK), W_G(JP_OF(IJPK)))
772 ENDIF
773 IF (.NOT. FPFOI ) &
774 MOM_HO = XSI_N(IJMK)*W_G(IJK)+ &
775 (1.0-XSI_N(IJMK))*W_G(IJMK)
776
777
778
779 IF(CUT_W_TREATMENT_AT(IJK)) THEN
780 Flux = (Theta_W_bn(IJMK) * Flux_gN(IJMK) + Theta_W_tn(IJMK) * Flux_gN(IJMKP))
781 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJMK),AW,HW,VELW)
782 Flux = Flux * AW
783 ELSE
784 = HALF * (Flux_gN(IJMK) + Flux_gN(IJMKP))
785 ENDIF
786
787
788
789 = Flux*(MOM_LO-MOM_HO)
790
791
792
793 = KM_OF(IJK)
794 KM = KM1(K)
795 IJKB = BOTTOM_OF(IJK)
796 IF(WW(IJK) >= ZERO)THEN
797 MOM_LO = W_G(IJKM)
798 IF ( FPFOI ) &
799 MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJKM), &
800 W_G(KM_OF(IJKM)), TMP4(KMMM4))
801 ELSE
802 MOM_LO = W_G(IJK)
803 IF ( FPFOI ) &
804 MOM_HO = FPFOI_OF(W_G(IJKM), W_G(IJK), &
805 W_G(IJKP), W_G(KP_OF(IJKP)))
806 ENDIF
807 IF (.NOT. FPFOI ) &
808 MOM_HO = XSI_T(IJKM)*W_G(IJK)+ &
809 (1.0-XSI_T(IJKM))*W_G(IJKM)
810
811
812
813 IF(CUT_W_TREATMENT_AT(IJK)) THEN
814 Flux = (Theta_Wt_bar(IJKM) * Flux_gT(IJKM) + Theta_Wt(IJKM) * Flux_gT(IJK))
815 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJKM),AW,HW,VELW)
816 Flux = Flux * AW
817 ELSE
818 = HALF * (Flux_gT(IJKM) + Flux_gT(IJK))
819 ENDIF
820
821
822
823 = Flux*(MOM_LO-MOM_HO)
824
825
826
827 (IJK) = B_M(IJK)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC&
828 +BOTTOM_DC-TOP_DC
829
830 ENDIF
831 END DO
832
833 call unlock_tmp4_array
834 call unlock_tmp_array
835 call unlock_xsi_array
836
837 RETURN
838 END SUBROUTINE STORE_A_W_GDC
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868 SUBROUTINE STORE_A_W_G1(A_W_G)
869
870
871
872
873
874
875
876
877 USE param
878 USE param1
879 USE parallel
880 USE matrix
881 USE geometry
882 USE indices
883 USE run
884 USE visc_g
885 USE toleranc
886 USE physprop
887 USE fldvar
888 USE output
889 USE vshear
890 USE xsi
891 USE xsi_array
892 Use tmp_array, U => Array1, V => Array2, WW => Array3
893 USE compar
894 USE mflux
895 USE fun_avg
896 USE functions
897
898
899
900 USE cutcell
901
902
903
904
905 IMPLICIT NONE
906
907
908
909
910
911
912
913
914 INTEGER I, J, K, IPJK, IJPK, IJKN, IJKC, KP, IJKE,&
915 IJKTE, IJKP, IJKT, IJKTN, IJK
916 INTEGER IMJK, IM, IJKW, IJKWT, IMJKP
917 INTEGER IJMK, JM, IJMKP, IJKS, IJKST
918 INTEGER IJKM, KM, IJKB
919
920
921 INTEGER incr
922
923
924
925 DOUBLE PRECISION Flux
926
927
928 DOUBLE PRECISION D_f
929
930
931 DOUBLE PRECISION A_W_g(DIMENSION_3, -3:3)
932
933
934
935
936
937
938
939
940
941 DOUBLE PRECISION :: AW,HW,VELW
942
943
944
945
946 call lock_tmp_array
947 call lock_xsi_array
948
949
950
951
952
953
954 DO IJK = ijkstart3, ijkend3
955 K = K_OF(IJK)
956 IJKT = TOP_OF(IJK)
957 IJKP = KP_OF(IJK)
958
959
960
961
962
963 IF(CUT_W_TREATMENT_AT(IJK)) THEN
964
965 (IJK) = (Theta_W_be(IJK) * U_G(IJK) + Theta_W_te(IJK) * U_G(IJKP))
966 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
967 U(IJK) = U(IJK) * AW
968
969 (IJK) = (Theta_W_bn(IJK) * V_G(IJK) + Theta_W_tn(IJK) * V_G(IJKP))
970 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
971 V(IJK) = V(IJK) * AW
972
973 (IJK) = (Theta_Wt_bar(IJK) * W_G(IJK) + Theta_Wt(IJK) * W_G(IJKP))
974 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
975 WW(IJK) = WW(IJK) * AW
976 ELSE
977 (IJK) = AVG_Z(U_G(IJK),U_G(IJKP),K)
978 V(IJK) = AVG_Z(V_G(IJK),V_G(IJKP),K)
979 WW(IJK) = AVG_Z_T(W_G(IJK),W_G(IJKP))
980 ENDIF
981
982
983
984
985 END DO
986
987
988 =0
989
990
991 CALL CALC_XSI (DISCRETIZE(5), W_G, U, V, WW, XSI_E, XSI_N,&
992 XSI_T,incr)
993
994
995
996 IF (SHEAR) THEN
997
998 DO IJK = ijkstart3, ijkend3
999 IF (FLUID_AT(IJK)) THEN
1000 V(IJK)=V(IJK)+VSH(IJK)
1001 END IF
1002 END DO
1003 END IF
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019 DO IJK = ijkstart3, ijkend3
1020
1021 IF (FLOW_AT_T(IJK)) THEN
1022 I = I_OF(IJK)
1023 J = J_OF(IJK)
1024 K = K_OF(IJK)
1025 IPJK = IP_OF(IJK)
1026 IJPK = JP_OF(IJK)
1027 IJKN = NORTH_OF(IJK)
1028 IJKT = TOP_OF(IJK)
1029 IF (WALL_AT(IJK)) THEN
1030 IJKC = IJKT
1031 ELSE
1032 IJKC = IJK
1033 ENDIF
1034 KP = KP1(K)
1035 IJKE = EAST_OF(IJK)
1036 IJKP = KP_OF(IJK)
1037 IJKTN = NORTH_OF(IJKT)
1038 IJKTE = EAST_OF(IJKT)
1039
1040
1041
1042
1043
1044 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1045 Flux = (Theta_W_be(IJK) * Flux_gE(IJK) + Theta_W_te(IJK) * Flux_gE(IJKP))
1046 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
1047 Flux = Flux * AW
1048 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
1049 ),MU_GT(IJKTE),I),K)*ONEoDX_E_W(IJK)*AYZ_W(IJK)
1050 ELSE
1051 = HALF * (Flux_gE(IJK) + Flux_gE(IJKP))
1052 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
1053 ),MU_GT(IJKTE),I),K)*ODX_E(I)*AYZ_W(IJK)
1054 ENDIF
1055
1056
1057
1058
1059 (IJK,E) = D_F - XSI_E(IJK)*Flux
1060
1061 (IPJK,W) = D_F + (ONE - XSI_E(IJK))*Flux
1062
1063
1064
1065
1066
1067 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1068 Flux = (Theta_W_bn(IJK) * Flux_gN(IJK) + Theta_W_tn(IJK) * Flux_gN(IJKP))
1069 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
1070 Flux = Flux * AW
1071 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
1072 ),MU_GT(IJKTN),J),K)*ONEoDY_N_W(IJK)*AXZ_W(IJK)
1073 ELSE
1074 = HALF * (Flux_gN(IJK) + Flux_gN(IJKP))
1075 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
1076 ),MU_GT(IJKTN),J),K)*ODY_N(J)*AXZ_W(IJK)
1077 ENDIF
1078
1079
1080
1081
1082 (IJK,N) = D_F - XSI_N(IJK)*Flux
1083
1084 (IJPK,S) = D_F + (ONE - XSI_N(IJK))*Flux
1085
1086
1087
1088
1089
1090 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1091 Flux = (Theta_Wt_bar(IJK) * Flux_gT(IJK) + Theta_Wt(IJK) * Flux_gT(IJKP))
1092 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
1093 Flux = Flux * AW
1094 D_F = MU_GT(IJKT)*ONEoDZ_T_W(IJK)*AXY_W(IJK)
1095 ELSE
1096 = HALF * (Flux_gT(IJK) + Flux_gT(IJKP))
1097 D_F = MU_GT(IJKT)*OX(I)*ODZ(KP)*AXY_W(IJK)
1098 ENDIF
1099
1100
1101
1102 (IJK,T) = D_F - XSI_T(IJK)*Flux
1103 A_W_G(IJKP,B) = D_F + (ONE - XSI_T(IJK))*Flux
1104
1105
1106 = IM_OF(IJK)
1107 IF (.NOT.FLOW_AT_T(IMJK)) THEN
1108 IM = IM1(I)
1109 IJKW = WEST_OF(IJK)
1110 IJKWT = TOP_OF(IJKW)
1111 IMJKP = KP_OF(IMJK)
1112
1113
1114
1115 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1116 Flux = (Theta_W_be(IMJK) * Flux_gE(IMJK) + Theta_W_te(IMJK) * Flux_gE(IMJKP))
1117 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IMJK),AW,HW,VELW)
1118 Flux = Flux * AW
1119 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
1120 IJKWT),MU_GT(IJKT),IM),K)*ONEoDX_E_W(IMJK)*AYZ_W(IMJK)
1121 ELSE
1122 = HALF * (Flux_gE(IMJK) + Flux_gE(IMJKP))
1123 D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
1124 IJKWT),MU_GT(IJKT),IM),K)*ODX_E(IM)*AYZ_W(IMJK)
1125 ENDIF
1126
1127
1128
1129 (IJK,W) = D_F + (ONE - XSI_E(IMJK))*Flux
1130 ENDIF
1131
1132
1133 = JM_OF(IJK)
1134 IF (.NOT.FLOW_AT_T(IJMK)) THEN
1135 JM = JM1(J)
1136 IJMKP = KP_OF(IJMK)
1137 IJKS = SOUTH_OF(IJK)
1138 IJKST = TOP_OF(IJKS)
1139
1140
1141
1142 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1143 Flux = (Theta_W_bn(IJMK) * Flux_gN(IJMK) + Theta_W_tn(IJMK) * Flux_gN(IJMKP))
1144 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJMK),AW,HW,VELW)
1145 Flux = Flux * AW
1146 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
1147 IJKST),MU_GT(IJKT),JM),K)*ONEoDY_N_W(IJMK)*AXZ_W(IJMK)
1148 ELSE
1149 = HALF * (Flux_gN(IJMK) + Flux_gN(IJMKP))
1150 D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
1151 IJKST),MU_GT(IJKT),JM),K)*ODY_N(JM)*AXZ_W(IJMK)
1152 ENDIF
1153
1154
1155
1156 (IJK,S) = D_F + (ONE - XSI_N(IJMK))*Flux
1157 ENDIF
1158
1159
1160 = KM_OF(IJK)
1161 IF (.NOT.FLOW_AT_T(IJKM)) THEN
1162 KM = KM1(K)
1163 IJKB = BOTTOM_OF(IJK)
1164
1165
1166
1167 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1168 Flux = (Theta_Wt_bar(IJKM) * Flux_gT(IJKM) + Theta_Wt(IJKM) * Flux_gT(IJK))
1169 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJKM),AW,HW,VELW)
1170 Flux = Flux * AW
1171 D_F = MU_GT(IJK)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
1172 ELSE
1173 = HALF * (Flux_gT(IJKM) + Flux_gT(IJK))
1174 D_F = MU_GT(IJK)*OX(I)*ODZ(K)*AXY_W(IJKM)
1175 ENDIF
1176
1177
1178
1179 (IJK,B) = D_F + (ONE - XSI_T(IJKM))*Flux
1180 ENDIF
1181 ENDIF
1182 END DO
1183
1184 call unlock_tmp_array
1185 call unlock_xsi_array
1186
1187 RETURN
1188 END SUBROUTINE STORE_A_W_G1
1189
1190
1191
1192
1193
1194