File: N:\mfix\model\conv_source_epp.f

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
614