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