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