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