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