MFIX  2016-1
conv_source_epp.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CONV_SOURCE_EPp C
4 ! Purpose: Determine convection & source terms for solids volume C
5 ! fraction correction equation. Master routine. C
6 ! C
7 ! Notes: MCP must be defined to call this routine. C
8 ! C
9 ! Author: M. Syamlal Date: 6-MAR-97 C
10 ! Reviewer: Date: C
11 ! C
12 ! C
13 ! C
14 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15 
16  SUBROUTINE conv_source_epp(A_M, B_M, B_mmax, M)
17 
18 !-----------------------------------------------
19 ! Modules
20 !-----------------------------------------------
21  USE compar
22  USE fldvar
23  USE geometry
24  USE param
25  USE param1
26  USE run
27  USE sendrecv
28  IMPLICIT NONE
29 !-----------------------------------------------
30 ! Dummy Arguments
31 !-----------------------------------------------
32 ! Septadiagonal matrix A_m
33  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
34 ! Vector b_m
35  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
36 ! Maximum term in b_m expression
37  DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(dimension_3, 0:dimension_m)
38 ! Lowest solids phase index of those solids phases that can
39 ! close packed (M=MCP)
40  INTEGER, INTENT(IN) :: M
41 !-----------------------------------------------
42 
43  IF (discretize(2) == 0) THEN ! 0 & 1 => first order upwinding
44  CALL conv_source_epp0 (a_m, b_m, b_mmax, m)
45  ELSE
46  CALL conv_source_epp1 (a_m, b_m, b_mmax, m)
47  ENDIF
48 
49  RETURN
50  END SUBROUTINE conv_source_epp
51 
52 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
53 ! C
54 ! Subroutine: CONV_SOURCE_EPp0 C
55 ! Purpose: Determine convection & source terms for solids volume C
56 ! fraction correction equation. First order upwinding. C
57 ! C
58 ! Author: M. Syamlal Date: 24-SEP-96 C
59 ! Reviewer: Date: C
60 ! C
61 ! C
62 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
63 
64  SUBROUTINE conv_source_epp0(A_M, B_M, B_MMAX, M)
65 
66 !-----------------------------------------------
67 ! Modules
68 !-----------------------------------------------
69  USE compar
70  USE constant
71  USE fldvar
72  USE functions
73  USE geometry
74  USE indices
75  USE parallel
76  USE param
77  USE param1
78  USE pgcor
79  USE physprop
80  USE pscor
81  USE run
82  USE rxns
83  USE sendrecv
84  USE solids_pressure
85  IMPLICIT NONE
86 !-----------------------------------------------
87 ! Dummy arguments
88 !-----------------------------------------------
89 ! Septadiagonal matrix A_m
90  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
91 ! Vector b_m
92  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
93 ! maximum term in b_m expression
94  DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(dimension_3, 0:dimension_m)
95 ! Lowest solids phase index of those solids phases that can
96 ! close packed (M=MCP)
97  INTEGER, INTENT(IN) :: M
98 !-----------------------------------------------
99 ! Local variables
100 !-----------------------------------------------
101 ! Indices
102  INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
103  INTEGER :: IMJK, IJMK, IJKM, IJKE, IJKW, IJKN, IJKS
104  INTEGER :: IJKB, IJKT
105 ! dPodEP_s(EP_s(IJK, M)); the derivative of the plastic or friction
106 ! pressure with respect to ep_s
107  DOUBLE PRECISION :: K_P
108 ! Mass source
109  DOUBLE PRECISION :: Src
110 ! error message
111  CHARACTER(LEN=80) :: LINE(1)
112 ! terms of bm expression
113  DOUBLE PRECISION :: bma, bme, bmw, bmn, bms, bmt, bmb, bmr
114 !-----------------------------------------------
115 
116 !!$omp parallel do &
117 !!$omp& private(I, J, K, IJK, IPJK, IJPK, IJKP, &
118 !!$omp& IMJK, IJMK, IJKM, &
119 !!$omp& IJKE, IJKW, IJKN, IJKS, IJKT, IJKB, &
120 !!$omp& K_P, SRC, bma, bme, bmw, bmn, bms, bmt, bmb, bmr )
121  DO ijk = ijkstart3, ijkend3
122 ! Determine whether IJK falls within 1 ghost layer........
123  i = i_of(ijk)
124  j = j_of(ijk)
125  k = k_of(ijk)
126 
127  IF (fluid_at(ijk)) THEN
128  ipjk = ip_of(ijk)
129  ijpk = jp_of(ijk)
130  ijkp = kp_of(ijk)
131  imjk = im_of(ijk)
132  ijmk = jm_of(ijk)
133  ijkm = km_of(ijk)
134 
135  ijke = east_of(ijk)
136  ijkw = west_of(ijk)
137  ijkn = north_of(ijk)
138  ijks = south_of(ijk)
139  ijkt = top_of(ijk)
140  ijkb = bottom_of(ijk)
141 
142 ! initializing
143  a_m(ijk,0,0) = zero
144  b_m(ijk,0) = zero
145  k_p = k_cp(ijk)
146 
147 ! Calculate convection-diffusion fluxes through each of the faces
148 
149 ! East face (i+1/2, j, k)
150  IF (u_s(ijk,m) < zero) THEN
151  a_m(ijk,east,0) = (rop_s(ijke,m)*e_e(ijk)*k_cp(ijke)-&
152  ro_s(ijke,m)*u_s(ijk,m))*ayz(ijk)
153  a_m(ijk,0,0) = a_m(ijk,0,0)+&
154  rop_s(ijke,m)*e_e(ijk)*k_p*ayz(ijk)
155  bme = (-rop_s(ijke,m)*u_s(ijk,m))*ayz(ijk)
156  b_m(ijk,0) = b_m(ijk,0) + bme
157  ELSE
158  a_m(ijk,east,0) = (rop_s(ijk,m)*&
159  e_e(ijk)*k_cp(ijke))*ayz(ijk)
160  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_s(ijk,m)*&
161  e_e(ijk)*k_p+ro_s(ijk,m)*u_s(ijk,m))*ayz(ijk)
162  bme = (-rop_s(ijk,m)*u_s(ijk,m))*ayz(ijk)
163  b_m(ijk,0) = b_m(ijk,0) + bme
164  ENDIF
165 
166 
167 ! West face (i-1/2, j, k)
168  IF (u_s(imjk,m) > zero) THEN
169  a_m(ijk,west,0) = (rop_s(ijkw,m)*e_e(imjk)*k_cp(ijkw)+&
170  ro_s(ijkw,m)*u_s(imjk,m))*ayz(imjk)
171  a_m(ijk,0,0) = a_m(ijk,0,0) + &
172  rop_s(ijkw,m)*e_e(imjk)*k_p*ayz(imjk)
173  bmw = (rop_s(ijkw,m)*u_s(imjk,m))*ayz(imjk)
174  b_m(ijk,0) = b_m(ijk,0) + bmw
175  ELSE
176  a_m(ijk,west,0) = (rop_s(ijk,m)*&
177  e_e(imjk)*k_cp(ijkw))*ayz(imjk)
178  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_s(ijk,m)*&
179  e_e(imjk)*k_p-ro_s(ijk,m)*u_s(imjk,m))*ayz(imjk)
180  bmw = (rop_s(ijk,m)*u_s(imjk,m))*ayz(imjk)
181  b_m(ijk,0) = b_m(ijk,0) + bmw
182  ENDIF
183 
184 
185 ! North face (i, j+1/2, k)
186  IF (v_s(ijk,m) < zero) THEN
187  a_m(ijk,north,0) = (rop_s(ijkn,m)*e_n(ijk)*k_cp(ijkn)-&
188  ro_s(ijkn,m)*v_s(ijk,m))*axz(ijk)
189  a_m(ijk,0,0)=a_m(ijk,0,0)+&
190  rop_s(ijkn,m)*e_n(ijk)*k_p*axz(ijk)
191  bmn = (-rop_s(ijkn,m)*v_s(ijk,m))*axz(ijk)
192  b_m(ijk,0) = b_m(ijk,0) + bmn
193  ELSE
194  a_m(ijk,north,0) = (rop_s(ijk,m)*&
195  e_n(ijk)*k_cp(ijkn))*axz(ijk)
196  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_s(ijk,m)*&
197  e_n(ijk)*k_p+ro_s(ijk,m)*v_s(ijk,m))*axz(ijk)
198  bmn = (-rop_s(ijk,m)*v_s(ijk,m))*axz(ijk)
199  b_m(ijk,0) = b_m(ijk,0) + bmn
200  ENDIF
201 
202 
203 ! South face (i, j-1/2, k)
204  IF (v_s(ijmk,m) > zero) THEN
205  a_m(ijk,south,0) = (rop_s(ijks,m)*e_n(ijmk)*k_cp(ijks)+&
206  ro_s(ijks,m)*v_s(ijmk,m))*axz(ijmk)
207  a_m(ijk,0,0) = a_m(ijk,0,0) + &
208  rop_s(ijks,m)*e_n(ijmk)*k_p*axz(ijmk)
209  bms = (rop_s(ijks,m)*v_s(ijmk,m))*axz(ijmk)
210  b_m(ijk,0) = b_m(ijk,0) + bms
211  ELSE
212  a_m(ijk,south,0) = (rop_s(ijk,m)*&
213  e_n(ijmk)*k_cp(ijks))*axz(ijmk)
214  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_s(ijk,m)*&
215  e_n(ijmk)*k_p-ro_s(ijk,m)*v_s(ijmk,m))*axz(ijmk)
216  bms = (rop_s(ijk,m)*v_s(ijmk,m))*axz(ijmk)
217  b_m(ijk,0) = b_m(ijk,0) + bms
218  ENDIF
219 
220 
221  IF (do_k) THEN
222 ! Top face (i, j, k+1/2)
223  IF (w_s(ijk,m) < zero) THEN
224  a_m(ijk,top,0) = (rop_s(ijkt,m)*e_t(ijk)*k_cp(ijkt)-&
225  ro_s(ijkt,m)*w_s(ijk,m))*axy(ijk)
226  a_m(ijk,0,0) = a_m(ijk,0,0) + &
227  rop_s(ijkt,m)*e_t(ijk)*k_p*axy(ijk)
228  bmt = (-rop_s(ijkt,m)*w_s(ijk,m))*axy(ijk)
229  b_m(ijk,0)=b_m(ijk,0) + bmt
230  ELSE
231  a_m(ijk,top,0) = (rop_s(ijk,m)*&
232  e_t(ijk)*k_cp(ijkt))*axy(ijk)
233  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_s(ijk,m)*&
234  e_t(ijk)*k_p+ro_s(ijk,m)*w_s(ijk,m))*axy(ijk)
235  bmt = (-rop_s(ijk,m)*w_s(ijk,m))*axy(ijk)
236  b_m(ijk,0) = b_m(ijk,0) + bmt
237  ENDIF
238 
239 
240 ! Bottom face (i, j, k-1/2)
241  IF (w_s(ijkm,m) > zero) THEN
242  a_m(ijk,bottom,0) = (rop_s(ijkb,m)*e_t(ijkm)*k_cp(ijkb)+&
243  ro_s(ijkb,m)*w_s(ijkm,m))*axy(ijkm)
244  a_m(ijk,0,0) = a_m(ijk,0,0) + &
245  rop_s(ijkb,m)*e_t(ijkm)*k_p*axy(ijkm)
246  bmb = (rop_s(ijkb,m)*w_s(ijkm,m))*axy(ijkm)
247  b_m(ijk,0) = b_m(ijk,0) + bmb
248  ELSE
249  a_m(ijk,bottom,0) = (rop_s(ijk,m)*&
250  e_t(ijkm)*k_cp(ijkb))*axy(ijkm)
251  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_s(ijk,m)*&
252  e_t(ijkm)*k_p-ro_s(ijk,m)*w_s(ijkm,m))*axy(ijkm)
253  bmb = (rop_s(ijk,m)*w_s(ijkm,m))*axy(ijkm)
254  b_m(ijk,0)=b_m(ijk,0) + bmb
255  ENDIF
256 
257  ELSE ! not(DO_K) branch
258  bmt = zero
259  bmb = zero
260  ENDIF ! end if/else (do_K)
261 
262  IF (rop_s(ijk,m)>zero .AND. sum_r_s(ijk,m)<zero) THEN
263  src = vol(ijk)*(-sum_r_s(ijk,m))/rop_s(ijk,m)
264  ELSE
265  src = zero
266  ENDIF
267 
268  a_m(ijk,0,0) = -(a_m(ijk,0,0)+vol(ijk)*odt*ro_s(ijk,m)+src*ro_s(ijk,m))
269 
270  bma = (rop_s(ijk,m)-rop_so(ijk,m))*vol(ijk)*odt
271  bmr = sum_r_s(ijk,m)*vol(ijk)
272  b_m(ijk,0) = -(b_m(ijk,0) - bma + bmr)
273  b_mmax(ijk,0) = max(abs(bma), abs(bme), abs(bmw), abs(bmn), abs(bms), abs(bmt), abs(bmb), abs(bmr) )
274 
275  IF ((-a_m(ijk,0,0)) < small_number) THEN
276  IF (abs(b_m(ijk,0)) < small_number) THEN
277  a_m(ijk,0,0) = -one ! Equation is undefined.
278  b_m(ijk,0) = zero ! Use existing value
279  ELSE
280 !!$omp critical
281  WRITE (line, '(A,I6,A,G12.5)') &
282  'Error: At IJK = ', ijk, &
283  ' A = 0 and b = ', b_m(ijk,0)
284  CALL write_error ('CONV_SOURCE_EPp0', line, 1)
285 !!$omp end critical
286  ENDIF
287  ENDIF
288  ELSE
289  a_m(ijk,0,0) = -one
290  b_m(ijk,0) = zero
291  ENDIF
292  ENDDO
293 
294  RETURN
295  END SUBROUTINE conv_source_epp0
296 
297 
298 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
299 ! C
300 ! Subroutine: CONV_SOURCE_EPp1 C
301 ! Purpose: Determine convection & source terms for solids volume C
302 ! fraction correction equation. Higher order scheme. C
303 ! C
304 ! Author: M. Syamlal Date: 24-SEP-96 C
305 ! Reviewer: Date: C
306 ! C
307 ! C
308 ! C
309 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
310 
311  SUBROUTINE conv_source_epp1(A_M, B_M, B_MMAX, M)
312 
313 !-----------------------------------------------
314 ! Modules
315 !-----------------------------------------------
316  USE compar
317  USE constant
318  USE fldvar
319  USE functions
320  USE geometry
321  USE indices
322  USE parallel
323  USE param
324  USE param1
325  USE pgcor
326  USE physprop
327  USE pscor
328  USE run
329  USE rxns
330  USE sendrecv
331  USE solids_pressure
332  USE vshear
333  USE xsi
334  IMPLICIT NONE
335 !-----------------------------------------------
336 ! Dummy arguments
337 !-----------------------------------------------
338 ! Septadiagonal matrix A_m
339  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
340 ! Vector b_m
341  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
342 ! maximum term in b_m expression
343  DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(dimension_3, 0:dimension_m)
344 ! Lowest solids phase index of those solids phases that can
345 ! close packed (M=MCP)
346  INTEGER, INTENT(IN) :: M
347 
348 !-----------------------------------------------
349 ! Local variables
350 !-----------------------------------------------
351 ! Indices
352  INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
353  INTEGER :: IMJK, IJMK, IJKM, IJKE, IJKW, IJKN, IJKS
354  INTEGER :: IJKB, IJKT
355 ! loezos : used for including shearing
356  INTEGER :: incr
357 ! dPodEP_s(EP_s(IJK, M)); the derivative of the plastic or friction
358 ! pressure with respect to ep_s
359  DOUBLE PRECISION :: K_P
360 ! Mass source
361  DOUBLE PRECISION :: Src
362 ! face value of ROP_s
363  DOUBLE PRECISION :: ROP_sf
364 ! terms of bm expression
365  DOUBLE PRECISION :: bma, bme, bmw, bmn, bms, bmt, bmb, bmr
366 ! error message
367  CHARACTER(LEN=80) :: LINE(1)
368 ! temporary use of global arrays:
369 ! xsi_array: convection weighting factors
370  DOUBLE PRECISION :: XSI_e(dimension_3), &
371  XSI_n(dimension_3),&
372  XSI_t(dimension_3)
373 !-----------------------------------------------
374 
375 ! Loezos:
376  incr=0
377 
378 ! Calculate convection factors
379  CALL calc_xsi (discretize(2), rop_s(1,m), u_s(1,m), v_s(1,m), w_s(1,m), &
380  xsi_e, xsi_n, xsi_t,incr)
381 
382 ! Loezos:
383 ! update to true velocity
384  IF (shear) THEN
385 !!$omp parallel do private(IJK)
386  DO ijk = ijkstart3, ijkend3
387  IF (fluid_at(ijk)) THEN
388  v_s(ijk,m)=v_s(ijk,m)+vsh(ijk)
389  ENDIF
390  ENDDO
391  ENDIF
392 
393 !!$omp parallel do &
394 !!$omp& private(I, J, K, IJK, IPJK, IJPK, IJKP, &
395 !!$omp& IMJK, IJMK, IJKM, IJKE, IJKW, IJKN, IJKS, IJKT, IJKB, &
396 !!$omp& K_P,ROP_SF,SRC, bma, bme, bmw, bmn, bms, bmt, bmb, bmr )
397  DO ijk = ijkstart3, ijkend3
398 ! Determine if IJK falls within 1 ghost layer........
399  i = i_of(ijk)
400  j = j_of(ijk)
401  k = k_of(ijk)
402 
403  IF (fluid_at(ijk)) THEN
404  ipjk = ip_of(ijk)
405  ijpk = jp_of(ijk)
406  ijkp = kp_of(ijk)
407  imjk = im_of(ijk)
408  ijmk = jm_of(ijk)
409  ijkm = km_of(ijk)
410 
411  ijke = east_of(ijk)
412  ijkw = west_of(ijk)
413  ijkn = north_of(ijk)
414  ijks = south_of(ijk)
415  ijkt = top_of(ijk)
416  ijkb = bottom_of(ijk)
417 
418 ! initializing
419  a_m(ijk,0,0) = zero
420  b_m(ijk,0) = zero
421  k_p = k_cp(ijk)
422 
423 ! Calculate convection-diffusion fluxes through each of the faces
424 
425 ! East face (i+1/2, j, k)
426  rop_sf = rop_s(ijke,m)*xsi_e(ijk) + rop_s(ijk,m)*(one - xsi_e(ijk))
427  a_m(ijk,east,0) = (rop_sf*e_e(ijk)*k_cp(ijke)-ro_s(ijk,m)*u_s(ijk,m)*xsi_e&
428  (ijk))*ayz(ijk)
429  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_sf*e_e(ijk)*k_p+ro_s(ijk,m)*u_s(ijk,&
430  m)*(one-xsi_e(ijk)))*ayz(ijk)
431  bme = (-rop_sf*u_s(ijk,m))*ayz(ijk)
432  b_m(ijk,0) = b_m(ijk,0) + bme
433 
434 ! West face (i-1/2, j, k)
435  rop_sf=rop_s(ijk,m)*xsi_e(imjk)+rop_s(ijkw,m)*(one-xsi_e(imjk))
436  a_m(ijk,west,0) = (rop_sf*e_e(imjk)*k_cp(ijkw)+ro_s(ijk,m)*u_s(imjk,m)*(&
437  one-xsi_e(imjk)))*ayz(imjk)
438  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_sf*e_e(imjk)*k_p-ro_s(ijk,m)*u_s(&
439  imjk,m)*xsi_e(imjk))*ayz(imjk)
440  bmw = (rop_sf*u_s(imjk,m))*ayz(imjk)
441  b_m(ijk,0) = b_m(ijk,0) + bmw
442 
443 
444 ! North face (i, j+1/2, k)
445  rop_sf = rop_s(ijkn,m)*xsi_n(ijk) + rop_s(ijk,m)*(one - xsi_n(ijk))
446  a_m(ijk,north,0) = (rop_sf*e_n(ijk)*k_cp(ijkn)-ro_s(ijk,m)*v_s(ijk,m)*xsi_n&
447  (ijk))*axz(ijk)
448  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_sf*e_n(ijk)*k_p+ro_s(ijk,m)*v_s(ijk,&
449  m)*(one-xsi_n(ijk)))*axz(ijk)
450  bmn = (-rop_sf*v_s(ijk,m))*axz(ijk)
451  b_m(ijk,0) = b_m(ijk,0) + bmn
452 
453 
454 ! South face (i, j-1/2, k)
455  rop_sf=rop_s(ijk,m)*xsi_n(ijmk)+rop_s(ijks,m)*(one-xsi_n(ijmk))
456  a_m(ijk,south,0) = (rop_sf*e_n(ijmk)*k_cp(ijks)+ro_s(ijk,m)*v_s(ijmk,m)*(&
457  one-xsi_n(ijmk)))*axz(ijmk)
458  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_sf*e_n(ijmk)*k_p-ro_s(ijk,m)*v_s(&
459  ijmk,m)*xsi_n(ijmk))*axz(ijmk)
460  bms = (rop_sf*v_s(ijmk,m))*axz(ijmk)
461  b_m(ijk,0) = b_m(ijk,0) + bms
462 
463 
464  IF (do_k) THEN
465 ! Top face (i, j, k+1/2)
466  rop_sf=rop_s(ijkt,m)*xsi_t(ijk)+rop_s(ijk,m)*(one-xsi_t(ijk))
467  a_m(ijk,top,0) = (rop_sf*e_t(ijk)*k_cp(ijkt)-ro_s(ijk,m)*w_s(ijk,m)*&
468  xsi_t(ijk))*axy(ijk)
469  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_sf*e_t(ijk)*k_p+ro_s(ijk,m)*w_s(&
470  ijk,m)*(one-xsi_t(ijk)))*axy(ijk)
471  bmt = (-rop_sf*w_s(ijk,m))*axy(ijk)
472  b_m(ijk,0) = b_m(ijk,0) + bmt
473 
474 ! Bottom face (i, j, k-1/2)
475  rop_sf = rop_s(ijk,m)*xsi_t(ijkm) + rop_s(ijkb,m)*(one - xsi_t(&
476  ijkm))
477  a_m(ijk,bottom,0) = (rop_sf*e_t(ijkm)*k_cp(ijkb)+ro_s(ijk,m)*w_s(ijkm,m)*&
478  (one-xsi_t(ijkm)))*axy(ijkm)
479  a_m(ijk,0,0) = a_m(ijk,0,0) + (rop_sf*e_t(ijkm)*k_p-ro_s(ijk,m)*w_s(&
480  ijkm,m)*xsi_t(ijkm))*axy(ijkm)
481  bmb = (rop_sf*w_s(ijkm,m))*axy(ijkm)
482  b_m(ijk,0) = b_m(ijk,0) + bmb
483  ELSE ! not(do_k) branch
484  bmt = zero
485  bmb = zero
486  ENDIF ! end if/else (do_k)
487 
488  IF (rop_s(ijk,m)>zero .AND. sum_r_s(ijk,m)<zero) THEN
489  src = vol(ijk)*(-sum_r_s(ijk,m))/rop_s(ijk,m)
490  ELSE
491  src = zero
492  ENDIF
493 
494  a_m(ijk,0,0) = -(a_m(ijk,0,0)+vol(ijk)*odt*ro_s(ijk,m)+src*ro_s(ijk,m))
495 
496  bma = (rop_s(ijk,m)-rop_so(ijk,m))*vol(ijk)*odt
497  bmr = sum_r_s(ijk,m)*vol(ijk)
498  b_m(ijk,0) = -(b_m(ijk,0)- bma + bmr)
499  b_mmax(ijk,0) = max(abs(bma), abs(bme), abs(bmw), abs(bmn), abs(bms), abs(bmt), abs(bmb), abs(bmr) ) !
500 
501  IF (abs(a_m(ijk,0,0)) < small_number) THEN
502  IF (abs(b_m(ijk,0)) < small_number) THEN
503  a_m(ijk,0,0) = -one ! Equation is undefined.
504  b_m(ijk,0) = zero ! Use existing value
505  ELSE
506 !!$omp critical
507  WRITE (line(1), '(A,I6,A,G12.5)') &
508  'Error: At IJK = ', ijk, &
509  ' A = 0 and b = ', b_m(ijk,0)
510 ! Having problem to compile this statement on SGI
511  CALL write_error ('CONV_SOURCE_EPp1', line, 1)
512 !!$omp end critical
513  ENDIF
514  ENDIF
515  ELSE
516  a_m(ijk,0,0) = -one
517  b_m(ijk,0) = zero
518  ENDIF
519  ENDDO
520 
521 ! Loezos:
522  IF (shear) THEN
523 !!$omp parallel do private(IJK)
524  DO ijk = ijkstart3, ijkend3
525  IF (fluid_at(ijk)) THEN
526  v_s(ijk,m)=v_s(ijk,m)-vsh(ijk)
527  ENDIF
528  ENDDO
529  ENDIF
530 
531  RETURN
532  END SUBROUTINE conv_source_epp1
533 
534 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
535 ! C
536 ! Subroutine: POINT_SOURCE_EPP C
537 ! Purpose: Adds point sources to the solids volume fraction C
538 ! correction equation. C
539 ! C
540 ! Notes: The off-diagonal coefficients are positive. The center C
541 ! coefficient and the source vector are negative. See C
542 ! conv_Pp_g C
543 ! C
544 ! Author: J. Musser Date: 10-JUN-13 C
545 ! Reviewer: Date: C
546 ! C
547 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
548  SUBROUTINE point_source_epp(B_M, B_MMAX, M)
549 
550 !-----------------------------------------------
551 ! Modules
552 !-----------------------------------------------
553  use compar
554  use constant
555  use geometry
556  use indices
557  use physprop
558  use ps
559  use pscor
560  use run
561  use functions
562 
563  IMPLICIT NONE
564 !-----------------------------------------------
565 ! Dummy arguments
566 !-----------------------------------------------
567 ! Vector b_m
568  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
569 ! maximum term in b_m expression
570  DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(dimension_3, 0:dimension_m)
571 ! Lowest solids phase index of those solids phases that can
572 ! close packed (M=MCP)
573  INTEGER, INTENT(IN) :: M
574 !-----------------------------------------------
575 ! Local Variables
576 !-----------------------------------------------
577 ! Indices
578  INTEGER :: IJK, I, J, K
579  INTEGER :: PSV
580 
581 ! terms of bm expression
582  DOUBLE PRECISION :: pSource
583 
584 !-----------------------------------------------
585 
586  ps_lp: do psv = 1, dimension_ps
587 
588  if(.NOT.ps_defined(psv)) cycle ps_lp
589  if(ps_massflow_s(psv,m) < small_number) cycle ps_lp
590 
591  do k = ps_k_b(psv), ps_k_t(psv)
592  do j = ps_j_s(psv), ps_j_n(psv)
593  do i = ps_i_w(psv), ps_i_e(psv)
594 
595  if(.NOT.is_on_mype_plus2layers(i,j,k)) cycle
596  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
597  ijk = funijk(i,j,k)
598  if(fluid_at(ijk)) then
599 
600  psource = ps_massflow_s(psv,m) * &
601  (vol(ijk)/ps_volume(psv))
602 
603  b_m(ijk,0) = b_m(ijk,0) - psource
604  b_mmax(ijk,0) = max(abs(b_mmax(ijk,0)), abs(b_m(ijk,0)))
605  endif
606  enddo
607  enddo
608  enddo
609 
610  enddo ps_lp
611 
612  RETURN
613  END SUBROUTINE point_source_epp
ps::ps_i_w
integer, dimension(dimension_ps) ps_i_w
Definition: ps_mod.f:40
conv_source_epp1
subroutine conv_source_epp1(A_M, B_M, B_MMAX, M)
Definition: conv_source_epp.f:312
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
pscor::e_n
double precision, dimension(:), allocatable e_n
Definition: pscor_mod.f:17
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
run::shear
logical shear
Definition: run_mod.f:175
pgcor
Definition: pgcor_mod.f:1
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
write_error
subroutine write_error(NAME, LINE, LMAX)
Definition: write_error.f:21
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
conv_source_epp
subroutine conv_source_epp(A_M, B_M, B_mmax, M)
Definition: conv_source_epp.f:17
rxns
Definition: rxns_mod.f:1
conv_source_epp0
subroutine conv_source_epp0(A_M, B_M, B_MMAX, M)
Definition: conv_source_epp.f:65
rxns::sum_r_s
double precision, dimension(:,:), allocatable sum_r_s
Definition: rxns_mod.f:35
pscor::e_t
double precision, dimension(:), allocatable e_t
Definition: pscor_mod.f:19
ps::ps_j_n
integer, dimension(dimension_ps) ps_j_n
Definition: ps_mod.f:43
xsi::calc_xsi
subroutine calc_xsi(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T, incr)
Definition: xsi_mod.f:23
param::east
integer east
Definition: param_mod.f:29
geometry::ayz
double precision, dimension(:), allocatable ayz
Definition: geometry_mod.f:206
ps::ps_defined
logical, dimension(dimension_ps) ps_defined
Definition: ps_mod.f:29
xsi
Definition: xsi_mod.f:15
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
vshear::vsh
double precision, dimension(:), allocatable vsh
Definition: vshear_mod.f:3
ps::ps_volume
double precision, dimension(dimension_ps) ps_volume
Definition: ps_mod.f:84
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
ps::ps_k_b
integer, dimension(dimension_ps) ps_k_b
Definition: ps_mod.f:44
solids_pressure
Definition: solids_pressure_mod.f:25
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
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
run::odt
double precision odt
Definition: run_mod.f:54
param::north
integer north
Definition: param_mod.f:37
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
pscor
Definition: pscor_mod.f:1
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
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
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
param::top
integer top
Definition: param_mod.f:45
run::discretize
integer, dimension(dim_eqs) discretize
Definition: run_mod.f:67
ps
Definition: ps_mod.f:22
fldvar::rop_s
double precision, dimension(:,:), allocatable rop_s
Definition: fldvar_mod.f:51
pscor::k_cp
double precision, dimension(:), allocatable k_cp
Definition: pscor_mod.f:22
pscor::e_e
double precision, dimension(:), allocatable e_e
Definition: pscor_mod.f:15
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
ps::ps_i_e
integer, dimension(dimension_ps) ps_i_e
Definition: ps_mod.f:41
param::bottom
integer bottom
Definition: param_mod.f:49
point_source_epp
subroutine point_source_epp(B_M, B_MMAX, M)
Definition: conv_source_epp.f:549
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27