MFIX  2016-1
source_v_s.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine name: SOURCE_V_s C
4 ! Purpose: Determine source terms for V_s 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: Allow for partial-slip boundary conditions proposed by C
16 ! by Johnson & Jackson (1987) if the Granular Temperature C
17 ! equation is used. C
18 ! Author: K. Agrawal, Princeton University Date: 24-JAN-98 C
19 ! Reviewer: Date: dd-mmm-yy C
20 ! C
21 ! Revision Number: 2 C
22 ! Purpose: To incorporate Cartesian grid modifications C
23 ! Author: Jeff Dietiker Date: 01-Jul-09 C
24 ! C
25 ! C
26 ! Literature/Document References: C
27 ! C
28 ! C
29 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
30 
31  SUBROUTINE source_v_s(A_M, B_M)
32 
33 !-----------------------------------------------
34 ! Modules
35 !-----------------------------------------------
36  USE param
37  USE param1
38  USE parallel
39  USE scales
40  USE constant
41  USE physprop
42  USE fldvar
43  USE visc_s
44  USE rxns
45  USE run
46  USE toleranc
47  USE geometry
48  USE indices
49  USE is
50  USE tau_s
51  USE tau_g, only: ctau_v_g
52  USE bc
53  USE vshear
54  USE compar
55  USE sendrecv
56  use kintheory
57  USE ghdtheory
58  USE drag
59  USE cutcell
60  USE quadric
61  USE mms
62  USE bodyforce
63  USE fun_avg
64  USE functions
65 
66  IMPLICIT NONE
67 !-----------------------------------------------
68 ! Dummy Arguments
69 !-----------------------------------------------
70 ! Septadiagonal matrix A_m
71  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
72 
73 ! Vector b_m
74  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
75 !-----------------------------------------------
76 ! Local Variables
77 !-----------------------------------------------
78 ! Indices
79  INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM, IJKN, &
80  IPJK, IJPK, IJKP, IMJPK, IJPKM
81 ! Phase index
82  INTEGER :: M, MM, L
83 ! Internal surface
84  INTEGER :: ISV
85 ! Pressure at north cell
86  DOUBLE PRECISION :: PgN
87 ! Average volume fraction
88  DOUBLE PRECISION :: EPSA, EPStmp, epse, epsw, epsn, epss, &
89  epst, epsb, epsMix, epsMixN
90  DOUBLE PRECISION :: SUM_EPS_CP
91 ! Average density
92  DOUBLE PRECISION :: ROPSA
93 ! Average density difference
94  DOUBLE PRECISION :: dro1, dro2, droa
95 ! Source terms (Surface)
96  DOUBLE PRECISION :: Sdp, Sdps
97 ! Source terms (Volumetric)
98  DOUBLE PRECISION :: V0, Vmt, Vbf, Vmttmp
99 ! Source terms (Volumetric) for GHD theory
100  DOUBLE PRECISION :: Ghd_drag, avgRop
101 ! Source terms for HYS drag relation
102  DOUBLE PRECISION :: HYS_drag, avgDrag
103 ! virtual (added) mass
104  DOUBLE PRECISION :: ROP_MA, Vgn, Vgs, Uge, Ugw, Vge,&
105  Vgw, Wgt, Wgb, Vgt, Vgb
106  DOUBLE PRECISION :: F_vir
107 ! terms for shear (loezos)
108  DOUBLE PRECISION :: VSH_n,VSH_s,VSH_e,VSH_w,VSH_p,Source_conv
109  DOUBLE PRECISION :: SRT
110 !-----------------------------------------------
111 
112 
113  DO m = 1, mmax
114  IF(kt_type_enum /= ghd_2007 .OR. &
115  (kt_type_enum == ghd_2007 .AND. m==mmax)) THEN
116 
117  IF (momentum_y_eq(m)) THEN
118 
119 
120 !$omp parallel do default(shared) &
121 !$omp private( I, J, K, IJK, IJKN, IMJK, IPJK, IJMK, IJPK, IMJPK, &
122 !$omp IJKM, IJPKM, IJKP, ISV, epsMix, epsMixN, EPStmp, &
123 !$omp EPSA, EPSw, EPSe, EPSn, EPSs, EPSt, EPSb, SUM_EPS_CP,&
124 !$omp PGN, Sdp, Sdps, MM, L, ROPSA, V0, ROP_MA, &
125 !$omp Vgn, Vgs, Uge, Ugw, Vge, Vgw, Wgt, Wgb, Vgt, Vgb, &
126 !$omp F_vir, VMT, VMTtmp, DRO1, DRO2, DROA, Vbf, &
127 !$omp avgRop, Ghd_drag, HYS_drag, avgDrag, &
128 !$omp VSH_n, VSH_s, VSH_e, VSH_w, VSH_p, Source_conv, SRT)
129  DO ijk = ijkstart3, ijkend3
130 
131 ! Skip walls where some values are undefined.
132  IF(wall_at(ijk)) cycle
133 
134  i = i_of(ijk)
135  j = j_of(ijk)
136  k = k_of(ijk)
137  imjk = im_of(ijk)
138  ipjk = ip_of(ijk)
139  ijmk = jm_of(ijk)
140  ijpk = jp_of(ijk)
141  imjpk = im_of(ijpk)
142  ijkm = km_of(ijk)
143  ijkn = north_of(ijk)
144  ijpkm = km_of(ijpk)
145  ijkp = kp_of(ijk)
146 
147  IF (kt_type_enum == ghd_2007) THEN
148  epstmp = zero
149  epsmix = zero
150  epsmixn= zero
151  DO l = 1, smax
152  epstmp = epstmp + avg_y(ep_s(ijk,l),ep_s(ijkn,l),j)
153  epsmix = epsmix + ep_s(ijk,l) ! epsMix, epsMixN to be used for model B
154  epsmixn = epsmixn + ep_s(ijkn,l)
155  IF(ip_at_n(ijk)) THEN
156  v_s(ijk,l) = zero
157  ELSEIF(sip_at_n(ijk)) THEN
158  isv = is_id_at_n(ijk)
159  v_s(ijk,l) = is_vel_s(isv,l)
160  ENDIF
161  ENDDO
162  epsa = epstmp
163  ELSE
164  epsa = avg_y(ep_s(ijk,m),ep_s(ijkn,m),j)
165  ENDIF
166 
167 ! Impermeable internal surface
168  IF (ip_at_n(ijk)) THEN
169  a_m(ijk,east,m) = zero
170  a_m(ijk,west,m) = zero
171  a_m(ijk,north,m) = zero
172  a_m(ijk,south,m) = zero
173  a_m(ijk,top,m) = zero
174  a_m(ijk,bottom,m) = zero
175  a_m(ijk,0,m) = -one
176  b_m(ijk,m) = zero
177 
178 ! Semi-permeable internal surface
179  ELSEIF (sip_at_n(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  isv = is_id_at_n(ijk)
188  b_m(ijk,m) = -is_vel_s(isv,m)
189 
190 ! Dilute flow
191  ELSEIF (epsa <= dil_ep_s) THEN
192  a_m(ijk,east,m) = zero
193  a_m(ijk,west,m) = zero
194  a_m(ijk,north,m) = zero
195  a_m(ijk,south,m) = zero
196  a_m(ijk,top,m) = zero
197  a_m(ijk,bottom,m) = zero
198  a_m(ijk,0,m) = -one
199  b_m(ijk,m) = zero
200  IF (kt_type_enum == ghd_2007) THEN
201  epsw = zero
202  epse = zero
203  epsn = zero
204  epss = zero
205  epst = zero
206  epsb = zero
207  DO l = 1, smax
208  epsw = epsw + ep_s(west_of(ijk),l)
209  epse = epse + ep_s(east_of(ijk),l)
210  epsn = epsn + ep_s(north_of(ijk),l)
211  epss = epss + ep_s(south_of(ijk),l)
212  IF(.NOT. no_k) THEN
213  epst = epst + ep_s(top_of(ijk),l)
214  epsb = epsb + ep_s(bottom_of(ijk),l)
215  ENDIF
216  ENDDO
217  ELSE
218  epsw = ep_s(west_of(ijk),m)
219  epse = ep_s(east_of(ijk),m)
220  epsn = ep_s(north_of(ijk),m)
221  epss = ep_s(south_of(ijk),m)
222  IF(.NOT. no_k) THEN
223  epst = ep_s(top_of(ijk),m)
224  epsb = ep_s(bottom_of(ijk),m)
225  ENDIF
226  ENDIF
227 ! using the average boundary cell values to compute V_s (sof, Aug 23 2005)
228  IF (epsw > dil_ep_s .AND. .NOT.is_at_e(imjk)) a_m(ijk,west,m) = one
229  IF (epse > dil_ep_s .AND. .NOT.is_at_e(ijk)) a_m(ijk,east,m) = one
230  IF (epss > dil_ep_s .AND. .NOT.is_at_n(ijmk)) a_m(ijk,south,m) = one
231  IF (epsn > dil_ep_s .AND. .NOT.is_at_n(ijk)) a_m(ijk,north,m) = one
232  IF(.NOT. no_k) THEN
233  IF (epsb > dil_ep_s .AND. .NOT.is_at_t(ijkm)) a_m(ijk,bottom,m) = one
234  IF (epst > dil_ep_s .AND. .NOT.is_at_t(ijk)) a_m(ijk,top,m) = one
235  ENDIF
236  IF((a_m(ijk,west,m)+a_m(ijk,east,m)+a_m(ijk,south,m)+a_m(ijk,north,m)+ &
237  a_m(ijk,bottom,m)+a_m(ijk,top,m)) == zero) THEN
238  b_m(ijk,m) = -v_s(ijk,m)
239  ELSE
240  a_m(ijk,0,m) = -(a_m(ijk,east,m)+a_m(ijk,west,m)+a_m(ijk,north,m)+ &
241  a_m(ijk,south,m)+a_m(ijk,top,m)+a_m(ijk,bottom,m))
242  ENDIF
243 
244 ! Cartesian grid implementation
245  ELSEIF (blocked_v_cell_at(ijk)) THEN
246  a_m(ijk,east,m) = zero
247  a_m(ijk,west,m) = zero
248  a_m(ijk,north,m) = zero
249  a_m(ijk,south,m) = zero
250  a_m(ijk,top,m) = zero
251  a_m(ijk,bottom,m) = zero
252  a_m(ijk,0,m) = -one
253  b_m(ijk,m) = zero
254 
255 ! Normal case
256  ELSE
257 
258 ! Surface forces:
259 ! Pressure term
260  pgn = p_g(ijkn)
261  IF (cyclic_y_pd) THEN
262 ! CYCLIC_AT_N Flag is not set correctly in DMP and causes issues. This
263 ! is avoided by using the DMP cyclic map. The flags need fixed.
264 ! IF (CYCLIC_AT_N(IJK)) PGN = P_G(IJKN) - DELP_Y
265  IF (jmap(j_of(ijk)).EQ.jmax1)pgn = p_g(ijkn) - delp_y
266  ENDIF
267  IF (model_b) THEN
268  sdp = zero
269  ELSE
270  IF(.NOT.cut_v_treatment_at(ijk)) THEN
271  sdp = -p_scale*epsa*(pgn - p_g(ijk))*axz(ijk)
272  ELSE
273  sdp = -p_scale*epsa*(pgn * a_vpg_n(ijk) - p_g(ijk) * a_vpg_s(ijk) )
274  ENDIF
275  ENDIF
276 
277  IF (close_packed(m)) THEN
278  IF(smax > 1 .AND. kt_type_enum /= ghd_2007) THEN
279  sum_eps_cp=0.0
280  DO mm=1,smax
281  IF (close_packed(mm))&
282  sum_eps_cp=sum_eps_cp+avg_y(ep_s(ijk,mm),ep_s(ijkn,mm),j)
283  ENDDO
284  sum_eps_cp = max(sum_eps_cp, small_number)
285  sdps = - ((p_s(ijkn,m)-p_s(ijk,m))+(epsa/sum_eps_cp)* &
286  (p_star(ijkn)-p_star(ijk)))*axz(ijk)
287  ELSE
288  IF(.NOT.cut_v_treatment_at(ijk)) THEN
289  sdps = - ((p_s(ijkn,m)-p_s(ijk,m))+(p_star(ijkn)-p_star(ijk)))*axz(ijk)
290  ELSE
291  sdps = - ((p_s(ijkn,m)* a_vpg_n(ijk)-p_s(ijk,m)* a_vpg_s(ijk)) &
292  +(p_star(ijkn)* a_vpg_n(ijk)-p_star(ijk)* a_vpg_s(ijk)))
293  ENDIF
294  ENDIF
295  ELSE
296  IF(.NOT.cut_v_treatment_at(ijk)) THEN
297  sdps = -(p_s(ijkn,m)-p_s(ijk,m))*axz(ijk)
298  ELSE
299  sdps = -(p_s(ijkn,m) * a_vpg_n(ijk)-p_s(ijk,m) * a_vpg_s(ijk))
300  ENDIF
301  ENDIF
302 
303 
304  IF(.NOT.cut_v_treatment_at(ijk)) THEN
305 ! Volumetric forces
306  ropsa = avg_y(rop_s(ijk,m),rop_s(ijkn,m),j)
307 ! Previous time step
308  v0 = avg_y(rop_so(ijk,m),rop_so(ijkn,m),j)*odt
309 ! Added mass implicit transient term {Cv eps rop_g dV/dt}
310  IF(added_mass .AND. m==m_am) THEN
311  rop_ma = avg_y(rop_g(ijk)*ep_s(ijk,m),rop_g(ijkn)*ep_s(ijkn,m),j)
312  v0 = v0 + cv * rop_ma * odt
313  ENDIF
314  ELSE
315 ! Volumetric forces
316  ropsa = (vol(ijk)*rop_s(ijk,m) +&
317  vol(ijkn)*rop_s(ijkn,m))/(vol(ijk) + vol(ijkn))
318 ! Previous time step
319  v0 = (vol(ijk)*rop_so(ijk,m) + &
320  vol(ijkn)*rop_so(ijkn,m))*odt/&
321  (vol(ijk) + vol(ijkn))
322 ! Added mass implicit transient term {Cv eps rop_g dV/dt}
323  IF(added_mass .AND. m==m_am) THEN
324  rop_ma = (vol(ijk)*rop_g(ijk)*ep_s(ijk,m) + &
325  vol(ijkn)*rop_g(ijkn)*ep_s(ijkn,m))/&
326  (vol(ijk) + vol(ijkn))
327  v0 = v0 + cv * rop_ma * odt
328  ENDIF
329  ENDIF
330 
331 ! VIRTUAL MASS SECTION (explicit terms)
332 ! adding transient term dvg/dt to virtual mass term
333  f_vir = zero
334  IF(added_mass .AND. m==m_am .AND.&
335  (.NOT.cut_v_treatment_at(ijk))) THEN
336  f_vir = ( (v_g(ijk) - v_go(ijk)) )*odt*vol_v(ijk)
337 
338 ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
339  vgs = avg_y_n(v_g(ijmk),v_g(ijk))
340  vgn = avg_y_n(v_g(ijk),v_g(ijpk))
341  uge = avg_y(u_g(ijk),u_g(ijpk),j)
342  ugw = avg_y(u_g(imjk),u_g(imjpk),j)
343  vge = avg_x(v_g(ijk),v_g(ipjk),ip1(i))
344  vgw = avg_x(v_g(imjk),v_g(ijk),i)
345  IF(do_k) THEN
346  wgt = avg_y(w_g(ijk),w_g(ijpk),j)
347  wgb = avg_y(w_g(ijkm),w_g(ijpkm),j)
348  vgt = avg_z(v_g(ijk),v_g(ijkp),kp1(k))
349  vgb = avg_z(v_g(ijkm),v_g(ijk),k)
350  f_vir = f_vir + avg_z_t(wgb,wgt)*&
351  ox(i) * (vgt - vgb)*axy(ijk)
352  ENDIF
353 ! adding convective terms (U dV/dx + V dV/dy) to virtual mass; W dV/dz added above.
354  f_vir = f_vir + v_g(ijk)*(vgn - vgs)*axz(ijk) + &
355  avg_x_e(ugw,uge,ip1(i))*(vge - vgw)*ayz(ijk)
356  f_vir = f_vir * cv * rop_ma
357  ENDIF
358 
359 ! Interphase mass transfer
360  IF (kt_type_enum == ghd_2007) THEN
361  vmttmp = zero
362  DO l = 1,smax
363  vmttmp = vmttmp + avg_y(sum_r_s(ijk,l),sum_r_s(ijkn,l),j)
364  ENDDO
365  vmt = vmttmp
366  ELSE
367  IF(.NOT.cut_v_treatment_at(ijk)) THEN
368  vmt = avg_y(sum_r_s(ijk,m),sum_r_s(ijkn,m),j)
369  ELSE
370  vmt = (vol(ijk)*sum_r_s(ijk,m) + vol(ijkn)*&
371  sum_r_s(ijkn,m))/(vol(ijk) + vol(ijkn))
372  ENDIF
373  ENDIF
374 
375 ! Body force
376  IF (model_b) THEN
377  IF (kt_type_enum == ghd_2007) THEN
378  dro1 = rop_s(ijk,m) - ro_g(ijk) *epsmix
379  dro2 = rop_s(ijkn,m) - ro_g(ijkn)*epsmixn
380  droa = avg_y(dro1,dro2,j)
381  vbf = droa*bfy_s(ijk,m)
382  ELSE
383  dro1 = (ro_s(ijk,m)-ro_g(ijk))*ep_s(ijk,m)
384  dro2 = (ro_s(ijk,m)-ro_g(ijkn))*ep_s(ijkn,m)
385  droa = avg_y(dro1,dro2,j)
386  vbf = droa*bfy_s(ijk,m)
387  ENDIF
388  ELSE
389  vbf = ropsa*bfy_s(ijk,m)
390  ENDIF
391 
392 ! Additional force for GHD from drag force sum(beta_ig * Joi/rhop_i)
393  ghd_drag = zero
394  IF (kt_type_enum == ghd_2007) THEN
395  DO l = 1,smax
396  avgrop = avg_y(rop_s(ijk,l),rop_s(ijkn,l),j)
397  if(avgrop > zero) ghd_drag = ghd_drag -&
398  avg_y(f_gs(ijk,l),f_gs(ijkn,l),j) * &
399  joiy(ijk,l) / avgrop
400  ENDDO
401  ENDIF
402 
403 ! Additional force for HYS drag force, do not use with mixture GHD theory
404  hys_drag = zero
405  avgdrag = zero
406  IF (drag_type_enum .EQ. hys .AND. kt_type_enum /= ghd_2007) THEN
407  DO l = 1,mmax
408  IF (l /= m) THEN
409  avgdrag = avg_y(beta_ij(ijk,m,l),beta_ij(ijkn,m,l),j)
410  hys_drag = hys_drag - avgdrag * (v_g(ijk) - v_s(ijk,l))
411  ENDIF
412  ENDDO
413  ENDIF
414 
415 ! Source terms from convective mom. flux (loezos)
416  IF (shear) THEN
417  srt=(2d0*v_sh/xlength)
418  vsh_p=vsh(ijk)
419  vsh_n=vsh_p
420  vsh_s=vsh_p
421  vsh_e=vsh(ijk)+srt*1d0/odx_e(i)
422  vsh_w=vsh(ijk)-srt*1d0/odx_e(im1(i))
423  source_conv=a_m(ijk,north,m)*vsh_n+a_m(ijk,south,m)*vsh_s&
424  +a_m(ijk,west,m)*vsh_w+a_m(ijk,east,m)*vsh_e&
425  -(a_m(ijk,north,m)+a_m(ijk,south,m)+a_m(ijk,west,m)+a_m(ijk,east,m))&
426  *vsh_p
427  ELSE
428  source_conv=0d0
429  ENDIF
430 
431 ! Collect the terms
432  a_m(ijk,0,m) = -(a_m(ijk,east,m)+a_m(ijk,west,m)+&
433  a_m(ijk,north,m)+a_m(ijk,south,m)+a_m(ijk,top,m)+&
434  a_m(ijk,bottom,m)+(v0+zmax(vmt))*vol_v(ijk))
435 
436  b_m(ijk,m) = b_m(ijk,m) - (sdp + sdps + &
437  tau_v_s(ijk,m) + epsa*ctau_v_g(ijk) + &
438  source_conv + f_vir + &
439  ( (v0+zmax((-vmt)))*v_so(ijk,m) + &
440  vbf + hys_drag)*vol_v(ijk) )
441 ! MMS Source term
442  IF(use_mms) b_m(ijk,m) = &
443  b_m(ijk,m) - mms_v_s_src(ijk)*vol_v(ijk)
444 
445  IF (kt_type_enum == ia_2005) THEN
446  b_m(ijk,m) = b_m(ijk,m) - ktmom_v_s(ijk,m)
447  ELSEIF (kt_type_enum == ghd_2007) THEN
448  b_m(ijk,m) = b_m(ijk,m) - ghd_drag*vol_v(ijk)
449  ENDIF
450 
451  ENDIF ! end branching on cell type (sip/ip/dilute/block/else branches)
452  ENDDO ! end do loop over ijk
453 !$omp end parallel do
454 
455 ! modifications for cartesian grid implementation
456  IF(cartesian_grid) CALL cg_source_v_s(a_m, b_m, m)
457 ! modifications for bc
458  CALL source_v_s_bc (a_m, b_m, m)
459  IF(cartesian_grid) CALL cg_source_v_s_bc(a_m, b_m, m)
460 
461  ENDIF ! end if (momentum_y_eq)
462  ENDIF ! end if for GHD Theory
463  ENDDO ! end do loop over mmax
464 
465  RETURN
466  END SUBROUTINE source_v_s
467 
468 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
469 ! C
470 ! Subroutine: SOURCE_V_s_BC C
471 ! Purpose: Determine source terms for V_s momentum eq. The terms C
472 ! appear in the center coefficient and RHS vector. The center C
473 ! coefficient and source vector are negative. The off-diagonal C
474 ! coefficients are positive. C
475 ! The drag terms are excluded from the source at this stage. C
476 ! C
477 ! Author: M. Syamlal Date: 7-JUN-96 C
478 ! Reviewer: Date: C
479 ! C
480 ! C
481 ! Comments: see source_v_g_bc for more detailed in-code comments C
482 ! C
483 ! C
484 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
485 
486  SUBROUTINE source_v_s_bc(A_M, B_M, M)
487 
488 !-----------------------------------------------
489 ! Modules
490 !-----------------------------------------------
491  USE param
492  USE param1
493  USE parallel
494  USE scales
495  USE constant
496  USE physprop
497  USE fldvar
498  USE visc_s
499  USE rxns
500  USE run
501  USE toleranc
502  USE geometry
503  USE indices
504  USE is
505  USE tau_s
506  USE bc
507  USE output
508  USE compar
509  USE fun_avg
510  USE functions
511  IMPLICIT NONE
512 !-----------------------------------------------
513 ! Dummy Arguments
514 !-----------------------------------------------
515 ! Septadiagonal matrix A_m
516  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
517 ! Vector b_m
518  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
519 ! Solids phase index
520  INTEGER, INTENT(IN) :: M
521 !-----------------------------------------------
522 ! Local Variables
523 !-----------------------------------------------
524 ! Boundary condition
525  INTEGER :: L
526 ! Indices
527  INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
528  IM, KM, IJKS, IJMK, IJPK
529 !-----------------------------------------------
530 
531 ! Setting the default boundary conditions
532  IF (do_k) THEN
533  k1 = 1
534  DO j1 = jmin3,jmax3
535  DO i1 = imin3, imax3
536  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
537  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
538  ijk = funijk(i1,j1,k1)
539  IF (ns_wall_at(ijk)) THEN
540 ! Setting the wall velocity to zero
541  a_m(ijk,east,m) = zero
542  a_m(ijk,west,m) = zero
543  a_m(ijk,north,m) = zero
544  a_m(ijk,south,m) = zero
545  a_m(ijk,top,m) = -one
546  a_m(ijk,bottom,m) = zero
547  a_m(ijk,0,m) = -one
548  b_m(ijk,m) = zero
549  ELSEIF (fs_wall_at(ijk)) THEN
550 ! Setting the wall velocity equal to the adjacent fluid velocity
551  a_m(ijk,east,m) = zero
552  a_m(ijk,west,m) = zero
553  a_m(ijk,north,m) = zero
554  a_m(ijk,south,m) = zero
555  a_m(ijk,top,m) = one
556  a_m(ijk,bottom,m) = zero
557  a_m(ijk,0,m) = -one
558  b_m(ijk,m) = zero
559  ENDIF
560  ENDDO
561  ENDDO
562 
563  k1 = kmax2
564  DO j1 = jmin3,jmax3
565  DO i1 = imin3, imax3
566  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
567  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
568  ijk = funijk(i1,j1,k1)
569  IF (ns_wall_at(ijk)) THEN
570  a_m(ijk,east,m) = zero
571  a_m(ijk,west,m) = zero
572  a_m(ijk,north,m) = zero
573  a_m(ijk,south,m) = zero
574  a_m(ijk,top,m) = zero
575  a_m(ijk,bottom,m) = -one
576  a_m(ijk,0,m) = -one
577  b_m(ijk,m) = zero
578  ELSEIF (fs_wall_at(ijk)) THEN
579  a_m(ijk,east,m) = zero
580  a_m(ijk,west,m) = zero
581  a_m(ijk,north,m) = zero
582  a_m(ijk,south,m) = zero
583  a_m(ijk,top,m) = zero
584  a_m(ijk,bottom,m) = one
585  a_m(ijk,0,m) = -one
586  b_m(ijk,m) = zero
587  ENDIF
588  ENDDO
589  ENDDO
590  ENDIF
591 
592  i1 = 1
593  DO k1 = kmin3, kmax3
594  DO j1 = jmin3, jmax3
595  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
596  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
597  ijk = funijk(i1,j1,k1)
598  IF (ns_wall_at(ijk)) THEN
599  a_m(ijk,east,m) = -one
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  ELSEIF (fs_wall_at(ijk)) THEN
608  a_m(ijk,east,m) = one
609  a_m(ijk,west,m) = zero
610  a_m(ijk,north,m) = zero
611  a_m(ijk,south,m) = zero
612  a_m(ijk,top,m) = zero
613  a_m(ijk,bottom,m) = zero
614  a_m(ijk,0,m) = -one
615  b_m(ijk,m) = zero
616  ENDIF
617  ENDDO
618  ENDDO
619 
620  i1 = imax2
621  DO k1 = kmin3, kmax3
622  DO j1 = jmin3, jmax3
623  IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
624  IF (dead_cell_at(i1,j1,k1)) cycle ! skip dead cells
625  ijk = funijk(i1,j1,k1)
626  IF (ns_wall_at(ijk)) THEN
627  a_m(ijk,east,m) = zero
628  a_m(ijk,west,m) = -one
629  a_m(ijk,north,m) = zero
630  a_m(ijk,south,m) = zero
631  a_m(ijk,top,m) = zero
632  a_m(ijk,bottom,m) = zero
633  a_m(ijk,0,m) = -one
634  b_m(ijk,m) = zero
635  ELSEIF (fs_wall_at(ijk)) THEN
636  a_m(ijk,east,m) = zero
637  a_m(ijk,west,m) = one
638  a_m(ijk,north,m) = zero
639  a_m(ijk,south,m) = zero
640  a_m(ijk,top,m) = zero
641  a_m(ijk,bottom,m) = zero
642  a_m(ijk,0,m) = -one
643  b_m(ijk,m) = zero
644  ENDIF
645  ENDDO
646  ENDDO
647 ! End setting the default boundary conditions
648 
649 ! Setting user specified boundary conditions
650  DO l = 1, dimension_bc
651  IF (bc_defined(l)) THEN
652 
653 ! Setting wall boundary conditions
654  IF (bc_type_enum(l) == no_slip_wall) THEN
655  i1 = bc_i_w(l)
656  i2 = bc_i_e(l)
657  j1 = bc_j_s(l)
658  j2 = bc_j_n(l)
659  k1 = bc_k_b(l)
660  k2 = bc_k_t(l)
661  IF (bc_jj_ps(l) == 0) THEN
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  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
670  a_m(ijk,east,m) = zero
671  a_m(ijk,west,m) = zero
672  a_m(ijk,north,m) = zero
673  a_m(ijk,south,m) = zero
674  a_m(ijk,top,m) = zero
675  a_m(ijk,bottom,m) = zero
676  a_m(ijk,0,m) = -one
677  b_m(ijk,m) = zero
678  IF (fluid_at(east_of(ijk))) THEN
679  a_m(ijk,east,m) = -one
680  ELSEIF (fluid_at(west_of(ijk))) THEN
681  a_m(ijk,west,m) = -one
682  ELSEIF (fluid_at(top_of(ijk))) THEN
683  a_m(ijk,top,m) = -one
684  ELSEIF (fluid_at(bottom_of(ijk))) THEN
685  a_m(ijk,bottom,m) = -one
686  ENDIF
687  ENDDO
688  ENDDO
689  ENDDO
690  ELSE ! Johnson and Jackson partial slip
691  CALL jj_bc_v_s (i1, i2, j1, j2, k1, k2, l, m, a_m, b_m)
692  ENDIF
693 
694  ELSEIF (bc_type_enum(l) == free_slip_wall) THEN
695  i1 = bc_i_w(l)
696  i2 = bc_i_e(l)
697  j1 = bc_j_s(l)
698  j2 = bc_j_n(l)
699  k1 = bc_k_b(l)
700  k2 = bc_k_t(l)
701  IF (bc_jj_ps(l) == 0) THEN
702  DO k = k1, k2
703  DO j = j1, j2
704  DO i = i1, i2
705  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
706  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
707  ijk = funijk(i,j,k)
708  IF (.NOT.wall_at(ijk)) cycle !skip redefined cells
709  a_m(ijk,east,m) = zero
710  a_m(ijk,west,m) = zero
711  a_m(ijk,north,m) = zero
712  a_m(ijk,south,m) = zero
713  a_m(ijk,top,m) = zero
714  a_m(ijk,bottom,m) = zero
715  a_m(ijk,0,m) = -one
716  b_m(ijk,m) = zero
717  IF (fluid_at(east_of(ijk))) THEN
718  a_m(ijk,east,m) = one
719  ELSEIF (fluid_at(west_of(ijk))) THEN
720  a_m(ijk,west,m) = one
721  ELSEIF (fluid_at(top_of(ijk))) THEN
722  a_m(ijk,top,m) = one
723  ELSEIF (fluid_at(bottom_of(ijk))) THEN
724  a_m(ijk,bottom,m) = one
725  ENDIF
726  ENDDO
727  ENDDO
728  ENDDO
729  ELSE ! Johnson and Jackson partial slip
730  CALL jj_bc_v_s (i1, i2, j1, j2, k1, k2, l, m, a_m, b_m)
731  ENDIF
732 
733  ELSEIF (bc_type_enum(l) == par_slip_wall) THEN
734  i1 = bc_i_w(l)
735  i2 = bc_i_e(l)
736  j1 = bc_j_s(l)
737  j2 = bc_j_n(l)
738  k1 = bc_k_b(l)
739  k2 = bc_k_t(l)
740  IF (bc_jj_ps(l) == 0) THEN
741  DO k = k1, k2
742  DO j = j1, j2
743  DO i = i1, i2
744  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
745  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
746  ijk = funijk(i,j,k)
747  IF (.NOT.wall_at(ijk)) cycle !skip redefined cells
748  im = im1(i)
749  km = km1(k)
750  a_m(ijk,east,m) = zero
751  a_m(ijk,west,m) = zero
752  a_m(ijk,north,m) = zero
753  a_m(ijk,south,m) = zero
754  a_m(ijk,top,m) = zero
755  a_m(ijk,bottom,m) = zero
756  a_m(ijk,0,m) = -one
757  b_m(ijk,m) = zero
758  IF (fluid_at(east_of(ijk))) THEN
759  IF (bc_hw_s(l,m) == undefined) THEN
760  a_m(ijk,east,m) = -half
761  a_m(ijk,0,m) = -half
762  b_m(ijk,m) = -bc_vw_s(l,m)
763  ELSE
764  a_m(ijk,0,m) = -(half*bc_hw_s(l,m)+odx_e(i))
765  a_m(ijk,east,m) = -(half*bc_hw_s(l,m)-odx_e(i))
766  b_m(ijk,m) = -bc_hw_s(l,m)*bc_vw_s(l,m)
767  ENDIF
768  ELSEIF (fluid_at(west_of(ijk))) THEN
769  IF (bc_hw_s(l,m) == undefined) THEN
770  a_m(ijk,west,m) = -half
771  a_m(ijk,0,m) = -half
772  b_m(ijk,m) = -bc_vw_s(l,m)
773  ELSE
774  a_m(ijk,west,m) = -(half*bc_hw_s(l,m)-odx_e(im))
775  a_m(ijk,0,m) = -(half*bc_hw_s(l,m)+odx_e(im))
776  b_m(ijk,m) = -bc_hw_s(l,m)*bc_vw_s(l,m)
777  ENDIF
778  ELSEIF (fluid_at(top_of(ijk))) THEN
779  IF (bc_hw_s(l,m) == undefined) THEN
780  a_m(ijk,top,m) = -half
781  a_m(ijk,0,m) = -half
782  b_m(ijk,m) = -bc_vw_s(l,m)
783  ELSE
784  a_m(ijk,0,m) = -(half*bc_hw_s(l,m)+odz_t(k)*ox&
785  (i))
786  a_m(ijk,top,m) = -(half*bc_hw_s(l,m)-odz_t(k)*ox&
787  (i))
788  b_m(ijk,m) = -bc_hw_s(l,m)*bc_vw_s(l,m)
789  ENDIF
790  ELSEIF (fluid_at(bottom_of(ijk))) THEN
791  IF (bc_hw_s(l,m) == undefined) THEN
792  a_m(ijk,bottom,m) = -half
793  a_m(ijk,0,m) = -half
794  b_m(ijk,m) = -bc_vw_s(l,m)
795  ELSE
796  a_m(ijk,bottom,m) = -(half*bc_hw_s(l,m)-odz_t(km)*&
797  ox(i))
798  a_m(ijk,0,m) = -(half*bc_hw_s(l,m)+odz_t(km)*&
799  ox(i))
800  b_m(ijk,m) = -bc_hw_s(l,m)*bc_vw_s(l,m)
801  ENDIF
802  ENDIF
803  ENDDO
804  ENDDO
805  ENDDO
806  ELSE ! Johnson and Jackson partial slip
807  CALL jj_bc_v_s (i1, i2, j1, j2, k1, k2, l, m, a_m, b_m)
808  ENDIF
809 
810 ! Setting flow boundary conditions
811  ELSEIF (bc_type_enum(l)==p_inflow .OR. bc_type_enum(l)==p_outflow) THEN
812  IF (bc_plane(l) == 'S') THEN
813  i1 = bc_i_w(l)
814  i2 = bc_i_e(l)
815  j1 = bc_j_s(l)
816  j2 = bc_j_n(l)
817  k1 = bc_k_b(l)
818  k2 = bc_k_t(l)
819  DO k = k1, k2
820  DO j = j1, j2
821  DO i = i1, i2
822  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
823  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
824  ijk = funijk(i,j,k)
825  a_m(ijk,east,m) = zero
826  a_m(ijk,west,m) = zero
827  a_m(ijk,north,m) = zero
828  a_m(ijk,south,m) = one
829  a_m(ijk,top,m) = zero
830  a_m(ijk,bottom,m) = zero
831  a_m(ijk,0,m) = -one
832  b_m(ijk,m) = zero
833  ENDDO
834  ENDDO
835  ENDDO
836  ENDIF
837 
838  ELSEIF (bc_type_enum(l) == outflow) THEN
839  IF (bc_plane(l) == 'S') THEN
840  i1 = bc_i_w(l)
841  i2 = bc_i_e(l)
842  j1 = bc_j_s(l)
843  j2 = bc_j_n(l)
844  k1 = bc_k_b(l)
845  k2 = bc_k_t(l)
846  DO k = k1, k2
847  DO j = j1, j2
848  DO i = i1, i2
849  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
850  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
851  ijk = funijk(i,j,k)
852  a_m(ijk,east,m) = zero
853  a_m(ijk,west,m) = zero
854  a_m(ijk,north,m) = zero
855  a_m(ijk,south,m) = one
856  a_m(ijk,top,m) = zero
857  a_m(ijk,bottom,m) = zero
858  a_m(ijk,0,m) = -one
859  b_m(ijk,m) = zero
860  ijmk = jm_of(ijk)
861  a_m(ijmk,east,m) = zero
862  a_m(ijmk,west,m) = zero
863  a_m(ijmk,north,m) = zero
864  a_m(ijmk,south,m) = one
865  a_m(ijmk,top,m) = zero
866  a_m(ijmk,bottom,m) = zero
867  a_m(ijmk,0,m) = -one
868  b_m(ijmk,m) = zero
869  END DO
870  END DO
871  END DO
872  ELSEIF (bc_plane(l) == 'N') THEN
873  i1 = bc_i_w(l)
874  i2 = bc_i_e(l)
875  j1 = bc_j_s(l)
876  j2 = bc_j_n(l)
877  k1 = bc_k_b(l)
878  k2 = bc_k_t(l)
879  DO k = k1, k2
880  DO j = j1, j2
881  DO i = i1, i2
882  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
883  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
884  ijk = funijk(i,j,k)
885  ijpk = jp_of(ijk)
886  a_m(ijpk,east,m) = zero
887  a_m(ijpk,west,m) = zero
888  a_m(ijpk,north,m) = one
889  a_m(ijpk,south,m) = zero
890  a_m(ijpk,top,m) = zero
891  a_m(ijpk,bottom,m) = zero
892  a_m(ijpk,0,m) = -one
893  b_m(ijpk,m) = zero
894  ENDDO
895  ENDDO
896  ENDDO
897  ENDIF
898 
899 ! Setting bc that are not ns, fs, psw, p_inflow, p_outflow, or outflow
900  ELSE
901  i1 = bc_i_w(l)
902  i2 = bc_i_e(l)
903  j1 = bc_j_s(l)
904  j2 = bc_j_n(l)
905  k1 = bc_k_b(l)
906  k2 = bc_k_t(l)
907  DO k = k1, k2
908  DO j = j1, j2
909  DO i = i1, i2
910  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
911  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
912  ijk = funijk(i,j,k)
913  a_m(ijk,east,m) = zero
914  a_m(ijk,west,m) = zero
915  a_m(ijk,north,m) = zero
916  a_m(ijk,south,m) = zero
917  a_m(ijk,top,m) = zero
918  a_m(ijk,bottom,m) = zero
919  a_m(ijk,0,m) = -one
920  b_m(ijk,m) = -v_s(ijk,m)
921  IF (bc_plane(l) == 'S') THEN
922  ijks = south_of(ijk)
923  a_m(ijks,east,m) = zero
924  a_m(ijks,west,m) = zero
925  a_m(ijks,north,m) = zero
926  a_m(ijks,south,m) = zero
927  a_m(ijks,top,m) = zero
928  a_m(ijks,bottom,m) = zero
929  a_m(ijks,0,m) = -one
930  b_m(ijks,m) = -v_s(ijks,m)
931  ENDIF
932  ENDDO
933  ENDDO
934  ENDDO
935 
936  ENDIF ! end if (bc_type)
937  ENDIF ! end if (bc_defined)
938  ENDDO ! end L do loop over dimension_bc
939 
940  RETURN
941  END SUBROUTINE source_v_s_bc
942 
943 
944 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
945 ! C
946 ! Module name: JJ_BC_V_s C
947 ! Purpose: Implement Johnson and Jackson boundary condition C
948 ! C
949 ! Author: K. Agrawal, A. Srivastava, Date: 14-APR-98 C
950 ! Princeton University C
951 ! Reviewer: Date: C
952 ! C
953 ! C
954 ! C
955 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
956 
957  SUBROUTINE jj_bc_v_s(I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
958 
959 !-----------------------------------------------
960 ! Modules
961 !-----------------------------------------------
962  USE bc
963  USE calc_gr_boundary
964  USE compar
965  USE constant
966  USE fldvar
967  USE functions
968  USE geometry
969  USE indices
970  USE is
971  USE output
972  USE parallel
973  USE param
974  USE param1
975  USE physprop
976  USE run
977  USE rxns
978  USE scales
979  USE tau_s
980  USE toleranc
981  USE visc_s
982  IMPLICIT NONE
983 !-----------------------------------------------
984 ! Dummy Arguments
985 !-----------------------------------------------
986 ! Boundary condition
987  INTEGER, INTENT(IN) :: L
988 ! Indices
989  INTEGER, INTENT(IN) :: I1, I2, J1, J2, K1, K2
990 ! Solids phase index
991  INTEGER, INTENT(IN) :: M
992 ! Septadiagonal matrix A_m
993  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
994 ! Vector b_m
995  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
996 !-----------------------------------------------
997 ! Local Variables
998 !-----------------------------------------------
999 ! Indices
1000  INTEGER :: I, J, K, IJK, IM, KM, IJPK
1001 ! coefficients for granular bc
1002  DOUBLE PRECISION :: hw, gw, cw
1003 !-----------------------------------------------
1004 
1005  DO k = k1, k2
1006  DO j = j1, j2
1007  DO i = i1, i2
1008  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
1009  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
1010 
1011  ijk = funijk(i,j,k)
1012  IF (.NOT.wall_at(ijk)) cycle ! skip redefined cells
1013  im = im1(i)
1014  km = km1(k)
1015  a_m(ijk,east,m) = zero
1016  a_m(ijk,west,m) = zero
1017  a_m(ijk,north,m) = zero
1018  a_m(ijk,south,m) = zero
1019  a_m(ijk,top,m) = zero
1020  a_m(ijk,bottom,m) = zero
1021  a_m(ijk,0,m) = -one
1022  b_m(ijk,m) = zero
1023 
1024  IF (fluid_at(east_of(ijk))) THEN
1025  ijpk = jp_of(east_of(ijk))
1026  IF (wall_at(ijpk)) cycle
1027  IF (ep_s(east_of(ijk),m) <= dil_ep_s) THEN
1028  a_m(ijk,east,m) = one
1029  ELSE
1030  IF (bc_jj_ps(l) == 1) THEN
1031 ! Setting the wall velocity based on Johnson and Jackson B.C which are
1032 ! also modified to include frictional effects if requested
1033  CALL calc_grbdry (ijk, east_of(ijk), 'E', 'V', &
1034  m, l, gw, hw, cw)
1035  ELSEIF (bc_jj_ps(l) == 2) THEN
1036 ! Setting the wall velocity equal to the user specified value
1037 ! (if wall velocity = 0 then it is equivalent to no slip wall)
1038  gw = 0d0
1039  hw = 1d0
1040  cw = bc_vw_s(l,m)
1041  ELSE
1042 ! Setting the wall velocity equal to the adjacent fluid cell
1043 ! velocity (i.e. zero flux/free slip)
1044  gw = 1d0
1045  cw = 0d0
1046  hw = 0d0
1047  ENDIF
1048  a_m(ijk,east,m) = -(half*hw - odx_e(i)*gw)
1049  a_m(ijk,0,m) = -(half*hw + odx_e(i)*gw)
1050  b_m(ijk,m) = -cw
1051  ENDIF
1052 
1053  ELSEIF (fluid_at(west_of(ijk))) THEN
1054  ijpk = jp_of(west_of(ijk))
1055  IF (wall_at(ijpk)) cycle
1056  IF (ep_s(west_of(ijk),m) <= dil_ep_s) THEN
1057  a_m(ijk,west,m) = one
1058  ELSE
1059  IF (bc_jj_ps(l) == 1) THEN
1060  CALL calc_grbdry (ijk, west_of(ijk), 'W', 'V',&
1061  m, l, gw, hw, cw)
1062  ELSEIF (bc_jj_ps(l) == 2) THEN
1063  gw = 0d0
1064  hw = 1d0
1065  cw = bc_vw_s(l,m)
1066  ELSE
1067  gw = 1d0
1068  cw = 0d0
1069  hw = 0d0
1070  ENDIF
1071  a_m(ijk,west,m) = -(half*hw - odx_e(im)*gw)
1072  a_m(ijk,0,m) = -(half*hw + odx_e(im)*gw)
1073  b_m(ijk,m) = -cw
1074  ENDIF
1075 
1076  ELSEIF (fluid_at(top_of(ijk))) THEN
1077  ijpk = jp_of(top_of(ijk))
1078  IF (wall_at(ijpk)) cycle
1079  IF (ep_s(top_of(ijk),m) <= dil_ep_s) THEN
1080  a_m(ijk,top,m) = one
1081  ELSE
1082  IF (bc_jj_ps(l) == 1) THEN
1083  CALL calc_grbdry (ijk, top_of(ijk), 'T', 'V',&
1084  m, l, gw, hw, cw)
1085  ELSEIF (bc_jj_ps(l) == 2) THEN
1086  gw = 0d0
1087  hw = 1d0
1088  cw = bc_vw_s(l,m)
1089  ELSE
1090  gw = 1d0
1091  cw = 0d0
1092  hw = 0d0
1093  ENDIF
1094  a_m(ijk,top,m) = -(half*hw - odz_t(k)*ox(i)*gw)
1095  a_m(ijk,0,m) = -(half*hw + odz_t(k)*ox(i)*gw)
1096  b_m(ijk,m) = -cw
1097  ENDIF
1098 
1099  ELSEIF (fluid_at(bottom_of(ijk))) THEN
1100  ijpk = jp_of(bottom_of(ijk))
1101  IF (wall_at(ijpk)) cycle
1102  IF (ep_s(bottom_of(ijk),m) <= dil_ep_s) THEN
1103  a_m(ijk,bottom,m) = one
1104  ELSE
1105  IF (bc_jj_ps(l) == 1) THEN
1106  CALL calc_grbdry (ijk, bottom_of(ijk), 'B', 'V',&
1107  m, l, gw, hw, cw)
1108  ELSEIF (bc_jj_ps(l) == 2) THEN
1109  gw = 0d0
1110  hw = 1d0
1111  cw = bc_vw_s(l,m)
1112  ELSE
1113  gw = 1d0
1114  cw = 0d0
1115  hw = 0d0
1116  ENDIF
1117  a_m(ijk,bottom,m) = -(half*hw - odz_t(km)*ox(i)*gw)
1118  a_m(ijk,0,m) = -(half*hw + odz_t(km)*ox(i)*gw)
1119  b_m(ijk,m) = -cw
1120  ENDIF
1121  ENDIF
1122 
1123  ENDDO
1124  ENDDO
1125  ENDDO
1126 
1127  RETURN
1128  END SUBROUTINE jj_bc_v_s
1129 
1130 
1131 
1132 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1133 ! C
1134 ! Subroutine: POINT_SOURCE_V_S C
1135 ! Purpose: Adds point sources to the solids V-Momentum equations. C
1136 ! C
1137 ! Author: J. Musser Date: 10-JUN-13 C
1138 ! Reviewer: Date: C
1139 ! C
1140 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1141  SUBROUTINE point_source_v_s(A_M, B_M)
1142 
1143 !-----------------------------------------------
1144 ! Modules
1145 !-----------------------------------------------
1146  use compar
1147  use constant
1148  use fldvar
1149  use geometry
1150  use indices
1151  use param
1152  use param1, only: one, small_number
1153  use physprop
1154  use ps
1155  use run
1156  use functions
1157  IMPLICIT NONE
1158 !-----------------------------------------------
1159 ! Dummy arguments
1160 !-----------------------------------------------
1161 ! Septadiagonal matrix A_m
1162  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
1163 ! Vector b_m
1164  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
1165 !-----------------------------------------------
1166 ! Local variables
1167 !-----------------------------------------------
1168 ! Indices
1169  INTEGER :: IJK, I, J, K
1170  INTEGER :: PSV, M
1171  INTEGER :: lJN, lJS
1172 ! terms of bm expression
1173  DOUBLE PRECISION :: pSource
1174 !-----------------------------------------------
1175 
1176  do m=1, mmax
1177 
1178  ps_lp: do psv = 1, dimension_ps
1179  if(.NOT.ps_defined(psv)) cycle ps_lp
1180  if(abs(ps_v_s(psv,m)) < small_number) cycle ps_lp
1181 
1182  if(ps_v_s(psv,m) < 0.0d0) then
1183  ljs = ps_j_s(psv) - 1
1184  ljn = ps_j_n(psv) - 1
1185  else
1186  ljs = ps_j_s(psv)
1187  ljn = ps_j_n(psv)
1188  endif
1189 
1190  do k = ps_k_b(psv), ps_k_t(psv)
1191  do j = ljs, ljn
1192  do i = ps_i_w(psv), ps_i_e(psv)
1193 
1194  if(.NOT.is_on_mype_plus2layers(i,j,k)) cycle
1195  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
1196 
1197  ijk = funijk(i,j,k)
1198  if(.NOT.fluid_at(ijk)) cycle
1199 
1200  if(a_m(ijk,0,m) == -one .AND. &
1201  b_m(ijk,m) == -v_s(ijk,m)) then
1202  b_m(ijk,m) = -ps_v_s(psv,m) * ps_vel_mag_s(psv,m)
1203  else
1204  psource = ps_massflow_s(psv,m) * &
1205  (vol(ijk)/ps_volume(psv))
1206 
1207  b_m(ijk,m) = b_m(ijk,m) - psource * &
1208  ps_v_s(psv,m) * ps_vel_mag_s(psv,m)
1209  endif
1210 
1211  enddo
1212  enddo
1213  enddo
1214 
1215 
1216  enddo ps_lp
1217 
1218  enddo ! do M=1, MMAX
1219 
1220  RETURN
1221  END SUBROUTINE point_source_v_s
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
geometry::imax2
integer imax2
Definition: geometry_mod.f:61
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
visc_s
Definition: visc_s_mod.f:1
ps::ps_v_s
double precision, dimension(dimension_ps, dim_m) ps_v_s
Definition: ps_mod.f:70
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
cg_source_v_s_bc
subroutine cg_source_v_s_bc(A_M, B_M, M)
Definition: CG_source_v_s.f:309
mms::mms_v_s_src
double precision, dimension(:), allocatable mms_v_s_src
Definition: mms_mod.f:63
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
cg_source_v_s
subroutine cg_source_v_s(A_M, B_M, M)
Definition: CG_source_v_s.f:16
point_source_v_s
subroutine point_source_v_s(A_M, B_M)
Definition: source_v_s.f:1142
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
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
tau_s::tau_v_s
double precision, dimension(:,:), allocatable tau_v_s
Definition: tau_s_mod.f:5
drag
Definition: drag_mod.f:11
ps::ps_vel_mag_s
double precision, dimension(dimension_ps, dim_m) ps_vel_mag_s
Definition: ps_mod.f:74
tau_s
Definition: tau_s_mod.f:1
rxns::sum_r_s
double precision, dimension(:,:), allocatable sum_r_s
Definition: rxns_mod.f:35
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
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
physprop::close_packed
logical, dimension(dim_m) close_packed
Definition: physprop_mod.f:56
source_v_s
subroutine source_v_s(A_M, B_M)
Definition: source_v_s.f:32
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
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
fun_avg
Definition: fun_avg_mod.f:1
calc_gr_boundary::calc_grbdry
subroutine calc_grbdry(IJK1, IJK2, FCELL, COM, M, L, Gw, Hw, Cw)
Definition: calc_grbdry.f:25
tau_g::ctau_v_g
double precision, dimension(:), allocatable ctau_v_g
Definition: tau_g_mod.f:11
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
source_v_s_bc
subroutine source_v_s_bc(A_M, B_M, M)
Definition: source_v_s.f:487
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
kintheory::ktmom_v_s
double precision, dimension(:,:), allocatable ktmom_v_s
Definition: kintheory_mod.f:76
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
is::is_vel_s
double precision, dimension(dimension_is, dim_m) is_vel_s
Definition: is_mod.f:88
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_s
double precision, dimension(dimension_bc, dim_m) bc_hw_s
Definition: bc_mod.f:310
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
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
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
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
fldvar::rop_so
double precision, dimension(:,:), allocatable rop_so
Definition: fldvar_mod.f:54
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
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
param::dimension_ps
integer, parameter dimension_ps
Definition: param_mod.f:65
geometry::no_k
logical no_k
Definition: geometry_mod.f:28
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
fldvar::p_s
double precision, dimension(:,:), allocatable p_s
Definition: fldvar_mod.f:123
fldvar::p_star
double precision, dimension(:), allocatable p_star
Definition: fldvar_mod.f:142
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
mms::use_mms
logical use_mms
Definition: mms_mod.f:15
param::west
integer west
Definition: param_mod.f:33
ps::ps_massflow_s
double precision, dimension(dimension_ps, dim_m) ps_massflow_s
Definition: ps_mod.f:66
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
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
jj_bc_v_s
subroutine jj_bc_v_s(I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
Definition: source_v_s.f:958
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
bc::bc_jj_ps
integer, dimension(dimension_bc) bc_jj_ps
Definition: bc_mod.f:212
cutcell::blocked_v_cell_at
logical, dimension(:), allocatable blocked_v_cell_at
Definition: cutcell_mod.f:365
calc_gr_boundary
Definition: calc_grbdry.f:1
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
bc::bc_vw_s
double precision, dimension(dimension_bc, dim_m) bc_vw_s
Definition: bc_mod.f:325
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
bodyforce::bfy_s
double precision function bfy_s(IJK, M)
Definition: bodyforce_mod.f:48
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
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