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