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