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

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