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