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