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