source_w

File: /nfs/home/0/users/jenkins/mfix.git/model/source_w_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_W_g                                              C
4     !  Purpose: Determine source terms for W_g 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: 17-JUN-96  C
12     !  Reviewer:                                          Date:            C
13     !                                                                      C
14     !  Revision Number: 1                                                  C
15     !  Purpose: To incorporate Cartesian grid modifications                C
16     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
17     !                                                                      C
18     !                                                                      C
19     !  Literature/Document References:                                     C
20     !                                                                      C
21     !                                                                      C
22     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
23     
24           SUBROUTINE SOURCE_W_G(A_M, B_M, IER)
25     
26     !-----------------------------------------------
27     ! Modules
28     !-----------------------------------------------
29           USE bc
30           USE bodyforce
31           USE compar
32           USE constant
33           USE cutcell
34           USE drag
35           USE fldvar
36           USE fun_avg
37           USE functions
38           USE geometry
39           USE ghdtheory
40           USE indices
41           USE is
42           USE matrix
43           USE mms
44           USE parallel
45           USE param
46           USE param1
47           USE physprop
48           USE quadric
49           USE run
50           USE rxns
51           USE scales
52           USE sendrecv
53           USE tau_g
54           USE toleranc
55           USE visc_g
56     
57           IMPLICIT NONE
58     !-----------------------------------------------
59     ! Dummy arguments
60     !-----------------------------------------------
61     ! Septadiagonal matrix A_m
62           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
63     ! Vector b_m
64           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
65     ! Error index
66           INTEGER, INTENT(INOUT) :: IER
67     !-----------------------------------------------
68     ! Local variables
69     !-----------------------------------------------
70     ! Indices
71           INTEGER :: I, J, K, IJK, IJKT, IMJK, IJKP, IMJKP,&
72                      IJKE, IJKW, IJKTE, IJKTW, IM, IPJK,   &
73                      IJKM, IJMK, IJMKP, IJPK
74     ! Phase index
75           INTEGER :: M, L, MM
76     ! Internal surface
77           INTEGER :: ISV
78     ! Pressure at top cell
79           DOUBLE PRECISION :: PgT
80     ! Average volume fraction
81           DOUBLE PRECISION :: EPGA
82     ! Average density
83           DOUBLE PRECISION :: ROPGA, ROGA
84     ! Average viscosity
85           DOUBLE PRECISION :: MUGA
86     ! Average coefficients
87           DOUBLE PRECISION :: Cte, Ctw, EPMUoX
88     ! Average U_g
89           DOUBLE PRECISION Ugt
90     ! Source terms (Surface)
91           DOUBLE PRECISION Sdp, Sxzb
92     ! Source terms (Volumetric)
93           DOUBLE PRECISION V0, Vpm, Vmt, Vbf, Vcoa, Vcob, Vxza, Vxzb
94     ! Source terms (Volumetric) for GHD theory
95           DOUBLE PRECISION Ghd_drag, avgRop
96     ! Source terms for HYS drag relation
97           DOUBLE PRECISION HYS_drag, avgDrag
98     ! virtual (added) mass
99           DOUBLE PRECISION :: ROP_MA, U_se, Usw, Ust, Vsb, Vst, &
100                               Wse, Wsw, Wsn, Wss, Wst, Wsb
101           DOUBLE PRECISION F_vir
102     !-----------------------------------------------
103     
104     ! Set reference phase to gas
105           M = 0
106     
107           IF (.NOT.MOMENTUM_Z_EQ(0)) RETURN
108     
109     !
110     ! trailing commas added to the first two lines in the
111     ! private list
112     ! Charles Crosby
113     ! CHPC, 26 September 2013
114     !
115     
116     !$omp  parallel do default(shared)                                   &
117     !$omp  private(I, J, K, IJK, IJKT, IJKM, IJKP, IMJK, IPJK, IJMK,     &
118     !$omp          IMJKP, IJPK, IJMKP, IJKTE, IJKTW, IM, IJKW, IJKE,     &
119     !$omp          EPGA, PGT, SDP, ROPGA, ROGA, V0, ISV, MUGA, Vpm,      &
120     !$omp          Vmt, Vbf, F_vir, Ghd_drag, avgRop, HYS_drag, avgDrag, &
121     !$omp          MM, L, VXZA, VXZB, VCOA, VCOB, CTE, CTW, UGT,         &
122     !$omp          SXZB, EPMUOX)
123           DO IJK = ijkstart3, ijkend3
124              I = I_OF(IJK)
125              J = J_OF(IJK)
126              K = K_OF(IJK)
127              IJKT = TOP_OF(IJK)
128              IJKM = KM_OF(IJK)
129              IJKP = KP_OF(IJK)
130              IMJK = IM_OF(IJK)
131              IPJK = IP_OF(IJK)
132              IMJKP = KP_OF(IMJK)
133              IJMK = JM_OF(IJK)
134              IJPK = JP_OF(IJK)
135              IJMKP = KP_OF(IJMK)
136     
137              EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
138     
139     ! Impermeable internal surface
140              IF (IP_AT_T(IJK)) THEN
141                 A_M(IJK,E,M) = ZERO
142                 A_M(IJK,W,M) = ZERO
143                 A_M(IJK,N,M) = ZERO
144                 A_M(IJK,S,M) = ZERO
145                 A_M(IJK,T,M) = ZERO
146                 A_M(IJK,B,M) = ZERO
147                 A_M(IJK,0,M) = -ONE
148                 B_M(IJK,M) = ZERO
149     
150     ! Dilute flow
151              ELSEIF (EPGA <= DIL_EP_S) THEN
152                 A_M(IJK,E,M) = ZERO
153                 A_M(IJK,W,M) = ZERO
154                 A_M(IJK,N,M) = ZERO
155                 A_M(IJK,S,M) = ZERO
156                 A_M(IJK,T,M) = ZERO
157                 A_M(IJK,B,M) = ZERO
158                 A_M(IJK,0,M) = -ONE
159                 B_M(IJK,M) = ZERO
160     ! set velocity equal to that of bottom or top cell if solids are
161     ! present in those cells else set velocity equal to known value
162                 IF (EP_G(BOTTOM_OF(IJK)) > DIL_EP_S) THEN
163                    A_M(IJK,B,M) = ONE
164                 ELSE IF (EP_G(TOP_OF(IJK)) > DIL_EP_S) THEN
165                    A_M(IJK,T,M) = ONE
166                 ELSE
167                    B_M(IJK,M) = -W_G(IJK)
168                 ENDIF
169     
170     ! Cartesian grid implementation
171              ELSEIF (BLOCKED_W_CELL_AT(IJK)) THEN
172                 A_M(IJK,E,M) = ZERO
173                 A_M(IJK,W,M) = ZERO
174                 A_M(IJK,N,M) = ZERO
175                 A_M(IJK,S,M) = ZERO
176                 A_M(IJK,T,M) = ZERO
177                 A_M(IJK,B,M) = ZERO
178                 A_M(IJK,0,M) = -ONE
179                 B_M(IJK,M) = ZERO
180     
181     ! Normal case
182              ELSE
183     
184     ! Surface forces
185     
186     ! Pressure term
187                 PGT = P_G(IJKT)
188                 IF (CYCLIC_Z_PD) THEN
189                    IF (KMAP(K_OF(IJK)).EQ.KMAX1) PGT = P_G(IJKT) - DELP_Z
190                 ENDIF
191                 IF (MODEL_B) THEN
192                    IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
193                        SDP = -P_SCALE*(PGT - P_G(IJK))*AXY(IJK)
194                    ELSE
195                        SDP = -P_SCALE*(PGT * A_WPG_T(IJK) - P_G(IJK) * A_WPG_B(IJK) )
196                    ENDIF
197                 ELSE
198                    IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
199                        SDP = -P_SCALE*EPGA*(PGT - P_G(IJK))*AXY(IJK)
200                    ELSE
201                        SDP = -P_SCALE*EPGA*(PGT * A_WPG_T(IJK) - P_G(IJK) * A_WPG_B(IJK) )
202                    ENDIF
203                 ENDIF
204     
205                 IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
206     ! Volumetric forces
207                    ROPGA = AVG_Z(ROP_G(IJK),ROP_G(IJKT),K)
208                    ROGA = AVG_Z(RO_G(IJK),RO_G(IJKT),K)
209     
210     ! Previous time step
211                    V0 = AVG_Z(ROP_GO(IJK),ROP_GO(IJKT),K)*ODT
212     
213     ! Added mass implicit transient term {Cv eps rop_g dW/dt}
214                    IF(Added_Mass) THEN
215                       ROP_MA = AVG_Z(ROP_g(IJK)*EP_s(IJK,M_AM),&
216                          ROP_g(IJKT)*EP_s(IJKT,M_AM),K)
217                       V0 = V0 + Cv * ROP_MA * ODT
218                    ENDIF
219                 ELSE
220     ! Volumetric forces
221                    ROPGA = (VOL(IJK)*ROP_G(IJK) + VOL(IJKT)*ROP_G(IJKT))/&
222                       (VOL(IJK) + VOL(IJKT))
223                    ROGA  = (VOL(IJK)*RO_G(IJK) + VOL(IJKT)*RO_G(IJKT))/&
224                       (VOL(IJK) + VOL(IJKT))
225     ! Previous time step
226                    V0 = (VOL(IJK)*ROP_GO(IJK) + VOL(IJKT)*ROP_GO(IJKT))*&
227                       ODT/(VOL(IJK) + VOL(IJKT))
228     ! Added mass implicit transient term {Cv eps rop_g dW/dt}
229                    IF(Added_Mass) THEN
230                       ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + &
231                          VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M_AM))/&
232                          (VOL(IJK) + VOL(IJKT))
233                       V0 = V0 + Cv * ROP_MA * ODT
234                    ENDIF
235                 ENDIF
236     
237     ! VIRTUAL MASS SECTION (explicit terms)
238     ! adding transient term dWs/dt to virtual mass term
239                 F_vir = ZERO
240                 IF(Added_Mass.AND.(.NOT.CUT_W_TREATMENT_AT(IJK))) THEN
241                    F_vir = ( (W_s(IJK,M_AM) - W_sO(IJK,M_AM)) )*ODT*VOL_W(IJK)
242     
243     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
244                    Wsb = AVG_Z_T(W_S(IJKM,M_AM),W_s(IJK,M_AM))
245                    Wst = AVG_Z_T(W_s(IJK,M_AM),W_s(IJKP,M_AM))
246                    U_se = AVG_Z(U_s(IJK,M_AM),U_s(IJKP,M_AM),K)
247                    Usw = AVG_Z(U_s(IMJK,M_AM),U_s(IMJKP,M_AM),K)
248                    Ust = AVG_X_E(Usw,U_se,IP1(I))
249                    Wse = AVG_X(W_s(IJK,M_AM),W_s(IPJK,M_AM),IP1(I))
250                    Wsw = AVG_X(W_s(IMJK,M_AM),W_s(IJK,M_AM),I)
251                    Vsb = AVG_Y_N(V_s(IJMK,M_AM),V_s(IJK,M_AM))
252                    Vst = AVG_Y_N(V_s(IJMKP,M_AM),V_s(IJKP,M_AM))
253                    Wss = AVG_Y(W_s(IJMK,M_AM),W_s(IJK,M_AM),JM1(J))
254                    Wsn = AVG_Y(W_s(IJK,M_AM),W_s(IJPK,M_AM),J)
255     
256     ! adding convective terms (U dW/dx + V dW/dy + W dW/dz) to virtual mass.
257                    F_vir = F_vir + W_s(IJK,M_AM)*OX(I)* &
258                       (Wst - Wsb)*AXY(IJK) + Ust*(Wse - Wsw)*AYZ(IJK) + &
259                       AVG_Z(Vsb,Vst,K)*(Wsn - Wss)*AXZ(IJK)
260     
261     ! Coriolis force
262                    IF (CYLINDRICAL) F_vir = F_vir + &
263                       Ust*W_s(IJK,M_AM)*OX(I)
264                    F_vir = F_vir * Cv * ROP_MA
265                 ENDIF
266     
267     ! pressure drop through porous media
268                 IF (SIP_AT_T(IJK)) THEN
269                    ISV = IS_ID_AT_T(IJK)
270                    MUGA = AVG_Z(MU_G(IJK),MU_G(IJKT),K)
271                    VPM = MUGA/IS_PC(ISV,1)
272                    IF (IS_PC(ISV,2) /= ZERO) VPM = VPM + &
273                       HALF*IS_PC(ISV,2)*ROPGA*ABS(W_G(IJK))
274                 ELSE
275                    VPM = ZERO
276                 ENDIF
277     
278     ! Interphase mass transfer
279                 IF(.NOT.CUT_W_TREATMENT_AT(IJK)) THEN
280                    VMT = AVG_Z(SUM_R_G(IJK),SUM_R_G(IJKT),K)
281                 ELSE
282                    VMT = (VOL(IJK)*SUM_R_G(IJK) + VOL(IJKT)*SUM_R_G(IJKT))/&
283                       (VOL(IJK) + VOL(IJKT))
284                 ENDIF
285     
286     ! Body force
287                 IF (MODEL_B) THEN
288                    VBF = ROGA*BFZ_G(IJK)
289                 ELSE
290                    VBF = ROPGA*BFZ_G(IJK)
291                 ENDIF
292     
293     ! Additional force for GHD from darg force sum(beta_ig * Joi/rhop_i)
294                 Ghd_drag = ZERO
295                 IF (KT_TYPE_ENUM .EQ. GHD_2007) THEN
296                    DO L = 1,SMAX
297                       avgRop = AVG_Z(ROP_S(IJK,L),ROP_S(IJKT,L),K)
298                       if(avgRop > ZERO) Ghd_drag = Ghd_drag +&
299                          AVG_Z(F_GS(IJK,L),F_GS(IJKT,L),K) * JoiZ(IJK,L) / avgRop
300                    ENDDO
301                 ENDIF
302     
303     ! Additional force for HYS drag force, do not use with mixture GHD theory
304                 avgDrag = ZERO
305                 HYS_drag = ZERO
306                 IF (DRAG_TYPE_ENUM .EQ. HYS .AND. KT_TYPE_ENUM .NE. GHD_2007) THEN
307                    DO MM=1,MMAX
308                       DO L = 1,MMAX
309                          IF (L /= MM) THEN
310                             avgDrag = AVG_Z(beta_ij(IJK,MM,L),beta_ij(IJKT,MM,L),K)
311                             HYS_drag = HYS_drag + avgDrag * (W_g(IJK) - W_s(IJK,L))
312                          ENDIF
313                       ENDDO
314                    ENDDO
315                 ENDIF
316     
317     ! Special terms for cylindrical coordinates
318                 VCOA = ZERO
319                 VCOB = ZERO
320                 SXZB = ZERO
321                 VXZA = ZERO
322                 VXZB = ZERO
323                 CTE  = ZERO
324                 CTW  = ZERO
325                 IF (CYLINDRICAL) THEN
326     ! Coriolis force
327                    IMJK = IM_OF(IJK)
328                    IJKP = KP_OF(IJK)
329                    IMJKP = KP_OF(IMJK)
330                    UGT = AVG_Z(HALF*(U_G(IJK)+U_G(IMJK)),&
331                       HALF*(U_G(IJKP)+U_G(IMJKP)),K)
332                    IF (UGT > ZERO) THEN
333                       VCOA = ROPGA*UGT*OX(I)
334                       VCOB = ZERO
335     ! virtual mass contribution
336                       IF(Added_Mass) VCOA = VCOA + Cv*ROP_MA*UGT*OX(I)
337                    ELSE
338                       VCOA = ZERO
339                       VCOB = -ROPGA*UGT*W_G(IJK)*OX(I)
340     ! virtual mass contribution
341                       IF(Added_Mass) VCOB = VCOB - Cv*ROP_MA*UGT*W_G(IJK)*OX(I)
342                    ENDIF
343     
344     ! Term from tau_xz: integral of (1/x)*(d/dx)(x*mu*(-w/x))
345                    IJKE = EAST_OF(IJK)
346                    IJKW = WEST_OF(IJK)
347                    IJKTE = TOP_OF(IJKE)
348                    IJKTW = TOP_OF(IJKW)
349                    IM = IM1(I)
350                    IPJK = IP_OF(IJK)
351                    CTE = HALF*AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),&
352                                       AVG_X_H(MU_GT(IJKT),MU_GT(IJKTE),I),K)*&
353                                       OX_E(I)*AYZ_W(IJK)
354                    CTW = HALF*AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),&
355                                       AVG_X_H(MU_GT(IJKTW),MU_GT(IJKT),IM),K)*&
356                                       DY(J)*(HALF*(DZ(K)+DZ(KP1(K))))
357                         ! same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
358     
359     ! (mu/x)*(dw/dx) part of tau_xz/x
360                    EPMUOX = AVG_Z(MU_GT(IJK),MU_GT(IJKT),K)*OX(I)
361                    VXZB = ZERO
362                    A_M(IJK,E,M) = A_M(IJK,E,M) + HALF*EPMUOX*ODX_E(I)*VOL_W(IJK)
363                    A_M(IJK,W,M) = A_M(IJK,W,M) - HALF*EPMUOX*ODX_E(IM)*VOL_W(IJK)
364     
365     ! -(mu/x)*(w/x) part of tau_xz/x
366                    VXZA = EPMUOX*OX(I)
367                 ELSE
368                    VCOA = ZERO
369                    VCOB = ZERO
370                    SXZB = ZERO
371                    VXZA = ZERO
372                    VXZB = ZERO
373                 ENDIF
374     
375     ! Collect the terms
376                 A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+&
377                    A_M(IJK,N,M)+A_M(IJK,S,M)+A_M(IJK,T,M)+A_M(IJK,B,M)+&
378                    (V0+VPM+ZMAX(VMT)+VCOA + VXZA)*VOL_W(IJK) + CTE - CTW)
379                 A_M(IJK,E,M) = A_M(IJK,E,M) - CTE
380                 A_M(IJK,W,M) = A_M(IJK,W,M) + CTW
381                 B_M(IJK,M) = B_M(IJK,M) - ( SDP + TAU_W_G(IJK) + SXZB + &
382                    ( (V0+ZMAX((-VMT)))*W_GO(IJK) + VBF + VCOB + VXZB + &
383                    Ghd_drag+HYS_drag)*VOL_W(IJK) )
384     ! adding explicit term of virtual mass force
385                 B_M(IJK,M) = B_M(IJK,M) - F_vir
386     ! MMS Source term.
387                 IF(USE_MMS) B_M(IJK,M) = &
388                    B_M(IJK,M) - MMS_W_G_SRC(IJK)*VOL_W(IJK)
389              ENDIF   ! end branching on cell type (ip/dilute/block/else branches)
390           ENDDO   ! end do loop over ijk
391     !$omp end parallel do
392     
393     ! modifications for cartesian grid implementation
394           IF(CARTESIAN_GRID) CALL CG_SOURCE_W_G(A_M, B_M, IER)
395     ! modifications for bc
396           CALL SOURCE_W_G_BC (A_M, B_M)
397     ! modifications for cartesian grid implementation
398           IF(CARTESIAN_GRID) CALL CG_SOURCE_W_G_BC(A_M, B_M, IER)
399     
400           RETURN
401           END SUBROUTINE SOURCE_W_G
402     
403     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
404     !                                                                      C
405     !  Subroutine: SOURCE_W_g_BC                                           C
406     !  Purpose: Determine source terms for W_g momentum eq. The terms      C
407     !     appear in the center coefficient and RHS vector. The center      C
408     !     coefficient and source vector are negative.  The off-diagonal    C
409     !     coefficients are positive.                                       C
410     !     The drag terms are excluded from the source at this stage        C
411     !                                                                      C
412     !  Author: M. Syamlal                                 Date: 17-JUN-96  C
413     !  Reviewer:                                          Date:            C
414     !                                                                      C
415     !                                                                      C
416     !                                                                      C
417     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
418     
419           SUBROUTINE SOURCE_W_G_BC(A_M, B_M)
420     
421     !-----------------------------------------------
422     ! Modules
423     !-----------------------------------------------
424           USE param
425           USE param1
426           USE parallel
427           USE matrix
428           USE scales
429           USE constant
430           USE physprop
431           USE fldvar
432           USE visc_g
433           USE rxns
434           USE run
435           USE toleranc
436           USE geometry
437           USE indices
438           USE is
439           USE tau_g
440           USE bc
441           USE output
442           USE compar
443           USE fun_avg
444           USE functions
445     
446           IMPLICIT NONE
447     !-----------------------------------------------
448     ! Dummy arguments
449     !-----------------------------------------------
450     ! Septadiagonal matrix A_m
451           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
452     ! Vector b_m
453           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
454     !-----------------------------------------------
455     ! Local parameters
456     !-----------------------------------------------
457     ! C_mu is constant in turbulent viscosity
458           DOUBLE PRECISION, PARAMETER :: C_mu = 0.09D0
459     ! Kappa is Von Karmen constant
460           DOUBLE PRECISION, PARAMETER :: Kappa = 0.42D0
461     !-----------------------------------------------
462     ! Local variables
463     !-----------------------------------------------
464     ! Boundary condition
465           INTEGER :: L
466     ! Indices
467           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
468                      IM, JM, IJKB, IJKM, IJKP
469     ! Phase index
470           INTEGER :: M
471     ! Turbulent shear at walls
472           DOUBLE PRECISION W_F_Slip
473     
474     !-----------------------------------------------
475     
476     ! Set reference phase to gas
477           M = 0
478     
479     
480     ! Set the default boundary conditions
481     ! The NS default setting is the where bc_type='dummy' or any default
482     ! (i.e., bc_type=undefined) wall boundary regions are handled. Note that
483     ! the top and bottom xy planes do not have to be explictly addressed for
484     ! the w-momentum equation. In this direction the velocities are defined
485     ! at the wall (due staggered grid). They are defined as zero for a
486     ! no penetration condition (see zero_norm_vel subroutine and code under
487     ! the ip_at_t branch in the above source routine).
488     ! ---------------------------------------------------------------->>>
489     
490     ! south xz plane
491           J1 = 1
492           DO K1 = kmin3,kmax3
493              DO I1 = imin3,imax3
494                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
495                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
496                 IJK = FUNIJK(I1,J1,K1)
497                 IF (NS_WALL_AT(IJK)) THEN
498                    A_M(IJK,E,M) = ZERO
499                    A_M(IJK,W,M) = ZERO
500                    A_M(IJK,N,M) = -ONE
501                    A_M(IJK,S,M) = ZERO
502                    A_M(IJK,T,M) = ZERO
503                    A_M(IJK,B,M) = ZERO
504                    A_M(IJK,0,M) = -ONE
505                    B_M(IJK,M) = ZERO
506                 ELSEIF (FS_WALL_AT(IJK)) THEN
507                    A_M(IJK,E,M) = ZERO
508                    A_M(IJK,W,M) = ZERO
509                    A_M(IJK,N,M) = ONE
510                    A_M(IJK,S,M) = ZERO
511                    A_M(IJK,T,M) = ZERO
512                    A_M(IJK,B,M) = ZERO
513                    A_M(IJK,0,M) = -ONE
514                    B_M(IJK,M) = ZERO
515                 ENDIF
516              ENDDO
517           ENDDO
518     
519     ! north xz plane
520           J1 = JMAX2
521           DO K1 = kmin3, kmax3
522              DO I1 = imin3, imax3
523                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
524                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
525                 IJK = FUNIJK(I1,J1,K1)
526                 IF (NS_WALL_AT(IJK)) THEN
527     ! Setting the wall velocity to zero (set the boundary cell value equal
528     ! and oppostive to the adjacent fluid cell value)
529                    A_M(IJK,E,M) = ZERO
530                    A_M(IJK,W,M) = ZERO
531                    A_M(IJK,N,M) = ZERO
532                    A_M(IJK,S,M) = -ONE
533                    A_M(IJK,T,M) = ZERO
534                    A_M(IJK,B,M) = ZERO
535                    A_M(IJK,0,M) = -ONE
536                    B_M(IJK,M) = ZERO
537                 ELSEIF (FS_WALL_AT(IJK)) THEN
538     ! Setting the wall velocity equal to the adjacent fluid velocity (set
539     ! the boundary cell value equal to adjacent fluid cell value)
540                    A_M(IJK,E,M) = ZERO
541                    A_M(IJK,W,M) = ZERO
542                    A_M(IJK,N,M) = ZERO
543                    A_M(IJK,S,M) = ONE
544                    A_M(IJK,T,M) = ZERO
545                    A_M(IJK,B,M) = ZERO
546                    A_M(IJK,0,M) = -ONE
547                    B_M(IJK,M) = ZERO
548                 ENDIF
549              ENDDO
550           ENDDO
551     
552     ! west yz plane
553           I1 = 1
554           DO K1 = kmin3, kmax3
555              DO J1 = jmin3, jmax3
556                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
557                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
558                 IJK = FUNIJK(I1,J1,K1)
559                 IF (NS_WALL_AT(IJK)) THEN
560                    A_M(IJK,E,M) = -ONE
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) = ZERO
565                    A_M(IJK,B,M) = ZERO
566                    A_M(IJK,0,M) = -ONE
567                    B_M(IJK,M) = ZERO
568                 ELSEIF (FS_WALL_AT(IJK)) THEN
569                    A_M(IJK,E,M) = ONE
570                    A_M(IJK,W,M) = ZERO
571                    A_M(IJK,N,M) = ZERO
572                    A_M(IJK,S,M) = ZERO
573                    A_M(IJK,T,M) = ZERO
574                    A_M(IJK,B,M) = ZERO
575                    A_M(IJK,0,M) = -ONE
576                    B_M(IJK,M) = ZERO
577                 ENDIF
578              ENDDO
579           ENDDO
580     
581     ! east yz plane
582           I1 = IMAX2
583           DO K1 = kmin3,kmax3
584              DO J1 = jmin3,jmax3
585                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
586                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
587                 IJK = FUNIJK(I1,J1,K1)
588                 IF (NS_WALL_AT(IJK)) THEN
589                    A_M(IJK,E,M) = ZERO
590                    A_M(IJK,W,M) = -ONE
591                    A_M(IJK,N,M) = ZERO
592                    A_M(IJK,S,M) = ZERO
593                    A_M(IJK,T,M) = ZERO
594                    A_M(IJK,B,M) = ZERO
595                    A_M(IJK,0,M) = -ONE
596                    B_M(IJK,M) = ZERO
597                 ELSEIF (FS_WALL_AT(IJK)) THEN
598                    A_M(IJK,E,M) = ZERO
599                    A_M(IJK,W,M) = ONE
600                    A_M(IJK,N,M) = ZERO
601                    A_M(IJK,S,M) = ZERO
602                    A_M(IJK,T,M) = ZERO
603                    A_M(IJK,B,M) = ZERO
604                    A_M(IJK,0,M) = -ONE
605                    B_M(IJK,M) = ZERO
606                 ENDIF
607              ENDDO
608           ENDDO
609     ! End setting the default boundary conditions
610     ! ----------------------------------------------------------------<<<
611     
612     ! Setting user specified boundary conditions
613     
614           DO L = 1, DIMENSION_BC
615              IF (BC_DEFINED(L)) THEN
616     
617     ! Setting wall boundary conditions
618     ! ---------------------------------------------------------------->>>
619                 IF (BC_TYPE(L) == 'NO_SLIP_WALL' .AND. .NOT. K_Epsilon) THEN
620                    I1 = BC_I_W(L)
621                    I2 = BC_I_E(L)
622                    J1 = BC_J_S(L)
623                    J2 = BC_J_N(L)
624                    K1 = BC_K_B(L)
625                    K2 = BC_K_T(L)
626                    DO K = K1, K2
627                       DO J = J1, J2
628                          DO I = I1, I2
629                            IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
630                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
631                             IJK = FUNIJK(I,J,K)
632                             IF (.NOT.WALL_AT(IJK)) CYCLE  !skip redefined cells
633                             A_M(IJK,E,M) = ZERO
634                             A_M(IJK,W,M) = ZERO
635                             A_M(IJK,N,M) = ZERO
636                             A_M(IJK,S,M) = ZERO
637                             A_M(IJK,T,M) = ZERO
638                             A_M(IJK,B,M) = ZERO
639                             A_M(IJK,0,M) = -ONE
640                             B_M(IJK,M) = ZERO
641                             IF (FLUID_AT(EAST_OF(IJK))) THEN
642                                A_M(IJK,E,M) = -ONE
643                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
644                                A_M(IJK,W,M) = -ONE
645                             ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
646                                A_M(IJK,N,M) = -ONE
647                             ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
648                                A_M(IJK,S,M) = -ONE
649                             ENDIF
650                          ENDDO
651                       ENDDO
652                    ENDDO
653     
654                 ELSEIF (BC_TYPE(L) == 'FREE_SLIP_WALL' .AND. .NOT. K_Epsilon) THEN
655                    I1 = BC_I_W(L)
656                    I2 = BC_I_E(L)
657                    J1 = BC_J_S(L)
658                    J2 = BC_J_N(L)
659                    K1 = BC_K_B(L)
660                    K2 = BC_K_T(L)
661                    DO K = K1, K2
662                       DO J = J1, J2
663                          DO I = I1, I2
664                            IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
665                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
666                             IJK = FUNIJK(I,J,K)
667                             IF (.NOT.WALL_AT(IJK)) CYCLE  !skip redefined cells
668                             A_M(IJK,E,M) = ZERO
669                             A_M(IJK,W,M) = ZERO
670                             A_M(IJK,N,M) = ZERO
671                             A_M(IJK,S,M) = ZERO
672                             A_M(IJK,T,M) = ZERO
673                             A_M(IJK,B,M) = ZERO
674                             A_M(IJK,0,M) = -ONE
675                             B_M(IJK,M) = ZERO
676                             IF (FLUID_AT(EAST_OF(IJK))) THEN
677                                A_M(IJK,E,M) = ONE
678                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
679                                A_M(IJK,W,M) = ONE
680                             ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
681                                A_M(IJK,N,M) = ONE
682                             ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
683                                A_M(IJK,S,M) = ONE
684                             ENDIF
685                          ENDDO
686                       ENDDO
687                    ENDDO
688     
689                 ELSEIF (BC_TYPE(L) == 'PAR_SLIP_WALL' .AND. .NOT. K_Epsilon) THEN
690                    I1 = BC_I_W(L)
691                    I2 = BC_I_E(L)
692                    J1 = BC_J_S(L)
693                    J2 = BC_J_N(L)
694                    K1 = BC_K_B(L)
695                    K2 = BC_K_T(L)
696                    DO K = K1, K2
697                       DO J = J1, J2
698                          DO I = I1, I2
699                            IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
700                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
701                             IJK = FUNIJK(I,J,K)
702                             IF (.NOT.WALL_AT(IJK)) CYCLE  ! skip redefined cells
703                             IM = IM1(I)
704                             JM = JM1(J)
705                             A_M(IJK,E,M) = ZERO
706                             A_M(IJK,W,M) = ZERO
707                             A_M(IJK,N,M) = ZERO
708                             A_M(IJK,S,M) = ZERO
709                             A_M(IJK,T,M) = ZERO
710                             A_M(IJK,B,M) = ZERO
711                             A_M(IJK,0,M) = -ONE
712                             B_M(IJK,M) = ZERO
713                             IF (FLUID_AT(EAST_OF(IJK))) THEN
714                                IF (BC_HW_G(L) == UNDEFINED) THEN
715                                   A_M(IJK,E,M) = -HALF
716                                   A_M(IJK,0,M) = -HALF
717                                   B_M(IJK,M) = -BC_WW_G(L)
718                                ELSE
719                                   IF (CYLINDRICAL) THEN
720                                      A_M(IJK,0,M) = -( HALF*(BC_HW_G(L)-&
721                                         OX_E(I)) + ODX_E(I) )
722                                      A_M(IJK,E,M) = -(HALF*(BC_HW_G(L)-&
723                                         OX_E(I)) - ODX_E(I))
724                                      B_M(IJK,M) = -BC_HW_G(L)*BC_WW_G(L)
725                                   ELSE
726                                      A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODX_E(I))
727                                      A_M(IJK,E,M) = -(HALF*BC_HW_G(L)-ODX_E(I))
728                                      B_M(IJK,M) = -BC_HW_G(L)*BC_WW_G(L)
729                                   ENDIF
730                                ENDIF
731                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
732                                IF (BC_HW_G(L) == UNDEFINED) THEN
733                                   A_M(IJK,W,M) = -HALF
734                                   A_M(IJK,0,M) = -HALF
735                                   B_M(IJK,M) = -BC_WW_G(L)
736                                ELSE
737                                   IF (CYLINDRICAL) THEN
738                                      A_M(IJK,W,M) = -(HALF*(BC_HW_G(L)-&
739                                         OX_E(IM)) - ODX_E(IM))
740                                      A_M(IJK,0,M) = -(HALF*(BC_HW_G(L)-&
741                                         OX_E(IM)) + ODX_E(IM))
742                                      B_M(IJK,M) = -BC_HW_G(L)*BC_WW_G(L)
743                                   ELSE
744                                      A_M(IJK,W,M) = -(HALF*BC_HW_G(L)-ODX_E(IM))
745                                      A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODX_E(IM))
746                                      B_M(IJK,M) = -BC_HW_G(L)*BC_WW_G(L)
747                                   ENDIF
748                                ENDIF
749                             ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
750                                IF (BC_HW_G(L) == UNDEFINED) THEN
751                                   A_M(IJK,N,M) = -HALF
752                                   A_M(IJK,0,M) = -HALF
753                                   B_M(IJK,M) = -BC_WW_G(L)
754                                ELSE
755                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODY_N(J))
756                                   A_M(IJK,N,M) = -(HALF*BC_HW_G(L)-ODY_N(J))
757                                   B_M(IJK,M) = -BC_HW_G(L)*BC_WW_G(L)
758                                ENDIF
759                             ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
760                                IF (BC_HW_G(L) == UNDEFINED) THEN
761                                   A_M(IJK,S,M) = -HALF
762                                   A_M(IJK,0,M) = -HALF
763                                   B_M(IJK,M) = -BC_WW_G(L)
764                                ELSE
765                                   A_M(IJK,S,M) = -(HALF*BC_HW_G(L)-ODY_N(JM))
766                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODY_N(JM))
767                                   B_M(IJK,M) = -BC_HW_G(L)*BC_WW_G(L)
768                                ENDIF
769                             ENDIF
770                          ENDDO
771                       ENDDO
772                    ENDDO
773     
774     ! Setting wall boundary conditions when K_EPSILON
775     ! wall functions for V-momentum are specify in this section of the code
776                 ELSEIF (BC_TYPE(L) == 'PAR_SLIP_WALL'   .OR.  &
777                         BC_TYPE(L) == 'NO_SLIP_WALL'    .OR.  &
778                         BC_TYPE(L) == 'FREE_SLIP_WALL'  .AND. &
779                         K_Epsilon )THEN
780                    I1 = BC_I_W(L)
781                    I2 = BC_I_E(L)
782                    J1 = BC_J_S(L)
783                    J2 = BC_J_N(L)
784                    K1 = BC_K_B(L)
785                    K2 = BC_K_T(L)
786                    DO K = K1, K2
787                       DO J = J1, J2
788                          DO I = I1, I2
789                            IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
790                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
791                             IJK = FUNIJK(I,J,K)
792                             IF (.NOT.WALL_AT(IJK)) CYCLE  !skip redefined cells
793                             IM = IM1(I)
794                             JM = JM1(J)
795                             A_M(IJK,E,M) = ZERO
796                             A_M(IJK,W,M) = ZERO
797                             A_M(IJK,N,M) = ZERO
798                             A_M(IJK,S,M) = ZERO
799                             A_M(IJK,T,M) = ZERO
800                             A_M(IJK,B,M) = ZERO
801                             A_M(IJK,0,M) = -ONE
802                             B_M(IJK,M) = ZERO
803                             IF (FLUID_AT(EAST_OF(IJK))) THEN
804                                IF (CYLINDRICAL) THEN
805                                   W_F_Slip = ( ONE/&
806                                      (ODX_E(I)+HALF*OX_E(I)) )*          &
807                                      ( ODX_E(I) - OX_E(I) -              &
808                                        RO_g(EAST_OF(IJK))*C_mu**0.25*    &
809                                        SQRT(K_Turb_G((EAST_OF(IJK))))/   &
810                                        MU_gT(EAST_OF(IJK))*Kappa/        &
811                                        LOG(9.81D0/ODX_E(I)/(2.D0)*       &
812                                        RO_g(EAST_OF(IJK))*C_mu**0.25*    &
813                                        SQRT(K_Turb_G((EAST_OF(IJK))))/   &
814                                        MU_g(EAST_OF(IJK))) )
815                                ELSE
816                                   CALL Wall_Function(IJK,EAST_OF(IJK),&
817                                      ODX_E(I),W_F_Slip)
818                                ENDIF
819                                A_M(IJK,E,M) = W_F_Slip
820                                A_M(IJK,0,M) = -ONE
821                                IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_WW_G(L)
822                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
823                                IF (CYLINDRICAL) THEN
824                                   W_F_Slip =  (ONE/&
825                                      (ONE*ODX_E(IM) + HALF*OX_E(IM)))*    &
826                                      ( ONE*ODX_E(IM) - OX_E(IM) -         &
827                                        RO_g(WEST_OF(IJK))*C_mu**0.25*     &
828                                        SQRT(K_Turb_G((WEST_OF(IJK))))/    &
829                                        MU_gT(WEST_OF(IJK))*Kappa/         &
830                                        LOG(9.81D0/ODX_E(IM)/(2.D0)*       &
831                                        RO_g(WEST_OF(IJK))*C_mu**0.25*     &
832                                        SQRT(K_Turb_G((WEST_OF(IJK))))/    &
833                                        MU_g(WEST_OF(IJK))) )
834                                ELSE
835                                   CALL Wall_Function(IJK,WEST_OF(IJK),&
836                                      ODX_E(IM),W_F_Slip)
837                                ENDIF
838                                A_M(IJK,W,M) = W_F_Slip
839                                A_M(IJK,0,M) = -ONE
840                                IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_WW_G(L)
841                             ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
842                                CALL Wall_Function(IJK,NORTH_OF(IJK),ODY_N(J),W_F_Slip)
843                                A_M(IJK,N,M) = W_F_Slip
844                                A_M(IJK,0,M) = -ONE
845                                IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_WW_G(L)
846                             ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
847                                CALL Wall_Function(IJK,SOUTH_OF(IJK),ODY_N(JM),W_F_Slip)
848                                A_M(IJK,S,M) = W_F_Slip
849                                A_M(IJK,0,M) = -ONE
850                                IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_WW_G(L)
851                             ENDIF
852                          ENDDO
853                       ENDDO
854                    ENDDO
855     
856     ! end setting of wall boundary conditions
857     ! ----------------------------------------------------------------<<<
858     
859     ! Setting p_inflow or p_outflow flow boundary conditions
860     ! ---------------------------------------------------------------->>>
861                 ELSEIF (BC_TYPE(L)=='P_INFLOW' .OR. BC_TYPE(L)=='P_OUTFLOW') THEN
862                    IF (BC_PLANE(L) == 'B') THEN
863     ! if the fluid cell is on the bottom side of the outflow/inflow boundary
864     ! then set the velocity in the boundary cell equal to the velocity of
865     ! the adjacent fluid cell
866                       I1 = BC_I_W(L)
867                       I2 = BC_I_E(L)
868                       J1 = BC_J_S(L)
869                       J2 = BC_J_N(L)
870                       K1 = BC_K_B(L)
871                       K2 = BC_K_T(L)
872                       DO K = K1, K2
873                          DO J = J1, J2
874                             DO I = I1, I2
875                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
876                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
877                                IJK = FUNIJK(I,J,K)
878                                A_M(IJK,E,M) = ZERO
879                                A_M(IJK,W,M) = ZERO
880                                A_M(IJK,N,M) = ZERO
881                                A_M(IJK,S,M) = ZERO
882                                A_M(IJK,T,M) = ZERO
883                                A_M(IJK,B,M) = ONE
884                                A_M(IJK,0,M) = -ONE
885                                B_M(IJK,M) = ZERO
886                             ENDDO
887                          ENDDO
888                       ENDDO
889                    ENDIF
890     ! end setting of p_inflow or p_otuflow flow boundary conditions
891     ! ----------------------------------------------------------------<<<
892     
893     ! Setting outflow flow boundary conditions
894     ! ---------------------------------------------------------------->>>
895                 ELSEIF (BC_TYPE(L) == 'OUTFLOW') THEN
896                    IF (BC_PLANE(L) == 'B') THEN
897                       I1 = BC_I_W(L)
898                       I2 = BC_I_E(L)
899                       J1 = BC_J_S(L)
900                       J2 = BC_J_N(L)
901                       K1 = BC_K_B(L)
902                       K2 = BC_K_T(L)
903                       DO K = K1, K2
904                          DO J = J1, J2
905                             DO I = I1, I2
906                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
907                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
908                                IJK = FUNIJK(I,J,K)
909                                A_M(IJK,E,M) = ZERO
910                                A_M(IJK,W,M) = ZERO
911                                A_M(IJK,N,M) = ZERO
912                                A_M(IJK,S,M) = ZERO
913                                A_M(IJK,T,M) = ZERO
914                                A_M(IJK,B,M) = ONE
915                                A_M(IJK,0,M) = -ONE
916                                B_M(IJK,M) = ZERO
917                                IJKM = KM_OF(IJK)
918                                A_M(IJKM,E,M) = ZERO
919                                A_M(IJKM,W,M) = ZERO
920                                A_M(IJKM,N,M) = ZERO
921                                A_M(IJKM,S,M) = ZERO
922                                A_M(IJKM,T,M) = ZERO
923                                A_M(IJKM,B,M) = ONE
924                                A_M(IJKM,0,M) = -ONE
925                                B_M(IJKM,M) = ZERO
926                             ENDDO
927                          ENDDO
928                       ENDDO
929                    ELSEIF (BC_PLANE(L) == 'T') THEN
930                       I1 = BC_I_W(L)
931                       I2 = BC_I_E(L)
932                       J1 = BC_J_S(L)
933                       J2 = BC_J_N(L)
934                       K1 = BC_K_B(L)
935                       K2 = BC_K_T(L)
936                       DO K = K1, K2
937                          DO J = J1, J2
938                             DO I = I1, I2
939                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
940                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
941                                IJK = FUNIJK(I,J,K)
942                                IJKP = KP_OF(IJK)
943                                A_M(IJKP,E,M) = ZERO
944                                A_M(IJKP,W,M) = ZERO
945                                A_M(IJKP,N,M) = ZERO
946                                A_M(IJKP,S,M) = ZERO
947                                A_M(IJKP,T,M) = ONE
948                                A_M(IJKP,B,M) = ZERO
949                                A_M(IJKP,0,M) = -ONE
950                                B_M(IJKP,M) = ZERO
951                             ENDDO
952                          ENDDO
953                       ENDDO
954                    ENDIF
955     ! end setting of outflow flow boundary conditions
956     ! ----------------------------------------------------------------<<<
957     
958     ! Setting bc that are defined but not nsw, fsw, psw, p_inflow,
959     ! p_outflow, or outflow (at this time, this section addresses
960     ! mass_inflow and mass_outflow type boundaries)
961     ! ---------------------------------------------------------------->>>
962                 ELSE
963                    I1 = BC_I_W(L)
964                    I2 = BC_I_E(L)
965                    J1 = BC_J_S(L)
966                    J2 = BC_J_N(L)
967                    K1 = BC_K_B(L)
968                    K2 = BC_K_T(L)
969                    DO K = K1, K2
970                       DO J = J1, J2
971                          DO I = I1, I2
972                            IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
973                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
974                             IJK = FUNIJK(I,J,K)
975     ! setting the velocity in the boundary cell equal to what is known
976                             A_M(IJK,E,M) = ZERO
977                             A_M(IJK,W,M) = ZERO
978                             A_M(IJK,N,M) = ZERO
979                             A_M(IJK,S,M) = ZERO
980                             A_M(IJK,T,M) = ZERO
981                             A_M(IJK,B,M) = ZERO
982                             A_M(IJK,0,M) = -ONE
983                             B_M(IJK,M) = -W_G(IJK)
984                             IF (BC_PLANE(L) == 'B') THEN
985     ! if the fluid cell is on the bottom side of the outflow/inflow boundary
986     ! then set the velocity in the adjacent fluid cell equal to what is
987     ! known in that cell
988                                IJKB = BOTTOM_OF(IJK)
989                                A_M(IJKB,E,M) = ZERO
990                                A_M(IJKB,W,M) = ZERO
991                                A_M(IJKB,N,M) = ZERO
992                                A_M(IJKB,S,M) = ZERO
993                                A_M(IJKB,T,M) = ZERO
994                                A_M(IJKB,B,M) = ZERO
995                                A_M(IJKB,0,M) = -ONE
996                                B_M(IJKB,M) = -W_G(IJKB)
997                             ENDIF
998                          ENDDO
999                       ENDDO
1000                    ENDDO
1001                 ENDIF   ! end if/else (bc_type)
1002                         ! ns, fs, psw; else
1003                         ! p_inflow, p_outflow, or outflow; else
1004     ! end setting of 'else' flow boundary conditions
1005     ! (mass_inflow/mass_outflow)
1006     ! ----------------------------------------------------------------<<<
1007     
1008              ENDIF   ! end if (bc_defined)
1009           ENDDO   ! end L do loop over dimension_bc
1010     
1011           RETURN
1012           END SUBROUTINE SOURCE_W_G_BC
1013     
1014     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1015     !                                                                      C
1016     !  Subroutine: POINT_SOURCE_W_G                                        C
1017     !  Purpose: Adds point sources to the gas phase W-Momentum equation.   C
1018     !                                                                      C
1019     !  Author: J. Musser                                  Date: 10-JUN-13  C
1020     !  Reviewer:                                          Date:            C
1021     !                                                                      C
1022     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1023           SUBROUTINE POINT_SOURCE_W_G(A_M, B_M)
1024     
1025     !-----------------------------------------------
1026     ! Modules
1027     !-----------------------------------------------
1028           use compar
1029           use constant
1030           use geometry
1031           use indices
1032           use physprop
1033           use ps
1034           use run
1035           use functions
1036           IMPLICIT NONE
1037     !-----------------------------------------------
1038     ! Dummy arguments
1039     !-----------------------------------------------
1040     ! Septadiagonal matrix A_m
1041           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1042     ! Vector b_m
1043           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1044     !-----------------------------------------------
1045     ! Local variables
1046     !-----------------------------------------------
1047     ! Indices
1048           INTEGER :: IJK, I, J, K
1049           INTEGER :: PSV, M
1050           INTEGER :: lKT, lKB
1051     ! terms of bm expression
1052           DOUBLE PRECISION :: pSource
1053     !-----------------------------------------------
1054     
1055     ! Set reference phase to gas
1056           M = 0
1057     
1058     ! Calculate the mass going into each IJK cell. This is done for each
1059     ! call in case the point source is time dependent.
1060           PS_LP: do PSV = 1, DIMENSION_PS
1061              if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
1062              if(abs(PS_W_g(PSV)) < small_number) cycle PS_LP
1063     
1064              if(PS_W_g(PSV) < ZERO) then
1065                 lKB = PS_K_B(PSV)-1
1066                 lKT = PS_K_T(PSV)-1
1067              else
1068                 lKB = PS_K_B(PSV)
1069                 lKT = PS_K_T(PSV)
1070              endif
1071     
1072              do k = lKB, lKT
1073              do j = PS_J_S(PSV), PS_J_N(PSV)
1074              do i = PS_I_W(PSV), PS_I_E(PSV)
1075     
1076                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
1077                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1078     
1079                 ijk = funijk(i,j,k)
1080                 if(.NOT.fluid_at(ijk)) cycle
1081     
1082                 pSource =  PS_MASSFLOW_G(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
1083     
1084                 B_M(IJK,M) = B_M(IJK,M) - pSource * &
1085                    PS_W_g(PSV) * PS_VEL_MAG_g(PSV)
1086     
1087              enddo
1088              enddo
1089              enddo
1090     
1091           enddo PS_LP
1092     
1093           RETURN
1094           END SUBROUTINE POINT_SOURCE_W_G
1095