source_w

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