File: RELATIVE:/../../../mfix.git/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)
17     
18     !-----------------------------------------------
19     ! Modules
20     !-----------------------------------------------
21           USE param
22           USE param1
23           USE fldvar
24           USE run
25           USE geometry
26           USE compar
27           USE sendrecv
28           Use xsi_array
29           IMPLICIT NONE
30     !-----------------------------------------------
31     ! Dummy Arguments
32     !-----------------------------------------------
33     ! Septadiagonal matrix A_m
34           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
35     ! Vector b_m
36           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
37     ! Maximum term in b_m expression
38           DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M)
39     !-----------------------------------------------
40     
41           IF (DISCRETIZE(2) == 0) THEN               ! 0 & 1 => first order upwinding
42              CALL CONV_SOURCE_EPP0 (A_M, B_M, B_MMAX)
43           ELSE
44              CALL CONV_SOURCE_EPP1 (A_M, B_M, B_MMAX)
45           ENDIF
46     
47           RETURN
48           END SUBROUTINE CONV_SOURCE_EPP
49     
50     
51     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
52     !                                                                      C
53     !  Subroutine: CONV_SOURCE_EPp0                                        C
54     !  Purpose: Determine convection & source terms for solids volume      C
55     !           fraction correction equation.  First order upwinding.      C
56     !                                                                      C
57     !  Author: M. Syamlal                                 Date: 24-SEP-96  C
58     !  Reviewer:                                          Date:            C
59     !                                                                      C
60     !                                                                      C
61     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
62     
63           SUBROUTINE CONV_SOURCE_EPP0(A_M, B_M, B_MMAX)
64     
65     !-----------------------------------------------
66     ! Modules
67     !-----------------------------------------------
68           USE param
69           USE param1
70           USE fldvar
71           USE run
72           USE parallel
73           USE matrix
74           USE constant
75           USE physprop
76           USE rxns
77           USE geometry
78           USE indices
79           USE pgcor
80           USE pscor
81           USE compar
82           USE sendrecv
83           USE solids_pressure
84           USE functions
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     !-----------------------------------------------
96     ! Local variables
97     !-----------------------------------------------
98     ! solids phase index locally assigned to mcp
99           INTEGER :: M
100     ! Indices
101           INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
102           INTEGER :: IMJK, IJMK, IJKM, IJKE, IJKW, IJKN, IJKS
103           INTEGER :: IJKB, IJKT
104     ! dPodEP_s(EP_s(IJK, M)); the derivative of the plastic or friction
105     ! pressure with respect to ep_s
106           DOUBLE PRECISION :: K_P
107     ! Mass source
108           DOUBLE PRECISION :: Src
109     ! error message
110           CHARACTER(LEN=80) :: LINE(1)
111     ! terms of bm expression
112           DOUBLE PRECISION :: bma, bme, bmw, bmn, bms, bmt, bmb, bmr
113     !-----------------------------------------------
114     
115           IF (MCP == UNDEFINED_I) THEN
116     ! this error should be caught earlier in the routines so that this
117     ! branch should never be entered
118              RETURN
119           ELSE
120     ! the lowest solids phase index of those solids phases that can close
121     ! pack (i.e. close_packed=T) and the index of the solids phase that is
122     ! used to form the solids correction equation.
123              M = MCP
124           ENDIF
125     
126     !!$omp    parallel do                                            &
127     !!$omp&   private(I, J, K, IJK, IPJK, IJPK, IJKP,                &
128     !!$omp&           IMJK, IJMK, IJKM,                              &
129     !!$omp&           IJKE, IJKW, IJKN, IJKS, IJKT, IJKB,            &
130     !!$omp&           K_P, SRC, bma, bme, bmw, bmn, bms, bmt, bmb, bmr )
131           DO IJK = ijkstart3, ijkend3
132     ! Determine whether IJK falls within 1 ghost layer........
133              I = I_OF(IJK)
134              J = J_OF(IJK)
135              K = K_OF(IJK)
136     
137              IF (FLUID_AT(IJK)) THEN
138                 IPJK = IP_OF(IJK)
139                 IJPK = JP_OF(IJK)
140                 IJKP = KP_OF(IJK)
141                 IMJK = IM_OF(IJK)
142                 IJMK = JM_OF(IJK)
143                 IJKM = KM_OF(IJK)
144     
145                 IJKE = EAST_OF(IJK)
146                 IJKW = WEST_OF(IJK)
147                 IJKN = NORTH_OF(IJK)
148                 IJKS = SOUTH_OF(IJK)
149                 IJKT = TOP_OF(IJK)
150                 IJKB = BOTTOM_OF(IJK)
151     
152     ! initializing
153                 A_M(IJK,0,0) = ZERO
154                 B_M(IJK,0) = ZERO
155                 K_P = K_CP(IJK)
156     
157     ! Calculate convection-diffusion fluxes through each of the faces
158     
159     ! East face (i+1/2, j, k)
160                 IF (U_S(IJK,M) < ZERO) THEN
161                    A_M(IJK,E,0) = (ROP_S(IJKE,M)*E_E(IJK)*K_CP(IJKE)-&
162                      RO_S(IJKE,M)*U_S(IJK,M))*AYZ(IJK)
163                    A_M(IJK,0,0) = A_M(IJK,0,0)+&
164                       ROP_S(IJKE,M)*E_E(IJK)*K_P*AYZ(IJK)
165                    bme = (-ROP_S(IJKE,M)*U_S(IJK,M))*AYZ(IJK)
166                    B_M(IJK,0) = B_M(IJK,0) +  bme
167                 ELSE
168                    A_M(IJK,E,0) = (ROP_S(IJK,M)*&
169                       E_E(IJK)*K_CP(IJKE))*AYZ(IJK)
170                    A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_S(IJK,M)*&
171                       E_E(IJK)*K_P+RO_S(IJK,M)*U_S(IJK,M))*AYZ(IJK)
172                    bme = (-ROP_S(IJK,M)*U_S(IJK,M))*AYZ(IJK)
173                    B_M(IJK,0) = B_M(IJK,0) +  bme
174                 ENDIF
175     
176     
177     ! West face (i-1/2, j, k)
178                 IF (U_S(IMJK,M) > ZERO) THEN
179                    A_M(IJK,W,0) = (ROP_S(IJKW,M)*E_E(IMJK)*K_CP(IJKW)+&
180                       RO_S(IJKW,M)*U_S(IMJK,M))*AYZ(IMJK)
181                    A_M(IJK,0,0) = A_M(IJK,0,0) + &
182                       ROP_S(IJKW,M)*E_E(IMJK)*K_P*AYZ(IMJK)
183                    bmw = (ROP_S(IJKW,M)*U_S(IMJK,M))*AYZ(IMJK)
184                    B_M(IJK,0) = B_M(IJK,0) + bmw
185                 ELSE
186                    A_M(IJK,W,0) = (ROP_S(IJK,M)*&
187                       E_E(IMJK)*K_CP(IJKW))*AYZ(IMJK)
188                    A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_S(IJK,M)*&
189                       E_E(IMJK)*K_P-RO_S(IJK,M)*U_S(IMJK,M))*AYZ(IMJK)
190                    bmw = (ROP_S(IJK,M)*U_S(IMJK,M))*AYZ(IMJK)
191                    B_M(IJK,0) = B_M(IJK,0) + bmw
192                 ENDIF
193     
194     
195     ! North face (i, j+1/2, k)
196                 IF (V_S(IJK,M) < ZERO) THEN
197                    A_M(IJK,N,0) = (ROP_S(IJKN,M)*E_N(IJK)*K_CP(IJKN)-&
198                       RO_S(IJKN,M)*V_S(IJK,M))*AXZ(IJK)
199                    A_M(IJK,0,0)=A_M(IJK,0,0)+&
200                       ROP_S(IJKN,M)*E_N(IJK)*K_P*AXZ(IJK)
201                    bmn = (-ROP_S(IJKN,M)*V_S(IJK,M))*AXZ(IJK)
202                    B_M(IJK,0) = B_M(IJK,0) + bmn
203                 ELSE
204                    A_M(IJK,N,0) = (ROP_S(IJK,M)*&
205                       E_N(IJK)*K_CP(IJKN))*AXZ(IJK)
206                    A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_S(IJK,M)*&
207                       E_N(IJK)*K_P+RO_S(IJK,M)*V_S(IJK,M))*AXZ(IJK)
208                    bmn = (-ROP_S(IJK,M)*V_S(IJK,M))*AXZ(IJK)
209                    B_M(IJK,0) = B_M(IJK,0) + bmn
210                 ENDIF
211     
212     
213     ! South face (i, j-1/2, k)
214                 IF (V_S(IJMK,M) > ZERO) THEN
215                    A_M(IJK,S,0) = (ROP_S(IJKS,M)*E_N(IJMK)*K_CP(IJKS)+&
216                       RO_S(IJKS,M)*V_S(IJMK,M))*AXZ(IJMK)
217                    A_M(IJK,0,0) = A_M(IJK,0,0) + &
218                       ROP_S(IJKS,M)*E_N(IJMK)*K_P*AXZ(IJMK)
219                    bms = (ROP_S(IJKS,M)*V_S(IJMK,M))*AXZ(IJMK)
220                    B_M(IJK,0) = B_M(IJK,0) + bms
221                 ELSE
222                    A_M(IJK,S,0) = (ROP_S(IJK,M)*&
223                       E_N(IJMK)*K_CP(IJKS))*AXZ(IJMK)
224                    A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_S(IJK,M)*&
225                       E_N(IJMK)*K_P-RO_S(IJK,M)*V_S(IJMK,M))*AXZ(IJMK)
226                    bms = (ROP_S(IJK,M)*V_S(IJMK,M))*AXZ(IJMK)
227                    B_M(IJK,0) = B_M(IJK,0) + bms
228                 ENDIF
229     
230     
231                 IF (DO_K) THEN
232     ! Top face (i, j, k+1/2)
233                    IF (W_S(IJK,M) < ZERO) THEN
234                       A_M(IJK,T,0) = (ROP_S(IJKT,M)*E_T(IJK)*K_CP(IJKT)-&
235                          RO_S(IJKT,M)*W_S(IJK,M))*AXY(IJK)
236                       A_M(IJK,0,0) = A_M(IJK,0,0) + &
237                          ROP_S(IJKT,M)*E_T(IJK)*K_P*AXY(IJK)
238                       bmt = (-ROP_S(IJKT,M)*W_S(IJK,M))*AXY(IJK)
239                       B_M(IJK,0)=B_M(IJK,0) + bmt
240                    ELSE
241                       A_M(IJK,T,0) = (ROP_S(IJK,M)*&
242                          E_T(IJK)*K_CP(IJKT))*AXY(IJK)
243                       A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_S(IJK,M)*&
244                          E_T(IJK)*K_P+RO_S(IJK,M)*W_S(IJK,M))*AXY(IJK)
245                       bmt = (-ROP_S(IJK,M)*W_S(IJK,M))*AXY(IJK)
246                       B_M(IJK,0) = B_M(IJK,0) + bmt
247                    ENDIF
248     
249     
250     ! Bottom face (i, j, k-1/2)
251                    IF (W_S(IJKM,M) > ZERO) THEN
252                       A_M(IJK,B,0) = (ROP_S(IJKB,M)*E_T(IJKM)*K_CP(IJKB)+&
253                          RO_S(IJKB,M)*W_S(IJKM,M))*AXY(IJKM)
254                       A_M(IJK,0,0) = A_M(IJK,0,0) + &
255                          ROP_S(IJKB,M)*E_T(IJKM)*K_P*AXY(IJKM)
256                       bmb = (ROP_S(IJKB,M)*W_S(IJKM,M))*AXY(IJKM)
257                       B_M(IJK,0) = B_M(IJK,0) + bmb
258                    ELSE
259                       A_M(IJK,B,0) = (ROP_S(IJK,M)*&
260                          E_T(IJKM)*K_CP(IJKB))*AXY(IJKM)
261                       A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_S(IJK,M)*&
262                          E_T(IJKM)*K_P-RO_S(IJK,M)*W_S(IJKM,M))*AXY(IJKM)
263                       bmb = (ROP_S(IJK,M)*W_S(IJKM,M))*AXY(IJKM)
264                       B_M(IJK,0)=B_M(IJK,0) + bmb
265                    ENDIF
266     
267                 ELSE   ! not(DO_K) branch
268                    bmt = zero
269                    bmb = zero
270                 ENDIF   ! end if/else (do_K)
271     
272                 IF (ROP_S(IJK,M)>ZERO .AND. SUM_R_S(IJK,M)<ZERO) THEN
273                    SRC = VOL(IJK)*(-SUM_R_S(IJK,M))/ROP_S(IJK,M)
274                 ELSE
275                    SRC = ZERO
276                 ENDIF
277     
278                 A_M(IJK,0,0) = -(A_M(IJK,0,0)+VOL(IJK)*ODT*RO_S(IJK,M)+SRC*RO_S(IJK,M))
279     
280                 bma = (ROP_S(IJK,M)-ROP_SO(IJK,M))*VOL(IJK)*ODT
281                 bmr = SUM_R_S(IJK,M)*VOL(IJK)
282                 B_M(IJK,0) = -(B_M(IJK,0) - bma + bmr)
283                 B_MMAX(IJK,0) = max(abs(bma), abs(bme), abs(bmw), abs(bmn), abs(bms), abs(bmt), abs(bmb), abs(bmr) )
284     
285                 IF ((-A_M(IJK,0,0)) < SMALL_NUMBER) THEN
286                    IF (ABS(B_M(IJK,0)) < SMALL_NUMBER) THEN
287                       A_M(IJK,0,0) = -ONE            ! Equation is undefined.
288                       B_M(IJK,0) = ZERO              ! Use existing value
289                    ELSE
290     !!$omp             critical
291                       WRITE (LINE, '(A,I6,A,G12.5)') &
292                          'Error: At IJK = ', IJK, &
293                          ' A = 0 and b = ', B_M(IJK,0)
294                       CALL WRITE_ERROR ('CONV_SOURCE_EPp0', LINE, 1)
295     !!$omp             end critical
296                    ENDIF
297                 ENDIF
298              ELSE
299                 A_M(IJK,0,0) = -ONE
300                 B_M(IJK,0) = ZERO
301              ENDIF
302           ENDDO
303     
304           RETURN
305           END SUBROUTINE CONV_SOURCE_EPP0
306     
307     
308     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
309     !                                                                      C
310     !  Subroutine: CONV_SOURCE_EPp1                                        C
311     !  Purpose: Determine convection & source terms for solids volume      C
312     !           fraction correction equation.  Higher order scheme.        C
313     !                                                                      C
314     !  Author: M. Syamlal                                 Date: 24-SEP-96  C
315     !  Reviewer:                                          Date:            C
316     !                                                                      C
317     !                                                                      C
318     !                                                                      C
319     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
320     
321           SUBROUTINE CONV_SOURCE_EPP1(A_M, B_M, B_MMAX)
322     
323     !-----------------------------------------------
324     ! Modules
325     !-----------------------------------------------
326           USE param
327           USE param1
328           USE fldvar
329           USE run
330           USE parallel
331           USE matrix
332           USE constant
333           USE physprop
334           USE rxns
335           USE geometry
336           USE indices
337           USE pgcor
338           USE pscor
339           USE xsi
340           USE xsi_array
341           USE vshear
342           USE compar
343           USE sendrecv
344           USE solids_pressure
345           USE functions
346           IMPLICIT NONE
347     !-----------------------------------------------
348     ! Dummy arguments
349     !-----------------------------------------------
350     ! Septadiagonal matrix A_m
351           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
352     ! Vector b_m
353           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
354     ! maximum term in b_m expression
355           DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M)
356     !-----------------------------------------------
357     ! Local variables
358     !-----------------------------------------------
359     ! solids phase index locally assigned to mcp
360           INTEGER :: M
361     ! Indices
362           INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
363           INTEGER :: IMJK, IJMK, IJKM, IJKE, IJKW, IJKN, IJKS
364           INTEGER :: IJKB, IJKT
365     ! loezos : used for including shearing
366           INTEGER :: incr
367     ! dPodEP_s(EP_s(IJK, M)); the derivative of the plastic or friction
368     ! pressure with respect to ep_s
369           DOUBLE PRECISION :: K_P
370     ! Mass source
371           DOUBLE PRECISION :: Src
372     ! face value of ROP_s
373           DOUBLE PRECISION :: ROP_sf
374     ! terms of bm expression
375           DOUBLE PRECISION :: bma, bme, bmw, bmn, bms, bmt, bmb, bmr
376     ! error message
377           CHARACTER(LEN=80) :: LINE(1)
378     ! temporary use of global arrays:
379     ! xsi_array: convection weighting factors
380     !      DOUBLE PRECISION :: XSI_e(DIMENSION_3), &
381     !                          XSI_n(DIMENSION_3),&
382     !                          XSI_t(DIMENSION_3)
383     !-----------------------------------------------
384     
385           call lock_xsi_array
386     
387           IF (MCP == UNDEFINED_I) THEN
388     ! this error should be caught earlier in the routines so that this
389     ! branch should never be entered
390              RETURN
391           ELSE
392     ! the lowest solids phase index of those solids phases that can close
393     ! pack (i.e. close_packed=T) and the index of the solids phase that is
394     ! used to form the solids correction equation.
395              M = MCP
396           ENDIF
397     
398     ! Loezos:
399           incr=0
400     
401     ! Calculate convection factors
402           CALL CALC_XSI (DISCRETIZE(2), ROP_S(1,M), U_S(1,M), V_S(1,M), W_S(1,M), &
403              XSI_E, XSI_N, XSI_T,incr)
404     
405     ! Loezos:
406     ! update to true velocity
407           IF (SHEAR) THEN
408     !!$omp parallel do private(IJK)
409              DO IJK = ijkstart3, ijkend3
410                 IF (FLUID_AT(IJK)) THEN
411                    V_S(IJK,m)=V_s(IJK,m)+VSH(IJK)
412                 ENDIF
413              ENDDO
414           ENDIF
415     
416     !!$omp parallel do                                                      &
417     !!$omp&   private(I, J, K, IJK, IPJK, IJPK, IJKP,                &
418     !!$omp&           IMJK, IJMK, IJKM, IJKE, IJKW, IJKN, IJKS, IJKT, IJKB, &
419     !!$omp&           K_P,ROP_SF,SRC, bma, bme, bmw, bmn, bms, bmt, bmb, bmr )
420           DO IJK = ijkstart3, ijkend3
421     ! Determine if IJK falls within 1 ghost layer........
422              I = I_OF(IJK)
423              J = J_OF(IJK)
424              K = K_OF(IJK)
425     
426              IF (FLUID_AT(IJK)) THEN
427                 IPJK = IP_OF(IJK)
428                 IJPK = JP_OF(IJK)
429                 IJKP = KP_OF(IJK)
430                 IMJK = IM_OF(IJK)
431                 IJMK = JM_OF(IJK)
432                 IJKM = KM_OF(IJK)
433     
434                 IJKE = EAST_OF(IJK)
435                 IJKW = WEST_OF(IJK)
436                 IJKN = NORTH_OF(IJK)
437                 IJKS = SOUTH_OF(IJK)
438                 IJKT = TOP_OF(IJK)
439                 IJKB = BOTTOM_OF(IJK)
440     
441     ! initializing
442                 A_M(IJK,0,0) = ZERO
443                 B_M(IJK,0) = ZERO
444                 K_P = K_CP(IJK)
445     
446     ! Calculate convection-diffusion fluxes through each of the faces
447     
448     ! East face (i+1/2, j, k)
449                 ROP_SF = ROP_S(IJKE,M)*XSI_E(IJK) + ROP_S(IJK,M)*(ONE - XSI_E(IJK))
450                 A_M(IJK,E,0) = (ROP_SF*E_E(IJK)*K_CP(IJKE)-RO_S(IJK,M)*U_S(IJK,M)*XSI_E&
451                    (IJK))*AYZ(IJK)
452                 A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_SF*E_E(IJK)*K_P+RO_S(IJK,M)*U_S(IJK,&
453                    M)*(ONE-XSI_E(IJK)))*AYZ(IJK)
454                 bme = (-ROP_SF*U_S(IJK,M))*AYZ(IJK)
455                 B_M(IJK,0) = B_M(IJK,0) + bme
456     
457     ! West face (i-1/2, j, k)
458                 ROP_SF=ROP_S(IJK,M)*XSI_E(IMJK)+ROP_S(IJKW,M)*(ONE-XSI_E(IMJK))
459                 A_M(IJK,W,0) = (ROP_SF*E_E(IMJK)*K_CP(IJKW)+RO_S(IJK,M)*U_S(IMJK,M)*(&
460                    ONE-XSI_E(IMJK)))*AYZ(IMJK)
461                 A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_SF*E_E(IMJK)*K_P-RO_S(IJK,M)*U_S(&
462                    IMJK,M)*XSI_E(IMJK))*AYZ(IMJK)
463                 bmw = (ROP_SF*U_S(IMJK,M))*AYZ(IMJK)
464                 B_M(IJK,0) = B_M(IJK,0) +  bmw
465     
466     
467     ! North face (i, j+1/2, k)
468                 ROP_SF = ROP_S(IJKN,M)*XSI_N(IJK) + ROP_S(IJK,M)*(ONE - XSI_N(IJK))
469                 A_M(IJK,N,0) = (ROP_SF*E_N(IJK)*K_CP(IJKN)-RO_S(IJK,M)*V_S(IJK,M)*XSI_N&
470                    (IJK))*AXZ(IJK)
471                 A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_SF*E_N(IJK)*K_P+RO_S(IJK,M)*V_S(IJK,&
472                    M)*(ONE-XSI_N(IJK)))*AXZ(IJK)
473                 bmn = (-ROP_SF*V_S(IJK,M))*AXZ(IJK)
474                 B_M(IJK,0) = B_M(IJK,0) + bmn
475     
476     
477     ! South face (i, j-1/2, k)
478                 ROP_SF=ROP_S(IJK,M)*XSI_N(IJMK)+ROP_S(IJKS,M)*(ONE-XSI_N(IJMK))
479                 A_M(IJK,S,0) = (ROP_SF*E_N(IJMK)*K_CP(IJKS)+RO_S(IJK,M)*V_S(IJMK,M)*(&
480                    ONE-XSI_N(IJMK)))*AXZ(IJMK)
481                 A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_SF*E_N(IJMK)*K_P-RO_S(IJK,M)*V_S(&
482                    IJMK,M)*XSI_N(IJMK))*AXZ(IJMK)
483                 bms = (ROP_SF*V_S(IJMK,M))*AXZ(IJMK)
484                 B_M(IJK,0) = B_M(IJK,0) + bms
485     
486     
487                 IF (DO_K) THEN
488     ! Top face (i, j, k+1/2)
489                    ROP_SF=ROP_S(IJKT,M)*XSI_T(IJK)+ROP_S(IJK,M)*(ONE-XSI_T(IJK))
490                    A_M(IJK,T,0) = (ROP_SF*E_T(IJK)*K_CP(IJKT)-RO_S(IJK,M)*W_S(IJK,M)*&
491                       XSI_T(IJK))*AXY(IJK)
492                    A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_SF*E_T(IJK)*K_P+RO_S(IJK,M)*W_S(&
493                       IJK,M)*(ONE-XSI_T(IJK)))*AXY(IJK)
494                    bmt = (-ROP_SF*W_S(IJK,M))*AXY(IJK)
495                    B_M(IJK,0) = B_M(IJK,0) + bmt
496     
497     ! Bottom face (i, j, k-1/2)
498                    ROP_SF = ROP_S(IJK,M)*XSI_T(IJKM) + ROP_S(IJKB,M)*(ONE - XSI_T(&
499                       IJKM))
500                    A_M(IJK,B,0) = (ROP_SF*E_T(IJKM)*K_CP(IJKB)+RO_S(IJK,M)*W_S(IJKM,M)*&
501                       (ONE-XSI_T(IJKM)))*AXY(IJKM)
502                    A_M(IJK,0,0) = A_M(IJK,0,0) + (ROP_SF*E_T(IJKM)*K_P-RO_S(IJK,M)*W_S(&
503                       IJKM,M)*XSI_T(IJKM))*AXY(IJKM)
504                    bmb = (ROP_SF*W_S(IJKM,M))*AXY(IJKM)
505                    B_M(IJK,0) = B_M(IJK,0) + bmb
506                 ELSE   ! not(do_k) branch
507                   bmt = zero
508                   bmb = zero
509                ENDIF   ! end if/else (do_k)
510     
511                 IF (ROP_S(IJK,M)>ZERO .AND. SUM_R_S(IJK,M)<ZERO) THEN
512                    SRC = VOL(IJK)*(-SUM_R_S(IJK,M))/ROP_S(IJK,M)
513                 ELSE
514                    SRC = ZERO
515                 ENDIF
516     
517                 A_M(IJK,0,0) = -(A_M(IJK,0,0)+VOL(IJK)*ODT*RO_S(IJK,M)+SRC*RO_S(IJK,M))
518     
519                 bma = (ROP_S(IJK,M)-ROP_SO(IJK,M))*VOL(IJK)*ODT
520                 bmr = SUM_R_S(IJK,M)*VOL(IJK)
521                 B_M(IJK,0) = -(B_M(IJK,0)- bma + bmr)
522                 B_MMAX(IJK,0) = max(abs(bma), abs(bme), abs(bmw), abs(bmn), abs(bms), abs(bmt), abs(bmb), abs(bmr) ) !
523     
524                 IF (ABS(A_M(IJK,0,0)) < SMALL_NUMBER) THEN
525                    IF (ABS(B_M(IJK,0)) < SMALL_NUMBER) THEN
526                       A_M(IJK,0,0) = -ONE            ! Equation is undefined.
527                       B_M(IJK,0) = ZERO              ! Use existing value
528                    ELSE
529     !!$omp             critical
530                       WRITE (LINE(1), '(A,I6,A,G12.5)') &
531                          'Error: At IJK = ', IJK, &
532                          ' A = 0 and b = ', B_M(IJK,0)
533     ! Having problem to compile this statement on SGI
534                       CALL WRITE_ERROR ('CONV_SOURCE_EPp1', LINE, 1)
535     !!$omp             end critical
536                    ENDIF
537                 ENDIF
538              ELSE
539                 A_M(IJK,0,0) = -ONE
540                 B_M(IJK,0) = ZERO
541              ENDIF
542           ENDDO
543     
544     ! Loezos:
545           IF (SHEAR) THEN
546     !!$omp parallel do private(IJK)
547              DO IJK = ijkstart3, ijkend3
548                 IF (FLUID_AT(IJK)) THEN
549                    V_S(IJK,m)=V_s(IJK,m)-VSH(IJK)
550                 ENDIF
551              ENDDO
552           ENDIF
553     
554           call unlock_xsi_array
555     
556     
557           RETURN
558           END SUBROUTINE CONV_SOURCE_EPP1
559     
560     
561     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
562     !                                                                      C
563     !  Subroutine: POINT_SOURCE_EPP                                        C
564     !  Purpose: Adds point sources to the solids volume fraction           C
565     !           correction equation.                                       C
566     !                                                                      C
567     !  Notes: The off-diagonal coefficients are positive. The center       C
568     !         coefficient and the source vector are negative. See          C
569     !         conv_Pp_g                                                    C
570     !                                                                      C
571     !  Author: J. Musser                                  Date: 10-JUN-13  C
572     !  Reviewer:                                          Date:            C
573     !                                                                      C
574     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
575           SUBROUTINE POINT_SOURCE_EPP(B_M, B_MMAX)
576     
577     !-----------------------------------------------
578     ! Modules
579     !-----------------------------------------------
580           use compar
581           use constant
582           use geometry
583           use indices
584           use physprop
585           use ps
586           use pscor
587           use run
588           use functions
589     
590           IMPLICIT NONE
591     !-----------------------------------------------
592     ! Dummy arguments
593     !-----------------------------------------------
594     ! Vector b_m
595           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
596     ! maximum term in b_m expression
597           DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M)
598     
599     !-----------------------------------------------
600     ! Local Variables
601     !-----------------------------------------------
602     
603     ! Indices
604           INTEGER :: IJK, I, J, K
605           INTEGER :: PSV, M
606     
607     ! terms of bm expression
608           DOUBLE PRECISION :: pSource
609     
610     !-----------------------------------------------
611     
612           IF (MCP == UNDEFINED_I) THEN
613     ! this error should be caught earlier in the routines so that this
614     ! branch should never be entered
615              RETURN
616           ELSE
617     ! the lowest solids phase index of those solids phases that can close
618     ! pack (i.e. close_packed=T) and the index of the solids phase that is
619     ! used to form the solids correction equation.
620              M = MCP
621           ENDIF
622     
623     
624           PS_LP: do PSV = 1, DIMENSION_PS
625     
626              if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
627              if(PS_MASSFLOW_S(PSV,M) < small_number) cycle PS_LP
628     
629              do k = PS_K_B(PSV), PS_K_T(PSV)
630              do j = PS_J_S(PSV), PS_J_N(PSV)
631              do i = PS_I_W(PSV), PS_I_E(PSV)
632     
633                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
634                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
635                 ijk = funijk(i,j,k)
636                 if(fluid_at(ijk)) then
637     
638                    pSource =  PS_MASSFLOW_S(PSV,M) * &
639                       (VOL(IJK)/PS_VOLUME(PSV))
640     
641                    B_M(IJK,0) = B_M(IJK,0) - pSource
642                    B_MMAX(IJK,0) = max(abs(B_MMAX(IJK,0)), abs(B_M(IJK,0)))
643                 endif
644              enddo
645              enddo
646              enddo
647     
648           enddo PS_LP
649     
650           RETURN
651           END SUBROUTINE POINT_SOURCE_EPP
652