MFIX  2016-1
source_v_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: SOURCE_V_g C
4 ! Purpose: Determine source terms for V_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: 7-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 ! C
20 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21  SUBROUTINE source_v_g(A_M, B_M)
22 
23 ! Modules
24 !---------------------------------------------------------------------//
25  USE bodyforce, only: bfy_g
26  USE bc, only: delp_y
27 
28  USE compar, only: ijkstart3, ijkend3, jmap
29 
30  USE drag, only: f_gs, beta_ij
31  USE fldvar, only: p_g, ro_g, rop_g, rop_go, rop_s
32  USE fldvar, only: ep_g, ep_s, epg_jfac
33  USE fldvar, only: v_g, v_go, u_s, v_s, w_s, v_so
34 
35  USE fun_avg, only: avg_x, avg_z, avg_y
36  USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
37  USE functions, only: ip_at_n, sip_at_n, is_id_at_n
38  USE functions, only: ip_of, jp_of, kp_of, im_of, jm_of, km_of
39  USE functions, only: north_of, south_of
40  USE functions, only: zmax
41  USE geometry, only: jmax1, cyclic_y_pd, do_k, xlength
42  USE geometry, only: vol, vol_v
43  USE geometry, only: axy, ayz, axz
44  USE geometry, only: ox, odx_e
45 
46  USE ghdtheory, only: joiy
47 
48  USE indices, only: i_of, j_of, k_of
49  USE indices, only: ip1, im1, kp1
50  USE is, only: is_pc
51 
52  USE mms, only: use_mms, mms_v_g_src
53  USE param
54  USE param1, only: zero, one, half
55  USE physprop, only: mmax, smax
56  USE physprop, only: mu_g, cv
57  USE run, only: momentum_y_eq
58  USE run, only: model_b, added_mass, m_am
59  USE run, only: kt_type_enum, drag_type_enum
60  USE run, only: ghd_2007, hys
61  USE run, only: odt
62  USE run, only: v_sh, shear
63  USE rxns, only: sum_r_g
64  USE scales, only: p_scale
65  USE vshear, only: vsh
66  USE tau_g, only: tau_v_g
67  USE toleranc, only: dil_ep_s
68  USE cutcell, only: cartesian_grid, cut_v_treatment_at
69  USE cutcell, only: blocked_v_cell_at
70  USE cutcell, only: a_vpg_n, a_vpg_s
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, IJKN, &
84  IMJK, IPJK, IJMK, IJPK, IJKP, IJKM, IMJPK, IJPKM
85 ! Phase index
86  INTEGER :: M, L, MM
87 ! Internal surface
88  INTEGER :: ISV
89 ! Pressure at north cell
90  DOUBLE PRECISION :: PgN
91 ! Average volume fraction
92  DOUBLE PRECISION :: EPGA, EPGAJ
93 ! Average density
94  DOUBLE PRECISION :: ROPGA, ROGA
95 ! Average viscosity
96  DOUBLE PRECISION :: MUGA
97 ! Source terms (Surface)
98  DOUBLE PRECISION :: Sdp
99 ! Source terms (Volumetric)
100  DOUBLE PRECISION :: V0, Vpm, Vmt, Vbf
101 ! Source terms (Volumetric) for GHD theory
102  DOUBLE PRECISION :: Ghd_drag, avgRop
103 ! Source terms for HYS drag relation
104  DOUBLE PRECISION :: HYS_drag, avgDrag
105 ! virtual (added) mass
106  DOUBLE PRECISION :: ROP_MA, Vsn, Vss, U_se, Usw, Vse, Vsw, &
107  Wst, Wsb, Vst, Vsb
108  DOUBLE PRECISION :: F_vir
109 ! loezos: shear terms
110  DOUBLE PRECISION :: VSH_n,VSH_s,VSH_e,VSH_w,VSH_p,Source_conv
111  DOUBLE PRECISION :: SRT
112 ! jackson terms: local stress tensor quantity
113  DOUBLE PRECISION :: ltau_v_g
114 !---------------------------------------------------------------------//
115 
116 ! Set reference phase to gas
117  m = 0
118 
119  IF (.NOT.momentum_y_eq(0)) RETURN
120 
121 
122 !$omp parallel do default(shared) &
123 !$omp private(I, J, K, IJK, IJKN, IMJK, IPJK, IJMK, IJPK, IMJPK, &
124 !$omp IJKM, IJPKM, IJKP, EPGA, EPGAJ, PGN, SDP, ROPGA, &
125 !$omp ROGA, ROP_MA, V0, ISV, MUGA, Vpm, Vmt, Vbf, L, MM, &
126 !$omp Vsn, Vss, U_se, Usw, Vse, Vsw, Wst, Wsb, Vst, &
127 !$omp Vsb, F_vir, Ghd_drag, avgRop, avgDrag, HYS_drag, &
128 !$omp VSH_n, VSH_s, VSH_e, VSH_w, VSH_p, Source_conv, SRT, &
129 !$omp ltau_v_g)
130  DO ijk = ijkstart3, ijkend3
131  i = i_of(ijk)
132  j = j_of(ijk)
133  k = k_of(ijk)
134  ijkn = north_of(ijk)
135  imjk = im_of(ijk)
136  ipjk = ip_of(ijk)
137  ijmk = jm_of(ijk)
138  ijpk = jp_of(ijk)
139  imjpk = im_of(ijpk)
140  ijkm = km_of(ijk)
141  ijpkm = km_of(ijpk)
142  ijkp = kp_of(ijk)
143 
144  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
145 ! if jackson avg ep_g otherwise 1
146  epgaj = avg_y(epg_jfac(ijk),epg_jfac(ijkn),j)
147 
148 ! Impermeable internal surface
149  IF (ip_at_n(ijk)) THEN
150  a_m(ijk,east,m) = zero
151  a_m(ijk,west,m) = zero
152  a_m(ijk,north,m) = zero
153  a_m(ijk,south,m) = zero
154  a_m(ijk,top,m) = zero
155  a_m(ijk,bottom,m) = zero
156  a_m(ijk,0,m) = -one
157  b_m(ijk,m) = zero
158 
159 ! dilute flow
160  ELSEIF (epga <= dil_ep_s) THEN
161  a_m(ijk,east,m) = zero
162  a_m(ijk,west,m) = zero
163  a_m(ijk,north,m) = zero
164  a_m(ijk,south,m) = zero
165  a_m(ijk,top,m) = zero
166  a_m(ijk,bottom,m) = zero
167  a_m(ijk,0,m) = -one
168  b_m(ijk,m) = zero
169  IF (ep_g(south_of(ijk)) > dil_ep_s) THEN
170  a_m(ijk,south,m) = one
171  ELSE IF (ep_g(north_of(ijk)) > dil_ep_s) THEN
172  a_m(ijk,north,m) = one
173  ELSE
174  b_m(ijk,m) = -v_g(ijk)
175  ENDIF
176 
177 ! Cartesian grid implementation
178  ELSEIF (blocked_v_cell_at(ijk)) THEN
179  a_m(ijk,east,m) = zero
180  a_m(ijk,west,m) = zero
181  a_m(ijk,north,m) = zero
182  a_m(ijk,south,m) = zero
183  a_m(ijk,top,m) = zero
184  a_m(ijk,bottom,m) = zero
185  a_m(ijk,0,m) = -one
186  b_m(ijk,m) = zero
187 
188 ! Normal case
189  ELSE
190 
191 ! Surface forces
192 ! Pressure term
193  pgn = p_g(ijkn)
194  IF (cyclic_y_pd) THEN
195  IF (jmap(j_of(ijk)).EQ.jmax1)pgn = p_g(ijkn) - delp_y
196  ENDIF
197  IF (model_b) THEN
198  IF(.NOT.cut_v_treatment_at(ijk)) THEN
199  sdp = -p_scale*(pgn - p_g(ijk))*axz(ijk)
200  ELSE
201  sdp = -p_scale*(pgn * a_vpg_n(ijk) - &
202  p_g(ijk) * a_vpg_s(ijk) )
203  ENDIF
204  ELSE
205  IF(.NOT.cut_v_treatment_at(ijk)) THEN
206  sdp = -p_scale*epga*(pgn - p_g(ijk))*axz(ijk)
207  ELSE
208  sdp = -p_scale*epga*(pgn * a_vpg_n(ijk) - &
209  p_g(ijk) * a_vpg_s(ijk) )
210  ENDIF
211  ENDIF
212 
213  IF(.NOT.cut_v_treatment_at(ijk)) THEN
214 ! Volumetric forces
215  ropga = avg_y(rop_g(ijk),rop_g(ijkn),j)
216  roga = avg_y(ro_g(ijk),ro_g(ijkn),j)
217 ! Previous time step
218  v0 = avg_y(rop_go(ijk),rop_go(ijkn),j)*odt
219 ! Added mass implicit transient term {Cv eps rop_g dV/dt}
220  IF(added_mass) THEN
221  rop_ma = avg_y(rop_g(ijk)*ep_s(ijk,m_am),&
222  rop_g(ijkn)*ep_s(ijkn,m_am),j)
223  v0 = v0 + cv * rop_ma * odt
224  ENDIF
225  ELSE
226 ! Volumetric forces
227  ropga = (vol(ijk)*rop_g(ijk) + &
228  vol(ijkn)*rop_g(ijkn))/(vol(ijk) + vol(ijkn))
229  roga = (vol(ijk)*ro_g(ijk) + &
230  vol(ijkn)*ro_g(ijkn) )/(vol(ijk) + vol(ijkn))
231 ! Previous time step
232  v0 = (vol(ijk)*rop_go(ijk) + vol(ijkn)*rop_go(ijkn))*&
233  odt/(vol(ijk) + vol(ijkn))
234 ! Added mass implicit transient term {Cv eps rop_g dV/dt}
235  IF(added_mass) THEN
236  rop_ma = (vol(ijk)*rop_g(ijk)*ep_s(ijk,m_am) + &
237  vol(ijkn)*rop_g(ijkn)*ep_s(ijkn,m_am) )/&
238  (vol(ijk) + vol(ijkn))
239  v0 = v0 + cv * rop_ma * odt
240  ENDIF
241  ENDIF
242 
243 ! VIRTUAL MASS SECTION (explicit terms)
244 ! adding transient term dVs/dt to virtual mass term
245  f_vir = zero
246  IF(added_mass.AND.(.NOT.cut_v_treatment_at(ijk))) THEN
247  f_vir = ((v_s(ijk,m_am) - v_so(ijk,m_am)))*&
248  odt*vol_v(ijk)
249 
250 ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
251  vss = avg_y_n(v_s(ijmk,m_am),v_s(ijk,m_am))
252  vsn = avg_y_n(v_s(ijk,m_am),v_s(ijpk,m_am))
253  u_se = avg_y(u_s(ijk,m_am),u_s(ijpk,m_am),j)
254  usw = avg_y(u_s(imjk,m_am),u_s(imjpk,m_am),j)
255  vse = avg_x(v_s(ijk,m_am),v_s(ipjk,m_am),ip1(i))
256  vsw = avg_x(v_s(imjk,m_am),v_s(ijk,m_am),i)
257  IF(do_k) THEN
258  wst = avg_y(w_s(ijk,m_am),w_s(ijpk,m_am),j)
259  wsb = avg_y(w_s(ijkm,m_am),w_s(ijpkm,m_am),j)
260  vst = avg_z(v_s(ijk,m_am),v_s(ijkp,m_am),kp1(k))
261  vsb = avg_z(v_s(ijkm,m_am),v_s(ijk,m_am),k)
262  f_vir = f_vir + avg_z_t(wsb,wst)*ox(i) * &
263  (vst - vsb)*axy(ijk)
264  ENDIF
265 ! adding convective terms (U dV/dx + V dV/dy) to virtual mass; W dV/dz added above.
266  f_vir = f_vir + v_s(ijk,m_am)*(vsn - vss)*axz(ijk) + &
267  avg_x_e(usw,u_se,ip1(i))*(vse - vsw)*ayz(ijk)
268  f_vir = f_vir * cv * rop_ma
269  ENDIF
270 
271 ! pressure drop through porous media
272  IF (sip_at_n(ijk)) THEN
273  isv = is_id_at_n(ijk)
274  muga = avg_y(mu_g(ijk),mu_g(ijkn),j)
275  vpm = muga/is_pc(isv,1)
276  IF (is_pc(isv,2) /= zero) vpm = vpm + half*is_pc(isv,2)*&
277  ropga*abs(v_g(ijk))
278  ELSE
279  vpm = zero
280  ENDIF
281 
282 ! Interphase mass transfer
283  IF(.NOT.cut_v_treatment_at(ijk)) THEN
284  vmt = avg_y(sum_r_g(ijk),sum_r_g(ijkn),j)
285  ELSE
286  vmt = (vol(ijk)*sum_r_g(ijk) + vol(ijkn)*sum_r_g(ijkn))/&
287  (vol(ijk) + vol(ijkn))
288  ENDIF
289 
290 ! Body force
291  IF (model_b) THEN
292  vbf = roga*bfy_g(ijk)
293  ELSE ! Model A
294  vbf = ropga*bfy_g(ijk)
295  ENDIF
296 
297 ! Additional force for GHD from darg force sum(beta_ig * Joi/rhop_i)
298  ghd_drag = zero
299  IF (kt_type_enum .EQ. ghd_2007) THEN
300  DO l = 1,smax
301  avgrop = avg_y(rop_s(ijk,l),rop_s(ijkn,l),j)
302  if(avgrop > zero) ghd_drag = ghd_drag +&
303  avg_y(f_gs(ijk,l),f_gs(ijkn,l),j) * joiy(ijk,l) / avgrop
304  ENDDO
305  ENDIF
306 
307 ! Additional force for HYS drag force, do not use with mixture GHD theory
308  avgdrag = zero
309  hys_drag = zero
310  IF (drag_type_enum .EQ. hys .AND. kt_type_enum .NE. ghd_2007) THEN
311  DO mm=1,mmax
312  DO l = 1,mmax
313  IF (l /= mm) THEN
314  avgdrag = avg_y(beta_ij(ijk,mm,l),beta_ij(ijkn,mm,l),j)
315  hys_drag = hys_drag + avgdrag * (v_g(ijk) - v_s(ijk,l))
316  ENDIF
317  ENDDO
318  ENDDO
319  ENDIF
320 ! if jackson, implement jackson form of governing equations (ep_g dot
321 ! del tau_g): multiply by void fraction otherwise by 1
322  ltau_v_g = epgaj*tau_v_g(ijk)
323 
324 
325 ! loezos: Shear Source terms from convective mom. flux
326  IF (shear) THEN
327  srt=(2d0*v_sh/xlength)
328  vsh_p=vsh(ijk)
329  vsh_n=vsh_p
330  vsh_s=vsh_p
331  vsh_e=vsh(ijk)+srt*1d0/odx_e(i)
332  vsh_w=vsh(ijk)-srt*1d0/odx_e(im1(i))
333  source_conv=a_m(ijk,north,m)*vsh_n + a_m(ijk,south,m)*vsh_s +&
334  a_m(ijk,west,m)*vsh_w + a_m(ijk,east,m)*vsh_e - &
335  (a_m(ijk,north,m) + a_m(ijk,south,m) + a_m(ijk,west,m) + &
336  a_m(ijk,east,m))*vsh_p
337  ELSE
338  source_conv=0d0
339  ENDIF
340 
341 ! Collect the terms
342  a_m(ijk,0,m) = -(a_m(ijk,east,m)+a_m(ijk,west,m)+&
343  a_m(ijk,north,m)+a_m(ijk,south,m)+a_m(ijk,top,m)+a_m(ijk,bottom,m)+&
344  (v0+vpm+zmax(vmt))*vol_v(ijk))
345  b_m(ijk,m) = b_m(ijk,m) - (sdp + ltau_v_g + f_vir + &
346  source_conv + ((v0+zmax((-vmt)))*v_go(ijk) + vbf + &
347  ghd_drag + hys_drag)*vol_v(ijk) )
348 ! MMS Source term.
349  IF(use_mms) b_m(ijk,m) = &
350  b_m(ijk,m) - mms_v_g_src(ijk)*vol_v(ijk)
351 
352  ENDIF
353  ENDDO
354 !$omp end parallel do
355 
356 ! modifications for cartesian grid implementation
357  IF(cartesian_grid) CALL cg_source_v_g(a_m, b_m)
358 ! modifications for bc
359  CALL source_v_g_bc(a_m, b_m)
360 ! modifications for cartesian grid implementation
361  IF(cartesian_grid) CALL cg_source_v_g_bc(a_m, b_m)
362 
363  RETURN
364  END SUBROUTINE source_v_g
365 
366 
367 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
368 ! C
369 ! Subroutine: SOURCE_V_g_BC C
370 ! Purpose: Determine source terms for V_g momentum eq. The terms C
371 ! appear in the center coefficient and RHS vector. The center C
372 ! coefficient and source vector are negative. The off-diagonal C
373 ! coefficients are positive. C
374 ! The drag terms are excluded from the source at this stage C
375 ! C
376 ! Author: M. Syamlal Date: 7-JUN-96 C
377 ! Reviewer: Date: C
378 ! C
379 ! C
380 ! C
381 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
382 
383  SUBROUTINE source_v_g_bc(A_M, B_M)
384 
385 !-----------------------------------------------
386 ! Modules
387 !-----------------------------------------------
388  USE param
389  USE param1
390  USE parallel
391  USE scales
392  USE constant
393  USE physprop
394  USE fldvar
395  USE visc_g
396  USE rxns
397  USE run
398  USE toleranc
399  USE geometry
400  USE indices
401  USE is
402  USE tau_g
403  USE bc
404  USE output
405  USE compar
406  USE fun_avg
407  USE functions
408  IMPLICIT NONE
409 !-----------------------------------------------
410 ! Dummy Arguments
411 !-----------------------------------------------
412 ! Septadiagonal matrix A_m
413  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
414 ! Vector b_m
415  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
416 !-----------------------------------------------
417 ! Local Variables
418 !-----------------------------------------------
419 ! Boundary condition
420  INTEGER :: L
421 ! Indices
422  INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK, &
423  IM, KM, IJKS, IJMK, IJPK
424 ! Phase index
425  INTEGER :: M
426 ! Turbulent shear stress
427  DOUBLE PRECISION :: W_F_Slip
428 !-----------------------------------------------
429 
430 ! Set reference phase to gas
431  m = 0
432 
433 
434 ! Set the default boundary conditions
435 ! The NS default setting is the where bc_type='dummy' or any default
436 ! (i.e., bc_type=undefined) wall boundary regions are handled. Note that
437 ! the north and south xz planes do not have to be explicitly addressed for
438 ! the v-momentum equation. In this direction the velocities are defined
439 ! at the wall (due staggered grid). They are defined as zero for a
440 ! no penetration condition (see zero_norm_vel subroutine and code under
441 ! the ip_at_n branch in the above source routine).
442 ! ---------------------------------------------------------------->>>
443  IF (do_k) THEN
444 ! bottom xy plane
445  k1 = 1
446  DO j1 = jmin3,jmax3
447  DO i1 = imin3, imax3
448  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
449  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
450  ijk = funijk(i1,j1,k1)
451  IF (ns_wall_at(ijk)) THEN
452 ! Setting the wall velocity to zero (set the boundary cell value equal
453 ! and oppostive to the adjacent fluid cell value)
454  a_m(ijk,east,m) = zero
455  a_m(ijk,west,m) = zero
456  a_m(ijk,north,m) = zero
457  a_m(ijk,south,m) = zero
458  a_m(ijk,top,m) = -one
459  a_m(ijk,bottom,m) = zero
460  a_m(ijk,0,m) = -one
461  b_m(ijk,m) = zero
462  ELSEIF (fs_wall_at(ijk)) THEN
463 ! Setting the wall velocity equal to the adjacent fluid velocity (set
464 ! the boundary cell value equal to 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  ENDIF
474  ENDDO
475  ENDDO
476 
477 ! top xy plane
478  k1 = kmax2
479  DO j1 = jmin3,jmax3
480  DO i1 = imin3, imax3
481  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
482  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
483  ijk = funijk(i1,j1,k1)
484  IF (ns_wall_at(ijk)) THEN
485  a_m(ijk,east,m) = zero
486  a_m(ijk,west,m) = zero
487  a_m(ijk,north,m) = zero
488  a_m(ijk,south,m) = zero
489  a_m(ijk,top,m) = zero
490  a_m(ijk,bottom,m) = -one
491  a_m(ijk,0,m) = -one
492  b_m(ijk,m) = zero
493  ELSE IF (fs_wall_at(ijk)) THEN
494  a_m(ijk,east,m) = zero
495  a_m(ijk,west,m) = zero
496  a_m(ijk,north,m) = zero
497  a_m(ijk,south,m) = zero
498  a_m(ijk,top,m) = zero
499  a_m(ijk,bottom,m) = one
500  a_m(ijk,0,m) = -one
501  b_m(ijk,m) = zero
502  ENDIF
503  ENDDO
504  ENDDO
505  ENDIF ! end if (do_k)
506 
507 
508 ! west zy plane
509  i1 = 1
510  DO k1 = kmin3, kmax3
511  DO j1 = jmin3, jmax3
512  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
513  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
514  ijk = funijk(i1,j1,k1)
515  IF (ns_wall_at(ijk)) THEN
516  a_m(ijk,east,m) = -one
517  a_m(ijk,west,m) = zero
518  a_m(ijk,north,m) = zero
519  a_m(ijk,south,m) = zero
520  a_m(ijk,top,m) = zero
521  a_m(ijk,bottom,m) = zero
522  a_m(ijk,0,m) = -one
523  b_m(ijk,m) = zero
524  ELSEIF (fs_wall_at(ijk)) THEN
525  a_m(ijk,east,m) = one
526  a_m(ijk,west,m) = zero
527  a_m(ijk,north,m) = zero
528  a_m(ijk,south,m) = zero
529  a_m(ijk,top,m) = zero
530  a_m(ijk,bottom,m) = zero
531  a_m(ijk,0,m) = -one
532  b_m(ijk,m) = zero
533  ENDIF
534  ENDDO
535  ENDDO
536 
537 ! east zy plane
538  i1 = imax2
539  DO k1 = kmin3, kmax3
540  DO j1 = jmin3, jmax3
541  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
542  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
543  ijk = funijk(i1,j1,k1)
544  IF (ns_wall_at(ijk)) THEN
545  a_m(ijk,east,m) = zero
546  a_m(ijk,west,m) = -one
547  a_m(ijk,north,m) = zero
548  a_m(ijk,south,m) = zero
549  a_m(ijk,top,m) = zero
550  a_m(ijk,bottom,m) = zero
551  a_m(ijk,0,m) = -one
552  b_m(ijk,m) = zero
553  ELSEIF (fs_wall_at(ijk)) THEN
554  a_m(ijk,east,m) = zero
555  a_m(ijk,west,m) = one
556  a_m(ijk,north,m) = zero
557  a_m(ijk,south,m) = zero
558  a_m(ijk,top,m) = zero
559  a_m(ijk,bottom,m) = zero
560  a_m(ijk,0,m) = -one
561  b_m(ijk,m) = zero
562  ENDIF
563  ENDDO
564  ENDDO
565 ! End setting the default boundary conditions
566 ! ----------------------------------------------------------------<<<
567 
568 ! Setting user specified boundary conditions
569  DO l = 1, dimension_bc
570  IF (bc_defined(l)) THEN
571 
572 ! Setting wall boundary conditions
573 ! ---------------------------------------------------------------->>>
574  IF (bc_type_enum(l) == no_slip_wall .AND. .NOT. k_epsilon) THEN
575  i1 = bc_i_w(l)
576  i2 = bc_i_e(l)
577  j1 = bc_j_s(l)
578  j2 = bc_j_n(l)
579  k1 = bc_k_b(l)
580  k2 = bc_k_t(l)
581  DO k = k1, k2
582  DO j = j1, j2
583  DO i = i1, i2
584  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
585  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
586  ijk = funijk(i,j,k)
587  IF (.NOT.wall_at(ijk)) cycle ! skip redefined cells
588  a_m(ijk,east,m) = zero
589  a_m(ijk,west,m) = zero
590  a_m(ijk,north,m) = zero
591  a_m(ijk,south,m) = zero
592  a_m(ijk,top,m) = zero
593  a_m(ijk,bottom,m) = zero
594  a_m(ijk,0,m) = -one
595  b_m(ijk,m) = zero
596  IF (fluid_at(east_of(ijk))) THEN
597  a_m(ijk,east,m) = -one
598  ELSEIF (fluid_at(west_of(ijk))) THEN
599  a_m(ijk,west,m) = -one
600  ELSEIF (fluid_at(top_of(ijk))) THEN
601  a_m(ijk,top,m) = -one
602  ELSEIF (fluid_at(bottom_of(ijk))) THEN
603  a_m(ijk,bottom,m) = -one
604  ENDIF
605  ENDDO
606  ENDDO
607  ENDDO
608 
609  ELSEIF (bc_type_enum(l) == free_slip_wall .AND. .NOT. k_epsilon) THEN
610  i1 = bc_i_w(l)
611  i2 = bc_i_e(l)
612  j1 = bc_j_s(l)
613  j2 = bc_j_n(l)
614  k1 = bc_k_b(l)
615  k2 = bc_k_t(l)
616  DO k = k1, k2
617  DO j = j1, j2
618  DO i = i1, i2
619  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
620  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
621  ijk = funijk(i,j,k)
622  IF (.NOT.wall_at(ijk)) cycle ! skip redefined cells
623  a_m(ijk,east,m) = zero
624  a_m(ijk,west,m) = zero
625  a_m(ijk,north,m) = zero
626  a_m(ijk,south,m) = zero
627  a_m(ijk,top,m) = zero
628  a_m(ijk,bottom,m) = zero
629  a_m(ijk,0,m) = -one
630  b_m(ijk,m) = zero
631  IF (fluid_at(east_of(ijk))) THEN
632  a_m(ijk,east,m) = one
633  ELSEIF (fluid_at(west_of(ijk))) THEN
634  a_m(ijk,west,m) = one
635  ELSEIF (fluid_at(top_of(ijk))) THEN
636  a_m(ijk,top,m) = one
637  ELSEIF (fluid_at(bottom_of(ijk))) THEN
638  a_m(ijk,bottom,m) = one
639  ENDIF
640  ENDDO
641  ENDDO
642  ENDDO
643 
644  ELSEIF (bc_type_enum(l) == par_slip_wall .AND. .NOT. k_epsilon) THEN
645  i1 = bc_i_w(l)
646  i2 = bc_i_e(l)
647  j1 = bc_j_s(l)
648  j2 = bc_j_n(l)
649  k1 = bc_k_b(l)
650  k2 = bc_k_t(l)
651  DO k = k1, k2
652  DO j = j1, j2
653  DO i = i1, i2
654  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
655  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
656  ijk = funijk(i,j,k)
657  IF (.NOT.wall_at(ijk)) cycle ! skip redefined cells
658  im = im1(i)
659  km = km1(k)
660  a_m(ijk,east,m) = zero
661  a_m(ijk,west,m) = zero
662  a_m(ijk,north,m) = zero
663  a_m(ijk,south,m) = zero
664  a_m(ijk,top,m) = zero
665  a_m(ijk,bottom,m) = zero
666  a_m(ijk,0,m) = -one
667  b_m(ijk,m) = zero
668  IF (fluid_at(east_of(ijk))) THEN
669  IF (bc_hw_g(l) == undefined) THEN
670  a_m(ijk,east,m) = -half
671  a_m(ijk,0,m) = -half
672  b_m(ijk,m) = -bc_vw_g(l)
673  ELSE
674  a_m(ijk,0,m) = -(half*bc_hw_g(l)+odx_e(i))
675  a_m(ijk,east,m) = -(half*bc_hw_g(l)-odx_e(i))
676  b_m(ijk,m) = -bc_hw_g(l)*bc_vw_g(l)
677  ENDIF
678  ELSEIF (fluid_at(west_of(ijk))) THEN
679  IF (bc_hw_g(l) == undefined) THEN
680  a_m(ijk,west,m) = -half
681  a_m(ijk,0,m) = -half
682  b_m(ijk,m) = -bc_vw_g(l)
683  ELSE
684  a_m(ijk,west,m) = -(half*bc_hw_g(l)-odx_e(im))
685  a_m(ijk,0,m) = -(half*bc_hw_g(l)+odx_e(im))
686  b_m(ijk,m) = -bc_hw_g(l)*bc_vw_g(l)
687  ENDIF
688  ELSEIF (fluid_at(top_of(ijk))) THEN
689  IF (bc_hw_g(l) == undefined) THEN
690  a_m(ijk,top,m) = -half
691  a_m(ijk,0,m) = -half
692  b_m(ijk,m) = -bc_vw_g(l)
693  ELSE
694  a_m(ijk,0,m) = -(half*bc_hw_g(l)+odz_t(k)*ox(i))
695  a_m(ijk,top,m) = -(half*bc_hw_g(l)-odz_t(k)*ox(i))
696  b_m(ijk,m) = -bc_hw_g(l)*bc_vw_g(l)
697  ENDIF
698  ELSEIF (fluid_at(bottom_of(ijk))) THEN
699  IF (bc_hw_g(l) == undefined) THEN
700  a_m(ijk,bottom,m) = -half
701  a_m(ijk,0,m) = -half
702  b_m(ijk,m) = -bc_vw_g(l)
703  ELSE
704  a_m(ijk,bottom,m) = -(half*bc_hw_g(l)-odz_t(km)*ox(i))
705  a_m(ijk,0,m) = -(half*bc_hw_g(l)+odz_t(km)*ox(i))
706  b_m(ijk,m) = -bc_hw_g(l)*bc_vw_g(l)
707  ENDIF
708  ENDIF
709  ENDDO
710  ENDDO
711  ENDDO
712 
713 ! Setting wall boundary conditions when K_EPSILON
714  ELSEIF (bc_type_enum(l) == par_slip_wall .OR. &
715  bc_type_enum(l) == no_slip_wall .OR. &
716  bc_type_enum(l) == free_slip_wall .AND. &
717  k_epsilon )THEN
718  i1 = bc_i_w(l)
719  i2 = bc_i_e(l)
720  j1 = bc_j_s(l)
721  j2 = bc_j_n(l)
722  k1 = bc_k_b(l)
723  k2 = bc_k_t(l)
724  DO k = k1, k2
725  DO j = j1, j2
726  DO i = i1, i2
727  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
728  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
729  ijk = funijk(i,j,k)
730  IF (.NOT.wall_at(ijk)) cycle ! skip redefined cells
731  im = im1(i)
732  km = km1(k)
733  a_m(ijk,east,m) = zero
734  a_m(ijk,west,m) = zero
735  a_m(ijk,north,m) = zero
736  a_m(ijk,south,m) = zero
737  a_m(ijk,top,m) = zero
738  a_m(ijk,bottom,m) = zero
739  a_m(ijk,0,m) = -one
740  b_m(ijk,m) = zero
741  IF (fluid_at(east_of(ijk))) THEN
742  CALL wall_function(ijk,east_of(ijk),odx_e(i),w_f_slip)
743  a_m(ijk,east,m) = w_f_slip
744  a_m(ijk,0,m) = -one
745  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_vw_g(l)
746  ELSEIF (fluid_at(west_of(ijk))) THEN
747  CALL wall_function(ijk,west_of(ijk),odx_e(im),w_f_slip)
748  a_m(ijk,west,m) = w_f_slip
749  a_m(ijk,0,m) = -one
750  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_vw_g(l)
751  ELSEIF (fluid_at(top_of(ijk))) THEN
752  CALL wall_function(ijk,top_of(ijk),odz_t(k)*ox(i),w_f_slip)
753  a_m(ijk,top,m) = w_f_slip
754  a_m(ijk,0,m) = -one
755  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_vw_g(l)
756  ELSEIF (fluid_at(bottom_of(ijk))) THEN
757  CALL wall_function(ijk,bottom_of(ijk),odz_t(km)*ox(i),w_f_slip)
758  a_m(ijk,bottom,m) = w_f_slip
759  a_m(ijk,0,m) = -one
760  IF (bc_type_enum(l) == par_slip_wall) b_m(ijk,m) = -bc_vw_g(l)
761  ENDIF
762  ENDDO
763  ENDDO
764  ENDDO
765 ! end setting of wall boundary conditions
766 ! ----------------------------------------------------------------<<<
767 
768 ! Setting p_inflow or p_outflow flow boundary conditions
769 ! ---------------------------------------------------------------->>>
770  ELSE IF (bc_type_enum(l)==p_inflow .OR. bc_type_enum(l)==p_outflow) THEN
771  IF (bc_plane(l) == 'S') THEN
772 ! if the fluid cell is on the south side of the outflow/inflow boundary
773 ! then set the velocity in the boundary cell equal to the velocity of
774 ! the adjacent fluid cell
775  i1 = bc_i_w(l)
776  i2 = bc_i_e(l)
777  j1 = bc_j_s(l)
778  j2 = bc_j_n(l)
779  k1 = bc_k_b(l)
780  k2 = bc_k_t(l)
781  DO k = k1, k2
782  DO j = j1, j2
783  DO i = i1, i2
784  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
785  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
786  ijk = funijk(i,j,k)
787  a_m(ijk,east,m) = zero
788  a_m(ijk,west,m) = zero
789  a_m(ijk,north,m) = zero
790  a_m(ijk,south,m) = one
791  a_m(ijk,top,m) = zero
792  a_m(ijk,bottom,m) = zero
793  a_m(ijk,0,m) = -one
794  b_m(ijk,m) = zero
795  ENDDO
796  ENDDO
797  ENDDO
798  ENDIF
799 ! end setting of p_inflow or p_otuflow flow boundary conditions
800 ! ----------------------------------------------------------------<<<
801 
802 ! Setting outflow flow boundary conditions
803 ! ---------------------------------------------------------------->>>
804  ELSEIF (bc_type_enum(l) == outflow) THEN
805  IF (bc_plane(l) == 'S') 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  a_m(ijk,east,m) = zero
819  a_m(ijk,west,m) = zero
820  a_m(ijk,north,m) = zero
821  a_m(ijk,south,m) = one
822  a_m(ijk,top,m) = zero
823  a_m(ijk,bottom,m) = zero
824  a_m(ijk,0,m) = -one
825  b_m(ijk,m) = zero
826  ijmk = jm_of(ijk)
827  a_m(ijmk,east,m) = zero
828  a_m(ijmk,west,m) = zero
829  a_m(ijmk,north,m) = zero
830  a_m(ijmk,south,m) = one
831  a_m(ijmk,top,m) = zero
832  a_m(ijmk,bottom,m) = zero
833  a_m(ijmk,0,m) = -one
834  b_m(ijmk,m) = zero
835  ENDDO
836  ENDDO
837  ENDDO
838  ELSEIF (bc_plane(l) == 'N') THEN
839  i1 = bc_i_w(l)
840  i2 = bc_i_e(l)
841  j1 = bc_j_s(l)
842  j2 = bc_j_n(l)
843  k1 = bc_k_b(l)
844  k2 = bc_k_t(l)
845  DO k = k1, k2
846  DO j = j1, j2
847  DO i = i1, i2
848  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
849  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
850  ijk = funijk(i,j,k)
851  ijpk = jp_of(ijk)
852  a_m(ijpk,east,m) = zero
853  a_m(ijpk,west,m) = zero
854  a_m(ijpk,north,m) = one
855  a_m(ijpk,south,m) = zero
856  a_m(ijpk,top,m) = zero
857  a_m(ijpk,bottom,m) = zero
858  a_m(ijpk,0,m) = -one
859  b_m(ijpk,m) = zero
860  ENDDO
861  ENDDO
862  ENDDO
863  ENDIF
864 ! end setting of outflow flow boundary conditions
865 ! ----------------------------------------------------------------<<<
866 
867 ! Setting bc that are defined but not nsw, fsw, psw, p_inflow,
868 ! p_outflow, or outflow (at this time, this section addresses
869 ! mass_inflow and mass_outflow type boundaries)
870 ! ---------------------------------------------------------------->>>
871  ELSE
872  i1 = bc_i_w(l)
873  i2 = bc_i_e(l)
874  j1 = bc_j_s(l)
875  j2 = bc_j_n(l)
876  k1 = bc_k_b(l)
877  k2 = bc_k_t(l)
878  DO k = k1, k2
879  DO j = j1, j2
880  DO i = i1, i2
881  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
882  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
883  ijk = funijk(i,j,k)
884 ! setting the velocity in the boundary cell equal to what is known
885  a_m(ijk,east,m) = zero
886  a_m(ijk,west,m) = zero
887  a_m(ijk,north,m) = zero
888  a_m(ijk,south,m) = zero
889  a_m(ijk,top,m) = zero
890  a_m(ijk,bottom,m) = zero
891  a_m(ijk,0,m) = -one
892  b_m(ijk,m) = -v_g(ijk)
893  IF (bc_plane(l) == 'S') THEN
894 ! if the fluid cell is on the south side of the outflow/inflow boundary
895 ! then set the velocity in the adjacent fluid cell equal to what is
896 ! known in that cell
897  ijks = south_of(ijk)
898  a_m(ijks,east,m) = zero
899  a_m(ijks,west,m) = zero
900  a_m(ijks,north,m) = zero
901  a_m(ijks,south,m) = zero
902  a_m(ijks,top,m) = zero
903  a_m(ijks,bottom,m) = zero
904  a_m(ijks,0,m) = -one
905  b_m(ijks,m) = -v_g(ijks)
906  ENDIF
907  ENDDO
908  ENDDO
909  ENDDO
910  ENDIF ! end if/else (bc_type)
911  ! ns, fs, psw; else
912  ! p_inflow, p_outflow, or outflow; else
913 ! end setting of 'else' flow boundary conditions
914 ! (mass_inflow/mass_outflow)
915 ! ----------------------------------------------------------------<<<
916 
917  ENDIF ! end if (bc_defined)
918  ENDDO ! end L do loop over dimension_bc
919 
920  RETURN
921  END SUBROUTINE source_v_g_bc
922 
923 
924 
925 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
926 ! C
927 ! Subroutine: POINT_SOURCE_V_G C
928 ! Purpose: Adds point sources to the gas phase V-Momentum equation. C
929 ! C
930 ! Author: J. Musser Date: 10-JUN-13 C
931 ! Reviewer: Date: C
932 ! C
933 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
934  SUBROUTINE point_source_v_g(A_M, B_M)
935 
936  use compar
937  use constant
938  use geometry
939  use indices
940  use param
941  use param1, only: small_number
942  use physprop
943  use ps
944  use run
945  use functions
946 
947  IMPLICIT NONE
948 !-----------------------------------------------
949 ! Dummy arguments
950 !-----------------------------------------------
951 ! Septadiagonal matrix A_m
952  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
953 ! Vector b_m
954  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
955 !-----------------------------------------------
956 ! Local Variables
957 !-----------------------------------------------
958 ! Indices
959  INTEGER :: IJK, I, J, K
960  INTEGER :: PSV, M
961  INTEGER :: lJN, lJS
962 ! terms of bm expression
963  DOUBLE PRECISION :: pSource
964 !-----------------------------------------------
965 
966 ! Set reference phase to gas
967  m = 0
968 
969 ! Calculate the mass going into each IJK cell. This is done for each
970 ! call in case the point source is time dependent.
971  ps_lp: do psv = 1, dimension_ps
972  if(.NOT.ps_defined(psv)) cycle ps_lp
973  if(abs(ps_v_g(psv)) < small_number) cycle ps_lp
974 
975  if(ps_v_g(psv) < 0.0d0) then
976  ljs = ps_j_s(psv) - 1
977  ljn = ps_j_n(psv) - 1
978  else
979  ljs = ps_j_s(psv)
980  ljn = ps_j_n(psv)
981  endif
982 
983  do k = ps_k_b(psv), ps_k_t(psv)
984  do j = ljs, ljn
985  do i = ps_i_w(psv), ps_i_e(psv)
986 
987  if(.NOT.is_on_mype_plus2layers(i,j,k)) cycle
988  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
989 
990  ijk = funijk(i,j,k)
991  if(.NOT.fluid_at(ijk)) cycle
992 
993  psource = ps_massflow_g(psv) * (vol(ijk)/ps_volume(psv))
994 
995  b_m(ijk,m) = b_m(ijk,m) - psource * &
996  ps_v_g(psv) * ps_vel_mag_g(psv)
997 
998  enddo
999  enddo
1000  enddo
1001 
1002  enddo ps_lp
1003 
1004  RETURN
1005  END SUBROUTINE point_source_v_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
run::momentum_y_eq
logical, dimension(0:dim_m) momentum_y_eq
Definition: run_mod.f:77
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
fldvar::v_so
double precision, dimension(:,:), allocatable v_so
Definition: fldvar_mod.f:108
ps::ps_v_g
double precision, dimension(dimension_ps) ps_v_g
Definition: ps_mod.f:52
geometry::imax2
integer imax2
Definition: geometry_mod.f:61
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
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
run::shear
logical shear
Definition: run_mod.f:175
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
cutcell::a_vpg_s
double precision, dimension(:), allocatable a_vpg_s
Definition: cutcell_mod.f:273
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
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
drag
Definition: drag_mod.f:11
rxns::sum_r_g
double precision, dimension(:), allocatable sum_r_g
Definition: rxns_mod.f:28
source_v_g
subroutine source_v_g(A_M, B_M)
Definition: source_v_g.f:22
point_source_v_g
subroutine point_source_v_g(A_M, B_M)
Definition: source_v_g.f:935
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
mms::mms_v_g_src
double precision, dimension(:), allocatable mms_v_g_src
Definition: mms_mod.f:52
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
fldvar::v_go
double precision, dimension(:), allocatable v_go
Definition: fldvar_mod.f:102
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
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
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
bc::bc_plane
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
vshear::vsh
double precision, dimension(:), allocatable vsh
Definition: vshear_mod.f:3
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
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
geometry::odx_e
double precision, dimension(:), allocatable odx_e
Definition: geometry_mod.f:121
tau_g
Definition: tau_g_mod.f:1
geometry::cyclic_y_pd
logical cyclic_y_pd
Definition: geometry_mod.f:157
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
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
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
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
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
ghdtheory::joiy
double precision, dimension(:,:), allocatable joiy
Definition: ghdtheory_mod.f:33
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
indices::kp1
integer, dimension(:), allocatable kp1
Definition: indices_mod.f:52
geometry::xlength
double precision xlength
Definition: geometry_mod.f:33
param1::half
double precision, parameter half
Definition: param1_mod.f:28
geometry::jmax1
integer jmax1
Definition: geometry_mod.f:56
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
cg_source_v_g_bc
subroutine cg_source_v_g_bc(A_M, B_M)
Definition: CG_source_v_g.f:282
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
param::dimension_ps
integer, parameter dimension_ps
Definition: param_mod.f:65
bodyforce::bfy_g
double precision function bfy_g(IJK)
Definition: bodyforce_mod.f:27
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
tau_g::tau_v_g
double precision, dimension(:), allocatable tau_v_g
Definition: tau_g_mod.f:5
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
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
bc::bc_vw_g
double precision, dimension(dimension_bc) bc_vw_g
Definition: bc_mod.f:316
cg_source_v_g
subroutine cg_source_v_g(A_M, B_M)
Definition: CG_source_v_g.f:13
cutcell::cut_v_treatment_at
logical, dimension(:), allocatable cut_v_treatment_at
Definition: cutcell_mod.f:351
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::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
param::top
integer top
Definition: param_mod.f:45
bc::delp_y
double precision delp_y
Definition: bc_mod.f:275
drag::f_gs
double precision, dimension(:,:), allocatable f_gs
Definition: drag_mod.f:14
compar::jmap
integer, dimension(:), allocatable jmap
Definition: compar_mod.f:77
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
cutcell::blocked_v_cell_at
logical, dimension(:), allocatable blocked_v_cell_at
Definition: cutcell_mod.f:365
cutcell::a_vpg_n
double precision, dimension(:), allocatable a_vpg_n
Definition: cutcell_mod.f:272
drag::beta_ij
double precision, dimension(:,:,:), allocatable beta_ij
Definition: drag_mod.f:20
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
source_v_g_bc
subroutine source_v_g_bc(A_M, B_M)
Definition: source_v_g.f:384
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
bc
Definition: bc_mod.f:23
run::v_sh
double precision v_sh
Definition: run_mod.f:177
geometry::vol_v
double precision, dimension(:), allocatable vol_v
Definition: geometry_mod.f:233