File: RELATIVE:/../../../mfix.git/model/source_u_s.f

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