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