File: /nfs/home/0/users/jenkins/mfix.git/model/conv_dif_w_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CONV_DIF_W_s                                            C
4     !  Purpose: Determine convection diffusion terms for W_s momentum eqs  C
5     !           The off-diagonal coefficients calculated here must be      C
6     !           positive. The center coefficient and the source vector     C
7     !           are negative;                                              C
8     !           See source_w_s                                             C
9     !                                                                      C
10     !  Author: M. Syamlal                                 Date: 24-DEC-96  C
11     !  Reviewer:                                          Date:            C
12     !                                                                      C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^CC
15     
16           SUBROUTINE CONV_DIF_W_S(A_M, B_M, IER)
17     
18     !-----------------------------------------------
19     ! Modules
20     !-----------------------------------------------
21     ! maximum number of computational cells, number of solids phases
22           USE param, only: dimension_3, dimension_m
23     
24     ! kinetic theories
25           USE run, only: kt_type_enum
26           USE run, only: ghd_2007
27     ! run time flag for deferred correction
28           USE run, only: def_cor
29     ! run time flag to solve z momentum equation
30           USE run, only: momentum_z_eq
31     ! discretization scheme for indicated equation
32           USE run, only: discretize
33     
34     ! number of solids phases
35           USE physprop, only: mmax
36     
37     ! solids phase viscosity
38           USE visc_s, only: mu_s
39     
40           IMPLICIT NONE
41     !-----------------------------------------------
42     ! Dummy Arguments
43     !-----------------------------------------------
44     ! Septadiagonal matrix A_m
45           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
46     
47     ! Vector b_m
48           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
49     ! Error index
50           INTEGER, INTENT(INOUT) :: IER
51     !-----------------------------------------------
52     ! Local Variables
53     !-----------------------------------------------
54     ! Solids phase index
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     ! IF DEFERRED CORRECTION IS TO BE USED TO SOLVE W_S
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     ! NO DEFERRED CORRECTION IS TO BE USED TO SOLVE FOR W_S
72                     IF (DISCRETIZE(5) == 0) THEN         ! 0 & 1 => FOUP
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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
89     !                                                                      C
90     !  Subroutine: STORE_A_W_s0                                            C
91     !  Purpose: Determine convection diffusion terms for W_s momentum eqs  C
92     !           The off-diagonal coefficients calculated here must be      C
93     !           positive. The center coefficient and the source vector     C
94     !           are negative. FOUP                                         C
95     !           See source_w_s                                             C
96     !                                                                      C
97     !  Author: M. Syamlal                                 Date: 7-Jun-96   C
98     !  Reviewer:                                          Date:            C
99     !                                                                      C
100     !  Revision Number: 1                                                  C
101     !  Purpose: To incorporate Cartesian grid modifications                C
102     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
103     !                                                                      C
104     !                                                                      C
105     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
106     
107           SUBROUTINE STORE_A_W_S0(A_W_S, M, IER)
108     
109     !-----------------------------------------------
110     ! Modules
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     ! Dummy arguments
132     !-----------------------------------------------
133     ! Solids phase index
134           INTEGER, INTENT(IN) :: M
135     ! Septadiagonal matrix A_W_s
136           DOUBLE PRECISION, INTENT(INOUT) :: A_W_s(DIMENSION_3, -3:3, M:M)
137     ! Error index
138           INTEGER, INTENT(INOUT) :: IER
139     !-----------------------------------------------
140     ! Local variables
141     !-----------------------------------------------
142     ! Indices
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     ! Face mass flux
149           DOUBLE PRECISION :: Flux
150     ! Diffusion parameter
151           DOUBLE PRECISION :: D_f
152     ! for cartesian grid:
153           DOUBLE PRECISION :: AW,HW,VELW
154     !-----------------------------------------------
155     
156     ! Calculate convection-diffusion fluxes through each of the faces
157     
158     !!$omp      parallel do         &
159     !!$omp&     private(IJK,  I,  J, K, IPJK, IJPK, IJKN, IJKC, KP, &
160     !!$omp&             IJKE, IJKTE, IJKP, IJKT, IJKTN, D_f,        &
161     !!$omp&             IMJK, IM, IJKW, IJKWT, IMJKP,       &
162     !!$omp&             IJMK, JM, IJMKP, IJKS, IJKST,       &
163     !!$omp&             IJKM, KM, IJKB)
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     ! East face (i+1/2, j, k+1/2)
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     ! North face (i, j+1/2, k+1/2)
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     ! Top face (i, j, k+1)
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     ! West face (i-1/2, j, k+1/2)
254                 IMJK = 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     ! South face (i, j-1/2, k+1/2)
283                 IJMK = 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     ! Bottom face (i, j, k)
312                 IJKM = 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   ! end if (do_k)
333     
334              ENDIF   ! end if (flow_at_t)
335           ENDDO   ! end do ijk
336     
337     
338           RETURN
339           END SUBROUTINE STORE_A_W_S0
340     
341     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
342     !                                                                      C
343     !  Subroutine: STORE_A_W_sdc                                           C
344     !  Purpose: To use deferred correction method to solve the w-momentum  C
345     !           equation. This method combines first order upwind and a    C
346     !           user specified higher order method                         C
347     !                                                                      C
348     !  Author: C. GUENTHER                                Date: 8-APR-99   C
349     !  Reviewer:                                          Date:            C
350     !                                                                      C
351     !  Revision Number: 1                                                  C
352     !  Purpose: To incorporate Cartesian grid modifications                C
353     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
354     !                                                                      C
355     !                                                                      C
356     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
357     
358           SUBROUTINE STORE_A_W_SDC(A_W_S, M, B_M, IER)
359     
360     !-----------------------------------------------
361     ! Modules
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     ! Dummy arguments
390     !-----------------------------------------------
391     ! Solids phase index
392           INTEGER, INTENT(IN) :: M
393     ! Septadiagonal matrix A_W_s
394           DOUBLE PRECISION, INTENT(INOUT) :: A_W_s(DIMENSION_3, -3:3, M:M)
395     ! Vector b_m
396           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
397     ! Error index
398           INTEGER, INTENT(INOUT) :: IER
399     !-----------------------------------------------
400     ! Local variables
401     !-----------------------------------------------
402     ! Indices
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     ! indicator for shear
411           INTEGER :: incr
412     ! Diffusion parameter
413           DOUBLE PRECISION :: D_f
414     ! Deferred correction contribution from high order method
415           DOUBLE PRECISION :: MOM_HO
416     ! low order approximation
417           DOUBLE PRECISION :: MOM_LO
418     ! convection factor at the face
419           DOUBLE PRECISION :: Flux
420     ! deferred correction contributions from each face
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     ! for cartesian grid:
429           DOUBLE PRECISION :: AW,HW,VELW
430     
431     ! temporary use of global arrays:
432     ! array1 (locally u)
433     ! the x directional velocity
434     !      DOUBLE PRECISION :: U(DIMENSION_3)
435     ! array2 (locally v)
436     ! the y directional velocity
437     !      DOUBLE PRECISION :: V(DIMENSION_3)
438     ! array3 (locally ww)
439     ! the z directional velocity
440     !      DOUBLE PRECISION :: WW(DIMENSION_3)
441     
442     !-----------------------------------------------
443     
444           call lock_tmp4_array
445           call lock_tmp_array   ! locks array1, array2, array3 (locally u, v, ww)
446           call lock_xsi_array
447     
448     ! Calculate convection factors
449     ! ---------------------------------------------------------------->>>
450     ! Send recv the third ghost layer
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     !!$omp parallel do private(IJK,K,IJKT,IJKP )
464           DO IJK = ijkstart3, ijkend3
465              K = K_OF(IJK)
466              IJKT = TOP_OF(IJK)
467              IJKP = KP_OF(IJK)
468     
469     ! East face (i+1/2, j, k+1/2)
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     ! North face (i, j+1/2, k+1/2)
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     ! Top face (i, j, k+1)
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     ! shear
504           incr=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     ! Calculate convection-diffusion fluxes through each of the faces
511     ! ---------------------------------------------------------------->>>
512     
513     !!!$omp      parallel do        &
514     !!!$omp&     private( I,  J, K, IPJK, IJPK, IJKN, IJKC, KP,     &
515     !!!$omp&             IJKE, IJKTE, IJKP, IJKT, IJKTN, IJK,  D_f, &
516     !!!$omp&             IMJK, IM, IJKW, IJKWT, IMJKP,      &
517     !!!$omp&             IJMK, JM, IJMKP, IJKS, IJKST,      &
518     !!!$omp&             IJKM, KM, IJKB, &
519     !!!$omp&              MOM_HO, MOM_LO, EAST_DC,WEST_DC,NORTH_DC,&
520     !!!$omp&              SOUTH_DC, TOP_DC,BOTTOM_DC)
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     ! Third Ghost layer information
545                 IPPP  = 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     ! DEFERRED CORRECTION CONTRIBUTION AT THE East face (i+1/2, j, k+1/2)
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     ! DEFERRED CORRECTION CONTRIBUTION AT THE North face (i, j+1/2, k+1/2)
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     ! DEFERRED CORRECTION CONTRIBUTION AT THE Top face (i, j, k+1)
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     ! DEFERRED CORRECTION CONTRIBUTION AT THE West face (i-1/2, j, k+1/2)
638                 IMJK = 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     ! CORRECTION CONTRIBUTION AT THE South face (i, j-1/2, k+1/2)
669                 IJMK = 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     ! DEFERRED CORRECTION CONTRIBUTION AT THE Bottom face (i, j, k)
700                 IJKM = 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     ! CONTRIBUTION DUE TO DEFERRED CORRECTION
729                 B_M(IJK,M) = B_M(IJK,M)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC+&
730                              BOTTOM_DC-TOP_DC
731     
732              ENDIF   ! end if flow_at_t
733           ENDDO   ! end do ijk
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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
744     !                                                                      C
745     !  Subroutine: STORE_A_W_s1                                            C
746     !  Purpose: Determine convection diffusion terms for W_s momentum eqs  C
747     !           The off-diagonal coefficients calculated here must be      C
748     !           positive. The center coefficient and the source vector     C
749     !           are negative. Higher order                                 C
750     !  See source_w_s                                                      C
751     !                                                                      C
752     !  Author: M. Syamlal                                 Date: 19-MAR-97  C
753     !  Reviewer:                                          Date:            C
754     !                                                                      C
755     !  Revision Number: 1                                                  C
756     !  Purpose: To incorporate Cartesian grid modifications                C
757     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
758     !                                                                      C
759     !                                                                      C
760     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
761     
762           SUBROUTINE STORE_A_W_S1(A_W_S, M)
763     
764     !-----------------------------------------------
765     ! Modules
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     ! Dummy arguments
791     !-----------------------------------------------
792     ! Solids phase
793           INTEGER, INTENT(IN) :: M
794     ! Septadiagonal matrix A_W_s
795           DOUBLE PRECISION, INTENT(INOUT) :: A_W_s(DIMENSION_3, -3:3, M:M)
796     !-----------------------------------------------
797     ! Local variables
798     !-----------------------------------------------
799     ! Indices
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     ! indicator for shear
806           INTEGER :: incr
807     ! Face mass flux
808           DOUBLE PRECISION :: Flux
809     ! Diffusion parameter
810           DOUBLE PRECISION :: D_f
811     ! for cartesian grid:
812           DOUBLE PRECISION :: AW,HW,VELW
813     
814     ! temporary use of global arrays:
815     ! array1 (locally u)
816     ! the x directional velocity
817     !      DOUBLE PRECISION :: U(DIMENSION_3)
818     ! array2 (locally v)
819     ! the y directional velocity
820     !      DOUBLE PRECISION :: V(DIMENSION_3)
821     ! array3 (locally ww)
822     ! the z directional velocity
823     !      DOUBLE PRECISION :: WW(DIMENSION_3)
824     
825     !-----------------------------------------------'
826     
827           call lock_tmp_array
828           call lock_xsi_array
829     
830     ! Calculate convection factors
831     ! ---------------------------------------------------------------->>>
832     !!$omp parallel do private(IJK,K,IJKT,IJKP )
833           DO IJK = ijkstart3, ijkend3
834              K = K_OF(IJK)
835              IJKT = TOP_OF(IJK)
836              IJKP = KP_OF(IJK)
837     
838     ! East face (i+1/2, j, k+1/2)
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     ! North face (i, j+1/2, k+1/2)
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     ! Top face (i, j, k+1)
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   ! end do ijk
871     
872     ! shear indicator
873           incr=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     ! shear: update to true velocity
880           IF (SHEAR) THEN
881     !!!$omp      parallel do private(IJK)
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     ! Calculate convection-diffusion fluxes through each of the faces
891     ! ---------------------------------------------------------------->>>
892     !!$omp      parallel do         &
893     !!$omp&     private( I,  J, K, IPJK, IJPK, IJKN, IJKC, KP,      &
894     !!$omp&             IJKE, IJKTE, IJKP, IJKT, IJKTN, IJK,  D_f,  &
895     !!$omp&             IMJK, IM, IJKW, IJKWT, IMJKP,       &
896     !!$omp&             IJMK, JM, IJMKP, IJKS, IJKST,       &
897     !!$omp&             IJKM, KM, IJKB)
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     ! East face (i+1/2, j, k+1/2)
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     ! North face (i, j+1/2, k+1/2)
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     ! Top face (i, j, k+1)
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     ! West face (i-1/2, j, k+1/2)
976                 IMJK = 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     ! South face (i, j-1/2, k+1/2)
1002                 IJMK = 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     ! Bottom face (i, j, k)
1028                 IJKM = 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   ! end if flow_at_t
1047           ENDDO   ! end do ijk
1048     
1049     
1050           call unlock_tmp_array
1051           call unlock_xsi_array
1052     
1053           RETURN
1054           END SUBROUTINE STORE_A_W_S1
1055