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