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