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