MFIX  2016-1
source_u_g.f
Go to the documentation of this file.
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
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
compar::imap
integer, dimension(:), allocatable imap
Definition: compar_mod.f:77
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
mpi_utility
Definition: mpi_utility_mod.f:6
ghdtheory::joix
double precision, dimension(:,:), allocatable joix
Definition: ghdtheory_mod.f:30
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
bc::bc_uw_g
double precision, dimension(dimension_bc) bc_uw_g
Definition: bc_mod.f:313
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
visc_g::mu_gt
double precision, dimension(:), allocatable mu_gt
Definition: visc_g_mod.f:8
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
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
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
geometry::x_e
double precision, dimension(:), allocatable x_e
Definition: geometry_mod.f:134
run::momentum_x_eq
logical, dimension(0:dim_m) momentum_x_eq
Definition: run_mod.f:74
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
bc::delp_x
double precision delp_x
Definition: bc_mod.f:272
drag
Definition: drag_mod.f:11
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
bodyforce::bfx_g
double precision function bfx_g(IJK)
Definition: bodyforce_mod.f:20
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
param::east
integer east
Definition: param_mod.f:29
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
cutcell::cut_u_treatment_at
logical, dimension(:), allocatable cut_u_treatment_at
Definition: cutcell_mod.f:350
geometry::imin3
integer imin3
Definition: geometry_mod.f:90
ps::ps_defined
logical, dimension(dimension_ps) ps_defined
Definition: ps_mod.f:29
fldvar::u_go
double precision, dimension(:), allocatable u_go
Definition: fldvar_mod.f:90
is
Definition: is_mod.f:11
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
point_source_u_g
subroutine point_source_u_g(A_M, B_M)
Definition: source_u_g.f:1023
bc::bc_plane
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
turb
Definition: turb_mod.f:2
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
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
tau_g
Definition: tau_g_mod.f:1
indices::jm1
integer, dimension(:), allocatable jm1
Definition: indices_mod.f:51
geometry::imax1
integer imax1
Definition: geometry_mod.f:54
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
geometry::jmax2
integer jmax2
Definition: geometry_mod.f:63
source_u_g
subroutine source_u_g(A_M, B_M)
Definition: source_u_g.f:24
tau_g::tau_u_g
double precision, dimension(:), allocatable tau_u_g
Definition: tau_g_mod.f:4
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
cg_source_u_g_bc
subroutine cg_source_u_g_bc(A_M, B_M)
Definition: CG_source_u_g.f:289
cg_source_u_g
subroutine cg_source_u_g(A_M, B_M)
Definition: CG_source_u_g.f:13
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
cutcell::blocked_u_cell_at
logical, dimension(:), allocatable blocked_u_cell_at
Definition: cutcell_mod.f:364
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
mms::mms_u_g_src
double precision, dimension(:), allocatable mms_u_g_src
Definition: mms_mod.f:51
param::south
integer south
Definition: param_mod.f:41
geometry::kmax2
integer kmax2
Definition: geometry_mod.f:65
ps::ps_k_t
integer, dimension(dimension_ps) ps_k_t
Definition: ps_mod.f:45
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
fldvar::u_so
double precision, dimension(:,:), allocatable u_so
Definition: fldvar_mod.f:96
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
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
fldvar::rop_go
double precision, dimension(:), allocatable rop_go
Definition: fldvar_mod.f:41
geometry::kmax3
integer kmax3
Definition: geometry_mod.f:91
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
run::k_epsilon
logical k_epsilon
Definition: run_mod.f:97
geometry::cyclic_x_pd
logical cyclic_x_pd
Definition: geometry_mod.f:155
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
cutcell::a_upg_w
double precision, dimension(:), allocatable a_upg_w
Definition: cutcell_mod.f:235
drag::f_gs
double precision, dimension(:,:), allocatable f_gs
Definition: drag_mod.f:14
geometry::vol_u
double precision, dimension(:), allocatable vol_u
Definition: geometry_mod.f:224
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
source_u_g_bc
subroutine source_u_g_bc(A_M, B_M)
Definition: source_u_g.f:395
drag::beta_ij
double precision, dimension(:,:,:), allocatable beta_ij
Definition: drag_mod.f:20
ps::ps_u_g
double precision, dimension(dimension_ps) ps_u_g
Definition: ps_mod.f:51
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
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
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
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
bc
Definition: bc_mod.f:23
cutcell::a_upg_e
double precision, dimension(:), allocatable a_upg_e
Definition: cutcell_mod.f:234