File: N:\mfix\model\source_v_s.f

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