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