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

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