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