File: N:\mfix\model\source_v_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_V_g                                              C
4     !  Purpose: Determine source terms for V_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: 7-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     !                                                                      C
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21           SUBROUTINE SOURCE_V_G(A_M, B_M)
22     
23     ! Modules
24     !---------------------------------------------------------------------//
25           USE bodyforce, only: bfy_g
26           USE bc, only: delp_y
27     
28           USE compar, only: ijkstart3, ijkend3, jmap
29     
30           USE drag, only: f_gs, beta_ij
31           USE fldvar, only: p_g, ro_g, rop_g, rop_go, rop_s
32           USE fldvar, only: ep_g, ep_s, epg_jfac
33           USE fldvar, only: v_g, v_go, u_s, v_s, w_s, v_so
34     
35           USE fun_avg, only: avg_x, avg_z, avg_y
36           USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
37           USE functions, only: ip_at_n, sip_at_n, is_id_at_n
38           USE functions, only: ip_of, jp_of, kp_of, im_of, jm_of, km_of
39           USE functions, only: north_of, south_of
40           USE functions, only: zmax
41           USE geometry, only: jmax1, cyclic_y_pd, do_k, xlength
42           USE geometry, only: vol, vol_v
43           USE geometry, only: axy, ayz, axz
44           USE geometry, only: ox, odx_e
45     
46           USE ghdtheory, only: joiy
47     
48           USE indices, only: i_of, j_of, k_of
49           USE indices, only: ip1, im1, kp1
50           USE is, only: is_pc
51     
52           USE mms, only: use_mms, mms_v_g_src
53           USE param
54           USE param1, only: zero, one, half
55           USE physprop, only: mmax, smax
56           USE physprop, only: mu_g, cv
57           USE run, only: momentum_y_eq
58           USE run, only: model_b, added_mass, m_am
59           USE run, only: kt_type_enum, drag_type_enum
60           USE run, only: ghd_2007, hys
61           USE run, only: odt
62           USE run, only: v_sh, shear
63           USE rxns, only: sum_r_g
64           USE scales, only: p_scale
65           USE vshear, only: vsh
66           USE tau_g, only: tau_v_g
67           USE toleranc, only: dil_ep_s
68           USE cutcell, only: cartesian_grid, cut_v_treatment_at
69           USE cutcell, only: blocked_v_cell_at
70           USE cutcell, only: a_vpg_n, a_vpg_s
71           IMPLICIT NONE
72     
73     ! Dummy Arguments
74     !---------------------------------------------------------------------//
75     ! Septadiagonal matrix A_m
76           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
77     ! Vector b_m
78           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
79     
80     ! Local variables
81     !---------------------------------------------------------------------//
82     ! Indices
83           INTEGER :: I, J, K, IJK, IJKN, &
84                      IMJK, IPJK, IJMK, IJPK, IJKP, IJKM, IMJPK, IJPKM
85     ! Phase index
86           INTEGER :: M, L, MM
87     ! Internal surface
88           INTEGER :: ISV
89     ! Pressure at north cell
90           DOUBLE PRECISION :: PgN
91     ! Average volume fraction
92           DOUBLE PRECISION :: EPGA, EPGAJ
93     ! Average density
94           DOUBLE PRECISION :: ROPGA, ROGA
95     ! Average viscosity
96           DOUBLE PRECISION :: MUGA
97     ! Source terms (Surface)
98           DOUBLE PRECISION :: Sdp
99     ! Source terms (Volumetric)
100           DOUBLE PRECISION :: V0, Vpm, Vmt, Vbf
101     ! Source terms (Volumetric) for GHD theory
102           DOUBLE PRECISION :: Ghd_drag, avgRop
103     ! Source terms for HYS drag relation
104           DOUBLE PRECISION :: HYS_drag, avgDrag
105     ! virtual (added) mass
106           DOUBLE PRECISION :: ROP_MA, Vsn, Vss, U_se, Usw, Vse, Vsw, &
107                               Wst, Wsb, Vst, Vsb
108           DOUBLE PRECISION :: F_vir
109     ! loezos: shear terms
110           DOUBLE PRECISION :: VSH_n,VSH_s,VSH_e,VSH_w,VSH_p,Source_conv
111           DOUBLE PRECISION :: SRT
112     ! jackson terms: local stress tensor quantity
113           DOUBLE PRECISION :: ltau_v_g
114     !---------------------------------------------------------------------//
115     
116     ! Set reference phase to gas
117           M = 0
118     
119           IF (.NOT.MOMENTUM_Y_EQ(0)) RETURN
120     
121     
122     !$omp  parallel do default(shared)                                   &
123     !$omp  private(I, J, K, IJK, IJKN, IMJK, IPJK, IJMK, IJPK, IMJPK,    &
124     !$omp          IJKM, IJPKM, IJKP, EPGA, EPGAJ, PGN, SDP, ROPGA,      &
125     !$omp          ROGA, ROP_MA, V0, ISV, MUGA, Vpm, Vmt, Vbf, L, MM,    &
126     !$omp          Vsn, Vss, U_se, Usw, Vse, Vsw, Wst, Wsb, Vst,         &
127     !$omp          Vsb, F_vir, Ghd_drag, avgRop, avgDrag, HYS_drag,      &
128     !$omp          VSH_n, VSH_s, VSH_e, VSH_w, VSH_p, Source_conv, SRT,  &
129     !$omp          ltau_v_g)
130           DO IJK = ijkstart3, ijkend3
131              I = I_OF(IJK)
132              J = J_OF(IJK)
133              K = K_OF(IJK)
134              IJKN = NORTH_OF(IJK)
135              IMJK = IM_OF(IJK)
136              IPJK = IP_OF(IJK)
137              IJMK = JM_OF(IJK)
138              IJPK = JP_OF(IJK)
139              IMJPK = IM_OF(IJPK)
140              IJKM = KM_OF(IJK)
141              IJPKM = KM_OF(IJPK)
142              IJKP = KP_OF(IJK)
143     
144              EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
145     ! if jackson avg ep_g otherwise 1
146              EPGAJ = AVG_Y(EPG_jfac(IJK),EPG_jfac(IJKN),J)
147     
148     ! Impermeable internal surface
149              IF (IP_AT_N(IJK)) THEN
150                 A_M(IJK,east,M) = ZERO
151                 A_M(IJK,west,M) = ZERO
152                 A_M(IJK,north,M) = ZERO
153                 A_M(IJK,south,M) = ZERO
154                 A_M(IJK,top,M) = ZERO
155                 A_M(IJK,bottom,M) = ZERO
156                 A_M(IJK,0,M) = -ONE
157                 B_M(IJK,M) = ZERO
158     
159     ! dilute flow
160              ELSEIF (EPGA <= DIL_EP_S) THEN
161                 A_M(IJK,east,M) = ZERO
162                 A_M(IJK,west,M) = ZERO
163                 A_M(IJK,north,M) = ZERO
164                 A_M(IJK,south,M) = ZERO
165                 A_M(IJK,top,M) = ZERO
166                 A_M(IJK,bottom,M) = ZERO
167                 A_M(IJK,0,M) = -ONE
168                 B_M(IJK,M) = ZERO
169                 IF (EP_G(SOUTH_OF(IJK)) > DIL_EP_S) THEN
170                    A_M(IJK,south,M) = ONE
171                 ELSE IF (EP_G(NORTH_OF(IJK)) > DIL_EP_S) THEN
172                    A_M(IJK,north,M) = ONE
173                 ELSE
174                    B_M(IJK,M) = -V_G(IJK)
175                 ENDIF
176     
177     ! Cartesian grid implementation
178              ELSEIF (BLOCKED_V_CELL_AT(IJK)) THEN
179                 A_M(IJK,east,M) = ZERO
180                 A_M(IJK,west,M) = ZERO
181                 A_M(IJK,north,M) = ZERO
182                 A_M(IJK,south,M) = ZERO
183                 A_M(IJK,top,M) = ZERO
184                 A_M(IJK,bottom,M) = ZERO
185                 A_M(IJK,0,M) = -ONE
186                 B_M(IJK,M) = ZERO
187     
188     ! Normal case
189              ELSE
190     
191     ! Surface forces
192     ! Pressure term
193                 PGN = P_G(IJKN)
194                 IF (CYCLIC_Y_PD) THEN
195                    IF (JMAP(J_OF(IJK)).EQ.JMAX1)PGN = P_G(IJKN) - DELP_Y
196                 ENDIF
197                 IF (MODEL_B) THEN
198                    IF(.NOT.CUT_V_TREATMENT_AT(IJK)) THEN
199                       SDP = -P_SCALE*(PGN - P_G(IJK))*AXZ(IJK)
200                    ELSE
201                       SDP = -P_SCALE*(PGN * A_VPG_N(IJK) - &
202                                       P_G(IJK) * A_VPG_S(IJK) )
203                    ENDIF
204                 ELSE
205                    IF(.NOT.CUT_V_TREATMENT_AT(IJK)) THEN
206                       SDP = -P_SCALE*EPGA*(PGN - P_G(IJK))*AXZ(IJK)
207                    ELSE
208                       SDP = -P_SCALE*EPGA*(PGN * A_VPG_N(IJK) - &
209                                            P_G(IJK) * A_VPG_S(IJK) )
210                    ENDIF
211                 ENDIF
212     
213                 IF(.NOT.CUT_V_TREATMENT_AT(IJK)) THEN
214     ! Volumetric forces
215                    ROPGA = AVG_Y(ROP_G(IJK),ROP_G(IJKN),J)
216                    ROGA = AVG_Y(RO_G(IJK),RO_G(IJKN),J)
217     ! Previous time step
218                    V0 = AVG_Y(ROP_GO(IJK),ROP_GO(IJKN),J)*ODT
219     ! Added mass implicit transient term {Cv eps rop_g dV/dt}
220                    IF(Added_Mass) THEN
221                      ROP_MA = AVG_Y(ROP_g(IJK)*EP_s(IJK,M_AM),&
222                                     ROP_g(IJKN)*EP_s(IJKN,M_AM),J)
223                      V0 = V0 + Cv * ROP_MA * ODT
224                    ENDIF
225                 ELSE
226     ! Volumetric forces
227                    ROPGA = (VOL(IJK)*ROP_G(IJK) + &
228                       VOL(IJKN)*ROP_G(IJKN))/(VOL(IJK) + VOL(IJKN))
229                    ROGA  = (VOL(IJK)*RO_G(IJK)  + &
230                       VOL(IJKN)*RO_G(IJKN) )/(VOL(IJK) + VOL(IJKN))
231     ! Previous time step
232                    V0 = (VOL(IJK)*ROP_GO(IJK) + VOL(IJKN)*ROP_GO(IJKN))*&
233                       ODT/(VOL(IJK) + VOL(IJKN))
234     ! Added mass implicit transient term {Cv eps rop_g dV/dt}
235                    IF(Added_Mass) THEN
236                      ROP_MA  = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM)  + &
237                         VOL(IJKN)*ROP_g(IJKN)*EP_s(IJKN,M_AM) )/&
238                         (VOL(IJK) + VOL(IJKN))
239                      V0 = V0 + Cv * ROP_MA * ODT
240                    ENDIF
241                 ENDIF
242     
243     ! VIRTUAL MASS SECTION (explicit terms)
244     ! adding transient term dVs/dt to virtual mass term
245                 F_vir = ZERO
246                 IF(Added_Mass.AND.(.NOT.CUT_V_TREATMENT_AT(IJK))) THEN
247                    F_vir = ((V_s(IJK,M_AM) - V_sO(IJK,M_AM)))*&
248                       ODT*VOL_V(IJK)
249     
250     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
251                    Vss = AVG_Y_N(V_S(IJMK,M_AM),V_s(IJK,M_AM))
252                    Vsn = AVG_Y_N(V_s(IJK,M_AM),V_s(IJPK,M_AM))
253                    U_se = AVG_Y(U_s(IJK,M_AM),U_s(IJPK,M_AM),J)
254                    Usw = AVG_Y(U_s(IMJK,M_AM),U_s(IMJPK,M_AM),J)
255                    Vse = AVG_X(V_s(IJK,M_AM),V_s(IPJK,M_AM),IP1(I))
256                    Vsw = AVG_X(V_s(IMJK,M_AM),V_s(IJK,M_AM),I)
257                    IF(DO_K) THEN
258                       Wst = AVG_Y(W_s(IJK,M_AM),W_s(IJPK,M_AM),J)
259                       Wsb = AVG_Y(W_s(IJKM,M_AM),W_s(IJPKM,M_AM),J)
260                       Vst = AVG_Z(V_s(IJK,M_AM),V_s(IJKP,M_AM),KP1(K))
261                       Vsb = AVG_Z(V_s(IJKM,M_AM),V_s(IJK,M_AM),K)
262                       F_vir = F_vir + AVG_Z_T(Wsb,Wst)*OX(I) * &
263                          (Vst - Vsb)*AXY(IJK)
264                    ENDIF
265     ! adding convective terms (U dV/dx + V dV/dy) to virtual mass; W dV/dz added above.
266                    F_vir = F_vir + V_s(IJK,M_AM)*(Vsn - Vss)*AXZ(IJK) + &
267                       AVG_X_E(Usw,U_se,IP1(I))*(Vse - Vsw)*AYZ(IJK)
268                    F_vir = F_vir * Cv * ROP_MA
269                 ENDIF
270     
271     ! pressure drop through porous media
272                 IF (SIP_AT_N(IJK)) THEN
273                    ISV = IS_ID_AT_N(IJK)
274                    MUGA = AVG_Y(MU_G(IJK),MU_G(IJKN),J)
275                    VPM = MUGA/IS_PC(ISV,1)
276                    IF (IS_PC(ISV,2) /= ZERO) VPM = VPM + HALF*IS_PC(ISV,2)*&
277                        ROPGA*ABS(V_G(IJK))
278                 ELSE
279                    VPM = ZERO
280                 ENDIF
281     
282     ! Interphase mass transfer
283                 IF(.NOT.CUT_V_TREATMENT_AT(IJK)) THEN
284                    VMT = AVG_Y(SUM_R_G(IJK),SUM_R_G(IJKN),J)
285                 ELSE
286                    VMT = (VOL(IJK)*SUM_R_G(IJK) + VOL(IJKN)*SUM_R_G(IJKN))/&
287                       (VOL(IJK) + VOL(IJKN))
288                 ENDIF
289     
290     ! Body force
291                 IF (MODEL_B) THEN
292                    VBF = ROGA*BFY_G(IJK)
293                 ELSE   ! Model A
294                    VBF = ROPGA*BFY_G(IJK)
295                 ENDIF
296     
297     ! Additional force for GHD from darg force sum(beta_ig * Joi/rhop_i)
298                 Ghd_drag = ZERO
299                 IF (KT_TYPE_ENUM .EQ. GHD_2007) THEN
300                   DO L = 1,SMAX
301                     avgRop = AVG_Y(ROP_S(IJK,L),ROP_S(IJKN,L),J)
302                     if(avgRop > ZERO) Ghd_drag = Ghd_drag +&
303                          AVG_Y(F_GS(IJK,L),F_GS(IJKN,L),J) * JoiY(IJK,L) / avgRop
304                   ENDDO
305                 ENDIF
306     
307     ! Additional force for HYS drag force, do not use with mixture GHD theory
308                 avgDrag = ZERO
309                 HYS_drag = ZERO
310                 IF (DRAG_TYPE_ENUM .EQ. HYS .AND. KT_TYPE_ENUM .NE. GHD_2007) THEN
311                    DO MM=1,MMAX
312                       DO L = 1,MMAX
313                          IF (L /= MM) THEN
314                             avgDrag = AVG_Y(beta_ij(IJK,MM,L),beta_ij(IJKN,MM,L),J)
315                             HYS_drag = HYS_drag + avgDrag * (V_g(IJK) - V_s(IJK,L))
316                          ENDIF
317                       ENDDO
318                    ENDDO
319                 ENDIF
320     ! if jackson, implement jackson form of governing equations (ep_g dot
321     ! del tau_g): multiply by void fraction otherwise by 1
322                 ltau_v_g = epgaj*tau_v_g(ijk)
323     
324     
325     ! loezos: Shear Source terms from convective mom. flux
326                 IF (SHEAR) THEN
327                     SRT=(2d0*V_sh/XLENGTH)
328                     VSH_p=VSH(IJK)
329                     VSH_n=VSH_p
330                     VSH_s=VSH_p
331                     VSH_e=VSH(IJK)+SRT*1d0/oDX_E(I)
332                     VSH_w=VSH(IJK)-SRT*1d0/oDX_E(IM1(I))
333                     Source_conv=A_M(IJK,north,m)*VSH_n + A_M(IJK,south,m)*VSH_s +&
334                        A_M(IJK,west,m)*VSH_w + A_M(IJK,east,m)*VSH_e - &
335                        (A_M(IJK,north,m) + A_M(IJK,south,m) + A_M(IJK,west,m) + &
336                         A_M(IJK,east,m))*VSH_p
337                  ELSE
338                     Source_conv=0d0
339                  ENDIF
340     
341     ! Collect the terms
342                 A_M(IJK,0,M) = -(A_M(IJK,east,M)+A_M(IJK,west,M)+&
343                    A_M(IJK,north,M)+A_M(IJK,south,M)+A_M(IJK,top,M)+A_M(IJK,bottom,M)+&
344                    (V0+VPM+ZMAX(VMT))*VOL_V(IJK))
345                 B_M(IJK,M) = B_M(IJK,M) - (SDP + lTAU_V_G + F_VIR + &
346                    Source_conv + ((V0+ZMAX((-VMT)))*V_GO(IJK) + VBF + &
347                    Ghd_drag + HYS_drag)*VOL_V(IJK) )
348     ! MMS Source term.
349                 IF(USE_MMS) B_M(IJK,M) = &
350                    B_M(IJK,M) - MMS_V_G_SRC(IJK)*VOL_V(IJK)
351     
352              ENDIF
353           ENDDO
354     !$omp end parallel do
355     
356     ! modifications for cartesian grid implementation
357           IF(CARTESIAN_GRID) CALL CG_SOURCE_V_G(A_M, B_M)
358     ! modifications for bc
359           CALL SOURCE_V_G_BC(A_M, B_M)
360     ! modifications for cartesian grid implementation
361           IF(CARTESIAN_GRID) CALL CG_SOURCE_V_G_BC(A_M, B_M)
362     
363           RETURN
364           END SUBROUTINE SOURCE_V_G
365     
366     
367     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
368     !                                                                      C
369     !  Subroutine: SOURCE_V_g_BC                                           C
370     !  Purpose: Determine source terms for V_g momentum eq. The terms      C
371     !     appear in the center coefficient and RHS vector. The center      C
372     !     coefficient and source vector are negative.  The off-diagonal    C
373     !     coefficients are positive.                                       C
374     !     The drag terms are excluded from the source at this stage        C
375     !                                                                      C
376     !  Author: M. Syamlal                                 Date: 7-JUN-96   C
377     !  Reviewer:                                          Date:            C
378     !                                                                      C
379     !                                                                      C
380     !                                                                      C
381     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
382     
383           SUBROUTINE SOURCE_V_G_BC(A_M, B_M)
384     
385     !-----------------------------------------------
386     ! Modules
387     !-----------------------------------------------
388           USE param
389           USE param1
390           USE parallel
391           USE scales
392           USE constant
393           USE physprop
394           USE fldvar
395           USE visc_g
396           USE rxns
397           USE run
398           USE toleranc
399           USE geometry
400           USE indices
401           USE is
402           USE tau_g
403           USE bc
404           USE output
405           USE compar
406           USE fun_avg
407           USE functions
408           IMPLICIT NONE
409     !-----------------------------------------------
410     ! Dummy Arguments
411     !-----------------------------------------------
412     ! Septadiagonal matrix A_m
413           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
414     ! Vector b_m
415           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
416     !-----------------------------------------------
417     ! Local Variables
418     !-----------------------------------------------
419     ! Boundary condition
420           INTEGER :: L
421     ! Indices
422           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK, &
423                      IM, KM, IJKS, IJMK, IJPK
424     ! Phase index
425           INTEGER :: M
426     ! Turbulent shear stress
427           DOUBLE PRECISION :: W_F_Slip
428     !-----------------------------------------------
429     
430     ! Set reference phase to gas
431           M = 0
432     
433     
434     ! Set the default boundary conditions
435     ! The NS default setting is the where bc_type='dummy' or any default
436     ! (i.e., bc_type=undefined) wall boundary regions are handled. Note that
437     ! the north and south xz planes do not have to be explicitly addressed for
438     ! the v-momentum equation. In this direction the velocities are defined
439     ! at the wall (due staggered grid). They are defined as zero for a
440     ! no penetration condition (see zero_norm_vel subroutine and code under
441     ! the ip_at_n branch in the above source routine).
442     ! ---------------------------------------------------------------->>>
443           IF (DO_K) THEN
444     ! bottom xy plane
445              K1 = 1
446              DO J1 = jmin3,jmax3
447                 DO I1 = imin3, imax3
448                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
449                    IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
450                    IJK = FUNIJK(I1,J1,K1)
451                    IF (NS_WALL_AT(IJK)) THEN
452     ! Setting the wall velocity to zero (set the boundary cell value equal
453     ! and oppostive to the adjacent fluid cell value)
454                       A_M(IJK,east,M) = ZERO
455                       A_M(IJK,west,M) = ZERO
456                       A_M(IJK,north,M) = ZERO
457                       A_M(IJK,south,M) = ZERO
458                       A_M(IJK,top,M) = -ONE
459                       A_M(IJK,bottom,M) = ZERO
460                       A_M(IJK,0,M) = -ONE
461                       B_M(IJK,M) = ZERO
462                    ELSEIF (FS_WALL_AT(IJK)) THEN
463     ! Setting the wall velocity equal to the adjacent fluid velocity (set
464     ! the boundary cell value equal to adjacent fluid cell value)
465                       A_M(IJK,east,M) = ZERO
466                       A_M(IJK,west,M) = ZERO
467                       A_M(IJK,north,M) = ZERO
468                       A_M(IJK,south,M) = ZERO
469                       A_M(IJK,top,M) = ONE
470                       A_M(IJK,bottom,M) = ZERO
471                       A_M(IJK,0,M) = -ONE
472                       B_M(IJK,M) = ZERO
473                    ENDIF
474                 ENDDO
475              ENDDO
476     
477     ! top xy plane
478              K1 = KMAX2
479              DO J1 = jmin3,jmax3
480                 DO I1 = imin3, imax3
481                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
482                    IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
483                    IJK = FUNIJK(I1,J1,K1)
484                    IF (NS_WALL_AT(IJK)) THEN
485                       A_M(IJK,east,M) = ZERO
486                       A_M(IJK,west,M) = ZERO
487                       A_M(IJK,north,M) = ZERO
488                       A_M(IJK,south,M) = ZERO
489                       A_M(IJK,top,M) = ZERO
490                       A_M(IJK,bottom,M) = -ONE
491                       A_M(IJK,0,M) = -ONE
492                       B_M(IJK,M) = ZERO
493                    ELSE IF (FS_WALL_AT(IJK)) THEN
494                       A_M(IJK,east,M) = ZERO
495                       A_M(IJK,west,M) = ZERO
496                       A_M(IJK,north,M) = ZERO
497                       A_M(IJK,south,M) = ZERO
498                       A_M(IJK,top,M) = ZERO
499                       A_M(IJK,bottom,M) = ONE
500                       A_M(IJK,0,M) = -ONE
501                       B_M(IJK,M) = ZERO
502                    ENDIF
503                 ENDDO
504              ENDDO
505           ENDIF   ! end if (do_k)
506     
507     
508     ! west zy plane
509           I1 = 1
510           DO K1 = kmin3, kmax3
511              DO J1 = jmin3, jmax3
512                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
513                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
514                 IJK = FUNIJK(I1,J1,K1)
515                 IF (NS_WALL_AT(IJK)) THEN
516                    A_M(IJK,east,M) = -ONE
517                    A_M(IJK,west,M) = ZERO
518                    A_M(IJK,north,M) = ZERO
519                    A_M(IJK,south,M) = ZERO
520                    A_M(IJK,top,M) = ZERO
521                    A_M(IJK,bottom,M) = ZERO
522                    A_M(IJK,0,M) = -ONE
523                    B_M(IJK,M) = ZERO
524                 ELSEIF (FS_WALL_AT(IJK)) THEN
525                    A_M(IJK,east,M) = ONE
526                    A_M(IJK,west,M) = ZERO
527                    A_M(IJK,north,M) = ZERO
528                    A_M(IJK,south,M) = ZERO
529                    A_M(IJK,top,M) = ZERO
530                    A_M(IJK,bottom,M) = ZERO
531                    A_M(IJK,0,M) = -ONE
532                    B_M(IJK,M) = ZERO
533                 ENDIF
534              ENDDO
535           ENDDO
536     
537     ! east zy plane
538           I1 = IMAX2
539           DO K1 = kmin3, kmax3
540              DO J1 = jmin3, jmax3
541                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
542                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
543                 IJK = FUNIJK(I1,J1,K1)
544                 IF (NS_WALL_AT(IJK)) THEN
545                    A_M(IJK,east,M) = ZERO
546                    A_M(IJK,west,M) = -ONE
547                    A_M(IJK,north,M) = ZERO
548                    A_M(IJK,south,M) = ZERO
549                    A_M(IJK,top,M) = ZERO
550                    A_M(IJK,bottom,M) = ZERO
551                    A_M(IJK,0,M) = -ONE
552                    B_M(IJK,M) = ZERO
553                 ELSEIF (FS_WALL_AT(IJK)) THEN
554                    A_M(IJK,east,M) = ZERO
555                    A_M(IJK,west,M) = ONE
556                    A_M(IJK,north,M) = ZERO
557                    A_M(IJK,south,M) = ZERO
558                    A_M(IJK,top,M) = ZERO
559                    A_M(IJK,bottom,M) = ZERO
560                    A_M(IJK,0,M) = -ONE
561                    B_M(IJK,M) = ZERO
562                 ENDIF
563              ENDDO
564           ENDDO
565     ! End setting the default boundary conditions
566     ! ----------------------------------------------------------------<<<
567     
568     ! Setting user specified boundary conditions
569           DO L = 1, DIMENSION_BC
570              IF (BC_DEFINED(L)) THEN
571     
572     ! Setting wall boundary conditions
573     ! ---------------------------------------------------------------->>>
574                 IF (BC_TYPE_ENUM(L) == NO_SLIP_WALL .AND. .NOT. K_Epsilon) THEN
575                    I1 = BC_I_W(L)
576                    I2 = BC_I_E(L)
577                    J1 = BC_J_S(L)
578                    J2 = BC_J_N(L)
579                    K1 = BC_K_B(L)
580                    K2 = BC_K_T(L)
581                    DO K = K1, K2
582                       DO J = J1, J2
583                          DO I = I1, I2
584                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
585                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
586                             IJK = FUNIJK(I,J,K)
587                             IF (.NOT.WALL_AT(IJK)) CYCLE  ! skip redefined cells
588                             A_M(IJK,east,M) = ZERO
589                             A_M(IJK,west,M) = ZERO
590                             A_M(IJK,north,M) = ZERO
591                             A_M(IJK,south,M) = ZERO
592                             A_M(IJK,top,M) = ZERO
593                             A_M(IJK,bottom,M) = ZERO
594                             A_M(IJK,0,M) = -ONE
595                             B_M(IJK,M) = ZERO
596                             IF (FLUID_AT(EAST_OF(IJK))) THEN
597                                A_M(IJK,east,M) = -ONE
598                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
599                                A_M(IJK,west,M) = -ONE
600                             ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
601                                A_M(IJK,top,M) = -ONE
602                             ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
603                                A_M(IJK,bottom,M) = -ONE
604                             ENDIF
605                          ENDDO
606                       ENDDO
607                    ENDDO
608     
609                 ELSEIF (BC_TYPE_ENUM(L) == FREE_SLIP_WALL .AND. .NOT. K_Epsilon) THEN
610                    I1 = BC_I_W(L)
611                    I2 = BC_I_E(L)
612                    J1 = BC_J_S(L)
613                    J2 = BC_J_N(L)
614                    K1 = BC_K_B(L)
615                    K2 = BC_K_T(L)
616                    DO K = K1, K2
617                       DO J = J1, J2
618                          DO I = I1, I2
619                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
620                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
621                             IJK = FUNIJK(I,J,K)
622                             IF (.NOT.WALL_AT(IJK)) CYCLE  ! skip redefined cells
623                             A_M(IJK,east,M) = ZERO
624                             A_M(IJK,west,M) = ZERO
625                             A_M(IJK,north,M) = ZERO
626                             A_M(IJK,south,M) = ZERO
627                             A_M(IJK,top,M) = ZERO
628                             A_M(IJK,bottom,M) = ZERO
629                             A_M(IJK,0,M) = -ONE
630                             B_M(IJK,M) = ZERO
631                             IF (FLUID_AT(EAST_OF(IJK))) THEN
632                                A_M(IJK,east,M) = ONE
633                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
634                                A_M(IJK,west,M) = ONE
635                             ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
636                                A_M(IJK,top,M) = ONE
637                             ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
638                                A_M(IJK,bottom,M) = ONE
639                             ENDIF
640                          ENDDO
641                       ENDDO
642                    ENDDO
643     
644                 ELSEIF (BC_TYPE_ENUM(L) == PAR_SLIP_WALL .AND. .NOT. K_Epsilon) THEN
645                    I1 = BC_I_W(L)
646                    I2 = BC_I_E(L)
647                    J1 = BC_J_S(L)
648                    J2 = BC_J_N(L)
649                    K1 = BC_K_B(L)
650                    K2 = BC_K_T(L)
651                    DO K = K1, K2
652                       DO J = J1, J2
653                          DO I = I1, I2
654                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
655                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
656                             IJK = FUNIJK(I,J,K)
657                             IF (.NOT.WALL_AT(IJK)) CYCLE  ! skip redefined cells
658                             IM = IM1(I)
659                             KM = KM1(K)
660                             A_M(IJK,east,M) = ZERO
661                             A_M(IJK,west,M) = ZERO
662                             A_M(IJK,north,M) = ZERO
663                             A_M(IJK,south,M) = ZERO
664                             A_M(IJK,top,M) = ZERO
665                             A_M(IJK,bottom,M) = ZERO
666                             A_M(IJK,0,M) = -ONE
667                             B_M(IJK,M) = ZERO
668                             IF (FLUID_AT(EAST_OF(IJK))) THEN
669                                IF (BC_HW_G(L) == UNDEFINED) THEN
670                                   A_M(IJK,east,M) = -HALF
671                                   A_M(IJK,0,M) = -HALF
672                                   B_M(IJK,M) = -BC_VW_G(L)
673                                ELSE
674                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODX_E(I))
675                                   A_M(IJK,east,M) = -(HALF*BC_HW_G(L)-ODX_E(I))
676                                   B_M(IJK,M) = -BC_HW_G(L)*BC_VW_G(L)
677                                ENDIF
678                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
679                                IF (BC_HW_G(L) == UNDEFINED) THEN
680                                   A_M(IJK,west,M) = -HALF
681                                   A_M(IJK,0,M) = -HALF
682                                   B_M(IJK,M) = -BC_VW_G(L)
683                                ELSE
684                                   A_M(IJK,west,M) = -(HALF*BC_HW_G(L)-ODX_E(IM))
685                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODX_E(IM))
686                                   B_M(IJK,M) = -BC_HW_G(L)*BC_VW_G(L)
687                                ENDIF
688                             ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
689                                IF (BC_HW_G(L) == UNDEFINED) THEN
690                                   A_M(IJK,top,M) = -HALF
691                                   A_M(IJK,0,M) = -HALF
692                                   B_M(IJK,M) = -BC_VW_G(L)
693                                ELSE
694                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODZ_T(K)*OX(I))
695                                   A_M(IJK,top,M) = -(HALF*BC_HW_G(L)-ODZ_T(K)*OX(I))
696                                   B_M(IJK,M) = -BC_HW_G(L)*BC_VW_G(L)
697                                ENDIF
698                             ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
699                                IF (BC_HW_G(L) == UNDEFINED) THEN
700                                   A_M(IJK,bottom,M) = -HALF
701                                   A_M(IJK,0,M) = -HALF
702                                   B_M(IJK,M) = -BC_VW_G(L)
703                                ELSE
704                                   A_M(IJK,bottom,M) = -(HALF*BC_HW_G(L)-ODZ_T(KM)*OX(I))
705                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODZ_T(KM)*OX(I))
706                                   B_M(IJK,M) = -BC_HW_G(L)*BC_VW_G(L)
707                                ENDIF
708                             ENDIF
709                          ENDDO
710                       ENDDO
711                    ENDDO
712     
713     ! Setting wall boundary conditions when K_EPSILON
714                 ELSEIF (BC_TYPE_ENUM(L) == PAR_SLIP_WALL   .OR.  &
715                         BC_TYPE_ENUM(L) == NO_SLIP_WALL    .OR.  &
716                         BC_TYPE_ENUM(L) == FREE_SLIP_WALL  .AND. &
717                         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                             KM = KM1(K)
733                             A_M(IJK,east,M) = ZERO
734                             A_M(IJK,west,M) = ZERO
735                             A_M(IJK,north,M) = ZERO
736                             A_M(IJK,south,M) = ZERO
737                             A_M(IJK,top,M) = ZERO
738                             A_M(IJK,bottom,M) = ZERO
739                             A_M(IJK,0,M) = -ONE
740                             B_M(IJK,M) = ZERO
741                             IF (FLUID_AT(EAST_OF(IJK))) THEN
742                                  CALL Wall_Function(IJK,EAST_OF(IJK),ODX_E(I),W_F_Slip)
743                                  A_M(IJK,east,M) = W_F_Slip
744                                  A_M(IJK,0,M) = -ONE
745                                  IF (BC_TYPE_ENUM(L) == PAR_SLIP_WALL) B_M(IJK,M) = -BC_VW_G(L)
746                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
747                                  CALL Wall_Function(IJK,WEST_OF(IJK),ODX_E(IM),W_F_Slip)
748                                  A_M(IJK,west,M) = W_F_Slip
749                                  A_M(IJK,0,M) = -ONE
750                                  IF (BC_TYPE_ENUM(L) == PAR_SLIP_WALL) B_M(IJK,M) = -BC_VW_G(L)
751                             ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
752                                  CALL Wall_Function(IJK,TOP_OF(IJK),ODZ_T(K)*OX(I),W_F_Slip)
753                                  A_M(IJK,top,M) = W_F_Slip
754                                  A_M(IJK,0,M) = -ONE
755                                  IF (BC_TYPE_ENUM(L) == PAR_SLIP_WALL) B_M(IJK,M) = -BC_VW_G(L)
756                             ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
757                                  CALL Wall_Function(IJK,BOTTOM_OF(IJK),ODZ_T(KM)*OX(I),W_F_Slip)
758                                  A_M(IJK,bottom,M) = W_F_Slip
759                                  A_M(IJK,0,M) = -ONE
760                                  IF (BC_TYPE_ENUM(L) == PAR_SLIP_WALL) B_M(IJK,M) = -BC_VW_G(L)
761                             ENDIF
762                          ENDDO
763                       ENDDO
764                    ENDDO
765     ! end setting of wall boundary conditions
766     ! ----------------------------------------------------------------<<<
767     
768     ! Setting p_inflow or p_outflow flow boundary conditions
769     ! ---------------------------------------------------------------->>>
770                 ELSE IF (BC_TYPE_ENUM(L)==P_INFLOW .OR. BC_TYPE_ENUM(L)==P_OUTFLOW) THEN
771                    IF (BC_PLANE(L) == 'S') THEN
772     ! if the fluid cell is on the south side of the outflow/inflow boundary
773     ! then set the velocity in the boundary cell equal to the velocity of
774     ! the adjacent fluid cell
775                       I1 = BC_I_W(L)
776                       I2 = BC_I_E(L)
777                       J1 = BC_J_S(L)
778                       J2 = BC_J_N(L)
779                       K1 = BC_K_B(L)
780                       K2 = BC_K_T(L)
781                       DO K = K1, K2
782                          DO J = J1, J2
783                             DO I = I1, I2
784                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
785                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
786                                IJK = FUNIJK(I,J,K)
787                                A_M(IJK,east,M) = ZERO
788                                A_M(IJK,west,M) = ZERO
789                                A_M(IJK,north,M) = ZERO
790                                A_M(IJK,south,M) = ONE
791                                A_M(IJK,top,M) = ZERO
792                                A_M(IJK,bottom,M) = ZERO
793                                A_M(IJK,0,M) = -ONE
794                                B_M(IJK,M) = ZERO
795                             ENDDO
796                          ENDDO
797                       ENDDO
798                    ENDIF
799     ! end setting of p_inflow or p_otuflow flow boundary conditions
800     ! ----------------------------------------------------------------<<<
801     
802     ! Setting outflow flow boundary conditions
803     ! ---------------------------------------------------------------->>>
804                 ELSEIF (BC_TYPE_ENUM(L) == OUTFLOW) THEN
805                    IF (BC_PLANE(L) == 'S') 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                                A_M(IJK,east,M) = ZERO
819                                A_M(IJK,west,M) = ZERO
820                                A_M(IJK,north,M) = ZERO
821                                A_M(IJK,south,M) = ONE
822                                A_M(IJK,top,M) = ZERO
823                                A_M(IJK,bottom,M) = ZERO
824                                A_M(IJK,0,M) = -ONE
825                                B_M(IJK,M) = ZERO
826                                IJMK = JM_OF(IJK)
827                                A_M(IJMK,east,M) = ZERO
828                                A_M(IJMK,west,M) = ZERO
829                                A_M(IJMK,north,M) = ZERO
830                                A_M(IJMK,south,M) = ONE
831                                A_M(IJMK,top,M) = ZERO
832                                A_M(IJMK,bottom,M) = ZERO
833                                A_M(IJMK,0,M) = -ONE
834                                B_M(IJMK,M) = ZERO
835                             ENDDO
836                          ENDDO
837                       ENDDO
838                    ELSEIF (BC_PLANE(L) == 'N') THEN
839                       I1 = BC_I_W(L)
840                       I2 = BC_I_E(L)
841                       J1 = BC_J_S(L)
842                       J2 = BC_J_N(L)
843                       K1 = BC_K_B(L)
844                       K2 = BC_K_T(L)
845                       DO K = K1, K2
846                          DO J = J1, J2
847                             DO I = I1, I2
848                                IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
849                                IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
850                                IJK = FUNIJK(I,J,K)
851                                IJPK = JP_OF(IJK)
852                                A_M(IJPK,east,M) = ZERO
853                                A_M(IJPK,west,M) = ZERO
854                                A_M(IJPK,north,M) = ONE
855                                A_M(IJPK,south,M) = ZERO
856                                A_M(IJPK,top,M) = ZERO
857                                A_M(IJPK,bottom,M) = ZERO
858                                A_M(IJPK,0,M) = -ONE
859                                B_M(IJPK,M) = ZERO
860                             ENDDO
861                          ENDDO
862                       ENDDO
863                    ENDIF
864     ! end setting of outflow flow boundary conditions
865     ! ----------------------------------------------------------------<<<
866     
867     ! Setting bc that are defined but not nsw, fsw, psw, p_inflow,
868     ! p_outflow, or outflow (at this time, this section addresses
869     ! mass_inflow and mass_outflow type boundaries)
870     ! ---------------------------------------------------------------->>>
871                 ELSE
872                    I1 = BC_I_W(L)
873                    I2 = BC_I_E(L)
874                    J1 = BC_J_S(L)
875                    J2 = BC_J_N(L)
876                    K1 = BC_K_B(L)
877                    K2 = BC_K_T(L)
878                    DO K = K1, K2
879                       DO J = J1, J2
880                          DO I = I1, I2
881                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
882                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
883                             IJK = FUNIJK(I,J,K)
884     ! setting the velocity in the boundary cell equal to what is known
885                             A_M(IJK,east,M) = ZERO
886                             A_M(IJK,west,M) = ZERO
887                             A_M(IJK,north,M) = ZERO
888                             A_M(IJK,south,M) = ZERO
889                             A_M(IJK,top,M) = ZERO
890                             A_M(IJK,bottom,M) = ZERO
891                             A_M(IJK,0,M) = -ONE
892                             B_M(IJK,M) = -V_G(IJK)
893                             IF (BC_PLANE(L) == 'S') THEN
894     ! if the fluid cell is on the south side of the outflow/inflow boundary
895     ! then set the velocity in the adjacent fluid cell equal to what is
896     ! known in that cell
897                                IJKS = SOUTH_OF(IJK)
898                                A_M(IJKS,east,M) = ZERO
899                                A_M(IJKS,west,M) = ZERO
900                                A_M(IJKS,north,M) = ZERO
901                                A_M(IJKS,south,M) = ZERO
902                                A_M(IJKS,top,M) = ZERO
903                                A_M(IJKS,bottom,M) = ZERO
904                                A_M(IJKS,0,M) = -ONE
905                                B_M(IJKS,M) = -V_G(IJKS)
906                             ENDIF
907                          ENDDO
908                       ENDDO
909                    ENDDO
910                 ENDIF   ! end if/else (bc_type)
911                         ! ns, fs, psw; else
912                         ! p_inflow, p_outflow, or outflow; else
913     ! end setting of 'else' flow boundary conditions
914     ! (mass_inflow/mass_outflow)
915     ! ----------------------------------------------------------------<<<
916     
917              ENDIF   ! end if (bc_defined)
918           ENDDO   ! end L do loop over dimension_bc
919     
920           RETURN
921           END SUBROUTINE SOURCE_V_G_BC
922     
923     
924     
925     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
926     !                                                                      C
927     !  Subroutine: POINT_SOURCE_V_G                                        C
928     !  Purpose: Adds point sources to the gas phase V-Momentum equation.   C
929     !                                                                      C
930     !  Author: J. Musser                                  Date: 10-JUN-13  C
931     !  Reviewer:                                          Date:            C
932     !                                                                      C
933     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
934           SUBROUTINE POINT_SOURCE_V_G(A_M, B_M)
935     
936           use compar
937           use constant
938           use geometry
939           use indices
940           use param
941           use param1, only: small_number
942           use physprop
943           use ps
944           use run
945           use functions
946     
947           IMPLICIT NONE
948     !-----------------------------------------------
949     ! Dummy arguments
950     !-----------------------------------------------
951     ! Septadiagonal matrix A_m
952           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
953     ! Vector b_m
954           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
955     !-----------------------------------------------
956     ! Local Variables
957     !-----------------------------------------------
958     ! Indices
959           INTEGER :: IJK, I, J, K
960           INTEGER :: PSV, M
961           INTEGER :: lJN, lJS
962     ! terms of bm expression
963           DOUBLE PRECISION :: pSource
964     !-----------------------------------------------
965     
966     ! Set reference phase to gas
967           M = 0
968     
969     ! Calculate the mass going into each IJK cell. This is done for each
970     ! call in case the point source is time dependent.
971           PS_LP: do PSV = 1, DIMENSION_PS
972              if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
973              if(abs(PS_V_g(PSV)) < small_number) cycle PS_LP
974     
975              if(PS_V_g(PSV) < 0.0d0) then
976                 lJS = PS_J_S(PSV) - 1
977                 lJN = PS_J_N(PSV) - 1
978              else
979                 lJS = PS_J_S(PSV)
980                 lJN = PS_J_N(PSV)
981              endif
982     
983              do k = PS_K_B(PSV), PS_K_T(PSV)
984              do j = lJS, lJN
985              do i = PS_I_W(PSV), PS_I_E(PSV)
986     
987                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
988                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
989     
990                 ijk = funijk(i,j,k)
991                 if(.NOT.fluid_at(ijk)) cycle
992     
993                 pSource =  PS_MASSFLOW_G(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
994     
995                 B_M(IJK,M) = B_M(IJK,M) - pSource * &
996                    PS_V_g(PSV) * PS_VEL_MAG_g(PSV)
997     
998              enddo
999              enddo
1000              enddo
1001     
1002           enddo PS_LP
1003     
1004           RETURN
1005           END SUBROUTINE POINT_SOURCE_V_G
1006