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