File: RELATIVE:/../../../mfix.git/model/calc_d_e.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CALC_d_e                                                !
4     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
5     !                                                                      !
6     !  Purpose: calculate coefficients linking velocity correction to      !
7     !  pressure correction -- East                                         !
8     !                                                                      !
9     !  Note that the D_E coefficients for phases M>0 are generally not     !
10     !  used unless the solids phase has close_packed=F, in which case the  !
11     !  D_E coefficients for that phase are employed in a mixture pressure  !
12     !  correction equation and for correcting velocities.                  !
13     !                                                                      !
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
15           SUBROUTINE CALC_D_E(A_M, VXF_GS, VXF_SS, D_E, IER)
16     
17     ! Global Variables:
18     !---------------------------------------------------------------------//
19     ! Number of solids phases.
20           use physprop, only: MMAX
21     ! Flag: Solve the X-Momentum Equations
22           use run, only: MOMENTUM_X_EQ
23     ! Flag: Coupled DEM, PIC, or Hybrid simulation
24           use discretelement, only: DES_CONTINUUM_COUPLED
25     ! Flag: TMF-DEM solids hybrid model
26           use discretelement, only: DES_CONTINUUM_HYBRID
27     ! Volume x average at momentum cell center drag for DEM/PIC
28           use discretelement, only: VXF_GDS, VXF_SDS
29     ! Number of solids phases.
30           use physprop, only: MMAX
31     
32     ! Global Parameters:
33     !---------------------------------------------------------------------//
34     ! Size of IJK arrays and size of solids phase arrays.
35           use param, only: DIMENSION_3, DIMENSION_M
36     ! Size of MxM upper triangular matrix.
37           use param1, only: DIMENSION_LM
38     
39           IMPLICIT NONE
40     
41     ! Dummy arguments
42     !---------------------------------------------------------------------//
43     ! Error index
44           INTEGER, INTENT(INOUT) :: IER
45     ! Septadiagonal matrix A_m.  The center coefficient is negative.
46           DOUBLE PRECISION, INTENT(IN):: A_m(DIMENSION_3,-3:3,0:DIMENSION_M)
47     ! Volume x average at momentum cell centers
48           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
49     ! Volume x average at momentum cell centers Solid-Solid Drag
50           DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
51     ! Coefficients for pressure correction equation. These coefficients are
52     ! initialized to zero in the subroutine time_march before the time loop.
53           DOUBLE PRECISION, INTENT(INOUT) :: d_e(DIMENSION_3, 0:DIMENSION_M)
54     
55     
56     ! Local variables:
57     !---------------------------------------------------------------------//
58     ! Numerator and demoninator needed for evaluation of pressure correction
59     ! coefficient for GAS phase. A temporary variable is used for local
60     ! manipulations.
61           DOUBLE PRECISION, dimension(:,:), allocatable :: AM0
62     ! Flag: One or more solids momentum equations are solved.
63           LOGICAL :: ANY_SOLIDS_X_MOMENTUM
64     
65     !......................................................................!
66     
67           allocate(AM0(DIMENSION_3, 0:DIMENSION_M))
68     
69     ! Initialize the error flag.
70           IER = 0
71     
72     ! Copy the gas phase A_M coefficients into temp array.
73           AM0(:,:) = A_M(:,0,:)
74     
75     ! Add DEM temp A_M so they are included in the pressure correciton eq.
76           IF(DES_CONTINUUM_COUPLED) THEN
77              AM0(:,0) = AM0(:,0) - VXF_GDS(:)
78              IF (DES_CONTINUUM_HYBRID) &
79                 AM0(:,1:MMAX) = AM0(:,1:MMAX) - VXF_SDS(:,1:MMAX)
80           ENDIF
81     
82           ANY_SOLIDS_X_MOMENTUM = any(MOMENTUM_X_EQ(1:MMAX))
83     
84           IF(MOMENTUM_X_EQ(0)) THEN
85              IF(ANY_SOLIDS_X_MOMENTUM) THEN
86                 CALL CALC_D_E_GAS_AND_SOLIDS(AM0, VXF_GS, VXF_SS, D_E)
87              ELSE
88                 CALL CALC_D_E_GAS_ONLY(AM0, VXF_GS, D_E)
89              ENDIF
90     
91           ELSEIF(ANY_SOLIDS_X_MOMENTUM) THEN
92              CALL CALC_D_E_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_E)
93           ENDIF
94     
95           deallocate(AM0)
96     
97           RETURN
98           END SUBROUTINE CALC_D_E
99     
100     
101     
102     
103     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
104     !                                                                      !
105     !  Subroutine: CALC_D_E_GAS_AND_SOLIDS                                 !
106     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
107     !                                                                      !
108     !  Purpose: calculate coefficients linking velocity correction to      !
109     !  pressure correction -- East                                         !
110     !                                                                      !
111     !  CASE: The gas phase X-momentum equation is solved and at least one  !
112     !    solids phase (m>0) X-momentum equation is solved.                 !
113     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
114           SUBROUTINE CALC_D_E_GAS_AND_SOLIDS(AM0, VXF_GS, VXF_SS, D_E)
115     
116     ! Global Variables:
117     !---------------------------------------------------------------------//
118     ! Volume fractions of gas and solids phases.
119           use fldvar, only: EP_G, EP_s
120     ! Number of solids phases.
121           use physprop, only: MMAX
122     ! Flag: Solve the X-Momentum Equations
123           use run, only: MOMENTUM_X_EQ
124     ! Flag: Use MODEL_B
125           use run, only: MODEL_B
126     ! Pressure scale factor
127           use scales, only: P_SCALE
128     ! Flag: Cartesian grid simulation
129           use cutcell, only: CARTESIAN_GRID
130     ! Area of cell's XZ face.
131           use geometry, only: AYZ
132     ! Function to average across X face.
133           use fun_avg, only: AVG_X
134     ! Fluid grid loop bounds.
135           use compar, only: IJKStart3, IJKEnd3
136     ! Function to convert to phases to single array, IJK of cell to east
137           use functions, only: FUNLM, EAST_OF
138     ! Flags: Impermeable surface and mass flow at east face
139           use functions, only: IP_AT_E, MFLOW_AT_E
140     ! Indices: I index of cell
141           use indices, only: I_OF
142     
143     ! Global Parameters:
144     !---------------------------------------------------------------------//
145     ! Size of IJK arrays and size of solids phase arrays.
146           use param, only: DIMENSION_3, DIMENSION_M
147     ! Size of MxM upper triangular matrix.
148           use param1, only: DIMENSION_LM
149     ! Double precision parameters.
150           use param1, only: ZERO, SMALL_NUMBER, ONE
151     
152           IMPLICIT NONE
153     
154     ! Dummy arguments
155     !---------------------------------------------------------------------//
156     ! Temporary variable to store A matrix for local manipulation
157           DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
158     ! Volume x average at momentum cell centers
159           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
160     ! Volume x average at momentum cell centers Solid-Solid Drag
161           DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
162     ! Coefficients for pressure correction equation. These coefficients are
163     ! initialized to zero in the subroutine time_march before the time loop.
164           DOUBLE PRECISION, INTENT(INOUT) :: D_E(DIMENSION_3, 0:DIMENSION_M)
165     
166     ! Local variables:
167     !---------------------------------------------------------------------//
168     ! Usual Indices
169           INTEGER :: LM, M, L, LPL, LP
170           INTEGER :: I, IJK, IJKE
171     ! Average solid volume fraction at momentum cell centers
172           DOUBLE PRECISION :: EPSA(DIMENSION_M)
173     ! Average volume fraction at momentum cell centers
174           DOUBLE PRECISION :: EPGA
175     ! local alias for face area
176           DOUBLE PRECISION :: AREA_FACE
177     ! sum of All Solid-Gas VolxDrag
178           DOUBLE PRECISION :: SUM_VXF_GS
179     ! sum of Solid M - All other Solid drag
180           DOUBLE PRECISION :: SUM_VXF_SS(DIMENSION_M)
181     ! sum of Solid L - All other Solid VolxDrag but M
182           DOUBLE PRECISION :: SUM_VXF_SS_wt_M
183     ! numerator and demoninator needed for evaluation of pressure correction
184     ! coefficient for GAS phase
185           DOUBLE PRECISION :: DEN_MGas, NUM_MGas
186     ! component of the total numerator needed for evaluation of pressure
187     ! correction coefficient for SOLIDS phase M. specifically, for when the
188     ! sum L is over solids phase only
189           DOUBLE PRECISION :: NUM_MSol_LSol(DIMENSION_M)
190     ! component of the total denominator needed for evaluation of pressure
191     ! correction coefficient for SOLIDS phase M.  specifically, for when the
192     ! sum L is over solids phases only
193           DOUBLE PRECISION :: DEN_MSol_LSol(DIMENSION_M)
194     ! component of the total numerator needed for evaluation of pressure
195     ! correction coefficient for SOLIDS phase M.  specifically, for when the
196     ! sum L is over gas phase only
197           DOUBLE PRECISION :: NUM_MSol_LGas
198     ! component of the total denominator needed for evaluation of pressure
199     ! correction coefficient for SOLIDS phase M.  specifically, for when the
200     ! sum L is over gas phase only
201           DOUBLE PRECISION :: DEN_MSol_LGas
202     ! Temp variable for double precision values.
203           DOUBLE PRECISION :: TMPdp
204     
205     !......................................................................!
206     
207     !$omp parallel default(none) &
208     !$omp          private(ijk,area_face,i,ijke,epga,sum_vxf_gs,epsa,lm,sum_vxf_ss,den_mgas,num_mgas,tmpdp, &
209     !$omp                  num_msol_lgas,den_msol_lgas,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
210     !$omp          shared(ijkstart3,ijkend3,d_e,mmax,cartesian_grid,ayz,i_of,ep_g,vxf_gs,vxf_ss,momentum_x_eq,am0,model_b,p_scale)
211     !$omp do
212           DO IJK = ijkstart3, ijkend3
213     
214              IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN
215                 DO M= 0, MMAX
216                    D_E(IJK,M) = ZERO
217                 ENDDO
218                 CYCLE
219              ENDIF
220     
221              AREA_FACE = merge(ONE, AYZ(IJK), CARTESIAN_GRID)
222              I = I_OF(IJK)
223              IJKE = EAST_OF(IJK)
224              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
225     
226              SUM_VXF_GS = ZERO
227              DO M= 1, MMAX
228                 EPSA(M) = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
229     ! Gas - All Solids VolxDrag summation
230                 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
231                 SUM_VXF_SS(M) = ZERO
232                 DO L = 1, MMAX
233                    IF (L .NE. M) THEN
234                       LM = FUNLM(L,M)
235     ! Solids M - All other Solids VolxDrag summation
236                       SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
237                    ENDIF
238                 ENDDO
239              ENDDO
240     
241     ! calculating DEN_MGas and NUM_MGas
242              DEN_MGas  = ZERO
243              NUM_MGas = ZERO
244              DO M= 1, MMAX
245                 IF (MOMENTUM_X_EQ(M)) THEN
246                    NUM_MGas = NUM_MGas + (EPSA(M)*VXF_GS(IJK,M)/           &
247                      ((-AM0(IJK,M)) + VXF_GS(IJK,M) + SUM_VXF_SS(M) +      &
248                      SMALL_NUMBER))
249                    DEN_MGas = DEN_MGas + (VXF_GS(IJK,M)*                   &
250                      ((-AM0(IJK,M))+SUM_VXF_SS(M))/                        &
251                      ((-AM0(IJK,M))+VXF_GS(IJK,M)+SUM_VXF_SS(M)+           &
252                      SMALL_NUMBER))
253                 ELSE
254                    DEN_MGas = DEN_MGas + VXF_GS(IJK,M)
255                 ENDIF
256              ENDDO
257     
258     ! Model B
259     ! -------------------------------------------------------------------
260              IF (MODEL_B) THEN
261     ! Linking velocity correction coefficient to pressure - GAS Phase
262                 TMPdp = -AM0(IJK,0) + DEN_MGas
263                 IF(abs(TMPdp) > SMALL_NUMBER) THEN
264                    D_E(IJK,0) = P_SCALE*AREA_FACE/TMPdp
265                 ELSE
266                    D_E(IJK,0) = ZERO
267                 ENDIF
268     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
269                 DO M = 1, MMAX
270                    IF(MOMENTUM_X_EQ(M)) THEN
271                       TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M)
272                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
273                          D_E(IJK,M) = D_E(IJK,0)*VXF_GS(IJK,M)/TMPdp
274                       ELSE
275                          D_E(IJK,M) = ZERO
276                       ENDIF
277                    ELSE
278                       D_E(IJK,M) = ZERO
279                    ENDIF
280                 ENDDO
281     
282     ! Model A
283     ! -------------------------------------------------------------------
284              ELSE
285     
286     ! Linking velocity correction coefficient to pressure - GAS Phase
287                 TMPdp = -AM0(IJK,0) + DEN_MGas
288                 IF(abs(TMPdp) >SMALL_NUMBER) THEN
289                    D_E(IJK,0) = P_SCALE*AREA_FACE*(EPGA+NUM_MGas)/TMPdp
290                 ELSE
291                    D_E(IJK,0) = ZERO
292                 ENDIF
293     
294                 DO M= 1, MMAX
295     ! calculating NUM_MSol_LGas and DEN_MSol_LGas
296                    NUM_MSol_LGas = VXF_GS(IJK,M)*EPGA/                     &
297                       ((-AM0(IJK,0))+SUM_VXF_GS+SMALL_NUMBER)
298                    DEN_MSol_LGas = VXF_GS(IJK,M)*(                         &
299                       ((-AM0(IJK,0)) + (SUM_VXF_GS - VXF_GS(IJK,M)))/      &
300                       ((-AM0(IJK,0)) + SUM_VXF_GS + SMALL_NUMBER) )
301     
302     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
303                    NUM_MSol_LSol(M) = ZERO
304                    DEN_MSol_LSol(M)  = ZERO
305                    DO L = 1, MMAX
306                       IF (L /= M) THEN
307                          LM = FUNLM(L,M)
308                          SUM_VXF_SS_wt_M = ZERO
309                          DO Lp = 1, MMAX
310                             IF((Lp /= L) .AND. (Lp /= M) ) THEN
311                                LpL = FUNLM(Lp,L)
312     ! Solids L - All other Solids VolxDrag but M summation
313                                SUM_VXF_SS_wt_M =                           &
314                                   SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
315                             ENDIF
316                          ENDDO
317     
318                          IF(MOMENTUM_X_EQ(L)) THEN
319                             NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +          &
320                                ( VXF_SS(IJK,LM)*EPSA(L)/                   &
321                                ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
322                                SMALL_NUMBER) )
323     
324                             DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +          &
325                                VXF_SS(IJK,LM)*(((-AM0(IJK,L))+             &
326                                VXF_GS(IJK,L)+SUM_VXF_SS_wt_M)/             &
327                                ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
328                                SMALL_NUMBER) )
329                          ELSE
330                             DEN_MSol_LSol(M) =                             &
331                                DEN_MSol_LSol(M) + VXF_SS(IJK,LM)
332                          ENDIF
333                       ENDIF  ! end if (l.ne.m)
334                    ENDDO  ! end do (l=1,mmax)
335     
336     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
337                    IF(MOMENTUM_X_EQ(M)) THEN
338                       TMPdp = -AM0(IJK,M) + DEN_MSol_LGas + DEN_MSol_LSol(M)
339                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
340                          D_E(IJK,M) = P_SCALE*AREA_FACE * (EPSA(M) +       &
341                             NUM_MSol_LSol(M) + NUM_MSol_LGas)/TMPdp
342                       ELSE
343                          D_E(IJK,M) = ZERO
344                       ENDIF
345                    ELSE
346                       D_E(IJK,M) = ZERO
347                    ENDIF
348                 ENDDO   ! end do (m=1,mmax)
349     
350              ENDIF    !end if/else branch Model_B/Model_A
351     
352           ENDDO  ! end do loop (ijk=ijkstart3,ijkend3)
353     !$omp end parallel
354     
355           RETURN
356           END SUBROUTINE CALC_D_E_GAS_AND_SOLIDS
357     
358     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
359     !                                                                      !
360     !  Subroutine: CALC_D_E_GAS_ONLY                                       !
361     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
362     !                                                                      !
363     !  Purpose: calculate coefficients linking velocity correction to      !
364     !  pressure correction -- East                                         !
365     !                                                                      !
366     !  Notes: The gas phase x-momentum equation is solved; NO solids phase !
367     !  x-momentum equations are solved.                                    !
368     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
369     
370           SUBROUTINE CALC_D_E_GAS_ONLY(AM0, VXF_GS, D_E)
371     
372     ! Global Variables:
373     !---------------------------------------------------------------------//
374     ! Volume fraction of gas phase.
375           use fldvar, only: EP_G
376     ! Number of solids phases.
377           use physprop, only: MMAX
378     ! Flag: Use MODEL_B
379           use run, only: MODEL_B
380     ! Pressure scale factor
381           use scales, only: P_SCALE
382     ! Flag: Cartesian grid simulation
383           use cutcell, only: CARTESIAN_GRID
384     ! Area of cell's YZ face.
385           use geometry, only: AYZ
386     ! Volume of V-momentum cell.
387           use geometry, only: VOL_U
388     ! Function to average across Y face.
389           use fun_avg, only: AVG_X
390     ! Fluid grid loop bounds.
391           use compar, only: IJKStart3, IJKEnd3
392     ! Flags: Impermeable surface and mass flow at east face, IJK of cell to east
393           use functions, only: IP_AT_E, MFLOW_AT_E, EAST_OF
394     ! Indices: J index of cell
395           use indices, only: I_OF
396     ! Flag and variables for QMOM implementation.
397           USE qmom_kinetic_equation, only: QMOMK, QMOMK_NN
398           USE qmom_kinetic_equation, only: QMOMK_F_GS, QMOMK_F_GS
399     
400     ! Global Parameters:
401     !---------------------------------------------------------------------//
402     ! Size of IJK arrays and size of solids phase arrays.
403           use param, only: DIMENSION_3, DIMENSION_M
404     ! Double precision parameters.
405           use param1, only: ZERO, SMALL_NUMBER, ONE
406     
407           IMPLICIT NONE
408     
409     ! Dummy arguments
410     !---------------------------------------------------------------------//
411     ! Temporary variable to store A matrix for local manipulation
412           DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
413     ! Volume x average at momentum cell centers
414           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
415     ! Coefficients for pressure correction equation. These coefficients are
416     ! initialized to zero in the subroutine time_march before the time loop.
417           DOUBLE PRECISION, INTENT(INOUT) :: D_E(DIMENSION_3, 0:DIMENSION_M)
418     
419     ! Local variables:
420     !---------------------------------------------------------------------//
421     ! Usual Indices
422           INTEGER :: M
423           INTEGER :: INN
424           INTEGER :: I, IJK, IJKE
425     ! Average solid volume fraction at momentum cell centers
426           !DOUBLE PRECISION :: EPSA(DIMENSION_M)
427     ! Average volume fraction at momentum cell centers
428           DOUBLE PRECISION :: EPGA
429     ! local alias for face area
430           DOUBLE PRECISION :: AREA_FACE
431     ! sum of All Solid-Gas VolxDrag
432           DOUBLE PRECISION :: SUM_VXF_GS
433     ! Temp variable for double precision values.
434           DOUBLE PRECISION :: TMPdp
435     
436     !......................................................................!
437     
438     !$omp parallel default(none) &
439     !$omp          private(ijk,area_face,i,ijke,epga,sum_vxf_gs,tmpdp) &
440     !$omp          shared(ijkstart3,ijkend3,d_e,mmax,cartesian_grid,ayz,i_of,ep_g,vxf_gs,am0,model_b,p_scale,qmomk,vol_u,qmomk_f_gs)
441     !$omp do
442           DO IJK = IJKSTART3, IJKEND3
443     
444              IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN   !impermeable
445                 D_E(IJK,0) = ZERO
446                 CYCLE
447              ENDIF
448     
449              AREA_FACE = merge(ONE, AYZ(IJK), CARTESIAN_GRID)
450     
451              I = I_OF(IJK)
452              IJKE = EAST_OF(IJK)
453              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
454     
455              SUM_VXF_GS = ZERO
456              IF (.NOT. QMOMK) THEN
457                 DO M= 1, MMAX
458     ! Gas - All Solids VolxDrag summation
459                    SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
460                 ENDDO
461              ELSE
462                 DO INN = 1, QMOMK_NN
463                    DO M = 1, MMAX
464                       SUM_VXF_GS = SUM_VXF_GS + VOL_U(IJK)* &
465                          AVG_X(QMOMK_F_GS(INN,IJK,M),QMOMK_F_GS(INN,IJKE,M),I)
466                    ENDDO
467                 ENDDO
468              ENDIF
469     
470              TMPdp = -AM0(IJK,0) + SUM_VXF_GS
471              IF(abs(TMPdp) > SMALL_NUMBER) THEN
472                 IF (MODEL_B) THEN
473                    D_E(IJK,0) = P_SCALE*AREA_FACE/TMPdp
474                 ELSE
475                    D_E(IJK,0) = P_SCALE*AREA_FACE*EPGA/TMPdp
476                 ENDIF   !end if/else branch Model_B/Model_A
477              ELSE
478                 D_E(IJK,0) = ZERO
479              ENDIF
480     
481           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
482     !$omp end parallel
483     
484           RETURN
485           END SUBROUTINE CALC_D_E_GAS_ONLY
486     
487     
488     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
489     !                                                                      !
490     !  Subroutine: CALC_D_E_SOLIDS_ONLY                                    !
491     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
492     !                                                                      !
493     !  Purpose: Calculate coefficients linking velocity correction to      !
494     !  pressure correction -- East                                         !
495     !                                                                      !
496     !  Notes: The gas phase momentum equations are NOT solved but at least !
497     !  one solids phase momentum equation is solved.                       !
498     !                                                                      !
499     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
500           SUBROUTINE CALC_D_E_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_E)
501     
502     ! Global Variables:
503     !---------------------------------------------------------------------//
504     ! Volume fractions of gas and solids phases.
505          use fldvar, only: EP_s
506     ! Number of solids phases.
507          use physprop, only: MMAX
508     ! Flag: Solve the X-Momentum Equations
509          use run, only: MOMENTUM_X_EQ
510     ! Flag: Use MODEL_B
511          use run, only: MODEL_B
512     ! Pressure scale factor
513          use scales, only: P_SCALE
514     ! Flag: Cartesian grid simulation
515           use cutcell, only: CARTESIAN_GRID
516     ! Area of cell's YZ face.
517           use geometry, only: AYZ
518     ! Function to average across X face.
519           use fun_avg, only: AVG_X
520     ! Fluid grid loop bounds.
521           use compar, only: IJKStart3, IJKEnd3
522     ! Function to convert to phases to single array, IJK of cell to east
523           use functions, only: FUNLM, EAST_OF
524     ! Flags: Impermeable surface and mass flow at east face
525           use functions, only: IP_AT_E, MFLOW_AT_E
526     ! Indices: I index of cell
527           use indices, only: I_OF
528     
529     ! Global Parameters:
530     !---------------------------------------------------------------------//
531     ! Size of IJK arrays and size of solids phase arrays.
532           use param, only: DIMENSION_3, DIMENSION_M
533     ! Size of MxM upper triangular matrix.
534           use param1, only: DIMENSION_LM
535     ! Double precision parameters.
536           use param1, only: ZERO, SMALL_NUMBER, ONE
537     
538           IMPLICIT NONE
539     
540     ! Dummy arguments
541     !---------------------------------------------------------------------//
542     ! Error index
543     !      INTEGER, INTENT(INOUT) :: IER
544     ! Temporary variable to store A matrix for local manipulation
545           DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
546     ! Volume x average at momentum cell centers
547           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
548     ! Volume x average at momentum cell centers Solid-Solid Drag
549           DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
550     ! Coefficients for pressure correction equation. These coefficients are
551     ! initialized to zero in the subroutine time_march before the time loop.
552           DOUBLE PRECISION, INTENT(INOUT) :: D_E(DIMENSION_3, 0:DIMENSION_M)
553     
554     
555     ! Local variables:
556     !---------------------------------------------------------------------//
557     ! Usual Indices
558           INTEGER :: LM, M, L, LPL, LP
559           INTEGER :: I, IJK, IJKE
560     ! Average solid volume fraction at momentum cell centers
561           DOUBLE PRECISION :: EPSA(DIMENSION_M)
562     ! local alias for face area
563           DOUBLE PRECISION :: AREA_FACE
564     ! sum of Solid M - All other Solid drag
565           DOUBLE PRECISION :: SUM_VXF_SS(DIMENSION_M)
566     ! sum of Solid L - All other Solid VolxDrag but M
567           DOUBLE PRECISION :: SUM_VXF_SS_wt_M
568     ! component of the total numerator needed for evaluation of pressure
569     ! correction coefficient for SOLIDS phase M. specifically, for when the
570     ! sum L is over solids phase only
571           DOUBLE PRECISION :: NUM_MSol_LSol(DIMENSION_M)
572     ! component of the total denominator needed for evaluation of pressure
573     ! correction coefficient for SOLIDS phase M.  specifically, for when the
574     ! sum L is over solids phases only
575           DOUBLE PRECISION :: DEN_MSol_LSol(DIMENSION_M)
576     ! Temp variable for double precision values.
577           DOUBLE PRECISION :: TMPdp
578     
579     !......................................................................!
580     
581     !$omp parallel default(none) &
582     !$omp          private(ijk,area_face,i,ijke,epsa,lm,sum_vxf_ss,tmpdp,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
583     !$omp          shared(ijkstart3,ijkend3,d_e,mmax,cartesian_grid,ayz,i_of,vxf_gs,vxf_ss,momentum_x_eq,am0,model_b,p_scale)
584     !$omp do
585           DO IJK = IJKSTART3, IJKEND3
586              IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK) .OR. MODEL_B) THEN
587                 DO M= 1, MMAX
588                    D_E(IJK,M) = ZERO
589                 ENDDO
590                 CYCLE
591              ENDIF
592     
593              AREA_FACE = merge(ONE, AYZ(IJK), CARTESIAN_GRID)
594              I = I_OF(IJK)
595              IJKE = EAST_OF(IJK)
596     
597              DO M= 1, MMAX
598                 EPSA(M) = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
599                 SUM_VXF_SS(M) = ZERO
600                 DO L = 1, MMAX
601                   IF (L .NE. M) THEN
602                      LM = FUNLM(L,M)
603     ! Solids M - All other Solids VolxDrag summation
604                      SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
605                   ENDIF
606                ENDDO
607              ENDDO
608     
609              DO M= 1, MMAX
610     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
611                 NUM_MSol_LSol(M) = ZERO
612                 DEN_MSol_LSol(M)  = ZERO
613                 DO L = 1, MMAX
614                    IF (L .NE. M) THEN
615                       LM = FUNLM(L,M)
616                       SUM_VXF_SS_wt_M = ZERO
617                       DO Lp = 1, MMAX
618                          IF ( (Lp .NE. L) .AND. (Lp .NE. M) ) THEN
619                             LpL = FUNLM(Lp,L)
620     ! Solids L - All other Solids VolxDrag but M summation
621                             SUM_VXF_SS_wt_M =                              &
622                                SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
623                          ENDIF
624                       ENDDO
625                       IF (MOMENTUM_X_EQ(L)) THEN
626                          NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +             &
627                             (VXF_SS(IJK,LM)*EPSA(L) /                      &
628                             ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+    &
629                             SMALL_NUMBER))
630     
631                          DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +             &
632                             VXF_SS(IJK,LM)*(((-AM0(IJK,L))+VXF_GS(IJK,L) + &
633                             SUM_VXF_SS_wt_M)/((-AM0(IJK,L))+VXF_GS(IJK,L)+ &
634                             SUM_VXF_SS(L)+ SMALL_NUMBER ))
635                       ELSE
636                          DEN_MSol_LSol(M) = DEN_MSol_LSol(M)+VXF_SS(IJK,LM)
637                       ENDIF
638                    ENDIF   ! end if (L.ne.M)
639                 ENDDO  ! end do (l=1,mmax)
640     
641     ! Linking velocity correction coefficient to pressure - SOLIDs Phase (Model_A only)
642                 IF (MOMENTUM_X_EQ(M)) THEN
643                    TMPdp = -AM0(IJK,M)+VXF_GS(IJK,M)+DEN_MSol_LSol(M)
644                    IF(abs(TMPdp) > SMALL_NUMBER) THEN
645                       D_E(IJK,M) = P_SCALE*AREA_FACE*                      &
646                          (EPSA(M) + NUM_MSol_LSol(M))/TMPdp
647                    ELSE
648                       D_E(IJK,M) = ZERO
649                    ENDIF
650                 ELSE
651                    D_E(IJK,M) = ZERO
652                 ENDIF
653              ENDDO  ! end do (m=1,mmax)
654           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
655     !$omp end parallel
656     
657           RETURN
658           END SUBROUTINE CALC_D_E_SOLIDS_ONLY
659