source_w

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