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

1     ! TODO:
2     ! 1. Account for solids-solids transfer terms for scalars.  Will need
3     !    to pass the exchange coefficients to this routine (unlike momentum
4     !    eqs)
5     ! 2. repetitive code may be eliminated by defining and using functions
6     
7     
8     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
9     !                                                                      C
10     !  Subroutine: PARTIAL_ELIM_S                                          C
11     !  Purpose: Do partial elimination for scalar quantities               C
12     !                                                                      C
13     !                                                                      C
14     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
15     !  Reviewer:                                          Date:            C
16     !                                                                      C
17     !  Comments:                                                           C
18     !     Currently does not evaluate partial elimination terms between    C
19     !     unlike solids phases.  The exchange coefficient must be hard     C
20     !     coded explicitly for this case - it is not a passed dummy        C
21     !     argument                                                         C
22     !                                                                      C
23     !  Literature/Document References:                                     C
24     !     Spalding, D.B., 1980, Numerical Computation of Multi-phase       C
25     !        fluid flow and heat transfer, Recent Advances in Numerical    C
26     !        Methods in Fluids, C. Taylor et al., eds, Pineridge Press     C
27     !     Syamlal, M., 1998, MFIX Documentation: Numerical Technique,      C
28     !        Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G      C
29     !        Technical Services of West Virginia, Inc., Morgantown, WV.    C
30     !                                                                      C
31     !                                                                      C
32     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
33     
34           SUBROUTINE PARTIAL_ELIM_S(VAR_G, VAR_S, VXF, A_M, B_M)
35     
36     !-----------------------------------------------
37     ! Modules
38     !-----------------------------------------------
39           USE param
40           USE param1
41           USE parallel
42           USE geometry
43           USE physprop
44           USE indices
45           USE compar
46           USE drag
47           USE fldvar
48           USE functions
49           IMPLICIT NONE
50     !-----------------------------------------------
51     ! Dummy arguments
52     !-----------------------------------------------
53     ! gas phase variable
54           DOUBLE PRECISION, INTENT(IN) :: Var_g(DIMENSION_3)
55     ! solids phase variable
56           DOUBLE PRECISION, INTENT(IN) :: Var_s(DIMENSION_3, DIMENSION_M)
57     ! Volume x gas-solids transfer coefficient
58           DOUBLE PRECISION, INTENT(IN) :: VxF(DIMENSION_3, DIMENSION_M)
59     ! Septadiagonal matrix A_m
60           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
61     ! Vector b_m
62           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
63     !-----------------------------------------------
64     ! Local variables
65     !-----------------------------------------------
66     ! Indices
67           INTEGER :: IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT
68           INTEGER :: L, M, LP
69     !
70           DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
71     ! a0, b0 etc.
72           DOUBLE PRECISION :: a(0:DIMENSION_M), BB(0:DIMENSION_M),&
73                               F(0:DIMENSION_M,0:DIMENSION_M),&
74                               Saxf(0:DIMENSION_M)
75     !-----------------------------------------------
76     
77     !!$omp  parallel do private( IJKW, IJKS, IJKB, IJKE, IJKN, IJKT,  &
78     !!$omp&  a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, den) &
79     !!$omp&  schedule(static)
80           DO IJK = ijkstart3, ijkend3
81              IF (FLUID_AT(IJK)) THEN
82                 IJKW = WEST_OF(IJK)
83                 IJKS = SOUTH_OF(IJK)
84                 IJKE = EAST_OF(IJK)
85                 IJKN = NORTH_OF(IJK)
86                 DO M=0, MMAX
87                    A(M)=A_M(IJK,0,M)
88                    BB(M)=B_M(IJK,M)
89     
90                    if (m .ne. 0) then
91                       F(M,0)=-VXF(IJK,M)
92                       F(0,M)=F(M,0)
93                    else
94                       F(0,0) = ZERO
95                    endif
96     
97                    DO L =1, MMAX
98                       IF ((L .NE. M) .AND. (M .NE. 0)) THEN
99                          F(M,L) = ZERO  !insert solids-solids exchange coefficients here
100                       ELSE
101                          F(M,L) = ZERO
102                       ENDIF
103                       F(L,M) = ZERO
104                    ENDDO
105     
106                    IF (M == 0 ) THEN
107                       SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IJKE)+A_M(IJK,west,M)*VAR_G(IJKW&
108                            )+A_M(IJK,north,M)*VAR_G(IJKN)+A_M(IJK,south,M)*VAR_G(IJKS))
109                    ELSE
110                       SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IJKE,M)+A_M(IJK,west,M)*VAR_S(&
111                            IJKW,M)+A_M(IJK,north,M)*VAR_S(IJKN,M)+A_M(IJK,south,M)*VAR_S(&
112                            IJKS,M))
113                    ENDIF
114     
115                    IF (DO_K) THEN
116                       IJKB = BOTTOM_OF(IJK)
117                       IJKT = TOP_OF(IJK)
118                       IF ( M ==0) THEN
119                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKT)+A_M(IJK,bottom,M)*&
120                                VAR_G(IJKB))
121                       ELSE
122                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKT,M)+A_M(IJK,bottom,M)*&
123                                VAR_S(IJKB,M))
124                       ENDIF
125                    ENDIF
126                 ENDDO
127     
128                 DO M=0,MMAX
129                    SUM_A = ZERO
130                    SUM_B = ZERO
131     
132                    DO L =0,MMAX
133                       SUM_A_LPRIME = ZERO
134                       SUM_B_LPRIME = ZERO
135     
136                       DO LP=0,MMAX
137                          IF ( LP .NE. M) THEN
138                             SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
139                             IF (LP == 0) THEN
140                                SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
141                             ELSE
142                                SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
143                             ENDIF
144                          ENDIF
145                       ENDDO
146     
147                       DEN = A(L) + SUM_A_LPRIME + F(L,M)
148                       IF ( DEN .NE. ZERO) THEN
149                          SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
150                          SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
151                                  )/DEN
152                       ENDIF
153                    ENDDO
154     
155                    A_M(IJK,0,M)= SUM_A+A(M)
156                    B_M(IJK,M)  = SUM_B+BB(M)
157                 ENDDO
158              ENDIF
159           ENDDO
160     
161           RETURN
162           END SUBROUTINE PARTIAL_ELIM_S
163     
164     
165     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
166     !                                                                      C
167     !  Subroutine: PARTIAL_ELIM_IA                                         C
168     !  Purpose: Do partial elimination for granular temperature coupling   C
169     !                                                                      C
170     !                                                                      C
171     !  Author: J. Galvin                                  Date: 21-MAY-05  C
172     !  Reviewer: S. Benyahia                              Date:            C
173     !                                                                      C
174     !  Literature/Document References:                                     C
175     !    Iddir, Y.H., Modeling of the multiphase mixture of particles      C
176     !       using the kinetic theory approach, PhD Thesis, Illinois        C
177     !       Institute of Technology, Chicago, Illinois, 2004               C
178     !    Iddir, Y.H., & H. Arastoopour, Modeling of multitype particle     C
179     !       flow using the kinetic theory approach, AIChE J., Vol 51,      C
180     !       No 6, June 2005                                                C
181     !                                                                      C
182     !                                                                      C
183     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
184     
185           SUBROUTINE PARTIAL_ELIM_IA(VAR_S, VXTCSS, A_M, B_M)
186     
187     !-----------------------------------------------
188     ! Modules
189     !-----------------------------------------------
190           USE param
191           USE param1
192           USE parallel
193           USE geometry
194           USE physprop
195           USE indices
196           USE compar
197           USE drag
198           USE fldvar
199           USE functions
200           IMPLICIT NONE
201     !-----------------------------------------------
202     ! Dummy arguments
203     !-----------------------------------------------
204     ! solids phase variable
205           DOUBLE PRECISION, INTENT(IN) :: Var_s(DIMENSION_3, DIMENSION_M)
206     ! Volume x solids-solids transfer coefficient
207           DOUBLE PRECISION, INTENT(IN) :: VxTCss (DIMENSION_3, DIMENSION_LM)
208     ! Septadiagonal matrix A_m
209           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
210     ! Vector b_m
211           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
212     !-----------------------------------------------
213     ! Local variables
214     !-----------------------------------------------
215     ! Indices
216           INTEGER :: IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT
217           INTEGER :: L, M, LP, LM
218     !
219           DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
220     ! a0, b0 etc.
221           DOUBLE PRECISION :: a(DIMENSION_M), BB(DIMENSION_M),&
222                               F(DIMENSION_M,DIMENSION_M),&
223                               Saxf(DIMENSION_M)
224     !-----------------------------------------------
225     
226     !!$omp  parallel do private( IJKW, IJKS, IJKB, IJKE, IJKN, IJKT,  &
227     !!$omp&  a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, den) &
228     !!$omp&  schedule(static)
229     
230           DO IJK = ijkstart3, ijkend3
231     
232              IF (FLUID_AT(IJK)) THEN
233                 IJKW = WEST_OF(IJK)
234                 IJKS = SOUTH_OF(IJK)
235                 IJKE = EAST_OF(IJK)
236                 IJKN = NORTH_OF(IJK)
237     
238                 F(:,:) = ZERO
239     
240                 DO M = 1, MMAX
241                    A(M)=A_M(IJK,0,M)
242                    BB(M)=B_M(IJK,M)
243     
244                    DO L = 1, MMAX
245                       IF ( L .NE. M ) THEN
246                          LM = FUNLM(L,M)
247                          F(M,L)=-VxTCSS(IJK,LM)
248                       ENDIF
249                       F(L,M) = F(M,L)
250                    ENDDO
251     
252                    SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IJKE,M)+&
253                       A_M(IJK,west,M)*VAR_S(IJKW,M)+&
254                       A_M(IJK,north,M)*VAR_S(IJKN,M)+&
255                       A_M(IJK,south,M)*VAR_S(IJKS,M))
256     
257                    IF (DO_K) THEN
258                       IJKB = BOTTOM_OF(IJK)
259                       IJKT = TOP_OF(IJK)
260                       SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKT,M)+&
261                          A_M(IJK,bottom,M)*VAR_S(IJKB,M))
262                    ENDIF
263                 ENDDO
264     
265                 DO M = 1, MMAX
266                    SUM_A = ZERO
267                    SUM_B = ZERO
268     
269                    DO L = 1, MMAX
270                       SUM_A_LPRIME = ZERO
271                       SUM_B_LPRIME = ZERO
272     
273                       DO LP = 1, MMAX
274                          IF ( LP .NE. M) THEN
275                             SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
276                             SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
277                          ENDIF
278                       ENDDO
279     
280                       DEN = A(L) + SUM_A_LPRIME + F(L,M)
281     
282                       IF ( DEN .NE. ZERO) THEN
283                          SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
284                          SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+&
285                             SUM_B_LPRIME)/DEN
286                       ENDIF
287                    ENDDO
288     
289                    A_M(IJK,0,M)= SUM_A+A(M)
290                    B_M(IJK,M)  = SUM_B+BB(M)
291                 ENDDO
292     
293              ENDIF
294           ENDDO
295     
296           RETURN
297           END SUBROUTINE PARTIAL_ELIM_IA
298     
299     
300     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
301     !                                                                      C
302     !  Subroutine: PARTIAL_ELIM_U                                          C
303     !  Purpose: Do partial elimination for X vector quantities             C
304     !                                                                      C
305     !                                                                      C
306     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
307     !  Reviewer:                                          Date:            C
308     !                                                                      C
309     !  Comments:                                                           C
310     !     Unlike gas-solids exchange coefficient (multiplied by volume),   C
311     !     the solids-solids exchange coefficient is not as a passed dummy  C
312     !     argument to this routine but simply re-evaluated here.           C
313     !                                                                      C
314     !  Literature/Document References:                                     C
315     !     Spalding, D.B., 1980, Numerical Computation of Multi-phase       C
316     !        fluid flow and heat transfer, Recent Advances in Numerical    C
317     !        Methods in Fluids, C. Taylor et al., eds, Pineridge Press     C
318     !     Syamlal, M., 1998, MFIX Documentation: Numerical Technique,      C
319     !        Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G      C
320     !        Technical Services of West Virginia, Inc., Morgantown, WV.    C
321     !                                                                      C
322     !                                                                      C
323     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
324     
325           SUBROUTINE PARTIAL_ELIM_U(VAR_G, VAR_S, VXF, A_M, B_M)
326     
327     !-----------------------------------------------
328     ! Modules
329     !-----------------------------------------------
330           USE param
331           USE param1
332           USE parallel
333           USE geometry
334           USE physprop
335           USE indices
336           USE run
337           USE compar
338           USE drag
339           USE fldvar
340           USE fun_avg
341           USE functions
342           IMPLICIT NONE
343     !-----------------------------------------------
344     ! Dummy arguments
345     !-----------------------------------------------
346     ! gas phase variable
347           DOUBLE PRECISION, INTENT(IN) :: Var_g(DIMENSION_3)
348     ! solids phase variable
349           DOUBLE PRECISION, INTENT(IN) :: Var_s(DIMENSION_3, DIMENSION_M)
350     ! Volume x gas-solids transfer coefficient
351           DOUBLE PRECISION, INTENT(IN) :: VxF(DIMENSION_3, DIMENSION_M)
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     !-----------------------------------------------
357     ! Local variables
358     !-----------------------------------------------
359     ! Indices
360           INTEGER :: IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, I, IJKE
361           INTEGER :: L, M, LP, LM
362     !
363           DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
364     ! a0, b0 etc.
365           DOUBLE PRECISION :: a(0:DIMENSION_M), BB(0:DIMENSION_M), &
366                               F(0:DIMENSION_M,0:DIMENSION_M),&
367                               Saxf(0:DIMENSION_M)
368     !-----------------------------------------------
369     
370     !!$omp  parallel do private( IMJK, IJMK, IJKM, IPJK, IJPK, IJKP,  &
371     !!$omp&  a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, den) &
372     !!$omp&  schedule(static)
373           DO IJK = ijkstart3, ijkend3
374              IF (FLOW_AT_E(IJK)) THEN
375                 IMJK = IM_OF(IJK)
376                 IJMK = JM_OF(IJK)
377                 IPJK = IP_OF(IJK)
378                 IJPK = JP_OF(IJK)
379     
380                 F = ZERO
381                 DO M=0, MMAX
382                    A(M) = A_M(IJK,0,M)
383                    BB(M)= B_M(IJK,M)
384     
385     ! locally storing vol x avg gas-solids drag coefficient at momentum cell
386     ! center
387                    IF (M .NE. 0) THEN
388                       F(M,0)=-VXF(IJK,M)
389                       F(0,M) = F(M,0)
390                    ENDIF
391     
392     ! locally storing vol x avg solids-solids drag coefficient at momentum cell
393     ! center
394                    DO L =1, MMAX
395                       IF ((L .NE. M) .AND. (M .NE. 0)) THEN
396                          LM = FUNLM(L,M)
397                          IF (.NOT.IP_AT_E(IJK)) THEN
398                             I = I_OF(IJK)
399                             IJKE = EAST_OF(IJK)
400                             F(M,L) = -AVG_X(F_SS(IJK,LM),F_SS(IJKE,LM),I)*VOL_U(IJK)
401                          ENDIF
402                       ENDIF
403                       F(L,M)=F(M,L)
404                    ENDDO
405     
406                    IF (M == 0) THEN
407                       SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IPJK)+A_M(IJK,west,M)*VAR_G(IMJK&
408                                 )+A_M(IJK,north,M)*VAR_G(IJPK)+A_M(IJK,south,M)*VAR_G(IJMK))
409                    ELSE
410                       SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IPJK,M)+A_M(IJK,west,M)*VAR_S(&
411                            IMJK,M)+A_M(IJK,north,M)*VAR_S(IJPK,M)+A_M(IJK,south,M)*VAR_S(&
412                            IJMK,M))
413                    ENDIF
414     
415                    IF (DO_K) THEN
416                       IJKM = KM_OF(IJK)
417                       IJKP = KP_OF(IJK)
418                       IF (M ==0) THEN
419                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKP)+A_M(IJK,bottom,M)*&
420                             VAR_G(IJKM))
421                       ELSE
422                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKP,M)+A_M(IJK,bottom,M)*&
423                             VAR_S(IJKM,M))
424                       ENDIF
425                    ENDIF
426                ENDDO   ! end do loop for M over 0 to MMAX
427     
428                DO M=0,MMAX
429                    IF (MOMENTUM_X_EQ(M)) THEN
430                       SUM_A = ZERO
431                       SUM_B = ZERO
432     
433                       DO L =0,MMAX
434     ! - should only need to loop for L/=M. since the resulting term should
435     !   evaluate to zero anyway, don't bother to add an explicit conditional.
436     ! - if phase L's momentum equation is not solved then do not evaluate the
437     !   outer sum for phase L (as phase L has no momentum equation). however,
438     !   the inner sum may include terms from phase L as all phases still see
439     !   phase L and therefore experience drag with phase L
440                          IF (MOMENTUM_X_EQ(L)) THEN
441                             SUM_A_LPRIME = ZERO
442                             SUM_B_LPRIME = ZERO
443     
444                             DO LP=0,MMAX    ! only need to loop for LP/=L; but
445                                             ! evaluates to zero so don't bother
446                                             ! with an explicit conditional
447                                IF ( LP .NE. M) THEN
448                                   SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
449                                   IF (LP == 0) THEN
450                                      SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
451                                   ELSE
452                                       SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
453                                   ENDIF
454                                ENDIF
455                             ENDDO
456     
457                             DEN = A(L) + SUM_A_LPRIME + F(L,M)
458                             IF ( DEN .NE. ZERO) THEN
459                                SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
460                                SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
461                                      )/DEN
462                             ENDIF
463                          ELSE
464                             SUM_A = SUM_A + F(L,M)
465                             IF (L == 0) THEN
466                                SUM_B = SUM_B + F(L,M)*VAR_G(IJK)
467                             ELSE
468                                SUM_B = SUM_B + F(L,M)*VAR_S(IJK,L)
469                             ENDIF
470                          ENDIF   ! end if momentum_x_eq(L)
471                       ENDDO   !end do loop for L over 0 to mmax
472     
473                       A_M(IJK,0,M) = SUM_A+A(M)
474                       B_M(IJK,M)  =  SUM_B+BB(M)
475                    ENDIF    ! end if momentum_x_eq(M)
476                 ENDDO       ! end do loop for M over 0 to MMAX
477     
478              ENDIF      ! end if flow_at_e(ijk)
479           ENDDO         ! end do loop for ijk over ijkstart3,ijkend3
480     
481           RETURN
482           END SUBROUTINE PARTIAL_ELIM_U
483     
484     
485     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
486     !                                                                      C
487     !  Module name: PARTIAL_ELIM_V                                         C
488     !  Purpose: Do partial elimination for Y vector quantities             C
489     !                                                                      C
490     !                                                                      C
491     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
492     !  Reviewer:                                          Date:            C
493     !                                                                      C
494     !  Comments:                                                           C
495     !     Unlike gas-solids exchange coefficient (multiplied by volume),   C
496     !     the solids-solids exchange coefficient is not as a passed dummy  C
497     !     argument to this routine but simply re-evaluated here.           C
498     !                                                                      C
499     !  Literature/Document References:                                     C
500     !     Spalding, D.B., 1980, Numerical Computation of Multi-phase       C
501     !        fluid flow and heat transfer, Recent Advances in Numerical    C
502     !        Methods in Fluids, C. Taylor et al., eds, Pineridge Press     C
503     !     Syamlal, M., 1998, MFIX Documentation: Numerical Technique,      C
504     !        Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G      C
505     !        Technical Services of West Virginia, Inc., Morgantown, WV.    C
506     !                                                                      C
507     !                                                                      C
508     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
509     
510           SUBROUTINE PARTIAL_ELIM_V(VAR_G, VAR_S, VXF, A_M, B_M)
511     
512     !-----------------------------------------------
513     ! Modules
514     !-----------------------------------------------
515           USE param
516           USE param1
517           USE parallel
518           USE geometry
519           USE physprop
520           USE indices
521           USE run
522           USE compar
523           USE drag
524           USE fldvar
525           USE fun_avg
526           USE functions
527           IMPLICIT NONE
528     !-----------------------------------------------
529     ! Dummy arguments
530     !-----------------------------------------------
531     ! gas phase variable
532           DOUBLE PRECISION, INTENT(IN) :: Var_g(DIMENSION_3)
533     ! solids phase variable
534           DOUBLE PRECISION, INTENT(IN) :: Var_s(DIMENSION_3, DIMENSION_M)
535     ! Volume x gas-solids transfer coefficient
536           DOUBLE PRECISION, INTENT(IN) :: VxF(DIMENSION_3, DIMENSION_M)
537     ! Septadiagonal matrix A_m
538           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
539     ! Vector b_m
540           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
541     !-----------------------------------------------
542     ! Local variables
543     !-----------------------------------------------
544     ! Indices
545           INTEGER :: IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, IJKN, J
546           INTEGER ::  L, M, LP, LM
547     !
548           DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
549     ! a0, b0 etc.
550           DOUBLE PRECISION :: a(0:DIMENSION_M), BB(0:DIMENSION_M),&
551                               F(0:DIMENSION_M,0:DIMENSION_M),&
552                               Saxf(0:DIMENSION_M)
553     !-----------------------------------------------
554     
555     !!$omp  parallel do private( IMJK, IJMK, IJKM, IPJK, IJPK, IJKP,  &
556     !!$omp&  a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, DEN) &
557     !!$omp&  schedule(static)
558           DO IJK = ijkstart3, ijkend3
559              IF (FLOW_AT_N(IJK)) THEN
560                 IMJK = IM_OF(IJK)
561                 IJMK = JM_OF(IJK)
562                 IPJK = IP_OF(IJK)
563                 IJPK = JP_OF(IJK)
564     
565                 F = ZERO
566                 DO M = 0, MMAX
567                    A(M)=A_M(IJK,0,M)
568                    BB(M)=B_M(IJK,M)
569     
570     ! locally storing vol x avg gas-solids drag coefficient at momentum cell
571     ! center
572                    IF (M .NE. 0) THEN
573                       F(M,0)=-VXF(IJK,M)
574                       F(0,M)=F(M,0)
575                    ENDIF
576     
577     ! locally storing vol x avg solids-solids drag coefficient at momentum cell
578     ! center
579                    DO L =1, MMAX
580                       IF ((L .NE. M) .AND. (M .NE. 0)) THEN
581                          LM = FUNLM(L,M)
582                          IF (.NOT.IP_AT_N(IJK)) THEN
583                             J = J_OF(IJK)
584                             IJKN = NORTH_OF(IJK)
585                             F(M,L)=-AVG_Y(F_SS(IJK,LM),F_SS(IJKN,LM),J)*VOL_V(IJK)
586                          ENDIF
587                       ENDIF
588                       F(L,M)=F(M,L)
589                    ENDDO
590     
591                    IF (M == 0) THEN
592                       SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IPJK)+A_M(IJK,west,M)*VAR_G(IMJK&
593                             )+A_M(IJK,north,M)*VAR_G(IJPK)+A_M(IJK,south,M)*VAR_G(IJMK))
594                    ELSE
595                       SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IPJK,M)+A_M(IJK,west,M)*VAR_S(&
596                            IMJK,M)+A_M(IJK,north,M)*VAR_S(IJPK,M)+A_M(IJK,south,M)*VAR_S(&
597                            IJMK,M))
598                    ENDIF
599     
600                    IF (DO_K) THEN
601                       IJKM = KM_OF(IJK)
602                       IJKP = KP_OF(IJK)
603                       IF (M == 0) THEN
604                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKP)+A_M(IJK,bottom,M)*&
605                             VAR_G(IJKM))
606                       ELSE
607                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKP,M)+A_M(IJK,bottom,M)*&
608                             VAR_S(IJKM,M))
609                       ENDIF
610                    ENDIF
611                 ENDDO     ! end do loop for M over 0 to MMAX
612     
613     
614                 DO M=0,MMAX
615                    IF (MOMENTUM_Y_EQ(M)) THEN
616                       SUM_A = ZERO
617                       SUM_B = ZERO
618     
619                       DO L =0,MMAX
620                          IF(MOMENTUM_Y_EQ(L)) THEN
621                             SUM_A_LPRIME = ZERO
622                             SUM_B_LPRIME = ZERO
623     
624                             DO LP=0,MMAX
625                                IF ( LP .NE. M) THEN
626                                   SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
627                                   IF (LP == 0) THEN
628                                      SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
629                                   ELSE
630                                      SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
631                                   ENDIF
632                                ENDIF
633                             ENDDO
634     
635                             DEN = A(L) + SUM_A_LPRIME + F(L,M)
636                             IF ( DEN .NE. ZERO) THEN
637                                SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
638                                SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
639                                        )/DEN
640                             ENDIF
641                          ELSE
642                             SUM_A = SUM_A + F(L,M)
643                             IF (L == 0) THEN
644                                SUM_B = SUM_B + F(L,M)*VAR_G(IJK)
645                             ELSE
646                                SUM_B = SUM_B + F(L,M)*VAR_S(IJK,L)
647                             ENDIF
648                          ENDIF   ! end if momentum_y_eq(l)
649                       ENDDO   !end do loop for L over 0 to mmax
650     
651                       A_M(IJK,0,M)=SUM_A+A(M)
652                       B_M(IJK,M) = SUM_B+BB(M)
653                    ENDIF    ! end if momentum_y_eq(M)
654                 ENDDO       ! end do loop for M over 0 to MMAX
655     
656              ENDIF      ! end if flow_at_n(ijk)
657           ENDDO         ! end do loop for ijk over ijkstart3,ijkend3
658     
659           RETURN
660           END SUBROUTINE PARTIAL_ELIM_V
661     
662     
663     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
664     !                                                                      C
665     !  Subroutine: PARTIAL_ELIM_W                                          C
666     !  Purpose: Do partial elimination for Z vector quantities
667     !                                                                      C
668     !                                                                      C
669     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
670     !  Reviewer:                                          Date:            C
671     !                                                                      C
672     !  Comments:                                                           C
673     !     Unlike gas-solids exchange coefficient (multiplied by volume),   C
674     !     the solids-solids exchange coefficient is not as a passed dummy  C
675     !     argument to this routine but simply re-evaluated here.           C
676     !                                                                      C
677     !  Literature/Document References:                                     C
678     !     Spalding, D.B., 1980, Numerical Computation of Multi-phase       C
679     !        fluid flow and heat transfer, Recent Advances in Numerical    C
680     !        Methods in Fluids, C. Taylor et al., eds, Pineridge Press     C
681     !     Syamlal, M., 1998, MFIX Documentation: Numerical Technique,      C
682     !        Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G      C
683     !        Technical Services of West Virginia, Inc., Morgantown, WV.    C
684     !                                                                      C
685     !                                                                      C
686     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
687     
688           SUBROUTINE PARTIAL_ELIM_W(VAR_G, VAR_S, VXF, A_M, B_M)
689     
690     !-----------------------------------------------
691     ! Modules
692     !-----------------------------------------------
693           USE param
694           USE param1
695           USE parallel
696           USE geometry
697           USE physprop
698           USE indices
699           USE run
700           USE compar
701           USE drag
702           USE fldvar
703           USE fun_avg
704           USE functions
705           IMPLICIT NONE
706     !-----------------------------------------------
707     ! Dummy arguments
708     !-----------------------------------------------
709     ! gas phase variable
710           DOUBLE PRECISION, INTENT(IN) :: Var_g(DIMENSION_3)
711     ! solids phase variable
712           DOUBLE PRECISION, INTENT(IN) :: Var_s(DIMENSION_3, DIMENSION_M)
713     ! Volume x gas-solids transfer coefficient
714           DOUBLE PRECISION, INTENT(IN) :: VxF(DIMENSION_3, DIMENSION_M)
715     ! Septadiagonal matrix A_m
716           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
717     ! Vector b_m
718           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
719     !-----------------------------------------------
720     ! Local variables
721     !-----------------------------------------------
722     ! Indices
723           INTEGER :: IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, K, IJKT
724           INTEGER :: L, M, LP, LM
725     !
726           DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
727     !                      a0, b0 etc.
728           DOUBLE PRECISION :: a(0:DIMENSION_M), BB(0:DIMENSION_M), &
729                               F(0:DIMENSION_M,0:DIMENSION_M),&
730                               Saxf(0:DIMENSION_M)
731     !-----------------------------------------------
732     
733     !!$omp  parallel do private( IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, &
734     !!$omp&  a, bb, F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, DEN) &
735     !!$omp&  schedule(static)
736           DO IJK = ijkstart3, ijkend3
737              IF (FLOW_AT_T(IJK)) THEN
738                 IMJK = IM_OF(IJK)
739                 IJMK = JM_OF(IJK)
740                 IPJK = IP_OF(IJK)
741                 IJPK = JP_OF(IJK)
742     
743                 F = ZERO
744                 DO M=0, MMAX
745                    A(M)=A_M(IJK,0,M)
746                    BB(M)=B_M(IJK,M)
747     
748     ! locally storing vol x avg gas-solids drag coefficient at momentum cell
749     ! center
750                    IF (M .NE. 0) THEN
751                       F(M,0)=-VXF(IJK,M)
752                       F(0,M)=F(M,0)
753                    ENDIF
754     
755     ! locally storing vol x avg solids-solids drag coefficient at momentum cell
756     ! center
757                    DO L =1, MMAX
758                       IF ((L .NE. M) .AND. (M .NE. 0)) THEN
759                          LM = FUNLM(L,M)
760                          IF (.NOT.IP_AT_T(IJK)) THEN
761                             K = K_OF(IJK)
762                             IJKT = TOP_OF(IJK)
763                             F(M,L) = -AVG_Z(F_SS(IJK,LM),F_SS(IJKT,LM),K)*VOL_W(IJK)
764                          ENDIF
765                       ENDIF
766                       F(L,M)=F(M,L)
767                    ENDDO
768     
769                    IF (M == 0) THEN
770                       SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IPJK)+A_M(IJK,west,M)*VAR_G(IMJK&
771                            )+A_M(IJK,north,M)*VAR_G(IJPK)+A_M(IJK,south,M)*VAR_G(IJMK))
772                    ELSE
773                       SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IPJK,M)+A_M(IJK,west,M)*VAR_S(&
774                            IMJK,M)+A_M(IJK,north,M)*VAR_S(IJPK,M)+A_M(IJK,south,M)*VAR_S(&
775                            IJMK,M))
776                    ENDIF
777     
778                    IF (DO_K) THEN
779                       IJKM = KM_OF(IJK)
780                       IJKP = KP_OF(IJK)
781                       IF (M == 0) THEN
782                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKP)+A_M(IJK,bottom,M)*&
783                                    VAR_G(IJKM))
784                       ELSE
785                          SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKP,M)+A_M(IJK,bottom,M)*&
786                                    VAR_S(IJKM,M))
787                       ENDIF
788                    ENDIF
789                 ENDDO     ! end do loop for M over 0 to MMAX
790     
791     
792                 DO M=0,MMAX
793                    IF (MOMENTUM_Z_EQ(M)) THEN
794                       SUM_A = ZERO
795                       SUM_B = ZERO
796     
797                       DO L =0,MMAX
798                          IF (MOMENTUM_Z_EQ(L)) THEN
799                             SUM_A_LPRIME = ZERO
800                             SUM_B_LPRIME = ZERO
801     
802                             DO LP=0,MMAX
803                                IF ( LP .NE. M) THEN
804                                   SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
805                                   IF (LP == 0) THEN
806                                      SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
807                                   ELSE
808                                      SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
809                                   ENDIF
810                                ENDIF
811                             ENDDO
812     
813                             DEN = A(L) + SUM_A_LPRIME + F(L,M)
814                             IF ( DEN .NE. ZERO) THEN
815                                SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
816                                SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
817                                )/DEN
818                             ENDIF
819                          ELSE
820                             SUM_A = SUM_A + F(L,M)
821                             IF (L == 0) THEN
822                                SUM_B = SUM_B + F(L,M)*VAR_G(IJK)
823                             ELSE
824                                SUM_B = SUM_B + F(L,M)*VAR_S(IJK,L)
825                             ENDIF
826                          ENDIF   ! end if momentum_z_eq(L)
827                       ENDDO   ! end do loop for L over 0 to mmax
828     
829                       A_M(IJK,0,M)=SUM_A+A(M)
830                       B_M(IJK,M) = SUM_B+BB(M)
831                    ENDIF    ! end if momentum_Z_eq(M)
832                 ENDDO       ! end do loop for M over 0 to MMAX
833     
834              ENDIF      ! end if flow_at_t(ijk)
835           ENDDO         ! end do loop for ijk over ijkstart3,ijkend3
836     
837           RETURN
838           END SUBROUTINE PARTIAL_ELIM_W
839     
840     
841