MFIX  2016-1
source_w_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: SOURCE_W_g C
4 ! Purpose: Determine source terms for W_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: 17-JUN-96 C
12 ! Reviewer: Date: C
13 ! C
14 ! Revision Number: 1 C
15 ! Purpose: To incorporate Cartesian grid modifications C
16 ! Author: Jeff Dietiker Date: 01-Jul-09 C
17 ! C
18 ! C
19 ! Literature/Document References: C
20 ! C
21 ! C
22 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
23  SUBROUTINE source_w_g(A_M, B_M)
24 
25 ! Modules
26 !---------------------------------------------------------------------//
27  USE bodyforce, only: bfz_g
28  USE bc, only: delp_z
29 
30  USE compar, only: ijkstart3, ijkend3, kmap
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, w_go, u_s, v_s, w_s, w_so
36 
37  USE fun_avg, only: avg_x, avg_z, avg_y
38  USE fun_avg, only: avg_x_h, avg_z_h
39  USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
40  USE functions, only: ip_at_t, sip_at_t, is_id_at_t
41  USE functions, only: ip_of, jp_of, kp_of, im_of, jm_of, km_of
42  USE functions, only: east_of, west_of, top_of, bottom_of
43  USE functions, only: zmax
44  USE geometry, only: kmax1, cyclic_z_pd, cylindrical
45  USE geometry, only: vol, vol_w
46  USE geometry, only: axy, ayz, axz, ayz_w
47  USE geometry, only: ox, ox_e, dy, dz, odx_e
48 
49  USE ghdtheory, only: joiz
50 
51  USE indices, only: i_of, j_of, k_of
52  USE indices, only: ip1, im1, jm1, kp1
53  USE is, only: is_pc
54 
55  USE mms, only: use_mms, mms_w_g_src
56  USE param
57  USE param1, only: zero, one, half
58  USE physprop, only: mmax, smax
59  USE physprop, only: mu_g, cv
60  USE run, only: momentum_z_eq
61  USE run, only: model_b, added_mass, m_am
62  USE run, only: kt_type_enum, drag_type_enum
63  USE run, only: ghd_2007, hys
64  USE run, only: odt
65  USE rxns, only: sum_r_g
66  USE scales, only: p_scale
67  USE tau_g, only: tau_w_g
68  USE toleranc, only: dil_ep_s
69  USE visc_g, only: epmu_gt
70  USE cutcell, only: cartesian_grid, cut_w_treatment_at
71  USE cutcell, only: blocked_w_cell_at
72  USE cutcell, only: a_wpg_t, a_wpg_b
73  IMPLICIT NONE
74 
75 ! Dummy arguments
76 !---------------------------------------------------------------------//
77 ! Septadiagonal matrix A_m
78  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
79 ! Vector b_m
80  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
81 
82 ! Local variables
83 !---------------------------------------------------------------------//
84 ! Indices
85  INTEGER :: I, J, K, IJK, IJKT, IMJK, IJKP, IMJKP,&
86  IJKE, IJKW, IJKTE, IJKTW, IM, IPJK, &
87  IJKM, IJMK, IJMKP, IJPK
88 ! Phase index
89  INTEGER :: M, L, MM
90 ! Internal surface
91  INTEGER :: ISV
92 ! Pressure at top cell
93  DOUBLE PRECISION :: PgT
94 ! Average volume fraction
95  DOUBLE PRECISION :: EPGA, EPGAJ
96 ! Average density
97  DOUBLE PRECISION :: ROPGA, ROGA
98 ! Average viscosity
99  DOUBLE PRECISION :: MUGA
100 ! Average coefficients
101  DOUBLE PRECISION :: Cte, Ctw, MUoX, cpe, cpw
102 ! Average U_g
103  DOUBLE PRECISION Ugt
104 ! Source terms (Surface)
105  DOUBLE PRECISION Sdp
106 ! Source terms (Volumetric)
107  DOUBLE PRECISION V0, Vpm, Vmt, Vbf, Vcoa, Vcob, Vxza
108 ! Source terms (Volumetric) for GHD theory
109  DOUBLE PRECISION Ghd_drag, avgRop
110 ! Source terms for HYS drag relation
111  DOUBLE PRECISION HYS_drag, avgDrag
112 ! virtual (added) mass
113  DOUBLE PRECISION :: ROP_MA, U_se, Usw, Ust, Vsb, Vst, &
114  Wse, Wsw, Wsn, Wss, Wst, Wsb
115  DOUBLE PRECISION F_vir
116 ! jackson terms: local stress tensor quantity
117  DOUBLE PRECISION :: ltau_w_g
118 !---------------------------------------------------------------------//
119 
120 ! Set reference phase to gas
121  m = 0
122 
123  IF (.NOT.momentum_z_eq(0)) RETURN
124 
125 !$omp parallel do default(shared) &
126 !$omp private(I, J, K, IJK, IJKT, IJKM, IJKP, IMJK, IPJK, IJMK, &
127 !$omp IMJKP, IJPK, IJMKP, IJKTE, IJKTW, IM, IJKW, IJKE, &
128 !$omp EPGA, epgaj, PGT, SDP, ROPGA, ROGA, V0, ISV, MUGA, &
129 !$omp vpm, Vmt, Vbf, F_vir, Ghd_drag, avgRop, HYS_drag, &
130 !$omp avgdrag, MM, L, VXZA, VCOA, VCOB, CTE, CTW, UGT, &
131 !$omp cpe, cpw, MUOX, ltau_w_g)
132  DO ijk = ijkstart3, ijkend3
133  i = i_of(ijk)
134  j = j_of(ijk)
135  k = k_of(ijk)
136  ijkt = top_of(ijk)
137  ijkm = km_of(ijk)
138  ijkp = kp_of(ijk)
139  imjk = im_of(ijk)
140  ipjk = ip_of(ijk)
141  imjkp = kp_of(imjk)
142  ijmk = jm_of(ijk)
143  ijpk = jp_of(ijk)
144  ijmkp = kp_of(ijmk)
145 
146  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
147 ! if jackson then avg ep_g otherwise 1
148  epgaj = avg_z(epg_jfac(ijk),epg_jfac(ijkt),k)
149 
150 ! Impermeable internal surface
151  IF (ip_at_t(ijk)) THEN
152  a_m(ijk,east,m) = zero
153  a_m(ijk,west,m) = zero
154  a_m(ijk,north,m) = zero
155  a_m(ijk,south,m) = zero
156  a_m(ijk,top,m) = zero
157  a_m(ijk,bottom,m) = zero
158  a_m(ijk,0,m) = -one
159  b_m(ijk,m) = zero
160 
161 ! Dilute flow
162  ELSEIF (epga <= dil_ep_s) THEN
163  a_m(ijk,east,m) = zero
164  a_m(ijk,west,m) = zero
165  a_m(ijk,north,m) = zero
166  a_m(ijk,south,m) = zero
167  a_m(ijk,top,m) = zero
168  a_m(ijk,bottom,m) = zero
169  a_m(ijk,0,m) = -one
170  b_m(ijk,m) = zero
171 ! set velocity equal to that of bottom or top cell if solids are
172 ! present in those cells else set velocity equal to known value
173  IF (ep_g(bottom_of(ijk)) > dil_ep_s) THEN
174  a_m(ijk,bottom,m) = one
175  ELSE IF (ep_g(top_of(ijk)) > dil_ep_s) THEN
176  a_m(ijk,top,m) = one
177  ELSE
178  b_m(ijk,m) = -w_g(ijk)
179  ENDIF
180 
181 ! Cartesian grid implementation
182  ELSEIF (blocked_w_cell_at(ijk)) THEN
183  a_m(ijk,east,m) = zero
184  a_m(ijk,west,m) = zero
185  a_m(ijk,north,m) = zero
186  a_m(ijk,south,m) = zero
187  a_m(ijk,top,m) = zero
188  a_m(ijk,bottom,m) = zero
189  a_m(ijk,0,m) = -one
190  b_m(ijk,m) = zero
191 
192 ! Normal case
193  ELSE
194 
195 ! Surface forces
196 
197 ! Pressure term
198  pgt = p_g(ijkt)
199  IF (cyclic_z_pd) THEN
200  IF (kmap(k_of(ijk)).EQ.kmax1) pgt = p_g(ijkt) - delp_z
201  ENDIF
202  IF (model_b) THEN
203  IF(.NOT.cut_w_treatment_at(ijk)) THEN
204  sdp = -p_scale*(pgt - p_g(ijk))*axy(ijk)
205  ELSE
206  sdp = -p_scale*(pgt * a_wpg_t(ijk) - p_g(ijk) * a_wpg_b(ijk) )
207  ENDIF
208  ELSE
209  IF(.NOT.cut_w_treatment_at(ijk)) THEN
210  sdp = -p_scale*epga*(pgt - p_g(ijk))*axy(ijk)
211  ELSE
212  sdp = -p_scale*epga*(pgt * a_wpg_t(ijk) - p_g(ijk) * a_wpg_b(ijk) )
213  ENDIF
214  ENDIF
215 
216  IF(.NOT.cut_w_treatment_at(ijk)) THEN
217 ! Volumetric forces
218  ropga = avg_z(rop_g(ijk),rop_g(ijkt),k)
219  roga = avg_z(ro_g(ijk),ro_g(ijkt),k)
220 
221 ! Previous time step
222  v0 = avg_z(rop_go(ijk),rop_go(ijkt),k)*odt
223 
224 ! Added mass implicit transient term {Cv eps rop_g dW/dt}
225  IF(added_mass) THEN
226  rop_ma = avg_z(rop_g(ijk)*ep_s(ijk,m_am),&
227  rop_g(ijkt)*ep_s(ijkt,m_am),k)
228  v0 = v0 + cv * rop_ma * odt
229  ENDIF
230  ELSE
231 ! Volumetric forces
232  ropga = (vol(ijk)*rop_g(ijk) + vol(ijkt)*rop_g(ijkt))/&
233  (vol(ijk) + vol(ijkt))
234  roga = (vol(ijk)*ro_g(ijk) + vol(ijkt)*ro_g(ijkt))/&
235  (vol(ijk) + vol(ijkt))
236 ! Previous time step
237  v0 = (vol(ijk)*rop_go(ijk) + vol(ijkt)*rop_go(ijkt))*&
238  odt/(vol(ijk) + vol(ijkt))
239 ! Added mass implicit transient term {Cv eps rop_g dW/dt}
240  IF(added_mass) THEN
241  rop_ma = (vol(ijk)*rop_g(ijk)*ep_s(ijk,m_am) + &
242  vol(ijkt)*rop_g(ijkt)*ep_s(ijkt,m_am))/&
243  (vol(ijk) + vol(ijkt))
244  v0 = v0 + cv * rop_ma * odt
245  ENDIF
246  ENDIF
247 
248 ! VIRTUAL MASS SECTION (explicit terms)
249 ! adding transient term dWs/dt to virtual mass term
250  f_vir = zero
251  IF(added_mass.AND.(.NOT.cut_w_treatment_at(ijk))) THEN
252  f_vir = ( (w_s(ijk,m_am) - w_so(ijk,m_am)) )*odt*vol_w(ijk)
253 
254 ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
255  wsb = avg_z_t(w_s(ijkm,m_am),w_s(ijk,m_am))
256  wst = avg_z_t(w_s(ijk,m_am),w_s(ijkp,m_am))
257  u_se = avg_z(u_s(ijk,m_am),u_s(ijkp,m_am),k)
258  usw = avg_z(u_s(imjk,m_am),u_s(imjkp,m_am),k)
259  ust = avg_x_e(usw,u_se,ip1(i))
260  wse = avg_x(w_s(ijk,m_am),w_s(ipjk,m_am),ip1(i))
261  wsw = avg_x(w_s(imjk,m_am),w_s(ijk,m_am),i)
262  vsb = avg_y_n(v_s(ijmk,m_am),v_s(ijk,m_am))
263  vst = avg_y_n(v_s(ijmkp,m_am),v_s(ijkp,m_am))
264  wss = avg_y(w_s(ijmk,m_am),w_s(ijk,m_am),jm1(j))
265  wsn = avg_y(w_s(ijk,m_am),w_s(ijpk,m_am),j)
266 
267 ! adding convective terms (U dW/dx + V dW/dy + W dW/dz) to virtual mass.
268  f_vir = f_vir + w_s(ijk,m_am)*ox(i)* &
269  (wst - wsb)*axy(ijk) + ust*(wse - wsw)*ayz(ijk) + &
270  avg_z(vsb,vst,k)*(wsn - wss)*axz(ijk)
271 
272 ! Coriolis force
273  IF (cylindrical) f_vir = f_vir + &
274  ust*w_s(ijk,m_am)*ox(i)
275  f_vir = f_vir * cv * rop_ma
276  ENDIF
277 
278 ! pressure drop through porous media
279  IF (sip_at_t(ijk)) THEN
280  isv = is_id_at_t(ijk)
281  muga = avg_z(mu_g(ijk),mu_g(ijkt),k)
282  vpm = muga/is_pc(isv,1)
283  IF (is_pc(isv,2) /= zero) vpm = vpm + &
284  half*is_pc(isv,2)*ropga*abs(w_g(ijk))
285  ELSE
286  vpm = zero
287  ENDIF
288 
289 ! Interphase mass transfer
290  IF(.NOT.cut_w_treatment_at(ijk)) THEN
291  vmt = avg_z(sum_r_g(ijk),sum_r_g(ijkt),k)
292  ELSE
293  vmt = (vol(ijk)*sum_r_g(ijk) + vol(ijkt)*sum_r_g(ijkt))/&
294  (vol(ijk) + vol(ijkt))
295  ENDIF
296 
297 ! Body force
298  IF (model_b) THEN
299  vbf = roga*bfz_g(ijk)
300  ELSE
301  vbf = ropga*bfz_g(ijk)
302  ENDIF
303 
304 ! Additional force for GHD from darg force sum(beta_ig * Joi/rhop_i)
305  ghd_drag = zero
306  IF (kt_type_enum .EQ. ghd_2007) THEN
307  DO l = 1,smax
308  avgrop = avg_z(rop_s(ijk,l),rop_s(ijkt,l),k)
309  if(avgrop > zero) ghd_drag = ghd_drag +&
310  avg_z(f_gs(ijk,l),f_gs(ijkt,l),k) * joiz(ijk,l) / avgrop
311  ENDDO
312  ENDIF
313 
314 ! Additional force for HYS drag force, do not use with mixture GHD theory
315  avgdrag = zero
316  hys_drag = zero
317  IF (drag_type_enum .EQ. hys .AND. kt_type_enum .NE. ghd_2007) THEN
318  DO mm=1,mmax
319  DO l = 1,mmax
320  IF (l /= mm) THEN
321  avgdrag = avg_z(beta_ij(ijk,mm,l),beta_ij(ijkt,mm,l),k)
322  hys_drag = hys_drag + avgdrag * (w_g(ijk) - w_s(ijk,l))
323  ENDIF
324  ENDDO
325  ENDDO
326  ENDIF
327 
328 ! Special terms for cylindrical coordinates
329  vcoa = zero
330  vcob = zero
331  vxza = zero
332  cte = zero
333  ctw = zero
334  cpe = zero
335  cpw = zero
336  IF (cylindrical) THEN
337 ! Coriolis force
338  imjk = im_of(ijk)
339  ijkp = kp_of(ijk)
340  imjkp = kp_of(imjk)
341  ugt = avg_z(half*(u_g(ijk)+u_g(imjk)),&
342  half*(u_g(ijkp)+u_g(imjkp)),k)
343  IF (ugt > zero) THEN
344  vcoa = ropga*ugt*ox(i)
345  vcob = zero
346  IF(added_mass) vcoa = vcoa + cv*rop_ma*ugt*ox(i)
347  ELSE
348  vcoa = zero
349  vcob = -ropga*ugt*w_g(ijk)*ox(i)
350  IF(added_mass) vcob = vcob - cv*rop_ma*ugt*w_g(ijk)*ox(i)
351  ENDIF
352 
353 ! if ishii, then viscosity is multiplied by void fraction otherwise by 1
354 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
355 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
356 ! 1/x d/dx (x.mu.(-w/x)) xdxdydz =>
357 ! delta (mu/x.(-w))Ayz |E-W : at (i+1/2 - i-1/2, j, k+1/2)
358  ijke = east_of(ijk)
359  ijkw = west_of(ijk)
360  ijkte = top_of(ijke)
361  ijktw = top_of(ijkw)
362  im = im1(i)
363  ipjk = ip_of(ijk)
364  cte = half*avg_z_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
365  avg_x_h(epmu_gt(ijkt),epmu_gt(ijkte),i),k)*&
366  ox_e(i)*ayz_w(ijk)
367  ctw = half*avg_z_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
368  avg_x_h(epmu_gt(ijktw),epmu_gt(ijkt),im),k)*&
369  dy(j)*(half*(dz(k)+dz(kp1(k))))
370 ! DY(J)*HALF(DZ(k)+DZ(kp)) = oX_E(IM)*AYZ_W(IMJK), but avoids singularity
371 
372 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
373 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
374 ! mu/x dw/dx xdxdydz =>
375 ! delta (mu/x.(dw/dx))Vp |p : at (i, j, k+1/2)
376  muox = avg_z(epmu_gt(ijk),epmu_gt(ijkt),k)*ox(i)
377  cpe = muox*half*vol_w(ijk)*odx_e(i)
378  cpw = muox*half*vol_w(ijk)*odx_e(im)
379 
380 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
381 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
382 ! 1/x d/dx (x.mu.(-w/x)) xdxdydz =>
383 ! delta (mu/x.(-w/x))Vp |p : at (i, j, k+1/2)
384  vxza = muox*ox(i)
385 
386  cte = epgaj*cte
387  ctw = epgaj*ctw
388  cpe = epgaj*cpe
389  cpw = epgaj*cpw
390  vxza = epgaj*vxza
391  ENDIF
392 
393 ! if jackson, implement jackson form of governing equations (ep_g dot
394 ! del tau_g): multiply by void fraction otherwise by 1
395  ltau_w_g = epgaj*tau_w_g(ijk)
396 
397 ! Collect the terms
398  a_m(ijk,east,m) = a_m(ijk,east,m) + cpe
399  a_m(ijk,west,m) = a_m(ijk,west,m) - cpw
400 
401  a_m(ijk,0,m) = -(a_m(ijk,east,m)+a_m(ijk,west,m)+&
402  a_m(ijk,north,m)+a_m(ijk,south,m)+a_m(ijk,top,m)+a_m(ijk,bottom,m)+&
403  (v0+vpm+zmax(vmt)+vcoa + vxza)*vol_w(ijk) + cte - ctw)
404 
405  a_m(ijk,east,m) = a_m(ijk,east,m) - cte
406  a_m(ijk,west,m) = a_m(ijk,west,m) + ctw
407 
408  b_m(ijk,m) = b_m(ijk,m) - ( sdp + ltau_w_g + f_vir + &
409  ( (v0+zmax((-vmt)))*w_go(ijk) + vbf + vcob + &
410  ghd_drag+hys_drag)*vol_w(ijk) )
411 
412 ! MMS Source term.
413  IF(use_mms) b_m(ijk,m) = &
414  b_m(ijk,m) - mms_w_g_src(ijk)*vol_w(ijk)
415  ENDIF ! end branching on cell type (ip/dilute/block/else branches)
416  ENDDO ! end do loop over ijk
417 !$omp end parallel do
418 
419 ! modifications for cartesian grid implementation
420  IF(cartesian_grid) CALL cg_source_w_g(a_m, b_m)
421 ! modifications for bc
422  CALL source_w_g_bc (a_m, b_m)
423 ! modifications for cartesian grid implementation
424  IF(cartesian_grid) CALL cg_source_w_g_bc(a_m, b_m)
425 
426  RETURN
427  END SUBROUTINE source_w_g
428 
429 
430 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
431 ! C
432 ! Subroutine: SOURCE_W_g_BC C
433 ! Purpose: Determine source terms for W_g momentum eq. The terms C
434 ! appear in the center coefficient and RHS vector. The center C
435 ! coefficient and source vector are negative. The off-diagonal C
436 ! coefficients are positive. C
437 ! The drag terms are excluded from the source at this stage C
438 ! C
439 ! Author: M. Syamlal Date: 17-JUN-96 C
440 ! Reviewer: Date: C
441 ! C
442 ! C
443 ! C
444 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
445 
446  SUBROUTINE source_w_g_bc(A_M, B_M)
447 
448 !-----------------------------------------------
449 ! Modules
450 !-----------------------------------------------
451  USE param
452  USE param1
453  USE parallel
454  USE scales
455  USE constant
456  USE physprop
457  USE fldvar
458  USE visc_g
459  USE rxns
460  USE run
461  USE toleranc
462  USE geometry
463  USE indices
464  USE is
465  USE tau_g
466  USE bc
467  USE output
468  USE compar
469  USE fun_avg
470  USE functions
471 
472  IMPLICIT NONE
473 !-----------------------------------------------
474 ! Dummy arguments
475 !-----------------------------------------------
476 ! Septadiagonal matrix A_m
477  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
478 ! Vector b_m
479  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
480 !-----------------------------------------------
481 ! Local parameters
482 !-----------------------------------------------
483 ! C_mu is constant in turbulent viscosity
484  DOUBLE PRECISION, PARAMETER :: C_mu = 0.09d0
485 ! Kappa is Von Karmen constant
486  DOUBLE PRECISION, PARAMETER :: Kappa = 0.42d0
487 !-----------------------------------------------
488 ! Local variables
489 !-----------------------------------------------
490 ! Boundary condition
491  INTEGER :: L
492 ! Indices
493  INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
494  IM, JM, IJKB, IJKM, IJKP
495 ! Phase index
496  INTEGER :: M
497 ! Turbulent shear at walls
498  DOUBLE PRECISION W_F_Slip
499 
500 !-----------------------------------------------
501 
502 ! Set reference phase to gas
503  m = 0
504 
505 
506 ! Set the default boundary conditions
507 ! The NS default setting is the where bc_type='dummy' or any default
508 ! (i.e., bc_type=undefined) wall boundary regions are handled. Note that
509 ! the top and bottom xy planes do not have to be explicitly addressed for
510 ! the w-momentum equation. In this direction the velocities are defined
511 ! at the wall (due staggered grid). They are defined as zero for a
512 ! no penetration condition (see zero_norm_vel subroutine and code under
513 ! the ip_at_t branch in the above source routine).
514 ! ---------------------------------------------------------------->>>
515 
516 ! south xz plane
517  j1 = 1
518  DO k1 = kmin3,kmax3
519  DO i1 = imin3,imax3
520  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
521  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
522  ijk = funijk(i1,j1,k1)
523  IF (ns_wall_at(ijk)) THEN
524  a_m(ijk,east,m) = zero
525  a_m(ijk,west,m) = zero
526  a_m(ijk,north,m) = -one
527  a_m(ijk,south,m) = zero
528  a_m(ijk,top,m) = zero
529  a_m(ijk,bottom,m) = zero
530  a_m(ijk,0,m) = -one
531  b_m(ijk,m) = zero
532  ELSEIF (fs_wall_at(ijk)) THEN
533  a_m(ijk,east,m) = zero
534  a_m(ijk,west,m) = zero
535  a_m(ijk,north,m) = one
536  a_m(ijk,south,m) = zero
537  a_m(ijk,top,m) = zero
538  a_m(ijk,bottom,m) = zero
539  a_m(ijk,0,m) = -one
540  b_m(ijk,m) = zero
541  ENDIF
542  ENDDO
543  ENDDO
544 
545 ! north xz plane
546  j1 = jmax2
547  DO k1 = kmin3, kmax3
548  DO i1 = imin3, imax3
549  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
550  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
551  ijk = funijk(i1,j1,k1)
552  IF (ns_wall_at(ijk)) THEN
553 ! Setting the wall velocity to zero (set the boundary cell value equal
554 ! and oppostive to the adjacent fluid cell value)
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 ! Setting the wall velocity equal to the adjacent fluid velocity (set
565 ! the boundary cell value equal to adjacent fluid cell value)
566  a_m(ijk,east,m) = zero
567  a_m(ijk,west,m) = zero
568  a_m(ijk,north,m) = zero
569  a_m(ijk,south,m) = one
570  a_m(ijk,top,m) = zero
571  a_m(ijk,bottom,m) = zero
572  a_m(ijk,0,m) = -one
573  b_m(ijk,m) = zero
574  ENDIF
575  ENDDO
576  ENDDO
577 
578 ! west yz plane
579  i1 = 1
580  DO k1 = kmin3, kmax3
581  DO j1 = jmin3, jmax3
582  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
583  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
584  ijk = funijk(i1,j1,k1)
585  IF (ns_wall_at(ijk)) THEN
586  a_m(ijk,east,m) = -one
587  a_m(ijk,west,m) = zero
588  a_m(ijk,north,m) = zero
589  a_m(ijk,south,m) = zero
590  a_m(ijk,top,m) = zero
591  a_m(ijk,bottom,m) = zero
592  a_m(ijk,0,m) = -one
593  b_m(ijk,m) = zero
594  ELSEIF (fs_wall_at(ijk)) THEN
595  a_m(ijk,east,m) = one
596  a_m(ijk,west,m) = zero
597  a_m(ijk,north,m) = zero
598  a_m(ijk,south,m) = zero
599  a_m(ijk,top,m) = zero
600  a_m(ijk,bottom,m) = zero
601  a_m(ijk,0,m) = -one
602  b_m(ijk,m) = zero
603  ENDIF
604  ENDDO
605  ENDDO
606 
607 ! east yz plane
608  i1 = imax2
609  DO k1 = kmin3,kmax3
610  DO j1 = jmin3,jmax3
611  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
612  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
613  ijk = funijk(i1,j1,k1)
614  IF (ns_wall_at(ijk)) THEN
615  a_m(ijk,east,m) = zero
616  a_m(ijk,west,m) = -one
617  a_m(ijk,north,m) = zero
618  a_m(ijk,south,m) = zero
619  a_m(ijk,top,m) = zero
620  a_m(ijk,bottom,m) = zero
621  a_m(ijk,0,m) = -one
622  b_m(ijk,m) = zero
623  ELSEIF (fs_wall_at(ijk)) THEN
624  a_m(ijk,east,m) = zero
625  a_m(ijk,west,m) = one
626  a_m(ijk,north,m) = zero
627  a_m(ijk,south,m) = zero
628  a_m(ijk,top,m) = zero
629  a_m(ijk,bottom,m) = zero
630  a_m(ijk,0,m) = -one
631  b_m(ijk,m) = zero
632  ENDIF
633  ENDDO
634  ENDDO
635 ! End setting the default boundary conditions
636 ! ----------------------------------------------------------------<<<
637 
638 ! Setting user specified boundary conditions
639 
640  DO l = 1, dimension_bc
641  IF (bc_defined(l)) THEN
642 
643 ! Setting wall boundary conditions
644 ! ---------------------------------------------------------------->>>
645  IF (bc_type_enum(l) == no_slip_wall .AND. .NOT. k_epsilon) THEN
646  i1 = bc_i_w(l)
647  i2 = bc_i_e(l)
648  j1 = bc_j_s(l)
649  j2 = bc_j_n(l)
650  k1 = bc_k_b(l)
651  k2 = bc_k_t(l)
652  DO k = k1, k2
653  DO j = j1, j2
654  DO i = i1, i2
655  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
656  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
657  ijk = funijk(i,j,k)
658  IF (.NOT.wall_at(ijk)) cycle !skip redefined cells
659  a_m(ijk,east,m) = zero
660  a_m(ijk,west,m) = zero
661  a_m(ijk,north,m) = zero
662  a_m(ijk,south,m) = zero
663  a_m(ijk,top,m) = zero
664  a_m(ijk,bottom,m) = zero
665  a_m(ijk,0,m) = -one
666  b_m(ijk,m) = zero
667  IF (fluid_at(east_of(ijk))) THEN
668  a_m(ijk,east,m) = -one
669  ELSEIF (fluid_at(west_of(ijk))) THEN
670  a_m(ijk,west,m) = -one
671  ELSEIF (fluid_at(north_of(ijk))) THEN
672  a_m(ijk,north,m) = -one
673  ELSEIF (fluid_at(south_of(ijk))) THEN
674  a_m(ijk,south,m) = -one
675  ENDIF
676  ENDDO
677  ENDDO
678  ENDDO
679 
680  ELSEIF (bc_type_enum(l) == free_slip_wall .AND. .NOT. k_epsilon) THEN
681  i1 = bc_i_w(l)
682  i2 = bc_i_e(l)
683  j1 = bc_j_s(l)
684  j2 = bc_j_n(l)
685  k1 = bc_k_b(l)
686  k2 = bc_k_t(l)
687  DO k = k1, k2
688  DO j = j1, j2
689  DO i = i1, i2
690  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
691  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
692  ijk = funijk(i,j,k)
693  IF (.NOT.wall_at(ijk)) cycle !skip redefined cells
694  a_m(ijk,east,m) = zero
695  a_m(ijk,west,m) = zero
696  a_m(ijk,north,m) = zero
697  a_m(ijk,south,m) = zero
698  a_m(ijk,top,m) = zero
699  a_m(ijk,bottom,m) = zero
700  a_m(ijk,0,m) = -one
701  b_m(ijk,m) = zero
702  IF (fluid_at(east_of(ijk))) THEN
703  a_m(ijk,east,m) = one
704  ELSEIF (fluid_at(west_of(ijk))) THEN
705  a_m(ijk,west,m) = one
706  ELSEIF (fluid_at(north_of(ijk))) THEN
707  a_m(ijk,north,m) = one
708  ELSEIF (fluid_at(south_of(ijk))) THEN
709  a_m(ijk,south,m) = one
710  ENDIF
711  ENDDO
712  ENDDO
713  ENDDO
714 
715  ELSEIF (bc_type_enum(l) == par_slip_wall .AND. .NOT. k_epsilon) THEN
716  i1 = bc_i_w(l)
717  i2 = bc_i_e(l)
718  j1 = bc_j_s(l)
719  j2 = bc_j_n(l)
720  k1 = bc_k_b(l)
721  k2 = bc_k_t(l)
722  DO k = k1, k2
723  DO j = j1, j2
724  DO i = i1, i2
725  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
726  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
727  ijk = funijk(i,j,k)
728  IF (.NOT.wall_at(ijk)) cycle ! skip redefined cells
729  im = im1(i)
730  jm = jm1(j)
731  a_m(ijk,east,m) = zero
732  a_m(ijk,west,m) = zero
733  a_m(ijk,north,m) = zero
734  a_m(ijk,south,m) = zero
735  a_m(ijk,top,m) = zero
736  a_m(ijk,bottom,m) = zero
737  a_m(ijk,0,m) = -one
738  b_m(ijk,m) = zero
739  IF (fluid_at(east_of(ijk))) THEN
740  IF (bc_hw_g(l) == undefined) THEN
741  a_m(ijk,east,m) = -half
742  a_m(ijk,0,m) = -half
743  b_m(ijk,m) = -bc_ww_g(l)
744  ELSE
745  IF (cylindrical) THEN
746  a_m(ijk,0,m) = -( half*(bc_hw_g(l)-&
747  ox_e(i)) + odx_e(i) )
748  a_m(ijk,east,m) = -(half*(bc_hw_g(l)-&
749  ox_e(i)) - odx_e(i))
750  b_m(ijk,m) = -bc_hw_g(l)*bc_ww_g(l)
751  ELSE
752  a_m(ijk,0,m) = -(half*bc_hw_g(l)+odx_e(i))
753  a_m(ijk,east,m) = -(half*bc_hw_g(l)-odx_e(i))
754  b_m(ijk,m) = -bc_hw_g(l)*bc_ww_g(l)
755  ENDIF
756  ENDIF
757  ELSEIF (fluid_at(west_of(ijk))) THEN
758  IF (bc_hw_g(l) == undefined) THEN
759  a_m(ijk,west,m) = -half
760  a_m(ijk,0,m) = -half
761  b_m(ijk,m) = -bc_ww_g(l)
762  ELSE
763  IF (cylindrical) THEN
764  a_m(ijk,west,m) = -(half*(bc_hw_g(l)-&
765  ox_e(im)) - odx_e(im))
766  a_m(ijk,0,m) = -(half*(bc_hw_g(l)-&
767  ox_e(im)) + odx_e(im))
768  b_m(ijk,m) = -bc_hw_g(l)*bc_ww_g(l)
769  ELSE
770  a_m(ijk,west,m) = -(half*bc_hw_g(l)-odx_e(im))
771  a_m(ijk,0,m) = -(half*bc_hw_g(l)+odx_e(im))
772  b_m(ijk,m) = -bc_hw_g(l)*bc_ww_g(l)
773  ENDIF
774  ENDIF
775  ELSEIF (fluid_at(north_of(ijk))) THEN
776  IF (bc_hw_g(l) == undefined) THEN
777  a_m(ijk,north,m) = -half
778  a_m(ijk,0,m) = -half
779  b_m(ijk,m) = -bc_ww_g(l)
780  ELSE
781  a_m(ijk,0,m) = -(half*bc_hw_g(l)+ody_n(j))
782  a_m(ijk,north,m) = -(half*bc_hw_g(l)-ody_n(j))
783  b_m(ijk,m) = -bc_hw_g(l)*bc_ww_g(l)
784  ENDIF
785  ELSEIF (fluid_at(south_of(ijk))) THEN
786  IF (bc_hw_g(l) == undefined) THEN
787  a_m(ijk,south,m) = -half
788  a_m(ijk,0,m) = -half
789  b_m(ijk,m) = -bc_ww_g(l)
790  ELSE
791  a_m(ijk,south,m) = -(half*bc_hw_g(l)-ody_n(jm))
792  a_m(ijk,0,m) = -(half*bc_hw_g(l)+ody_n(jm))
793  b_m(ijk,m) = -bc_hw_g(l)*bc_ww_g(l)
794  ENDIF
795  ENDIF
796  ENDDO
797  ENDDO
798  ENDDO
799 
800 ! Setting wall boundary conditions when K_EPSILON
801 ! wall functions for V-momentum are specify in this section of the code
802  ELSEIF (bc_type_enum(l) == par_slip_wall .OR. &
803  bc_type_enum(l) == no_slip_wall .OR. &
804  bc_type_enum(l) == free_slip_wall .AND. &
805  k_epsilon )THEN
806  i1 = bc_i_w(l)
807  i2 = bc_i_e(l)
808  j1 = bc_j_s(l)
809  j2 = bc_j_n(l)
810  k1 = bc_k_b(l)
811  k2 = bc_k_t(l)
812  DO k = k1, k2
813  DO j = j1, j2
814  DO i = i1, i2
815  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
816  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
817  ijk = funijk(i,j,k)
818  IF (.NOT.wall_at(ijk)) cycle !skip redefined cells
819  im = im1(i)
820  jm = jm1(j)
821  a_m(ijk,east,m) = zero
822  a_m(ijk,west,m) = zero
823  a_m(ijk,north,m) = zero
824  a_m(ijk,south,m) = zero
825  a_m(ijk,top,m) = zero
826  a_m(ijk,bottom,m) = zero
827  a_m(ijk,0,m) = -one
828  b_m(ijk,m) = zero
829  IF (fluid_at(east_of(ijk))) THEN
830  IF (cylindrical) THEN
831  w_f_slip = ( one/&
832  (odx_e(i)+half*ox_e(i)) )* &
833  ( odx_e(i) - ox_e(i) - &
834  ro_g(east_of(ijk))*c_mu**0.25* &
835  sqrt(k_turb_g((east_of(ijk))))/ &
836  mu_gt(east_of(ijk))*kappa/ &
837  log(9.81d0/odx_e(i)/(2.d0)* &
838  ro_g(east_of(ijk))*c_mu**0.25* &
839  sqrt(k_turb_g((east_of(ijk))))/ &
840  mu_g(east_of(ijk))) )
841  ELSE
842  CALL wall_function(ijk,east_of(ijk),&
843  odx_e(i),w_f_slip)
844  ENDIF
845  a_m(ijk,east,m) = w_f_slip
846  a_m(ijk,0,m) = -one
847  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_ww_g(l)
848  ELSEIF (fluid_at(west_of(ijk))) THEN
849  IF (cylindrical) THEN
850  w_f_slip = (one/&
851  (one*odx_e(im) + half*ox_e(im)))* &
852  ( one*odx_e(im) - ox_e(im) - &
853  ro_g(west_of(ijk))*c_mu**0.25* &
854  sqrt(k_turb_g((west_of(ijk))))/ &
855  mu_gt(west_of(ijk))*kappa/ &
856  log(9.81d0/odx_e(im)/(2.d0)* &
857  ro_g(west_of(ijk))*c_mu**0.25* &
858  sqrt(k_turb_g((west_of(ijk))))/ &
859  mu_g(west_of(ijk))) )
860  ELSE
861  CALL wall_function(ijk,west_of(ijk),&
862  odx_e(im),w_f_slip)
863  ENDIF
864  a_m(ijk,west,m) = w_f_slip
865  a_m(ijk,0,m) = -one
866  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_ww_g(l)
867  ELSEIF (fluid_at(north_of(ijk))) THEN
868  CALL wall_function(ijk,north_of(ijk),ody_n(j),w_f_slip)
869  a_m(ijk,north,m) = w_f_slip
870  a_m(ijk,0,m) = -one
871  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_ww_g(l)
872  ELSEIF (fluid_at(south_of(ijk))) THEN
873  CALL wall_function(ijk,south_of(ijk),ody_n(jm),w_f_slip)
874  a_m(ijk,south,m) = w_f_slip
875  a_m(ijk,0,m) = -one
876  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_ww_g(l)
877  ENDIF
878  ENDDO
879  ENDDO
880  ENDDO
881 
882 ! end setting of wall boundary conditions
883 ! ----------------------------------------------------------------<<<
884 
885 ! Setting p_inflow or p_outflow flow boundary conditions
886 ! ---------------------------------------------------------------->>>
887  ELSEIF (bc_type_enum(l)==p_inflow .OR. bc_type_enum(l)==p_outflow) THEN
888  IF (bc_plane(l) == 'B') THEN
889 ! if the fluid cell is on the bottom side of the outflow/inflow boundary
890 ! then set the velocity in the boundary cell equal to the velocity of
891 ! the adjacent fluid cell
892  i1 = bc_i_w(l)
893  i2 = bc_i_e(l)
894  j1 = bc_j_s(l)
895  j2 = bc_j_n(l)
896  k1 = bc_k_b(l)
897  k2 = bc_k_t(l)
898  DO k = k1, k2
899  DO j = j1, j2
900  DO i = i1, i2
901  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
902  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
903  ijk = funijk(i,j,k)
904  a_m(ijk,east,m) = zero
905  a_m(ijk,west,m) = zero
906  a_m(ijk,north,m) = zero
907  a_m(ijk,south,m) = zero
908  a_m(ijk,top,m) = zero
909  a_m(ijk,bottom,m) = one
910  a_m(ijk,0,m) = -one
911  b_m(ijk,m) = zero
912  ENDDO
913  ENDDO
914  ENDDO
915  ENDIF
916 ! end setting of p_inflow or p_otuflow flow boundary conditions
917 ! ----------------------------------------------------------------<<<
918 
919 ! Setting outflow flow boundary conditions
920 ! ---------------------------------------------------------------->>>
921  ELSEIF (bc_type_enum(l) == outflow) THEN
922  IF (bc_plane(l) == 'B') THEN
923  i1 = bc_i_w(l)
924  i2 = bc_i_e(l)
925  j1 = bc_j_s(l)
926  j2 = bc_j_n(l)
927  k1 = bc_k_b(l)
928  k2 = bc_k_t(l)
929  DO k = k1, k2
930  DO j = j1, j2
931  DO i = i1, i2
932  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
933  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
934  ijk = funijk(i,j,k)
935  a_m(ijk,east,m) = zero
936  a_m(ijk,west,m) = zero
937  a_m(ijk,north,m) = zero
938  a_m(ijk,south,m) = zero
939  a_m(ijk,top,m) = zero
940  a_m(ijk,bottom,m) = one
941  a_m(ijk,0,m) = -one
942  b_m(ijk,m) = zero
943  ijkm = km_of(ijk)
944  a_m(ijkm,east,m) = zero
945  a_m(ijkm,west,m) = zero
946  a_m(ijkm,north,m) = zero
947  a_m(ijkm,south,m) = zero
948  a_m(ijkm,top,m) = zero
949  a_m(ijkm,bottom,m) = one
950  a_m(ijkm,0,m) = -one
951  b_m(ijkm,m) = zero
952  ENDDO
953  ENDDO
954  ENDDO
955  ELSEIF (bc_plane(l) == 'T') THEN
956  i1 = bc_i_w(l)
957  i2 = bc_i_e(l)
958  j1 = bc_j_s(l)
959  j2 = bc_j_n(l)
960  k1 = bc_k_b(l)
961  k2 = bc_k_t(l)
962  DO k = k1, k2
963  DO j = j1, j2
964  DO i = i1, i2
965  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
966  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
967  ijk = funijk(i,j,k)
968  ijkp = kp_of(ijk)
969  a_m(ijkp,east,m) = zero
970  a_m(ijkp,west,m) = zero
971  a_m(ijkp,north,m) = zero
972  a_m(ijkp,south,m) = zero
973  a_m(ijkp,top,m) = one
974  a_m(ijkp,bottom,m) = zero
975  a_m(ijkp,0,m) = -one
976  b_m(ijkp,m) = zero
977  ENDDO
978  ENDDO
979  ENDDO
980  ENDIF
981 ! end setting of outflow flow boundary conditions
982 ! ----------------------------------------------------------------<<<
983 
984 ! Setting bc that are defined but not nsw, fsw, psw, p_inflow,
985 ! p_outflow, or outflow (at this time, this section addresses
986 ! mass_inflow and mass_outflow type boundaries)
987 ! ---------------------------------------------------------------->>>
988  ELSE
989  i1 = bc_i_w(l)
990  i2 = bc_i_e(l)
991  j1 = bc_j_s(l)
992  j2 = bc_j_n(l)
993  k1 = bc_k_b(l)
994  k2 = bc_k_t(l)
995  DO k = k1, k2
996  DO j = j1, j2
997  DO i = i1, i2
998  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
999  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
1000  ijk = funijk(i,j,k)
1001 ! setting the velocity in the boundary cell equal to what is known
1002  a_m(ijk,east,m) = zero
1003  a_m(ijk,west,m) = zero
1004  a_m(ijk,north,m) = zero
1005  a_m(ijk,south,m) = zero
1006  a_m(ijk,top,m) = zero
1007  a_m(ijk,bottom,m) = zero
1008  a_m(ijk,0,m) = -one
1009  b_m(ijk,m) = -w_g(ijk)
1010  IF (bc_plane(l) == 'B') THEN
1011 ! if the fluid cell is on the bottom side of the outflow/inflow boundary
1012 ! then set the velocity in the adjacent fluid cell equal to what is
1013 ! known in that cell
1014  ijkb = bottom_of(ijk)
1015  a_m(ijkb,east,m) = zero
1016  a_m(ijkb,west,m) = zero
1017  a_m(ijkb,north,m) = zero
1018  a_m(ijkb,south,m) = zero
1019  a_m(ijkb,top,m) = zero
1020  a_m(ijkb,bottom,m) = zero
1021  a_m(ijkb,0,m) = -one
1022  b_m(ijkb,m) = -w_g(ijkb)
1023  ENDIF
1024  ENDDO
1025  ENDDO
1026  ENDDO
1027  ENDIF ! end if/else (bc_type)
1028  ! ns, fs, psw; else
1029  ! p_inflow, p_outflow, or outflow; else
1030 ! end setting of 'else' flow boundary conditions
1031 ! (mass_inflow/mass_outflow)
1032 ! ----------------------------------------------------------------<<<
1033 
1034  ENDIF ! end if (bc_defined)
1035  ENDDO ! end L do loop over dimension_bc
1036 
1037  RETURN
1038  END SUBROUTINE source_w_g_bc
1039 
1040 
1041 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1042 ! C
1043 ! Subroutine: POINT_SOURCE_W_G C
1044 ! Purpose: Adds point sources to the gas phase W-Momentum equation. C
1045 ! C
1046 ! Author: J. Musser Date: 10-JUN-13 C
1047 ! Reviewer: Date: C
1048 ! C
1049 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1050  SUBROUTINE point_source_w_g(A_M, B_M)
1051 
1052 !-----------------------------------------------
1053 ! Modules
1054 !-----------------------------------------------
1055  use compar
1056  use constant
1057  use geometry
1058  use indices
1059  use param
1060  use param1, only: small_number, zero
1061  use physprop
1062  use ps
1063  use run
1064  use functions
1065  IMPLICIT NONE
1066 !-----------------------------------------------
1067 ! Dummy arguments
1068 !-----------------------------------------------
1069 ! Septadiagonal matrix A_m
1070  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
1071 ! Vector b_m
1072  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
1073 !-----------------------------------------------
1074 ! Local variables
1075 !-----------------------------------------------
1076 ! Indices
1077  INTEGER :: IJK, I, J, K
1078  INTEGER :: PSV, M
1079  INTEGER :: lKT, lKB
1080 ! terms of bm expression
1081  DOUBLE PRECISION :: pSource
1082 !-----------------------------------------------
1083 
1084 ! Set reference phase to gas
1085  m = 0
1086 
1087 ! Calculate the mass going into each IJK cell. This is done for each
1088 ! call in case the point source is time dependent.
1089  ps_lp: do psv = 1, dimension_ps
1090  if(.NOT.ps_defined(psv)) cycle ps_lp
1091  if(abs(ps_w_g(psv)) < small_number) cycle ps_lp
1092 
1093  if(ps_w_g(psv) < zero) then
1094  lkb = ps_k_b(psv)-1
1095  lkt = ps_k_t(psv)-1
1096  else
1097  lkb = ps_k_b(psv)
1098  lkt = ps_k_t(psv)
1099  endif
1100 
1101  do k = lkb, lkt
1102  do j = ps_j_s(psv), ps_j_n(psv)
1103  do i = ps_i_w(psv), ps_i_e(psv)
1104 
1105  if(.NOT.is_on_mype_plus2layers(i,j,k)) cycle
1106  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
1107 
1108  ijk = funijk(i,j,k)
1109  if(.NOT.fluid_at(ijk)) cycle
1110 
1111  psource = ps_massflow_g(psv) * (vol(ijk)/ps_volume(psv))
1112 
1113  b_m(ijk,m) = b_m(ijk,m) - psource * &
1114  ps_w_g(psv) * ps_vel_mag_g(psv)
1115 
1116  enddo
1117  enddo
1118  enddo
1119 
1120  enddo ps_lp
1121 
1122  RETURN
1123  END SUBROUTINE point_source_w_g
indices::ip1
integer, dimension(:), allocatable ip1
Definition: indices_mod.f:50
bc::bc_k_b
integer, dimension(dimension_bc) bc_k_b
Definition: bc_mod.f:70
physprop::cv
double precision cv
Definition: physprop_mod.f:65
ps::ps_i_w
integer, dimension(dimension_ps) ps_i_w
Definition: ps_mod.f:40
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
geometry::vol_w
double precision, dimension(:), allocatable vol_w
Definition: geometry_mod.f:242
geometry::imax2
integer imax2
Definition: geometry_mod.f:61
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
cg_source_w_g_bc
subroutine cg_source_w_g_bc(A_M, B_M)
Definition: CG_source_w_g.f:289
ps::ps_vel_mag_g
double precision, dimension(dimension_ps) ps_vel_mag_g
Definition: ps_mod.f:56
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
fldvar::k_turb_g
double precision, dimension(:), allocatable k_turb_g
Definition: fldvar_mod.f:161
geometry::ox_e
double precision, dimension(:), allocatable ox_e
Definition: geometry_mod.f:143
geometry::imax3
integer imax3
Definition: geometry_mod.f:91
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
bodyforce::bfz_g
double precision function bfz_g(IJK)
Definition: bodyforce_mod.f:34
visc_g::mu_gt
double precision, dimension(:), allocatable mu_gt
Definition: visc_g_mod.f:8
cutcell::a_wpg_t
double precision, dimension(:), allocatable a_wpg_t
Definition: cutcell_mod.f:311
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
mms::mms_w_g_src
double precision, dimension(:), allocatable mms_w_g_src
Definition: mms_mod.f:53
bc::bc_i_w
integer, dimension(dimension_bc) bc_i_w
Definition: bc_mod.f:54
run::added_mass
logical added_mass
Definition: run_mod.f:91
bc::delp_z
double precision delp_z
Definition: bc_mod.f:278
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
bc::bc_j_n
integer, dimension(dimension_bc) bc_j_n
Definition: bc_mod.f:66
rxns
Definition: rxns_mod.f:1
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
scales::p_scale
double precision p_scale
Definition: scales_mod.f:13
fldvar::epg_jfac
double precision, dimension(:), allocatable epg_jfac
Definition: fldvar_mod.f:18
cg_source_w_g
subroutine cg_source_w_g(A_M, B_M)
Definition: CG_source_w_g.f:13
drag
Definition: drag_mod.f:11
point_source_w_g
subroutine point_source_w_g(A_M, B_M)
Definition: source_w_g.f:1051
run::momentum_z_eq
logical, dimension(0:dim_m) momentum_z_eq
Definition: run_mod.f:80
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
rxns::sum_r_g
double precision, dimension(:), allocatable sum_r_g
Definition: rxns_mod.f:28
ps::ps_j_n
integer, dimension(dimension_ps) ps_j_n
Definition: ps_mod.f:43
fldvar::p_g
double precision, dimension(:), allocatable p_g
Definition: fldvar_mod.f:26
param::dimension_bc
integer, parameter dimension_bc
Definition: param_mod.f:61
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
param::east
integer east
Definition: param_mod.f:29
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
geometry::ayz
double precision, dimension(:), allocatable ayz
Definition: geometry_mod.f:206
ps::ps_massflow_g
double precision, dimension(dimension_ps) ps_massflow_g
Definition: ps_mod.f:48
geometry::imin3
integer imin3
Definition: geometry_mod.f:90
ps::ps_defined
logical, dimension(dimension_ps) ps_defined
Definition: ps_mod.f:29
is
Definition: is_mod.f:11
ps::ps_w_g
double precision, dimension(dimension_ps) ps_w_g
Definition: ps_mod.f:53
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
geometry::cyclic_z_pd
logical cyclic_z_pd
Definition: geometry_mod.f:159
source_w_g
subroutine source_w_g(A_M, B_M)
Definition: source_w_g.f:24
bc::bc_plane
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
ps::ps_volume
double precision, dimension(dimension_ps) ps_volume
Definition: ps_mod.f:84
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
ps::ps_k_b
integer, dimension(dimension_ps) ps_k_b
Definition: ps_mod.f:44
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
wall_function
subroutine wall_function(IJK1, IJK2, ODX_WF, W_F_Slip)
Definition: source_u_g.f:958
fun_avg
Definition: fun_avg_mod.f:1
cutcell::blocked_w_cell_at
logical, dimension(:), allocatable blocked_w_cell_at
Definition: cutcell_mod.f:366
geometry::kmax1
integer kmax1
Definition: geometry_mod.f:58
bc::bc_k_t
integer, dimension(dimension_bc) bc_k_t
Definition: bc_mod.f:74
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
compar::dead_cell_at
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
geometry::odx_e
double precision, dimension(:), allocatable odx_e
Definition: geometry_mod.f:121
tau_g
Definition: tau_g_mod.f:1
indices::jm1
integer, dimension(:), allocatable jm1
Definition: indices_mod.f:51
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
geometry::jmax2
integer jmax2
Definition: geometry_mod.f:63
bc::bc_j_s
integer, dimension(dimension_bc) bc_j_s
Definition: bc_mod.f:62
run::m_am
integer m_am
Definition: run_mod.f:94
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
run::odt
double precision odt
Definition: run_mod.f:54
visc_g::epmu_gt
double precision, dimension(:), allocatable epmu_gt
Definition: visc_g_mod.f:13
geometry::jmax3
integer jmax3
Definition: geometry_mod.f:91
mms
Definition: mms_mod.f:12
bc::bc_hw_g
double precision, dimension(dimension_bc) bc_hw_g
Definition: bc_mod.f:307
param::north
integer north
Definition: param_mod.f:37
bc::bc_defined
logical, dimension(dimension_bc) bc_defined
Definition: bc_mod.f:207
is::is_pc
double precision, dimension(dimension_is, 2) is_pc
Definition: is_mod.f:85
param::south
integer south
Definition: param_mod.f:41
ps::ps_k_t
integer, dimension(dimension_ps) ps_k_t
Definition: ps_mod.f:45
indices::kp1
integer, dimension(:), allocatable kp1
Definition: indices_mod.f:52
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
geometry::axz
double precision, dimension(:), allocatable axz
Definition: geometry_mod.f:208
geometry::ayz_w
double precision, dimension(:), allocatable ayz_w
Definition: geometry_mod.f:236
fldvar::w_so
double precision, dimension(:,:), allocatable w_so
Definition: fldvar_mod.f:120
fldvar::w_go
double precision, dimension(:), allocatable w_go
Definition: fldvar_mod.f:114
param
Definition: param_mod.f:2
toleranc::dil_ep_s
double precision, parameter dil_ep_s
Definition: toleranc_mod.f:24
geometry::jmin3
integer jmin3
Definition: geometry_mod.f:90
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
param::dimension_ps
integer, parameter dimension_ps
Definition: param_mod.f:65
cutcell::cut_w_treatment_at
logical, dimension(:), allocatable cut_w_treatment_at
Definition: cutcell_mod.f:352
fldvar::rop_go
double precision, dimension(:), allocatable rop_go
Definition: fldvar_mod.f:41
geometry::kmax3
integer kmax3
Definition: geometry_mod.f:91
run::k_epsilon
logical k_epsilon
Definition: run_mod.f:97
physprop::mu_g
double precision, dimension(:), allocatable mu_g
Definition: physprop_mod.f:68
geometry::cylindrical
logical cylindrical
Definition: geometry_mod.f:168
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
visc_g
Definition: visc_g_mod.f:1
mms::use_mms
logical use_mms
Definition: mms_mod.f:15
param::west
integer west
Definition: param_mod.f:33
geometry::kmin3
integer kmin3
Definition: geometry_mod.f:90
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
fldvar::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
param::top
integer top
Definition: param_mod.f:45
drag::f_gs
double precision, dimension(:,:), allocatable f_gs
Definition: drag_mod.f:14
tau_g::tau_w_g
double precision, dimension(:), allocatable tau_w_g
Definition: tau_g_mod.f:6
ps
Definition: ps_mod.f:22
run::model_b
logical model_b
Definition: run_mod.f:88
fldvar::rop_s
double precision, dimension(:,:), allocatable rop_s
Definition: fldvar_mod.f:51
physprop::smax
integer smax
Definition: physprop_mod.f:22
drag::beta_ij
double precision, dimension(:,:,:), allocatable beta_ij
Definition: drag_mod.f:20
source_w_g_bc
subroutine source_w_g_bc(A_M, B_M)
Definition: source_w_g.f:447
ps::ps_j_s
integer, dimension(dimension_ps) ps_j_s
Definition: ps_mod.f:42
geometry::vol
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
cutcell::a_wpg_b
double precision, dimension(:), allocatable a_wpg_b
Definition: cutcell_mod.f:312
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
bc::bc_ww_g
double precision, dimension(dimension_bc) bc_ww_g
Definition: bc_mod.f:319
fldvar::ro_g
double precision, dimension(:), allocatable ro_g
Definition: fldvar_mod.f:32
ps::ps_i_e
integer, dimension(dimension_ps) ps_i_e
Definition: ps_mod.f:41
fldvar::rop_g
double precision, dimension(:), allocatable rop_g
Definition: fldvar_mod.f:38
param::bottom
integer bottom
Definition: param_mod.f:49
compar::kmap
integer, dimension(:), allocatable kmap
Definition: compar_mod.f:77
bc::bc_i_e
integer, dimension(dimension_bc) bc_i_e
Definition: bc_mod.f:58
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
ghdtheory::joiz
double precision, dimension(:,:), allocatable joiz
Definition: ghdtheory_mod.f:36
bc
Definition: bc_mod.f:23