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