source_w

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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_W_s                                              C
4     !  Purpose: Determine source terms for W_g momentum eq. The terms      C
5     !     appear in the center coefficient and RHS vector. The center      C
6     !     coefficient and source vector are negative.  The off-diagonal    C
7     !     coefficients are positive.                                       C
8     !     The drag terms are excluded from the source at this stage.       C
9     !                                                                      C
10     !                                                                      C
11     !  Author: M. Syamlal                                 Date: 17-JUN-96  C
12     !  Reviewer:                                          Date:            C
13     !                                                                      C
14     !  Revision Number: 1                                                  C
15     !  Purpose: Allow for partial-slip boundary conditions proposed by     C
16     !           by Johnson & Jackson (1987) if the Granular Temperature    C
17     !           equation is used.                                          C
18     !  Author: K. Agrawal, Princeton University           Date: 24-JAN-98  C
19     !  Reviewer:                                          Date: dd-mmm-yy  C
20     !                                                                      C
21     !  Revision Number: 2                                                  C
22     !  Purpose: To incorporate Cartesian grid modifications                C
23     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
24     !                                                                      C
25     !                                                                      C
26     !  Literature/Document References:                                     C
27     !                                                                      C
28     !                                                                      C
29     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
30     
31           SUBROUTINE SOURCE_W_S(A_M, B_M, IER)
32     
33     !-----------------------------------------------
34     ! Modules
35     !-----------------------------------------------
36           USE param
37           USE param1
38           USE parallel
39           USE matrix
40           USE scales
41           USE constant
42           USE physprop
43           USE fldvar
44           USE visc_s
45           USE rxns
46           USE run
47           USE toleranc
48           USE geometry
49           USE indices
50           USE is
51           USE tau_s
52           USE bc
53           USE compar
54           USE sendrecv
55           use kintheory
56           USE ghdtheory
57           USE drag
58           USE cutcell
59           USE quadric
60           USE mms
61           USE bodyforce
62           USE fun_avg
63           USE functions
64     
65           IMPLICIT NONE
66     !-----------------------------------------------
67     ! Dummy Arguments
68     !-----------------------------------------------
69     ! Septadiagonal matrix A_m
70           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
71     ! Vector b_m
72           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
73     ! Error index
74           INTEGER, INTENT(INOUT) :: IER
75     !-----------------------------------------------
76     ! Local Variables
77     !-----------------------------------------------
78     ! Indices
79           INTEGER :: I, J, K, IJK, IJKT, IMJK, IJMK, IJKM, IJKP, IMJKP
80           INTEGER :: IJKE, IJKW, IJKTE, IJKTW, IM, IPJK, IJPK, IJMKP
81     ! Phase index
82           INTEGER ::  M, MM, L
83     ! Internal surface
84           INTEGER :: ISV
85     ! Pressure at top cell
86           DOUBLE PRECISION :: PgT
87     ! Average volume fraction
88           DOUBLE PRECISION :: EPSA, EPStmp, epse, epsw, epsn, epss, &
89                               epst, epsb, epsMix, epsMixT
90           DOUBLE PRECISION :: SUM_EPS_CP
91     ! Average density
92           DOUBLE PRECISION :: ROPSA
93     ! Average density differences
94           DOUBLE PRECISION :: dro1, dro2, droa
95     ! Average quantities
96           DOUBLE PRECISION :: ugt, Cte, Ctw, EPMUoX
97     ! Source terms (Surface)
98           DOUBLE PRECISION :: Sdp, Sdps, Sxzb, Vxza, Vxzb
99     ! Source terms (Volumetric)
100           DOUBLE PRECISION :: V0, Vmt, Vbf, Vcoa, Vcob, Vmttmp
101     ! Source terms (Volumetric) for GHD theory
102           DOUBLE PRECISION :: Ghd_drag, avgRop
103     ! Source terms for HYS drag relation
104           DOUBLE PRECISION :: HYS_drag, avgDrag
105     ! virtual (added) mass
106           DOUBLE PRECISION :: F_vir, ROP_MA, Uge, Ugw, Vgb, Vgt, Wge, &
107                               Wgw, Wgn, Wgs, Wgt, Wgb
108     !-----------------------------------------------
109     
110           DO M = 1, MMAX
111             IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
112                (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
113     
114               IF (MOMENTUM_Z_EQ(M)) THEN
115     
116     
117     !$omp  parallel do default(shared)                                   &
118     !$omp  private(I, J, K, IJK, IMJK, IJMK, IJKM, IJKP, IPJK, IMJKP,    &
119     !$omp         IJPK, IJMKP, IJKT, IJKE, IJKW, IJKTE, IJKTW, IM,       &
120     !$omp         ISV, EPStmp, epsMix, epsMixT, EPSA, EPSE, EPSW, EPSN,  &
121     !$omp         EPSS, EPST, EPSB, SUM_EPS_CP, PGT, SDP, SDPS, ROPSA,   &
122     !$omp         V0, ROP_MA, MM, L, F_vir,                              &
123     !$omp         Uge, Ugw, Vgb, Vgt, Wge, Wgw, Wgn, Wgs, Wgt, Wgb,      &
124     !$omp         Ugt, VMTtmp, VMT, DRO1, DRO2, DROA, VBF, Ghd_drag,     &
125     !$omp         avgRop, HYS_drag, avgDrag, VCOA, VCOB, CTE, CTW,       &
126     !$omp         SXZB, VXZA, VXZB, EPMUOX)
127                 DO IJK = ijkstart3, ijkend3
128     
129     ! Skip walls where some values are undefined.
130                     IF(WALL_AT(IJK)) cycle
131     
132                     I = I_OF(IJK)
133                     J = J_OF(IJK)
134                     K = K_OF(IJK)
135                     IMJK = IM_OF(IJK)
136                     IJMK = JM_OF(IJK)
137                     IJKM = KM_OF(IJK)
138                     IJKP = KP_OF(IJK)
139                     IPJK = IP_OF(IJK)
140                     IMJKP = KP_OF(IMJK)
141                     IJPK = JP_OF(IJK)
142                     IJMKP = KP_OF(IJMK)
143                     IJKT = TOP_OF(IJK)
144     
145                     IF (KT_TYPE_ENUM == GHD_2007) THEN
146     ! with ghd theory, m = mmax
147                        EPStmp = ZERO
148                        epsMix = ZERO
149                        epsMixT= ZERO
150                        DO L = 1, SMAX
151                           EPStmp = EPStmp + AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
152                           epsMix  = epsMix  + EP_S(IJK,L) ! epsMix, epsMixT to be used for model B
153                           epsMixT = epsMixT + EP_S(IJKT,L)
154                           IF(IP_AT_T(IJK)) THEN
155                              W_S(IJK,L) = ZERO
156                           ELSEIF(SIP_AT_T(IJK)) THEN
157                              ISV = IS_ID_AT_T(IJK)
158                              W_S(IJK,L) = IS_VEL_S(ISV,L)
159                           ENDIF
160                        ENDDO
161                        EPSA = EPStmp
162                     ELSE
163                        EPSA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
164                     ENDIF
165     
166     ! Impermeable internal surface
167                     IF (IP_AT_T(IJK)) THEN
168                       A_M(IJK,E,M) = ZERO
169                       A_M(IJK,W,M) = ZERO
170                       A_M(IJK,N,M) = ZERO
171                       A_M(IJK,S,M) = ZERO
172                       A_M(IJK,T,M) = ZERO
173                       A_M(IJK,B,M) = ZERO
174                       A_M(IJK,0,M) = -ONE
175                       B_M(IJK,M) = ZERO
176     
177     ! Semi-permeable internal surface
178                     ELSEIF (SIP_AT_T(IJK)) THEN
179                       A_M(IJK,E,M) = ZERO
180                       A_M(IJK,W,M) = ZERO
181                       A_M(IJK,N,M) = ZERO
182                       A_M(IJK,S,M) = ZERO
183                       A_M(IJK,T,M) = ZERO
184                       A_M(IJK,B,M) = ZERO
185                       A_M(IJK,0,M) = -ONE
186                       ISV = IS_ID_AT_T(IJK)
187                       B_M(IJK,M) = -IS_VEL_S(ISV,M)
188     
189     ! Dilute flow
190                     ELSEIF (EPSA <= DIL_EP_S) THEN
191                       A_M(IJK,E,M) = ZERO
192                       A_M(IJK,W,M) = ZERO
193                       A_M(IJK,N,M) = ZERO
194                       A_M(IJK,S,M) = ZERO
195                       A_M(IJK,T,M) = ZERO
196                       A_M(IJK,B,M) = ZERO
197                       A_M(IJK,0,M) = -ONE
198                       B_M(IJK,M) = ZERO
199                       IF (KT_TYPE_ENUM == GHD_2007) THEN
200                           EPSw = ZERO
201                           EPSe = ZERO
202                           EPSn = ZERO
203                           EPSs = ZERO
204                           EPSt = ZERO
205                           EPSb = ZERO
206                           DO L = 1, SMAX
207                             EPSw = EPSw + EP_S(WEST_OF(IJK),L)
208                             EPSe = EPSe + EP_S(EAST_OF(IJK),L)
209                             EPSn = EPSn + EP_S(NORTH_OF(IJK),L)
210                             EPSs = EPSs + EP_S(SOUTH_OF(IJK),L)
211                             EPSt = EPSt + EP_S(TOP_OF(IJK),L)
212                             EPSb = EPSb + EP_S(BOTTOM_OF(IJK),L)
213                           ENDDO
214                       ELSE
215                           EPSw = EP_S(WEST_OF(IJK),M)
216                           EPSe = EP_S(EAST_OF(IJK),M)
217                           EPSn = EP_S(NORTH_OF(IJK),M)
218                           EPSs = EP_S(SOUTH_OF(IJK),M)
219                           EPSt = EP_S(TOP_OF(IJK),M)
220                           EPSb = EP_S(BOTTOM_OF(IJK),M)
221                       ENDIF
222     ! using the average boundary cell values to compute U_s (sof, Aug 23 2005)
223                       IF (EPSw > DIL_EP_S .AND. .NOT.IS_AT_E(IMJK)) A_M(IJK,W,M) = ONE
224                       IF (EPSe > DIL_EP_S .AND. .NOT.IS_AT_E(IJK)) A_M(IJK,E,M) = ONE
225                       IF (EPSs > DIL_EP_S .AND. .NOT.IS_AT_N(IJMK)) A_M(IJK,S,M) = ONE
226                       IF (EPSn > DIL_EP_S .AND. .NOT.IS_AT_N(IJK)) A_M(IJK,N,M) = ONE
227                       IF (EPSb > DIL_EP_S .AND. .NOT.IS_AT_T(IJKM)) A_M(IJK,B,M) = ONE
228                       IF (EPSt > DIL_EP_S .AND. .NOT.IS_AT_T(IJK)) A_M(IJK,T,M) = ONE
229                       IF((A_M(IJK,W,M)+A_M(IJK,E,M)+A_M(IJK,S,M)+A_M(IJK,N,M)+ &
230                         A_M(IJK,B,M)+A_M(IJK,T,M)) == ZERO) THEN
231                         B_M(IJK,M) = -W_S(IJK,M)
232                       ELSE
233                         A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+A_M(IJK,N,M)+ &
234                                          A_M(IJK,S,M)+A_M(IJK,T,M)+A_M(IJK,B,M))
235                       ENDIF
236     
237     ! Cartesian grid implementation
238                    ELSEIF (BLOCKED_W_CELL_AT(IJK)) THEN
239                       A_M(IJK,E,M) = ZERO
240                       A_M(IJK,W,M) = ZERO
241                       A_M(IJK,N,M) = ZERO
242                       A_M(IJK,S,M) = ZERO
243                       A_M(IJK,T,M) = ZERO
244                       A_M(IJK,B,M) = ZERO
245                       A_M(IJK,0,M) = -ONE
246                       B_M(IJK,M) = ZERO
247     
248     ! Normal case
249                     ELSE
250     
251     ! Surface forces:
252     ! Pressure term
253                       PGT = P_G(IJKT)
254                       IF (CYCLIC_Z_PD) THEN
255     ! CYCLIC_AT_T Flag is not set correctly in DMP and causes issues. This
256     ! is avoided by using the DMP cyclic map. The flags need fixed.
257     !                     IF (CYCLIC_AT_T(IJK)) PGT = P_G(IJKT) - DELP_Z
258                          IF (KMAP(K_OF(IJK)).EQ.KMAX1) PGT = P_G(IJKT) - DELP_Z
259                       ENDIF
260                       IF (MODEL_B) THEN
261                          SDP = ZERO
262                       ELSE
263                          IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
264                             SDP = -P_SCALE*EPSA*(PGT - P_G(IJK))*AXY(IJK)
265                          ELSE
266                             SDP = -P_SCALE*EPSA*(PGT * A_WPG_T(IJK) - P_G(IJK) * A_WPG_B(IJK) )
267                          ENDIF
268                       ENDIF
269     
270                       IF (CLOSE_PACKED(M)) THEN
271                          IF(SMAX > 1 .AND. KT_TYPE_ENUM /= GHD_2007) THEN
272                             SUM_EPS_CP=0.0
273                             DO MM=1,SMAX
274                               IF (CLOSE_PACKED(MM))&
275                                 SUM_EPS_CP=SUM_EPS_CP+AVG_Z(EP_S(IJK,MM),EP_S(IJKT,MM),K)
276                             ENDDO
277                             SUM_EPS_CP = Max(SUM_EPS_CP, small_number)
278                             SDPS = -((P_S(IJKT,M)-P_S(IJK,M))+(EPSA/SUM_EPS_CP)*&
279                                (P_STAR(IJKT)-P_STAR(IJK)))*AXY(IJK)
280                          ELSE
281                             IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
282                                SDPS =-((P_S(IJKT,M)-P_S(IJK,M))+(P_STAR(IJKT)-P_STAR(IJK)))*AXY(IJK)
283                             ELSE
284                                SDPS =-((P_S(IJKT,M)* A_WPG_T(IJK)-P_S(IJK,M)* A_WPG_B(IJK)) &
285                                      +(P_STAR(IJKT)* A_WPG_T(IJK)-P_STAR(IJK)* A_WPG_B(IJK)))
286                             ENDIF
287                          ENDIF
288                       ELSE
289                          IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
290                             SDPS = -(P_S(IJKT,M)-P_S(IJK,M))*AXY(IJK)
291                          ELSE
292                             SDPS = -(P_S(IJKT,M) * A_WPG_T(IJK) - P_S(IJK,M) * A_WPG_B(IJK))
293                          ENDIF
294                       ENDIF
295     
296     
297                       IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
298     ! Volumetric forces
299                          ROPSA = AVG_Z(ROP_S(IJK,M),ROP_S(IJKT,M),K)
300     ! Previous time step
301                          V0 = AVG_Z(ROP_SO(IJK,M),ROP_SO(IJKT,M),K)*ODT
302     ! Added mass implicit transient term {Cv eps rop_g dW/dt}
303                          IF(Added_Mass .AND. M == M_AM) THEN
304                            ROP_MA = AVG_Z(ROP_g(IJK)*EP_s(IJK,M),&
305                                           ROP_g(IJKT)*EP_s(IJKT,M),K)
306                            V0 = V0 + Cv * ROP_MA * ODT
307                          ENDIF
308                       ELSE
309     ! Volumetric forces
310                          ROPSA = (VOL(IJK)*ROP_S(IJK,M) + &
311                             VOL(IJKT)*ROP_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
312     ! Previous time step
313                          V0 = (VOL(IJK)*ROP_SO(IJK,M) + &
314                             VOL(IJKT)*ROP_SO(IJKT,M))*ODT/&
315                             (VOL(IJK) + VOL(IJKT))
316     ! Added mass implicit transient term {Cv eps rop_g dW/dt}
317                          IF(Added_Mass .AND. M == M_AM) THEN
318                             ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M) +&
319                                VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M))/&
320                                (VOL(IJK) + VOL(IJKT))
321                             V0 = V0 + Cv * ROP_MA * ODT
322                          ENDIF
323                       ENDIF
324     
325     ! VIRTUAL MASS SECTION (explicit terms)
326     ! adding transient term dWg/dt to virtual mass term
327                       F_vir = ZERO
328                       IF(Added_Mass .AND. M == M_AM .AND.&
329                          (.NOT.CUT_W_TREATMENT_AT(IJK))) THEN
330                          F_vir = ( (W_g(IJK) - W_gO(IJK)) )*ODT*VOL_W(IJK)
331     
332     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
333                           Wgb = AVG_Z_T(W_g(IJKM),W_g(IJK))
334                           Wgt = AVG_Z_T(W_g(IJK),W_g(IJKP))
335                           Uge = AVG_Z(U_g(IJK),U_g(IJKP),K)
336                           Ugw = AVG_Z(U_g(IMJK),U_g(IMJKP),K)
337                           Ugt = AVG_X_E(Ugw,Uge,IP1(I))
338                           Wge = AVG_X(W_g(IJK),W_g(IPJK),IP1(I))
339                           Wgw = AVG_X(W_g(IMJK),W_g(IJK),I)
340                           Vgb = AVG_Y_N(V_g(IJMK),V_g(IJK))
341                           Vgt = AVG_Y_N(V_g(IJMKP),V_g(IJKP))
342                           Wgs = AVG_Y(W_g(IJMK),W_g(IJK),JM1(J))
343                           Wgn = AVG_Y(W_g(IJK),W_g(IJPK),J)
344     
345     ! Coriolis force
346                          IF (CYLINDRICAL) F_vir = F_vir + &
347                             Ugt*W_g(IJK)*OX(I)
348     
349     ! adding convective terms (U dW/dx + V dW/dy + W dW/dz) to virtual mass.
350                          F_vir = F_vir + W_g(IJK)*OX(I) * &
351                             (Wgt - Wgb)*AXY(IJK) + Ugt*(Wge - Wgw)*AYZ(IJK) + &
352                             AVG_Z(Vgb,Vgt,K) * (Wgn - Wgs)*AXZ(IJK)
353                          F_vir = F_vir * Cv * ROP_MA
354                       ENDIF
355     
356     ! Interphase mass transfer
357                       IF (KT_TYPE_ENUM == GHD_2007) THEN
358                          VMTtmp = ZERO
359                          DO L = 1,SMAX
360                             VMTtmp = VMTtmp + AVG_Z(SUM_R_S(IJK,L),SUM_R_S(IJKT,L),K)
361                          ENDDO
362                          VMT = VMTtmp
363                       ELSE
364                          IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
365                             VMT = AVG_Z(SUM_R_S(IJK,M),SUM_R_S(IJKT,M),K)
366                          ELSE
367                             VMT = (VOL(IJK)*SUM_R_S(IJK,M) + VOL(IJKT)*&
368                                SUM_R_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
369                          ENDIF
370                       ENDIF
371     
372     ! Body force
373                       IF (MODEL_B) THEN
374                          IF (KT_TYPE_ENUM == GHD_2007) THEN
375                            DRO1 = ROP_S(IJK,M)  - RO_G(IJK) *epsMix
376                            DRO2 = ROP_S(IJKT,M) - RO_G(IJKT)*epsMixT
377                            DROA = AVG_Z(DRO1,DRO2,K)
378                            VBF = DROA*BFZ_S(IJK,M)
379                          ELSE
380                            DRO1 = (RO_S(IJK,M)-RO_G(IJK))*EP_S(IJK,M)
381                            DRO2 = (RO_S(IJK,M)-RO_G(IJKT))*EP_S(IJKT,M)
382                            DROA = AVG_Z(DRO1,DRO2,K)
383                            VBF = DROA*BFZ_S(IJK,M)
384                          ENDIF
385                       ELSE
386                          VBF = ROPSA*BFZ_S(IJK,M)
387                       ENDIF
388     
389     ! Additional force for GHD from darg force sum(beta_ig * Joi/rhop_i)
390                       Ghd_drag = ZERO
391                       IF (KT_TYPE_ENUM == GHD_2007) THEN
392                         DO L = 1,SMAX
393                           avgRop = AVG_Z(ROP_S(IJK,L),ROP_S(IJKT,L),K)
394                           if(avgRop > ZERO) Ghd_drag = Ghd_drag -&
395                                AVG_Z(F_GS(IJK,L),F_GS(IJKT,L),K) * JoiZ(IJK,L) / avgRop
396                         ENDDO
397                       ENDIF
398     
399     ! Additional force for HYS drag force, do not use with mixture GHD theory
400                       HYS_drag = ZERO
401                       IF (DRAG_TYPE_ENUM .EQ. HYS .AND. &
402                           KT_TYPE_ENUM /= GHD_2007) THEN
403                          DO L = 1,MMAX
404                             IF (L /= M) THEN
405                                avgDrag = AVG_Z(beta_ij(IJK,M,L),beta_ij(IJKT,M,L),K)
406                                HYS_drag = HYS_drag - avgDrag * (W_g(ijk) - W_s(IJK,L))
407                             ENDIF
408                          ENDDO
409                       ENDIF
410     
411     ! Special terms for cylindrical coordinates
412                       VCOA = ZERO
413                       VCOB = ZERO
414                       SXZB = ZERO
415                       VXZA = ZERO
416                       VXZB = ZERO
417                       CTE  = ZERO
418                       CTW  = ZERO
419                       IF (CYLINDRICAL) THEN
420     ! Coriolis force
421                          IMJK = IM_OF(IJK)
422                          IJKP = KP_OF(IJK)
423                          IMJKP = KP_OF(IMJK)
424                          UGT = AVG_Z(HALF*(U_S(IJK,M)+U_S(IMJK,M)),&
425                                      HALF*(U_S(IJKP,M)+U_S(IMJKP,M)),K)
426                          IF (UGT > ZERO) THEN
427                             VCOA = ROPSA*UGT*OX(I)
428                             VCOB = ZERO
429                             IF(Added_Mass .AND. M == M_AM) &
430                                VCOA = VCOA + Cv*ROP_MA*UGT*OX(I)
431                          ELSE
432                             VCOA = ZERO
433                             VCOB = -ROPSA*UGT*W_S(IJK,M)*OX(I)
434                             IF(Added_Mass .AND. M == M_AM) &
435                                VCOB = VCOB - Cv*ROP_MA*UGT*W_S(IJK,M)*OX(I)
436                          ENDIF
437     
438     ! Term from tau_xz: intergral of (1/x)*(d/dx)(x*mu*(-w/x))
439                          IJKE = EAST_OF(IJK)
440                          IJKW = WEST_OF(IJK)
441                          IJKTE = TOP_OF(IJKE)
442                          IJKTW = TOP_OF(IJKW)
443                          IM = IM1(I)
444                          IPJK = IP_OF(IJK)
445     
446                          CTE = HALF*AVG_Z_H(AVG_X_H(MU_S(IJK,M),MU_S(IJKE,M),I),&
447                             AVG_X_H(MU_S(IJKT,M),MU_S(IJKTE,M),I),K)*OX_E(I)*&
448                             AYZ_W(IJK)
449                          CTW = HALF*AVG_Z_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJK,M),IM),&
450                             AVG_X_H(MU_S(IJKTW,M),MU_S(IJKT,M),IM),K)*DY(J)*&
451                             (HALF*(DZ(K)+DZ(KP1(K))))
452                             !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
453     
454     ! (mu/x)*(dw/dx) part of tau_xz/x
455                          EPMUOX = AVG_Z(MU_S(IJK,M),MU_S(IJKT,M),K)*OX(I)
456                          VXZB = ZERO
457                          A_M(IJK,E,M) = A_M(IJK,E,M) + HALF*EPMUOX*ODX_E(I)*&
458                                VOL_W(IJK)
459                          A_M(IJK,W,M) = A_M(IJK,W,M) - HALF*EPMUOX*ODX_E(IM)*&
460                                VOL_W(IJK)
461     
462     ! -(mu/x)*(w/x) part of tau_xz/x
463                          VXZA = EPMUOX*OX(I)
464                       ELSE
465     
466                          VCOA = ZERO
467                          VCOB = ZERO
468                          SXZB = ZERO
469                          VXZA = ZERO
470                          VXZB = ZERO
471                       ENDIF
472     
473     ! Collect the terms
474                       A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+&
475                          A_M(IJK,N,M)+A_M(IJK,S,M)+A_M(IJK,T,M)+&
476                          A_M(IJK,B,M)+(V0+ZMAX(VMT)+VCOA+VXZA)*&
477                          VOL_W(IJK)+ CTE - CTW)
478                       A_M(IJK,E,M) = A_M(IJK,E,M) - CTE
479                       A_M(IJK,W,M) = A_M(IJK,W,M) + CTW
480     
481                       B_M(IJK,M) = B_m(IJK, M) - (SDP + SDPS + &
482                          TAU_W_S(IJK,M) + SXZB + F_vir + &
483                          ( (V0+ZMAX((-VMT)))*W_SO(IJK,M)+ &
484                          VBF + VCOB + VXZB + HYS_drag)*VOL_W(IJK) )
485     ! MMS Source term.
486                       IF(USE_MMS) B_M(IJK,M) = &
487                          B_M(IJK,M) - MMS_W_S_SRC(IJK)*VOL_W(IJK)
488     
489                       IF (KT_TYPE_ENUM == IA_2005) THEN
490                         B_M(IJK,M) = B_M(IJK,M) - KTMOM_W_S(IJK,M)
491                       ELSEIF (KT_TYPE_ENUM == GHD_2007) THEN
492                         B_M(IJK,M) = B_M(IJK,M) - Ghd_drag*VOL_W(IJK)
493                       ENDIF
494     
495                     ENDIF   ! end branching on cell type (sip/ip/dilute/block/else branches)
496                 ENDDO   ! end do loop over ijk
497     !$omp end parallel do
498     
499     
500     ! modifications for cartesian grid implementation
501                 IF(CARTESIAN_GRID) CALL CG_SOURCE_W_S(A_M, B_M, M, IER)
502     ! modifications for bc
503                 CALL SOURCE_W_S_BC (A_M, B_M, M)
504                 IF(CARTESIAN_GRID) CALL CG_SOURCE_W_S_BC(A_M, B_M, M, IER)
505     
506               ENDIF   ! end if (momentum_z_eq)
507             ENDIF   ! end if for ghd theory
508           ENDDO   ! end do loop over mmax
509     
510           RETURN
511           END SUBROUTINE SOURCE_W_S
512     
513     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
514     !                                                                      C
515     !  Subroutine: SOURCE_W_s_BC                                           C
516     !  Purpose: Determine source terms for W_s momentum eq. The terms      C
517     !     appear in the center coefficient and RHS vector.  The center     C
518     !     coefficient and source vector are negative.  The off-diagonal    C
519     !     coefficients are positive.                                       C
520     !     The drag terms are excluded from the source at this stage.       C
521     !                                                                      C
522     !  Author: M. Syamlal                                 Date: 17-JUN-96  C
523     !  Reviewer:                                          Date:            C
524     !                                                                      C
525     !                                                                      C
526     !  Comments: see source_w_g_bc for more detailed in-code comments      C
527     !                                                                      C
528     !                                                                      C
529     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
530     
531           SUBROUTINE SOURCE_W_S_BC(A_M, B_M, M)
532     
533     !-----------------------------------------------
534     ! Modules
535     !-----------------------------------------------
536           USE param
537           USE param1
538           USE parallel
539           USE matrix
540           USE scales
541           USE constant
542           USE physprop
543           USE fldvar
544           USE visc_s
545           USE rxns
546           USE run
547           USE toleranc
548           USE geometry
549           USE indices
550           USE is
551           USE tau_s
552           USE bc
553           USE output
554           USE compar
555           USE fun_avg
556           USE functions
557     
558           IMPLICIT NONE
559     !-----------------------------------------------
560     ! Dummy Arguments
561     !-----------------------------------------------
562     ! Septadiagonal matrix A_m
563           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
564     ! Vector b_m
565           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
566     ! Solids phase index
567           INTEGER, INTENT(IN) :: M
568     !-----------------------------------------------
569     ! Local Variables
570     !-----------------------------------------------
571     !  Boundary condition
572           INTEGER :: L
573     ! Indices
574           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
575                      IM, JM, IJKB, IJKM, IJKP
576     !-----------------------------------------------
577     
578     ! Setting the default boundary conditions
579           J1 = 1
580           DO K1 = kmin3,kmax3
581              DO I1 = imin3,imax3
582                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
583                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
584                 IJK = FUNIJK(I1,J1,K1)
585                 IF (NS_WALL_AT(IJK)) THEN
586     ! Setting the wall velocity to zero
587                    A_M(IJK,E,M) = ZERO
588                    A_M(IJK,W,M) = ZERO
589                    A_M(IJK,N,M) = -ONE
590                    A_M(IJK,S,M) = ZERO
591                    A_M(IJK,T,M) = ZERO
592                    A_M(IJK,B,M) = ZERO
593                    A_M(IJK,0,M) = -ONE
594                    B_M(IJK,M) = ZERO
595                 ELSEIF (FS_WALL_AT(IJK)) THEN
596     ! Setting the wall velocity equal to the adjacent fluid velocity
597                    A_M(IJK,E,M) = ZERO
598                    A_M(IJK,W,M) = ZERO
599                    A_M(IJK,N,M) = ONE
600                    A_M(IJK,S,M) = ZERO
601                    A_M(IJK,T,M) = ZERO
602                    A_M(IJK,B,M) = ZERO
603                    A_M(IJK,0,M) = -ONE
604                    B_M(IJK,M) = ZERO
605                 ENDIF
606              ENDDO
607           ENDDO
608     
609           J1 = JMAX2
610           DO K1 = kmin3, kmax3
611              DO I1 = imin3, imax3
612                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
613                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
614                 IJK = FUNIJK(I1,J1,K1)
615                 IF (NS_WALL_AT(IJK)) THEN
616                    A_M(IJK,E,M) = ZERO
617                    A_M(IJK,W,M) = ZERO
618                    A_M(IJK,N,M) = ZERO
619                    A_M(IJK,S,M) = -ONE
620                    A_M(IJK,T,M) = ZERO
621                    A_M(IJK,B,M) = ZERO
622                    A_M(IJK,0,M) = -ONE
623                    B_M(IJK,M) = ZERO
624                 ELSEIF (FS_WALL_AT(IJK)) THEN
625                    A_M(IJK,E,M) = ZERO
626                    A_M(IJK,W,M) = ZERO
627                    A_M(IJK,N,M) = ZERO
628                    A_M(IJK,S,M) = ONE
629                    A_M(IJK,T,M) = ZERO
630                    A_M(IJK,B,M) = ZERO
631                    A_M(IJK,0,M) = -ONE
632                    B_M(IJK,M) = ZERO
633                 ENDIF
634              ENDDO
635           ENDDO
636     
637           I1 = 1
638           DO K1 = kmin3, kmax3
639              DO J1 = jmin3, jmax3
640                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
641                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
642                 IJK = FUNIJK(I1,J1,K1)
643                 IF (NS_WALL_AT(IJK)) THEN
644                    A_M(IJK,E,M) = -ONE
645                    A_M(IJK,W,M) = ZERO
646                    A_M(IJK,N,M) = ZERO
647                    A_M(IJK,S,M) = ZERO
648                    A_M(IJK,T,M) = ZERO
649                    A_M(IJK,B,M) = ZERO
650                    A_M(IJK,0,M) = -ONE
651                    B_M(IJK,M) = ZERO
652                 ELSEIF (FS_WALL_AT(IJK)) THEN
653                    A_M(IJK,E,M) = ONE
654                    A_M(IJK,W,M) = ZERO
655                    A_M(IJK,N,M) = ZERO
656                    A_M(IJK,S,M) = ZERO
657                    A_M(IJK,T,M) = ZERO
658                    A_M(IJK,B,M) = ZERO
659                    A_M(IJK,0,M) = -ONE
660                    B_M(IJK,M) = ZERO
661                 ENDIF
662              ENDDO
663           ENDDO
664     
665           I1 = IMAX2
666           DO K1 = kmin3,kmax3
667              DO J1 = jmin3,jmax3
668                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
669                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
670                 IJK = FUNIJK(I1,J1,K1)
671                 IF (NS_WALL_AT(IJK)) THEN
672                    A_M(IJK,E,M) = ZERO
673                    A_M(IJK,W,M) = -ONE
674                    A_M(IJK,N,M) = ZERO
675                    A_M(IJK,S,M) = ZERO
676                    A_M(IJK,T,M) = ZERO
677                    A_M(IJK,B,M) = ZERO
678                    A_M(IJK,0,M) = -ONE
679                    B_M(IJK,M) = ZERO
680                 ELSEIF (FS_WALL_AT(IJK)) THEN
681                    A_M(IJK,E,M) = ZERO
682                    A_M(IJK,W,M) = ONE
683                    A_M(IJK,N,M) = ZERO
684                    A_M(IJK,S,M) = ZERO
685                    A_M(IJK,T,M) = ZERO
686                    A_M(IJK,B,M) = ZERO
687                    A_M(IJK,0,M) = -ONE
688                    B_M(IJK,M) = ZERO
689                 ENDIF
690              ENDDO
691           ENDDO
692     ! End setting the default boundary conditions
693     
694     
695     ! Setting user specified boundary conditions
696           DO L = 1, DIMENSION_BC
697              IF (BC_DEFINED(L)) THEN
698     
699     ! Setting wall boundary conditions
700                 IF (BC_TYPE(L) == 'NO_SLIP_WALL') THEN
701                    I1 = BC_I_W(L)
702                    I2 = BC_I_E(L)
703                    J1 = BC_J_S(L)
704                    J2 = BC_J_N(L)
705                    K1 = BC_K_B(L)
706                    K2 = BC_K_T(L)
707                    IF (BC_JJ_PS(L) == 0) THEN
708                       DO K = K1, K2
709                          DO J = J1, J2
710                             DO I = I1, I2
711                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
712                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
713                                IJK = FUNIJK(I,J,K)
714                                IF (.NOT.WALL_AT(IJK)) CYCLE  !skip redefined cells
715                                A_M(IJK,E,M) = ZERO
716                                A_M(IJK,W,M) = ZERO
717                                A_M(IJK,N,M) = ZERO
718                                A_M(IJK,S,M) = ZERO
719                                A_M(IJK,T,M) = ZERO
720                                A_M(IJK,B,M) = ZERO
721                                A_M(IJK,0,M) = -ONE
722                                B_M(IJK,M) = ZERO
723                                IF (FLUID_AT(EAST_OF(IJK))) THEN
724                                   A_M(IJK,E,M) = -ONE
725                                ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
726                                   A_M(IJK,W,M) = -ONE
727                                ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
728                                   A_M(IJK,N,M) = -ONE
729                                ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
730                                   A_M(IJK,S,M) = -ONE
731                                ENDIF
732                             ENDDO
733                          ENDDO
734                       ENDDO
735                    ELSE   ! Johnson and Jackson partial slip
736                       CALL JJ_BC_W_S (I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
737                    ENDIF
738     
739                 ELSEIF (BC_TYPE(L) == 'FREE_SLIP_WALL') THEN
740                    I1 = BC_I_W(L)
741                    I2 = BC_I_E(L)
742                    J1 = BC_J_S(L)
743                    J2 = BC_J_N(L)
744                    K1 = BC_K_B(L)
745                    K2 = BC_K_T(L)
746                    IF (BC_JJ_PS(L) == 0) THEN
747                       DO K = K1, K2
748                          DO J = J1, J2
749                             DO I = I1, I2
750                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
751                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
752                                IJK = FUNIJK(I,J,K)
753                                IF (.NOT.WALL_AT(IJK)) CYCLE  !skip redefined cells
754                                A_M(IJK,E,M) = ZERO
755                                A_M(IJK,W,M) = ZERO
756                                A_M(IJK,N,M) = ZERO
757                                A_M(IJK,S,M) = ZERO
758                                A_M(IJK,T,M) = ZERO
759                                A_M(IJK,B,M) = ZERO
760                                A_M(IJK,0,M) = -ONE
761                                B_M(IJK,M) = ZERO
762                                IF (FLUID_AT(EAST_OF(IJK))) THEN
763                                   A_M(IJK,E,M) = ONE
764                                ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
765                                   A_M(IJK,W,M) = ONE
766                                ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
767                                   A_M(IJK,N,M) = ONE
768                                ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
769                                   A_M(IJK,S,M) = ONE
770                                ENDIF
771                             ENDDO
772                          ENDDO
773                       ENDDO
774                    ELSE   ! Johnson and Jackson partial slip
775                       CALL JJ_BC_W_S (I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
776                    ENDIF
777     
778                 ELSEIF (BC_TYPE(L) == 'PAR_SLIP_WALL') THEN
779                    I1 = BC_I_W(L)
780                    I2 = BC_I_E(L)
781                    J1 = BC_J_S(L)
782                    J2 = BC_J_N(L)
783                    K1 = BC_K_B(L)
784                    K2 = BC_K_T(L)
785                    IF (BC_JJ_PS(L) == 0) THEN
786                       DO K = K1, K2
787                          DO J = J1, J2
788                             DO I = I1, I2
789                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
790                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
791                                IJK = FUNIJK(I,J,K)
792                                IF (.NOT.WALL_AT(IJK)) CYCLE  !skip redefined cells
793                                IM = IM1(I)
794                                JM = JM1(J)
795                                A_M(IJK,E,M) = ZERO
796                                A_M(IJK,W,M) = ZERO
797                                A_M(IJK,N,M) = ZERO
798                                A_M(IJK,S,M) = ZERO
799                                A_M(IJK,T,M) = ZERO
800                                A_M(IJK,B,M) = ZERO
801                                A_M(IJK,0,M) = -ONE
802                                B_M(IJK,M) = ZERO
803                                IF (FLUID_AT(EAST_OF(IJK))) THEN
804                                   IF (BC_HW_S(L,M) == UNDEFINED) THEN
805                                      A_M(IJK,E,M) = -HALF
806                                      A_M(IJK,0,M) = -HALF
807                                      B_M(IJK,M) = -BC_WW_S(L,M)
808                                   ELSE
809                                      IF (CYLINDRICAL) THEN
810                                         A_M(IJK,0,M) = -(HALF*(BC_HW_S(L,M)-OX_E(I)&
811                                            )+ODX_E(I))
812                                         A_M(IJK,E,M) = -(HALF*(BC_HW_S(L,M)-OX_E(I)&
813                                            )-ODX_E(I))
814                                      ELSE
815                                         A_M(IJK,0,M)=-(HALF*BC_HW_S(L,M)+ODX_E(I))
816                                         A_M(IJK,E,M)=-(HALF*BC_HW_S(L,M)-ODX_E(I))
817                                      ENDIF
818                                      B_M(IJK,M) = -BC_HW_S(L,M)*BC_WW_S(L,M)
819                                   ENDIF
820                                ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
821                                   IF (BC_HW_S(L,M) == UNDEFINED) THEN
822                                      A_M(IJK,W,M) = -HALF
823                                      A_M(IJK,0,M) = -HALF
824                                      B_M(IJK,M) = -BC_WW_S(L,M)
825                                   ELSE
826                                      IF (CYLINDRICAL) THEN
827                                         A_M(IJK,W,M) = -(HALF*(BC_HW_S(L,M)-OX_E(IM&
828                                            ))-ODX_E(IM))
829                                         A_M(IJK,0,M) = -(HALF*(BC_HW_S(L,M)-OX_E(IM&
830                                            ))+ODX_E(IM))
831                                      ELSE
832                                         A_M(IJK,W,M) = -(HALF*BC_HW_S(L,M)-ODX_E(IM&
833                                            ))
834                                         A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODX_E(IM&
835                                            ))
836                                      ENDIF
837                                      B_M(IJK,M) = -BC_HW_S(L,M)*BC_WW_S(L,M)
838                                   ENDIF
839                                ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
840                                   IF (BC_HW_S(L,M) == UNDEFINED) THEN
841                                      A_M(IJK,N,M) = -HALF
842                                      A_M(IJK,0,M) = -HALF
843                                      B_M(IJK,M) = -BC_WW_S(L,M)
844                                   ELSE
845                                      A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODY_N(J))
846                                      A_M(IJK,N,M) = -(HALF*BC_HW_S(L,M)-ODY_N(J))
847                                      B_M(IJK,M) = -BC_HW_S(L,M)*BC_WW_S(L,M)
848                                   ENDIF
849                                ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
850                                   IF (BC_HW_S(L,M) == UNDEFINED) THEN
851                                      A_M(IJK,S,M) = -HALF
852                                      A_M(IJK,0,M) = -HALF
853                                      B_M(IJK,M) = -BC_WW_S(L,M)
854                                   ELSE
855                                      A_M(IJK,S,M) = -(HALF*BC_HW_S(L,M)-ODY_N(JM))
856                                      A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODY_N(JM))
857                                      B_M(IJK,M) = -BC_HW_S(L,M)*BC_WW_S(L,M)
858                                   ENDIF
859                                ENDIF
860                             ENDDO
861                          ENDDO
862                       ENDDO
863                    ELSE   !Johnson and Jackson partial slip
864                       CALL JJ_BC_W_S (I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
865                    ENDIF
866     
867     ! Setting flow boundary conditions
868                 ELSEIF (BC_TYPE(L)=='P_INFLOW' .OR. BC_TYPE(L)=='P_OUTFLOW') THEN
869                    IF (BC_PLANE(L) == 'B') THEN
870                       I1 = BC_I_W(L)
871                       I2 = BC_I_E(L)
872                       J1 = BC_J_S(L)
873                       J2 = BC_J_N(L)
874                       K1 = BC_K_B(L)
875                       K2 = BC_K_T(L)
876                       DO K = K1, K2
877                          DO J = J1, J2
878                             DO I = I1, I2
879                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
880                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
881                                IJK = FUNIJK(I,J,K)
882                                A_M(IJK,E,M) = ZERO
883                                A_M(IJK,W,M) = ZERO
884                                A_M(IJK,N,M) = ZERO
885                                A_M(IJK,S,M) = ZERO
886                                A_M(IJK,T,M) = ZERO
887                                A_M(IJK,B,M) = ONE
888                                A_M(IJK,0,M) = -ONE
889                                B_M(IJK,M) = ZERO
890                             ENDDO
891                          ENDDO
892                       ENDDO
893                    ENDIF
894     
895                 ELSEIF (BC_TYPE(L) == 'OUTFLOW') THEN
896                    IF (BC_PLANE(L) == 'B') THEN
897                       I1 = BC_I_W(L)
898                       I2 = BC_I_E(L)
899                       J1 = BC_J_S(L)
900                       J2 = BC_J_N(L)
901                       K1 = BC_K_B(L)
902                       K2 = BC_K_T(L)
903                       DO K = K1, K2
904                          DO J = J1, J2
905                             DO I = I1, I2
906                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
907                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
908                                IJK = FUNIJK(I,J,K)
909                                A_M(IJK,E,M) = ZERO
910                                A_M(IJK,W,M) = ZERO
911                                A_M(IJK,N,M) = ZERO
912                                A_M(IJK,S,M) = ZERO
913                                A_M(IJK,T,M) = ZERO
914                                A_M(IJK,B,M) = ONE
915                                A_M(IJK,0,M) = -ONE
916                                B_M(IJK,M) = ZERO
917                                IJKM = KM_OF(IJK)
918                                A_M(IJKM,E,M) = ZERO
919                                A_M(IJKM,W,M) = ZERO
920                                A_M(IJKM,N,M) = ZERO
921                                A_M(IJKM,S,M) = ZERO
922                                A_M(IJKM,T,M) = ZERO
923                                A_M(IJKM,B,M) = ONE
924                                A_M(IJKM,0,M) = -ONE
925                                B_M(IJKM,M) = ZERO
926                             ENDDO
927                          ENDDO
928                       ENDDO
929                    ELSEIF (BC_PLANE(L) == 'T') THEN
930                       I1 = BC_I_W(L)
931                       I2 = BC_I_E(L)
932                       J1 = BC_J_S(L)
933                       J2 = BC_J_N(L)
934                       K1 = BC_K_B(L)
935                       K2 = BC_K_T(L)
936                       DO K = K1, K2
937                          DO J = J1, J2
938                             DO I = I1, I2
939                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
940                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
941                                IJK = FUNIJK(I,J,K)
942                                IJKP = KP_OF(IJK)
943                                A_M(IJKP,E,M) = ZERO
944                                A_M(IJKP,W,M) = ZERO
945                                A_M(IJKP,N,M) = ZERO
946                                A_M(IJKP,S,M) = ZERO
947                                A_M(IJKP,T,M) = ONE
948                                A_M(IJKP,B,M) = ZERO
949                                A_M(IJKP,0,M) = -ONE
950                                B_M(IJKP,M) = ZERO
951                             ENDDO
952                          ENDDO
953                       ENDDO
954                    ENDIF
955     
956     ! Setting bc that are not ns, fs, psw, p_inflow, p_outflow, or outflow
957                 ELSE
958                    I1 = BC_I_W(L)
959                    I2 = BC_I_E(L)
960                    J1 = BC_J_S(L)
961                    J2 = BC_J_N(L)
962                    K1 = BC_K_B(L)
963                    K2 = BC_K_T(L)
964                    DO K = K1, K2
965                       DO J = J1, J2
966                          DO I = I1, I2
967                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
968                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
969                             IJK = FUNIJK(I,J,K)
970                             A_M(IJK,E,M) = ZERO
971                             A_M(IJK,W,M) = ZERO
972                             A_M(IJK,N,M) = ZERO
973                             A_M(IJK,S,M) = ZERO
974                             A_M(IJK,T,M) = ZERO
975                             A_M(IJK,B,M) = ZERO
976                             A_M(IJK,0,M) = -ONE
977                             B_M(IJK,M) = -W_S(IJK,M)
978                             IF (BC_PLANE(L) == 'B') THEN
979                                IJKB = BOTTOM_OF(IJK)
980                                A_M(IJKB,E,M) = ZERO
981                                A_M(IJKB,W,M) = ZERO
982                                A_M(IJKB,N,M) = ZERO
983                                A_M(IJKB,S,M) = ZERO
984                                A_M(IJKB,T,M) = ZERO
985                                A_M(IJKB,B,M) = ZERO
986                                A_M(IJKB,0,M) = -ONE
987                                B_M(IJKB,M) = -W_S(IJKB,M)
988                             ENDIF
989                          ENDDO
990                       ENDDO
991                    ENDDO
992     
993                 ENDIF   ! end if (bc_type)
994              ENDIF   ! end if (bc_defined)
995           ENDDO   ! end L do loop over dimension_bc
996     
997           RETURN
998           END SUBROUTINE SOURCE_W_S_BC
999     
1000     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1001     !                                                                      C
1002     !  Subroutine: JJ_BC_W_s                                               C
1003     !  Purpose: Implement Johnson and Jackson boundary condition           C
1004     !                                                                      C
1005     !  Author: K. Agrawal, A. Srivastava,                 Date: 14-APR-98  C
1006     !          Princeton University                                        C
1007     !  Reviewer:                                          Date:            C
1008     !                                                                      C
1009     !                                                                      C
1010     !                                                                      C
1011     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1012     
1013           SUBROUTINE JJ_BC_W_S(I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
1014     
1015     !-----------------------------------------------
1016     ! Modules
1017     !-----------------------------------------------
1018           USE param
1019           USE param1
1020           USE parallel
1021           USE matrix
1022           USE scales
1023           USE constant
1024           USE physprop
1025           USE fldvar
1026           USE visc_s
1027           USE rxns
1028           USE run
1029           USE toleranc
1030           USE geometry
1031           USE indices
1032           USE is
1033           USE tau_s
1034           USE bc
1035           USE output
1036           USE compar
1037           USE functions
1038     
1039           IMPLICIT NONE
1040     !-----------------------------------------------
1041     ! Dummy Arguments
1042     !-----------------------------------------------
1043     ! Boundary condition
1044           INTEGER, INTENT(IN) :: L
1045     ! Indices
1046           INTEGER, INTENT(IN) :: I1, I2, J1, J2, K1, K2
1047     ! Solids phase index
1048           INTEGER, INTENT(IN) ::  M
1049     ! Septadiagonal matrix A_m
1050           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1051     ! Vector b_m
1052           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1053     !-----------------------------------------------
1054     ! Local Variables
1055     !-----------------------------------------------
1056     ! Indices
1057           INTEGER :: I, J, K, IJK, JM, IM, IJKP
1058     ! coefficients for granular bc
1059           DOUBLE PRECISION :: hw, gw, cw
1060     !-----------------------------------------------
1061     
1062           DO K = K1, K2
1063              DO J = J1, J2
1064                 DO I = I1, I2
1065                    IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
1066                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1067                    IJK = FUNIJK(I,J,K)
1068                    IF (.NOT.WALL_AT(IJK)) CYCLE  ! skip redefined cells
1069                    IM = IM1(I)
1070                    JM = JM1(J)
1071                    A_M(IJK,E,M) = ZERO
1072                    A_M(IJK,W,M) = ZERO
1073                    A_M(IJK,N,M) = ZERO
1074                    A_M(IJK,S,M) = ZERO
1075                    A_M(IJK,T,M) = ZERO
1076                    A_M(IJK,B,M) = ZERO
1077                    A_M(IJK,0,M) = -ONE
1078                    B_M(IJK,M) = ZERO
1079     
1080                    IF (FLUID_AT(EAST_OF(IJK))) THEN
1081                       IJKP = KP_OF(EAST_OF(IJK))
1082                       IF (WALL_AT(IJKP)) CYCLE
1083                       IF (EP_S(EAST_OF(IJK),M) <= DIL_EP_S) THEN
1084                          A_M(IJK,E,M) = ONE
1085                       ELSE
1086                          IF (BC_JJ_PS(L) == 1) THEN
1087                             CALL CALC_GRBDRY (IJK, EAST_OF(IJK), 'E', 'W',&
1088                                M, L, GW, HW, CW)
1089                          ELSEIF (BC_JJ_PS(L) == 2) THEN
1090                             GW = 0D0
1091                             HW = 1D0
1092                             CW = BC_WW_S(L,M)
1093                          ELSE
1094                             GW = 1D0
1095                             CW = 0D0
1096                             HW = 0D0
1097                          ENDIF
1098                          IF (CYLINDRICAL) THEN
1099                             A_M(IJK,E,M) = -(HALF*(HW - OX_E(I)*GW)-ODX_E(I)*GW)
1100                             A_M(IJK,0,M) = -(HALF*(HW - OX_E(I)*GW)+ODX_E(I)*GW)
1101                          ELSE
1102                             A_M(IJK,E,M) = -(HALF*HW - ODX_E(I)*GW)
1103                             A_M(IJK,0,M) = -(HALF*HW + ODX_E(I)*GW)
1104                          ENDIF
1105                          B_M(IJK,M) = -CW
1106                       ENDIF
1107     
1108                    ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
1109                       IJKP = KP_OF(WEST_OF(IJK))
1110                       IF (WALL_AT(IJKP)) CYCLE
1111                       IF (EP_S(WEST_OF(IJK),M) <= DIL_EP_S) THEN
1112                          A_M(IJK,W,M) = ONE
1113                       ELSE
1114                          IF (BC_JJ_PS(L) == 1) THEN
1115                             CALL CALC_GRBDRY (IJK, WEST_OF(IJK), 'W', 'W',&
1116                                M, L, GW, HW, CW)
1117                          ELSEIF (BC_JJ_PS(L) == 2) THEN
1118                             GW = 0D0
1119                             HW = 1D0
1120                             CW = BC_WW_S(L,M)
1121                          ELSE
1122                             GW = 1D0
1123                             CW = 0D0
1124                             HW = 0D0
1125                          ENDIF
1126                          IF (CYLINDRICAL) THEN
1127                             A_M(IJK,W,M) = -(HALF*(HW - OX_E(IM)*GW)-ODX_E(IM)*GW)
1128                             A_M(IJK,0,M) = -(HALF*(HW - OX_E(IM)*GW)+ODX_E(IM)*GW)
1129                          ELSE
1130                             A_M(IJK,W,M) = -(HALF*HW - ODX_E(IM)*GW)
1131                             A_M(IJK,0,M) = -(HALF*HW + ODX_E(IM)*GW)
1132                          ENDIF
1133                          B_M(IJK,M) = -CW
1134                       ENDIF
1135     
1136                    ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
1137                       IJKP = KP_OF(NORTH_OF(IJK))
1138                       IF (WALL_AT(IJKP)) CYCLE
1139                       IF (EP_S(NORTH_OF(IJK),M) <= DIL_EP_S) THEN
1140                          A_M(IJK,N,M) = ONE
1141                       ELSE
1142                          IF (BC_JJ_PS(L) == 1) THEN
1143                             CALL CALC_GRBDRY (IJK, NORTH_OF(IJK), 'N', 'W',&
1144                                M, L, GW, HW, CW)
1145                          ELSEIF (BC_JJ_PS(L) == 2) THEN
1146                             GW = 0D0
1147                             HW = 1D0
1148                             CW = BC_WW_S(L,M)
1149                          ELSE
1150                             GW = 1D0
1151                             CW = 0D0
1152                             HW = 0D0
1153                          ENDIF
1154                          A_M(IJK,N,M) = -(HALF*HW - ODY_N(J)*GW)
1155                          A_M(IJK,0,M) = -(HALF*HW + ODY_N(J)*GW)
1156                          B_M(IJK,M) = -CW
1157                       ENDIF
1158     
1159                    ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
1160                       IJKP = KP_OF(SOUTH_OF(IJK))
1161                       IF (WALL_AT(IJKP)) CYCLE
1162                       IF (EP_S(SOUTH_OF(IJK),M) <= DIL_EP_S) THEN
1163                          A_M(IJK,S,M) = ONE
1164                       ELSE
1165                          IF (BC_JJ_PS(L) == 1) THEN
1166                             CALL CALC_GRBDRY (IJK, SOUTH_OF(IJK), 'S', 'W',&
1167                                M, L, GW, HW, CW)
1168                          ELSEIF (BC_JJ_PS(L) == 2) THEN
1169                             GW = 0D0
1170                             HW = 1D0
1171                             CW = BC_WW_S(L,M)
1172                          ELSE
1173                             GW = 1D0
1174                             CW = 0D0
1175                             HW = 0D0
1176                          ENDIF
1177                          A_M(IJK,S,M) = -(HALF*HW - ODY_N(JM)*GW)
1178                          A_M(IJK,0,M) = -(HALF*HW + ODY_N(JM)*GW)
1179                          B_M(IJK,M) = -CW
1180                       ENDIF
1181                    ENDIF
1182     
1183                 ENDDO
1184              ENDDO
1185           ENDDO
1186     
1187           RETURN
1188           END SUBROUTINE JJ_BC_W_S
1189     
1190     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1191     !                                                                      C
1192     !  Subroutine: POINT_SOURCE_W_S                                        C
1193     !  Purpose: Adds point sources to the solids W-Momentum equation.      C
1194     !                                                                      C
1195     !  Author: J. Musser                                  Date: 10-JUN-13  C
1196     !  Reviewer:                                          Date:            C
1197     !                                                                      C
1198     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1199           SUBROUTINE POINT_SOURCE_W_S(A_M, B_M)
1200     
1201     !-----------------------------------------------
1202     ! Modules
1203     !-----------------------------------------------
1204           use compar
1205           use constant
1206           use fldvar
1207           use geometry
1208           use indices
1209           use physprop
1210           use ps
1211           use run
1212           use functions
1213     
1214           IMPLICIT NONE
1215     !-----------------------------------------------
1216     ! Dummy arguments
1217     !-----------------------------------------------
1218     ! Septadiagonal matrix A_m
1219           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1220     ! Vector b_m
1221           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1222     !-----------------------------------------------
1223     ! Local variables
1224     !-----------------------------------------------
1225     ! Indices
1226           INTEGER :: IJK, I, J, K
1227           INTEGER :: PSV, M
1228           INTEGER :: lKT, lKB
1229     ! terms of bm expression
1230           DOUBLE PRECISION :: pSource
1231     !-----------------------------------------------
1232           do M=1, MMAX
1233     
1234           PS_LP: do PSV = 1, DIMENSION_PS
1235              if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
1236              if(abs(PS_W_s(PSV,M)) < small_number) cycle PS_LP
1237     
1238              if(PS_W_s(PSV,M) < ZERO) then
1239                 lKB = PS_K_B(PSV)-1
1240                 lKT = PS_K_T(PSV)-1
1241              else
1242                 lKB = PS_K_B(PSV)
1243                 lKT = PS_K_T(PSV)
1244              endif
1245     
1246              do k = lKB, lKT
1247              do j = PS_J_S(PSV), PS_J_N(PSV)
1248              do i = PS_I_W(PSV), PS_I_E(PSV)
1249     
1250                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
1251                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1252     
1253                 ijk = funijk(i,j,k)
1254                 if(.NOT.fluid_at(ijk)) cycle
1255     
1256     
1257                 if(A_M(IJK,0,M) == -ONE .AND.                           &
1258                    B_M(IJK,M) == -W_s(IJK,M)) then
1259                    B_M(IJK,M) = -PS_W_s(PSV,M) * PS_VEL_MAG_S(PSV,M)
1260                 else
1261                    pSource = PS_MASSFLOW_S(PSV,M) *                     &
1262                       (VOL(IJK)/PS_VOLUME(PSV))
1263     
1264                    B_M(IJK,M) = B_M(IJK,M) - pSource *                  &
1265                       PS_W_s(PSV,M) * PS_VEL_MAG_S(PSV,M)
1266                 endif
1267     
1268              enddo
1269              enddo
1270              enddo
1271     
1272              enddo PS_LP
1273     
1274           enddo ! do M=1, MMAX
1275     
1276           RETURN
1277           END SUBROUTINE POINT_SOURCE_W_S
1278