File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/calc_d_ghd.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_D_ghd_e(A_m, VxF_gs, d_e, IER)                    C
4     !  Purpose: calculte coefficients linking velocity correction to       C
5     !           pressure correction -- East                                C
6     !                                                                      C
7     !  Note:  MFIX convention: center coeff is negative, hence:            C
8     !                            (-A_M(IJK,0,M)) > or = 0                  C
9     !         The EAST face area is AYZ                                    C
10     !         Modifed for GHD theory                                       C
11     !                                                                      C
12     !  Author: M. Syamlal                                 Date: 21-JUN-96  C
13     !  Reviewer:                                          Date:            C
14     !                                                                      C
15     !  Revision Number: 2                                                  C
16     !  Purpose: allow multiparticle in D_E calculation in accounting for   C
17     !           the averaged Solid-Solid drag and multi-solid-particle dragC
18     !                                                                      C
19     !  Author: S. Dartevelle, LANL                        Date: 28-FEB-04  C
20     !  Reviewer:                                          Date: dd-mmm-yy  C
21     !                                                                      C
22     !  Revision Number: 3                                                  C
23     !  Purpose: more accurate formulas in the D_s(M) for Model_A           C
24     !                                                                      C
25     !  Author: S. Dartevelle, LANL                        Date: 17-MAR-04  C
26     !  Reviewer:                                          Date: dd-mmm-yy  C
27     !                                                                      C
28     !                                                                      C
29     !  Literature/Document References:                                     C
30     !                                                                      C
31     !  Variables referenced:                                               C
32     !  Variables modified:                                                 C
33     !                                                                      C
34     !  Local variables:                                                    C
35     !                                                                      C
36     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
37     !
38           SUBROUTINE CALC_D_ghd_E(A_M, VXF_GS, D_E)
39     !
40     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
41     !...Switches: -xf
42     !
43     !-----------------------------------------------
44     !   M o d u l e s
45     !-----------------------------------------------
46           USE param
47           USE param1
48           USE parallel
49           USE fldvar
50           USE geometry
51           USE indices
52           USE physprop
53           USE run
54           USE scales
55           USE compar
56           USE sendrecv
57           USE fun_avg
58           USE functions
59           IMPLICIT NONE
60     !-----------------------------------------------
61     !   G l o b a l   P a r a m e t e r s
62     !-----------------------------------------------
63     !-----------------------------------------------
64     !   D u m m y   A r g u m e n t s
65     !-----------------------------------------------
66     !
67     !                      Septadiagonal matrix A_m
68           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
69     !
70     !                      Volume x average at momentum cell centers
71           DOUBLE PRECISION VxF_gs(DIMENSION_3, DIMENSION_M)
72     !
73           DOUBLE PRECISION d_e(DIMENSION_3, 0:DIMENSION_M)
74     !
75     !                      Average volume fraction at momentum cell centers
76           DOUBLE PRECISION EPGA                                        !S. Dartevelle, LANL, Feb.2004
77     !          Usual Indices
78           INTEGER           M, L, I, IJK, IJKE                         !S. Dartevelle, LANL, Feb.2004
79     !          Average solid volume fraction at momentum cell centers
80           DOUBLE PRECISION  EPSA(DIMENSION_M)                          !S. Dartevelle, LANL, Feb.2004
81     !          ratio of drag and A0 or A_solid and other sum of Solid-Solid drag
82           DOUBLE PRECISION  other_ratio_1, FOA1                        !S. Dartevelle, LANL, Feb.2004
83     !          sum of All Solid-Gas VolxDrag
84           DOUBLE PRECISION  SUM_VXF_GS                                 !S. Dartevelle, LANL, Feb.2004
85     !          What phase momentum equation is activated?
86           LOGICAL  Pass1, Pass2                                        !S. Dartevelle, LANL, Feb.2004
87     !          numerator needed for solid phase M time EPSA
88           DOUBLE PRECISION  numeratorxEP(DIMENSION_M)                  !S. Dartevelle, LANL, Feb.2004
89     !          denominator needed for solid phase M
90           DOUBLE PRECISION  denominator(DIMENSION_M)                   !S. Dartevelle, LANL, Feb.2004
91     !          denominator needed for solid phase M
92           DOUBLE PRECISION  other_denominator(DIMENSION_M)             !S. Dartevelle, LANL, Feb.2004
93     !          total solids volume fraction
94           DOUBLE PRECISION  EPStmp
95     !-----------------------------------------------
96     
97        Pass1 = .FALSE.     !initialization
98        Pass2 = .FALSE.
99        M = MMAX
100              if (MOMENTUM_X_EQ(0) .AND. MOMENTUM_X_EQ(M)) then
101                Pass1 = .TRUE.    !we have at least one solid phase X-momentum equation
102                                  !with the gas phase X-momentum equation
103              elseif (MOMENTUM_X_EQ(M)) then
104                Pass2 = .TRUE.    !we have at least one solid phase X-momentum equation
105                                  !but the gas phase X-momentum is not solved
106              endif
107     
108     
109      IF (Pass1) THEN    !Gas and at least one solid phase X-momentum equation
110     !!!$omp   parallel do private(I,IJK, IJKE, EPGA, EPSA,EPStmp, numeratorxEP, &
111     !!!$omp&  M, L, Lp, LpL, LM, other_ratio_1, FOA1, &
112     !!!$omp&  SUM_VXF_GS, other_denominator, denominator ),&
113     !!!$omp&  schedule(static)
114       DO IJK = ijkstart3, ijkend3
115          IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN   !impermeable
116             DO M= 0, MMAX
117              D_E(IJK,M) = ZERO
118             END DO
119          ELSE
120             I = I_OF(IJK)
121             IJKE = EAST_OF(IJK)
122             EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
123     
124             SUM_VXF_GS = ZERO
125             EPSA(MMAX) = ZERO
126             DO M= 1, SMAX
127               EPSA(M) = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
128               EPSA(MMAX) = EPSA(MMAX) + EPSA(M)
129               SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)              !Gas - All Solids VolxDrag summation
130             END DO
131     
132             M = MMAX
133             other_ratio_1 = ( VXF_GS(IJK,M)* (-A_M(IJK,0,M)) /&
134                                           ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
135                             )
136     
137             IF (MODEL_B) THEN   !Model B
138               !Linking velocity correction coefficient to pressure - GAS Phase
139               if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR.  (other_ratio_1>SMALL_NUMBER) ) then
140                  D_E(IJK,0) = P_SCALE*AYZ(IJK)/( (-A_M(IJK,0,0))+other_ratio_1 )
141               else
142                  D_E(IJK,0) = ZERO
143               endif
144               !Linking velocity correction coefficient to pressure - SOLID Phase
145               M = MMAX
146                 if ( MOMENTUM_X_EQ(M) .AND. (  ((-A_M(IJK,0,M))>SMALL_NUMBER) .OR. &
147                                                          (VXF_GS(IJK,M)>SMALL_NUMBER) ) )   then
148                    D_E(IJK,M) = D_E(IJK,0)*(&
149                                     VXF_GS(IJK,M)/((-A_M(IJK,0,M))+VXF_GS(IJK,M))&
150                                 )
151                 else
152                    D_E(IJK,M) = ZERO
153                 endif
154             ELSE                !Model A
155               FOA1 = ZERO
156               M = MMAX
157                 FOA1 = FOA1 + (EPSA(M)*VXF_GS(IJK,M)/&
158                                  ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
159                               )
160                 other_denominator(M) = VXF_GS(IJK,M)*( (-A_M(IJK,0,0))/&
161                                                        ((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)&
162                                                      )
163                 numeratorxEP(M) = ZERO
164                 denominator(M)  = ZERO
165               !Linking velocity correction coefficient to pressure - GAS Phase
166               if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR.  (other_ratio_1>SMALL_NUMBER) ) then
167                 D_E(IJK,0) = P_SCALE*AYZ(IJK)*(EPGA+FOA1)/( (-A_M(IJK,0,0))+other_ratio_1 )
168               else
169                 D_E(IJK,0) = ZERO
170               endif
171               !Linking velocity correction coefficient to pressure - SOLID Phase
172               M = MMAX
173                 if ( MOMENTUM_X_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER)        .OR. &
174                                                     (other_denominator(M)>SMALL_NUMBER) ) )     then
175                   D_E(IJK,M) = P_SCALE*AYZ(IJK)*(&
176                             ( EPSA(M) + (VXF_GS(IJK,M)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)) )/&
177                             ( (-A_M(IJK,0,M))+other_denominator(M) )&
178                                                        )
179                 else
180                   D_E(IJK,M) = ZERO
181                 endif
182             ENDIF    !end of Model_B/Model_A if then condition
183          ENDIF
184       ENDDO
185     
186      ELSE IF (MOMENTUM_X_EQ(0)) THEN    !the solid X-momentum equations are not solved
187                                         !only gas phase X-momentum
188     !!!$omp   parallel do private(IJK, I, IJKE, EPGA, M, SUM_VXF_GS), &
189     !!!$omp&  schedule(static)
190        DO IJK = ijkstart3, ijkend3
191          IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN   !impermeable
192            D_E(IJK,0) = ZERO
193          ELSE
194            I = I_OF(IJK)
195            IJKE = EAST_OF(IJK)
196            EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
197     
198     
199            SUM_VXF_GS = VXF_GS(IJK,MMAX)              !Gas - All Solids VolxDrag summation
200     
201            IF ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (SUM_VXF_GS>SMALL_NUMBER) ) THEN
202              IF (MODEL_B) THEN
203                D_E(IJK,0) = P_SCALE*AYZ(IJK)/((-A_M(IJK,0,0))+SUM_VXF_GS)
204              ELSE
205                D_E(IJK,0) = P_SCALE*AYZ(IJK)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS)
206              ENDIF
207            ELSE
208              D_E(IJK,0) = ZERO
209            ENDIF
210          ENDIF
211        END DO
212     
213      ELSE IF (Pass2) THEN    !at least one solid phase X-momentum equation is solved
214                              !but the gas phase X-momentum is not solved
215     !!!$omp    parallel do private(IJK, I, IJKE, EPSA,EPStmp, L, Lp, M, &
216     !!!$omp&   numeratorxEP, denominator), &
217     !!!$omp&  schedule(static)
218        DO IJK = ijkstart3, ijkend3
219          IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK) .OR. MODEL_B) THEN
220            DO M= 1, MMAX
221              D_E(IJK,M) = ZERO
222            END DO
223          ELSE
224            I = I_OF(IJK)
225            IJKE = EAST_OF(IJK)
226     
227            M = MMAX
228            EPStmp = ZERO
229            DO L=1,SMAX
230               EPStmp = EPStmp+ AVG_X(EP_S(IJK,L),EP_S(IJKE,L),I)
231            ENDDO
232            EPSA(M) = EPStmp
233     
234            !Linking velocity correction coefficient to pressure - SOLID Phase (Model_A only)
235            M = MMAX
236              if ( MOMENTUM_X_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
237                                             (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
238                D_E(IJK,M) = P_SCALE*AYZ(IJK)*( EPSA(M) )/&
239                                               ( (-A_M(IJK,0,M))+VXF_GS(IJK,M) )
240              else
241                D_E(IJK,M) = ZERO
242              endif
243          ENDIF
244        END DO
245     
246      ENDIF
247     
248      RETURN
249      END SUBROUTINE CALC_D_ghd_E
250     
251     
252     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
253     !                                                                      C
254     !  Module name: CALC_D_ghd_n(A_m, VxF_gs, d_n, IER)                C
255     !  Purpose: calculte coefficients linking velocity correction to       C
256     !           pressure correction -- North                               C
257     !                                                                      C
258     !  Author: M. Syamlal                                 Date: 21-JUN-96  C
259     !  Reviewer:                                          Date:            C
260     !                                                                      C
261     !  Note:  MFIX convention: center coeff is negative, hence             C
262     !                            (-A_M(IJK,0,M)) > or = 0                  C
263     !         The NORTH face area is AXZ                                   C
264     !         Modifed for GHD theory                                       C
265     !                                                                      C
266     !  Revision Number: 2                                                  C
267     !  Purpose: allow multiparticle in D_N calculation in accounting for   C
268     !                       the averaged Solid-Solid drag                  C
269     !  Author: S. Dartevelle, LANL                        Date: 28-FEB-04  C
270     !  Reviewer:                                          Date: dd-mmm-yy  C
271     !                                                                      C
272     !  Revision Number: 3                                                  C
273     !  Purpose: more accurate formulas in the D_s(M) for Model_A           C
274     !                                                                      C
275     !  Author: S. Dartevelle, LANL                        Date: 17-MAR-04  C
276     !  Reviewer:                                          Date: dd-mmm-yy  C
277     !                                                                      C
278     !  Literature/Document References:                                     C
279     !                                                                      C
280     !  Variables referenced:                                               C
281     !  Variables modified:                                                 C
282     !                                                                      C
283     !  Local variables:                                                    C
284     !                                                                      C
285     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
286     !
287           SUBROUTINE CALC_D_ghd_N(A_M, VXF_GS, D_N)
288     !
289     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
290     !...Switches: -xf
291     !
292     !-----------------------------------------------
293     !   M o d u l e s
294     !-----------------------------------------------
295           USE param
296           USE param1
297           USE parallel
298           USE fldvar
299           USE geometry
300           USE indices
301           USE physprop
302           USE run
303           USE scales
304           USE compar
305           USE sendrecv
306           USE fun_avg
307           USE functions
308           IMPLICIT NONE
309     !-----------------------------------------------
310     !   G l o b a l   P a r a m e t e r s
311     !-----------------------------------------------
312     !-----------------------------------------------
313     !   D u m m y   A r g u m e n t s
314     !-----------------------------------------------
315     !
316     !                      Septadiagonal matrix A_m
317           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
318     !
319     !                      Volume x average at momentum cell centers
320           DOUBLE PRECISION VxF_gs(DIMENSION_3, DIMENSION_M)
321     !
322           DOUBLE PRECISION d_n(DIMENSION_3, 0:DIMENSION_M)
323     !
324     !                      Average volume fraction at momentum cell centers
325           DOUBLE PRECISION EPGA                                        !S. Dartevelle, LANL, Feb.2004
326     !          Usual Indices
327           INTEGER           M, L, J, IJK, IJKN                         !S. Dartevelle, LANL, Feb.2004
328     !          Average solid volume fraction at momentum cell centers
329           DOUBLE PRECISION  EPSA(DIMENSION_M)                          !S. Dartevelle, LANL, Feb.2004
330     !          ratio of drag and A0 or A_solid and other sum of Solid-Solid drag
331           DOUBLE PRECISION  other_ratio_1, FOA1                        !S. Dartevelle, LANL, Feb.2004
332     !          sum of All Solid-Gas VolxDrag
333           DOUBLE PRECISION  SUM_VXF_GS                                 !S. Dartevelle, LANL, Feb.2004
334     !          What phase momentum equation is activated?
335           LOGICAL  Pass1, Pass2                                        !S. Dartevelle, LANL, Feb.2004
336     !          numerator needed for solid phase M time EPSA
337           DOUBLE PRECISION  numeratorxEP(DIMENSION_M)                  !S. Dartevelle, LANL, Feb.2004
338     !          denominator needed for solid phase M
339           DOUBLE PRECISION  denominator(DIMENSION_M)                   !S. Dartevelle, LANL, Feb.2004
340     !          denominator needed for solid phase M
341           DOUBLE PRECISION  other_denominator(DIMENSION_M)             !S. Dartevelle, LANL, Feb.2004
342     !          total solids volume fraction
343           DOUBLE PRECISION  EPStmp
344     !-----------------------------------------------
345     
346        Pass1 = .FALSE.     !initialization
347        Pass2 = .FALSE.
348        M = MMAX
349          if (MOMENTUM_Y_EQ(0) .AND. MOMENTUM_Y_EQ(M)) then
350            Pass1 = .TRUE.    !we have at least one solid phase Y-momentum equation
351                              !with the gas phase Y-momentum equation
352          elseif (MOMENTUM_Y_EQ(M)) then
353            Pass2 = .TRUE.    !we have at least one solid phase Y-momentum equation
354                              !but the gas phase Y-momentum is not solved
355          endif
356     
357     
358      IF (Pass1) THEN    !Gas and at least one solid phases Y-momentum equation
359     !!!$omp   parallel do private(J,IJK, IJKN, EPGA, EPSA, EPStmp, numeratorxEP, &
360     !!!$omp&  M, L, Lp, LpL, LM, other_ratio_1, FOA1, &
361     !!!$omp&  SUM_VXF_GS, other_denominator, denominator ),&
362     !!!$omp&  schedule(static)
363        DO IJK = ijkstart3, ijkend3
364          IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN   !impermeable
365            DO M= 0, MMAX
366              D_N(IJK,M) = ZERO
367            END DO
368          ELSE
369            J = J_OF(IJK)
370            IJKN = NORTH_OF(IJK)
371            EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
372     
373            SUM_VXF_GS = ZERO
374            EPSA(MMAX) = ZERO
375            DO M= 1, SMAX
376              EPSA(M) = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
377              EPSA(MMAX) = EPSA(MMAX) + EPSA(M)
378              SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)              !Gas - All Solids VolxDrag summation
379            END DO
380     
381            M= MMAX
382            other_ratio_1 = ( VXF_GS(IJK,M)* (-A_M(IJK,0,M)) /&
383                                           ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
384                           )
385     
386            IF (MODEL_B) THEN   !Model B
387               !Linking velocity correction coefficient to pressure - GAS Phase
388              if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR.  (other_ratio_1>SMALL_NUMBER) ) then
389                D_N(IJK,0) = P_SCALE*AXZ(IJK)/( (-A_M(IJK,0,0))+other_ratio_1 )
390              else
391                D_N(IJK,0) = ZERO
392              endif
393              !Linking velocity correction coefficient to pressure - SOLID Phase
394              M = MMAX
395                if ( MOMENTUM_Y_EQ(M) .AND. (  ((-A_M(IJK,0,M))>SMALL_NUMBER) .OR. &
396                                                          (VXF_GS(IJK,M)>SMALL_NUMBER) ) )   then
397                  D_N(IJK,M) = D_N(IJK,0)*(&
398                                            VXF_GS(IJK,M)/((-A_M(IJK,0,M))+VXF_GS(IJK,M))&
399                                                    )
400                else
401                  D_N(IJK,M) = ZERO
402                endif
403            ELSE                !Model A
404              FOA1 = ZERO
405              M = MMAX
406                FOA1 = FOA1 + (EPSA(M)*VXF_GS(IJK,M)/&
407                                  ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
408                               )
409                other_denominator(M) = VXF_GS(IJK,M)*( (-A_M(IJK,0,0))/&
410                                                        ((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)&
411                                                      )
412                numeratorxEP(M) = ZERO
413                denominator(M)  = ZERO
414              !Linking velocity correction coefficient to pressure - GAS Phase
415              if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR.  (other_ratio_1>SMALL_NUMBER) ) then
416                D_N(IJK,0) = P_SCALE*AXZ(IJK)*(EPGA+FOA1)/( (-A_M(IJK,0,0))+other_ratio_1 )
417              else
418                D_N(IJK,0) = ZERO
419              endif
420              !Linking velocity correction coefficient to pressure - SOLID Phase
421              M = MMAX
422                if ( MOMENTUM_Y_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER)        .OR. &
423                                                     (other_denominator(M)>SMALL_NUMBER) ) )     then
424                  D_N(IJK,M) = P_SCALE*AXZ(IJK)*(&
425                             ( EPSA(M) + (VXF_GS(IJK,M)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)) )/&
426                             ( (-A_M(IJK,0,M))+other_denominator(M) )&
427                                                        )
428                else
429                  D_N(IJK,M) = ZERO
430                endif
431            ENDIF    !end of Model_B/Model_A if then condition
432          ENDIF
433        END DO
434     
435      ELSE IF (MOMENTUM_Y_EQ(0)) THEN    !the solid Y-momentum equations are not solved
436                                         !only the gas phase
437     !!!$omp   parallel do private(IJK, J, IJKN, EPGA, M, SUM_VXF_GS), &
438     !!!$omp&  schedule(static)
439        DO IJK = ijkstart3, ijkend3
440          IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN   !impermeable
441            D_N(IJK,0) = ZERO
442          ELSE
443            J = J_OF(IJK)
444            IJKN = NORTH_OF(IJK)
445            EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
446     
447     
448            SUM_VXF_GS = VXF_GS(IJK,MMAX)              !Gas - All Solids VolxDrag summation
449     
450            IF ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (SUM_VXF_GS>SMALL_NUMBER) ) THEN
451              IF (MODEL_B) THEN
452                D_N(IJK,0) = P_SCALE*AXZ(IJK)/((-A_M(IJK,0,0))+SUM_VXF_GS)
453              ELSE
454                D_N(IJK,0) = P_SCALE*AXZ(IJK)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS)
455              ENDIF
456            ELSE
457              D_N(IJK,0) = ZERO
458            ENDIF
459          ENDIF
460        END DO
461     
462      ELSE IF (Pass2) THEN    !at least one solid phase momentum Y-equation is solved
463                              !but the gas phase Y-momentum is not solved
464     !!!$omp    parallel do private(IJK, J, IJKN, EPSA, EPStmp, L, Lp, M, &
465     !!!$omp&   numeratorxEP, denominator), &
466     !!!$omp&  schedule(static)
467        DO IJK = ijkstart3, ijkend3
468          IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK) .OR. MODEL_B) THEN
469            DO M= 1, MMAX
470              D_N(IJK,M) = ZERO
471            END DO
472          ELSE
473            J = J_OF(IJK)
474            IJKN = NORTH_OF(IJK)
475     
476            M = MMAX
477            EPStmp = ZERO
478            DO L=1,SMAX
479              EPStmp = EPStmp+ AVG_Y(EP_S(IJK,L),EP_S(IJKN,L),J)
480            ENDDO
481            EPSA(M) = EPStmp
482     
483            !Linking velocity correction coefficient to pressure - SOLID Phase (Model_A only)
484            M = MMAX
485              if ( MOMENTUM_Y_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
486                                             (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
487                D_N(IJK,M) = P_SCALE*AXZ(IJK)*( EPSA(M) )/&
488                                               ( (-A_M(IJK,0,M))+VXF_GS(IJK,M) )
489              else
490                D_N(IJK,M) = ZERO
491              endif
492          ENDIF
493        END DO
494     
495      ENDIF
496     
497      RETURN
498      END SUBROUTINE CALC_D_ghd_N
499     
500     
501     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
502     !                                                                      C
503     !  Module name: CALC_D_ghd_t(A_m, VxF_gs, d_t, IER)                C
504     !  Purpose: calculte coefficients linking velocity correction to       C
505     !           pressure correction -- Top                                 C
506     !                                                                      C
507     !  Author: M. Syamlal                                 Date: 21-JUN-96  C
508     !  Reviewer:                                          Date:            C
509     !                                                                      C
510     !  Note:  MFIX convention: center coeff is negative, hence             C
511     !                            (-A_M(IJK,0,M)) > or = 0                  C
512     !         The TOP face area is AXY                                     C
513     !         Modifed for GHD theory                                       C
514     !                                                                      C
515     !  Revision Number: 2                                                  C
516     !  Purpose: allow multiparticle in D_T calculation in accounting for   C
517     !                       the averaged Solid-Solid drag                  C
518     !  Author: S. Dartevelle, LANL                        Date: 28-FEB-04  C
519     !  Reviewer:                                          Date: dd-mmm-yy  C
520     !                                                                      C
521     !  Revision Number: 3                                                  C
522     !  Purpose: more accurate formulas in the D_s(M) for Model_A           C
523     !                                                                      C
524     !  Author: S. Dartevelle, LANL                        Date: 17-MAR-04  C
525     !  Reviewer:                                          Date: dd-mmm-yy  C
526     !                                                                      C
527     !  Literature/Document References:                                     C
528     !                                                                      C
529     !  Variables referenced:                                               C
530     !  Variables modified:                                                 C
531     !                                                                      C
532     !  Local variables:                                                    C
533     !                                                                      C
534     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
535     !
536           SUBROUTINE CALC_D_ghd_T(A_M, VXF_GS, D_T)
537     !
538     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
539     !...Switches: -xf
540     !
541     !-----------------------------------------------
542     !   M o d u l e s
543     !-----------------------------------------------
544           USE param
545           USE param1
546           USE parallel
547           USE fldvar
548           USE geometry
549           USE indices
550           USE physprop
551           USE run
552           USE scales
553           USE compar
554           USE sendrecv
555           USE fun_avg
556           USE functions
557           IMPLICIT NONE
558     !-----------------------------------------------
559     !   G l o b a l   P a r a m e t e r s
560     !-----------------------------------------------
561     !-----------------------------------------------
562     !   D u m m y   A r g u m e n t s
563     !-----------------------------------------------
564     !
565     !                      Septadiagonal matrix A_m
566           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
567     !
568     !                      Volume x average at momentum cell centers
569           DOUBLE PRECISION VxF_gs(DIMENSION_3, DIMENSION_M)
570     !
571           DOUBLE PRECISION d_t(DIMENSION_3, 0:DIMENSION_M)
572     !
573     !                      Average volume fraction at momentum cell centers
574           DOUBLE PRECISION EPGA                                        !S. Dartevelle, LANL, Feb.2004
575     !          Usual Indices
576           INTEGER           M, L, K, IJK, IJKT                         !S. Dartevelle, LANL, Feb.2004
577     !          Average solid volume fraction at momentum cell centers
578           DOUBLE PRECISION  EPSA(DIMENSION_M)                          !S. Dartevelle, LANL, Feb.2004
579     !          ratio of drag and A0 or A_solid and other sum of Solid-Solid drag
580           DOUBLE PRECISION  other_ratio_1, FOA1                        !S. Dartevelle, LANL, Feb.2004
581     !          sum of All Solid-Gas VolxDrag
582           DOUBLE PRECISION  SUM_VXF_GS                                 !S. Dartevelle, LANL, Feb.2004
583     !          What phase momentum equation is activated?
584           LOGICAL  Pass1, Pass2                                        !S. Dartevelle, LANL, Feb.2004
585     !          numerator needed for solid phase M time EPSA
586           DOUBLE PRECISION  numeratorxEP(DIMENSION_M)                  !S. Dartevelle, LANL, Feb.2004
587     !          denominator needed for solid phase M
588           DOUBLE PRECISION  denominator(DIMENSION_M)                   !S. Dartevelle, LANL, Feb.2004
589     !          denominator needed for solid phase M
590           DOUBLE PRECISION  other_denominator(DIMENSION_M)             !S. Dartevelle, LANL, Feb.2004
591     !          tmp variable for total solids volume fraction
592           DOUBLE PRECISION  EPStmp
593     !-----------------------------------------------
594     
595        Pass1 = .FALSE.     !initialization
596        Pass2 = .FALSE.
597        M = MMAX
598          if (MOMENTUM_Z_EQ(0) .AND. MOMENTUM_Z_EQ(M)) then
599            Pass1 = .TRUE.    !we have at least one solid phase Z-momentum equation
600                              !with the gas phase Z-momentum equation
601          elseif (MOMENTUM_Z_EQ(M)) then
602            Pass2 = .TRUE.    !we have at least one solid phase Z-momentum equation
603                              !but the gas phase Z-momentum is not solved
604          endif
605     
606     
607      IF (Pass1) THEN    !Gas and at least one solid phases Z-momentum equation
608     !!!$omp   parallel do private(K,IJK, IJKT, EPGA, EPSA, EPStmp, numeratorxEP, &
609     !!!$omp&  M, L, Lp, LpL, LM, other_ratio_1, FOA1,  &
610     !!!$omp&  SUM_VXF_GS, other_denominator, denominator ),&
611     !!!$omp&  schedule(static)
612        DO IJK = ijkstart3, ijkend3
613          IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN   !impermeable
614            DO M= 0, MMAX
615              D_T(IJK,M) = ZERO
616            END DO
617          ELSE
618            K = K_OF(IJK)
619            IJKT = TOP_OF(IJK)
620            EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
621     
622            SUM_VXF_GS = ZERO
623            EPSA(MMAX) = ZERO
624            DO M= 1, SMAX
625              EPSA(M) = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
626              EPSA(MMAX) = EPSA(MMAX) + EPSA(M)
627              SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)              !Gas - All Solids VolxDrag summation
628            END DO
629     
630            M = MMAX
631            other_ratio_1 = ( VXF_GS(IJK,M)* (-A_M(IJK,0,M)) /&
632                                           ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
633                           )
634     
635            IF (MODEL_B) THEN   !Model B
636              !Linking velocity correction coefficient to pressure - GAS Phase
637              if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR.  (other_ratio_1>SMALL_NUMBER) ) then
638                D_T(IJK,0) = P_SCALE*AXY(IJK)/( (-A_M(IJK,0,0))+other_ratio_1 )
639              else
640                D_T(IJK,0) = ZERO
641              endif
642              !Linking velocity correction coefficient to pressure - SOLID Phase
643              M = MMAX
644                if ( MOMENTUM_Z_EQ(M) .AND. (  ((-A_M(IJK,0,M))>SMALL_NUMBER) .OR. &
645                                                          (VXF_GS(IJK,M)>SMALL_NUMBER) ) )   then
646                  D_T(IJK,M) = D_T(IJK,0)*(&
647                                            VXF_GS(IJK,M)/((-A_M(IJK,0,M))+VXF_GS(IJK,M))&
648                                                    )
649                else
650                  D_T(IJK,M) = ZERO
651                endif
652            ELSE                !Model A
653              FOA1 = ZERO
654              M = MMAX
655                FOA1 = FOA1 + (EPSA(M)*VXF_GS(IJK,M)/&
656                                  ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
657                               )
658                other_denominator(M) = VXF_GS(IJK,M)*( (-A_M(IJK,0,0))/&
659                                                        ((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)&
660                                                      )
661                numeratorxEP(M) = ZERO
662                denominator(M)  = ZERO
663              !Linking velocity correction coefficient to pressure - GAS Phase
664              if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR.  (other_ratio_1>SMALL_NUMBER) ) then
665                D_T(IJK,0) = P_SCALE*AXY(IJK)*(EPGA+FOA1)/( (-A_M(IJK,0,0))+other_ratio_1 )
666              else
667                D_T(IJK,0) = ZERO
668              endif
669              !Linking velocity correction coefficient to pressure - SOLID Phase
670              M = MMAX
671                if ( MOMENTUM_Z_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER)        .OR. &
672                                                     (other_denominator(M)>SMALL_NUMBER) ) )     then
673                  D_T(IJK,M) = P_SCALE*AXY(IJK)*(&
674                             ( EPSA(M) + (VXF_GS(IJK,M)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)) )/&
675                             ( (-A_M(IJK,0,M))+other_denominator(M) )&
676                                                        )
677                else
678                  D_T(IJK,M) = ZERO
679                endif
680            ENDIF    !end of Model_B/Model_A if then condition
681          ENDIF
682        END DO
683     
684      ELSE IF (MOMENTUM_Z_EQ(0)) THEN    !the solid Z-momentum equations are not solved
685                                         !only the gas phase Z-momentum is solved
686     !!!$omp   parallel do private(IJK, K, IJKT, EPGA, M, SUM_VXF_GS), &
687     !!!$omp&  schedule(static)
688        DO IJK = ijkstart3, ijkend3
689          IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN   !impermeable
690            D_T(IJK,0) = ZERO
691          ELSE
692            K = K_OF(IJK)
693            IJKT = TOP_OF(IJK)
694            EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
695     
696            SUM_VXF_GS = VXF_GS(IJK,MMAX)              !Gas - All Solids VolxDrag summation
697     
698            IF ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (SUM_VXF_GS>SMALL_NUMBER) ) THEN
699              IF (MODEL_B) THEN
700                D_T(IJK,0) = P_SCALE*AXY(IJK)/((-A_M(IJK,0,0))+SUM_VXF_GS)
701              ELSE
702                D_T(IJK,0) = P_SCALE*AXY(IJK)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS)
703              ENDIF
704            ELSE
705              D_T(IJK,0) = ZERO
706            ENDIF
707          ENDIF
708        END DO
709     
710      ELSE IF (Pass2) THEN    !at least one solid phase momentum Z-equation is solved
711                              !but the gas phase Z-momentum is not solved
712     !!!$omp    parallel do private(IJK, K, IJKT, EPSA, EPStmp, L, Lp, M, &
713     !!!$omp&   numeratorxEP, denominator), &
714     !!!$omp&  schedule(static)
715        DO IJK = ijkstart3, ijkend3
716          IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK) .OR. MODEL_B) THEN
717            DO M= 1, MMAX
718              D_T(IJK,M) = ZERO
719            END DO
720          ELSE
721            K = K_OF(IJK)
722            IJKT = TOP_OF(IJK)
723     
724            M = MMAX
725            EPStmp = ZERO
726            DO L=1,SMAX
727              EPStmp = EPStmp+ AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
728            ENDDO
729            EPSA(M) = EPStmp
730     
731            !Linking velocity correction coefficient to pressure - SOLID Phase (Model_A only)
732            M = MMAX
733              if ( MOMENTUM_Z_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
734                                             (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
735                 D_T(IJK,M) = P_SCALE*AXY(IJK)*( EPSA(M) )/&
736                                               ( (-A_M(IJK,0,M))+VXF_GS(IJK,M) )
737              else
738                D_T(IJK,M) = ZERO
739              endif
740          ENDIF
741        END DO
742     
743      ENDIF
744     
745      RETURN
746      END SUBROUTINE CALC_D_ghd_T
747     
748     !// Comments on the modifications for DMP version implementation
749     !// 001 Include header file and common declarations for parallelization
750     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
751