File: /nfs/home/0/users/jenkins/mfix.git/model/source_u_g.f

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