File: /nfs/home/0/users/jenkins/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,num_msol_lgas,den_msol_lgas,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
209     !$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)
210     !$omp do
211           DO IJK = ijkstart3, ijkend3
212     
213              IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN
214                 DO M= 0, MMAX
215                    D_E(IJK,M) = ZERO
216                 ENDDO
217                 CYCLE
218              ENDIF
219     
220              AREA_FACE = merge(ONE, AYZ(IJK), CARTESIAN_GRID)
221              I = I_OF(IJK)
222              IJKE = EAST_OF(IJK)
223              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
224     
225              SUM_VXF_GS = ZERO
226              DO M= 1, MMAX
227                 EPSA(M) = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
228     ! Gas - All Solids VolxDrag summation
229                 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
230                 SUM_VXF_SS(M) = ZERO
231                 DO L = 1, MMAX
232                    IF (L .NE. M) THEN
233                       LM = FUNLM(L,M)
234     ! Solids M - All other Solids VolxDrag summation
235                       SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
236                    ENDIF
237                 ENDDO
238              ENDDO
239     
240     ! calculating DEN_MGas and NUM_MGas
241              DEN_MGas  = ZERO
242              NUM_MGas = ZERO
243              DO M= 1, MMAX
244                 IF (MOMENTUM_X_EQ(M)) THEN
245                    NUM_MGas = NUM_MGas + (EPSA(M)*VXF_GS(IJK,M)/           &
246                      ((-AM0(IJK,M)) + VXF_GS(IJK,M) + SUM_VXF_SS(M) +      &
247                      SMALL_NUMBER))
248                    DEN_MGas = DEN_MGas + (VXF_GS(IJK,M)*                   &
249                      ((-AM0(IJK,M))+SUM_VXF_SS(M))/                        &
250                      ((-AM0(IJK,M))+VXF_GS(IJK,M)+SUM_VXF_SS(M)+           &
251                      SMALL_NUMBER))
252                 ELSE
253                    DEN_MGas = DEN_MGas + VXF_GS(IJK,M)
254                 ENDIF
255              ENDDO
256     
257     ! Model B
258     ! -------------------------------------------------------------------
259              IF (MODEL_B) THEN
260     ! Linking velocity correction coefficient to pressure - GAS Phase
261                 TMPdp = -AM0(IJK,0) + DEN_MGas
262                 IF(abs(TMPdp) > SMALL_NUMBER) THEN
263                    D_E(IJK,0) = P_SCALE*AREA_FACE/TMPdp
264                 ELSE
265                    D_E(IJK,0) = ZERO
266                 ENDIF
267     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
268                 DO M = 1, MMAX
269                    IF(MOMENTUM_X_EQ(M)) THEN
270                       TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M)
271                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
272                          D_E(IJK,M) = D_E(IJK,0)*VXF_GS(IJK,M)/TMPdp
273                       ELSE
274                          D_E(IJK,M) = ZERO
275                       ENDIF
276                    ELSE
277                       D_E(IJK,M) = ZERO
278                    ENDIF
279                 ENDDO
280     
281     ! Model A
282     ! -------------------------------------------------------------------
283              ELSE
284     
285     ! Linking velocity correction coefficient to pressure - GAS Phase
286                 TMPdp = -AM0(IJK,0) + DEN_MGas
287                 IF(abs(TMPdp) >SMALL_NUMBER) THEN
288                    D_E(IJK,0) = P_SCALE*AREA_FACE*(EPGA+NUM_MGas)/TMPdp
289                 ELSE
290                    D_E(IJK,0) = ZERO
291                 ENDIF
292     
293                 DO M= 1, MMAX
294     ! calculating NUM_MSol_LGas and DEN_MSol_LGas
295                    NUM_MSol_LGas = VXF_GS(IJK,M)*EPGA/                     &
296                       ((-AM0(IJK,0))+SUM_VXF_GS+SMALL_NUMBER)
297                    DEN_MSol_LGas = VXF_GS(IJK,M)*(                         &
298                       ((-AM0(IJK,0)) + (SUM_VXF_GS - VXF_GS(IJK,M)))/      &
299                       ((-AM0(IJK,0)) + SUM_VXF_GS + SMALL_NUMBER) )
300     
301     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
302                    NUM_MSol_LSol(M) = ZERO
303                    DEN_MSol_LSol(M)  = ZERO
304                    DO L = 1, MMAX
305                       IF (L /= M) THEN
306                          LM = FUNLM(L,M)
307                          SUM_VXF_SS_wt_M = ZERO
308                          DO Lp = 1, MMAX
309                             IF((Lp /= L) .AND. (Lp /= M) ) THEN
310                                LpL = FUNLM(Lp,L)
311     ! Solids L - All other Solids VolxDrag but M summation
312                                SUM_VXF_SS_wt_M =                           &
313                                   SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
314                             ENDIF
315                          ENDDO
316     
317                          IF(MOMENTUM_X_EQ(L)) THEN
318                             NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +          &
319                                ( VXF_SS(IJK,LM)*EPSA(L)/                   &
320                                ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
321                                SMALL_NUMBER) )
322     
323                             DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +          &
324                                VXF_SS(IJK,LM)*(((-AM0(IJK,L))+             &
325                                VXF_GS(IJK,L)+SUM_VXF_SS_wt_M)/             &
326                                ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
327                                SMALL_NUMBER) )
328                          ELSE
329                             DEN_MSol_LSol(M) =                             &
330                                DEN_MSol_LSol(M) + VXF_SS(IJK,LM)
331                          ENDIF
332                       ENDIF  ! end if (l.ne.m)
333                    ENDDO  ! end do (l=1,mmax)
334     
335     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
336                    IF(MOMENTUM_X_EQ(M)) THEN
337                       TMPdp = -AM0(IJK,M) + DEN_MSol_LGas + DEN_MSol_LSol(M)
338                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
339                          D_E(IJK,M) = P_SCALE*AREA_FACE * (EPSA(M) +       &
340                             NUM_MSol_LSol(M) + NUM_MSol_LGas)/TMPdp
341                       ELSE
342                          D_E(IJK,M) = ZERO
343                       ENDIF
344                    ELSE
345                       D_E(IJK,M) = ZERO
346                    ENDIF
347                 ENDDO   ! end do (m=1,mmax)
348     
349              ENDIF    !end if/else branch Model_B/Model_A
350     
351           ENDDO  ! end do loop (ijk=ijkstart3,ijkend3)
352     !$omp end parallel
353     
354           RETURN
355           END SUBROUTINE CALC_D_E_GAS_AND_SOLIDS
356     
357     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
358     !                                                                      !
359     !  Subroutine: CALC_D_E_GAS_ONLY                                       !
360     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
361     !                                                                      !
362     !  Purpose: calculate coefficients linking velocity correction to      !
363     !  pressure correction -- East                                         !
364     !                                                                      !
365     !  Notes: The gas phase x-momentum equation is solved; NO solids phase !
366     !  x-momentum equations are solved.                                    !
367     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
368     
369           SUBROUTINE CALC_D_E_GAS_ONLY(AM0, VXF_GS, D_E)
370     
371     ! Global Variables:
372     !---------------------------------------------------------------------//
373     ! Volume fraction of gas phase.
374           use fldvar, only: EP_G
375     ! Number of solids phases.
376           use physprop, only: MMAX
377     ! Flag: Use MODEL_B
378           use run, only: MODEL_B
379     ! Pressure scale factor
380           use scales, only: P_SCALE
381     ! Flag: Cartesian grid simulation
382           use cutcell, only: CARTESIAN_GRID
383     ! Area of cell's YZ face.
384           use geometry, only: AYZ
385     ! Volume of V-momentum cell.
386           use geometry, only: VOL_U
387     ! Function to average across Y face.
388           use fun_avg, only: AVG_X
389     ! Fluid grid loop bounds.
390           use compar, only: IJKStart3, IJKEnd3
391     ! Flags: Impermeable surface and mass flow at east face, IJK of cell to east
392           use functions, only: IP_AT_E, MFLOW_AT_E, EAST_OF
393     ! Indices: J index of cell
394           use indices, only: I_OF
395     ! Flag and variables for QMOM implementation.
396           USE qmom_kinetic_equation, only: QMOMK, QMOMK_NN
397           USE qmom_kinetic_equation, only: QMOMK_F_GS, QMOMK_F_GS
398     
399     ! Global Parameters:
400     !---------------------------------------------------------------------//
401     ! Size of IJK arrays and size of solids phase arrays.
402           use param, only: DIMENSION_3, DIMENSION_M
403     ! Size of MxM upper triangular matrix.
404           use param1, only: DIMENSION_LM
405     ! Double precision parameters.
406           use param1, only: ZERO, SMALL_NUMBER, ONE
407     
408           IMPLICIT NONE
409     
410     ! Dummy arguments
411     !---------------------------------------------------------------------//
412     ! Temporary variable to store A matrix for local manipulation
413           DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
414     ! Volume x average at momentum cell centers
415           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
416     ! Coefficients for pressure correction equation. These coefficients are
417     ! initialized to zero in the subroutine time_march before the time loop.
418           DOUBLE PRECISION, INTENT(INOUT) :: D_E(DIMENSION_3, 0:DIMENSION_M)
419     
420     ! Local variables:
421     !---------------------------------------------------------------------//
422     ! Usual Indices
423           INTEGER :: M
424           INTEGER :: INN
425           INTEGER :: I, IJK, IJKE
426     ! Average solid volume fraction at momentum cell centers
427           !DOUBLE PRECISION :: EPSA(DIMENSION_M)
428     ! Average volume fraction at momentum cell centers
429           DOUBLE PRECISION :: EPGA
430     ! local alias for face area
431           DOUBLE PRECISION :: AREA_FACE
432     ! sum of All Solid-Gas VolxDrag
433           DOUBLE PRECISION :: SUM_VXF_GS
434     ! Temp variable for double precision values.
435           DOUBLE PRECISION :: TMPdp
436     
437     !......................................................................!
438     
439     !$omp parallel default(none) &
440     !$omp          private(ijk,area_face,i,ijke,epga,sum_vxf_gs,tmpdp) &
441     !$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)
442     !$omp do
443           DO IJK = IJKSTART3, IJKEND3
444     
445              IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN   !impermeable
446                 D_E(IJK,0) = ZERO
447                 CYCLE
448              ENDIF
449     
450              AREA_FACE = merge(ONE, AYZ(IJK), CARTESIAN_GRID)
451     
452              I = I_OF(IJK)
453              IJKE = EAST_OF(IJK)
454              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
455     
456              SUM_VXF_GS = ZERO
457              IF (.NOT. QMOMK) THEN
458                 DO M= 1, MMAX
459     ! Gas - All Solids VolxDrag summation
460                    SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
461                 ENDDO
462              ELSE
463                 DO INN = 1, QMOMK_NN
464                    DO M = 1, MMAX
465                       SUM_VXF_GS = SUM_VXF_GS + VOL_U(IJK)* &
466                          AVG_X(QMOMK_F_GS(INN,IJK,M),QMOMK_F_GS(INN,IJKE,M),I)
467                    ENDDO
468                 ENDDO
469              ENDIF
470     
471              TMPdp = -AM0(IJK,0) + SUM_VXF_GS
472              IF(abs(TMPdp) > SMALL_NUMBER) THEN
473                 IF (MODEL_B) THEN
474                    D_E(IJK,0) = P_SCALE*AREA_FACE/TMPdp
475                 ELSE
476                    D_E(IJK,0) = P_SCALE*AREA_FACE*EPGA/TMPdp
477                 ENDIF   !end if/else branch Model_B/Model_A
478              ELSE
479                 D_E(IJK,0) = ZERO
480              ENDIF
481     
482           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
483     !$omp end parallel
484     
485           RETURN
486           END SUBROUTINE CALC_D_E_GAS_ONLY
487     
488     
489     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
490     !                                                                      !
491     !  Subroutine: CALC_D_E_SOLIDS_ONLY                                    !
492     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
493     !                                                                      !
494     !  Purpose: Calculate coefficients linking velocity correction to      !
495     !  pressure correction -- East                                         !
496     !                                                                      !
497     !  Notes: The gas phase momentum equations are NOT solved but at least !
498     !  one solids phase momentum equation is solved.                       !
499     !                                                                      !
500     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
501           SUBROUTINE CALC_D_E_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_E)
502     
503     ! Global Variables:
504     !---------------------------------------------------------------------//
505     ! Volume fractions of gas and solids phases.
506          use fldvar, only: EP_G, EP_s
507     ! Number of solids phases.
508          use physprop, only: MMAX
509     ! Flag: Solve the X-Momentum Equations
510          use run, only: MOMENTUM_X_EQ
511     ! Flag: Use MODEL_B
512          use run, only: MODEL_B
513     ! Pressure scale factor
514          use scales, only: P_SCALE
515     ! Flag: Cartesian grid simulation
516           use cutcell, only: CARTESIAN_GRID
517     ! Area of cell's YZ face.
518           use geometry, only: AYZ
519     ! Function to average across X face.
520           use fun_avg, only: AVG_X
521     ! Fluid grid loop bounds.
522           use compar, only: IJKStart3, IJKEnd3
523     ! Function to convert to phases to single array, IJK of cell to east
524           use functions, only: FUNLM, EAST_OF
525     ! Flags: Impermeable surface and mass flow at east face
526           use functions, only: IP_AT_E, MFLOW_AT_E
527     ! Indices: I index of cell
528           use indices, only: I_OF
529     
530     ! Global Parameters:
531     !---------------------------------------------------------------------//
532     ! Size of IJK arrays and size of solids phase arrays.
533           use param, only: DIMENSION_3, DIMENSION_M
534     ! Size of MxM upper triangular matrix.
535           use param1, only: DIMENSION_LM
536     ! Double precision parameters.
537           use param1, only: ZERO, SMALL_NUMBER, ONE
538     
539           IMPLICIT NONE
540     
541     ! Dummy arguments
542     !---------------------------------------------------------------------//
543     ! Error index
544     !      INTEGER, INTENT(INOUT) :: IER
545     ! Temporary variable to store A matrix for local manipulation
546           DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
547     ! Volume x average at momentum cell centers
548           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
549     ! Volume x average at momentum cell centers Solid-Solid Drag
550           DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
551     ! Coefficients for pressure correction equation. These coefficients are
552     ! initialized to zero in the subroutine time_march before the time loop.
553           DOUBLE PRECISION, INTENT(INOUT) :: D_E(DIMENSION_3, 0:DIMENSION_M)
554     
555     
556     ! Local variables:
557     !---------------------------------------------------------------------//
558     ! Usual Indices
559           INTEGER :: LM, M, L, LPL, LP
560           INTEGER :: I, IJK, IJKE
561     ! Average solid volume fraction at momentum cell centers
562           DOUBLE PRECISION :: EPSA(DIMENSION_M)
563     ! local alias for face area
564           DOUBLE PRECISION :: AREA_FACE
565     ! sum of Solid M - All other Solid drag
566           DOUBLE PRECISION :: SUM_VXF_SS(DIMENSION_M)
567     ! sum of Solid L - All other Solid VolxDrag but M
568           DOUBLE PRECISION :: SUM_VXF_SS_wt_M
569     ! component of the total numerator needed for evaluation of pressure
570     ! correction coefficient for SOLIDS phase M. specifically, for when the
571     ! sum L is over solids phase only
572           DOUBLE PRECISION :: NUM_MSol_LSol(DIMENSION_M)
573     ! component of the total denominator needed for evaluation of pressure
574     ! correction coefficient for SOLIDS phase M.  specifically, for when the
575     ! sum L is over solids phases only
576           DOUBLE PRECISION :: DEN_MSol_LSol(DIMENSION_M)
577     ! Temp variable for double precision values.
578           DOUBLE PRECISION :: TMPdp
579     
580     !......................................................................!
581     
582     !$omp parallel default(none) &
583     !$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) &
584     !$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)
585     !$omp do
586           DO IJK = IJKSTART3, IJKEND3
587              IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK) .OR. MODEL_B) THEN
588                 DO M= 1, MMAX
589                    D_E(IJK,M) = ZERO
590                 ENDDO
591                 CYCLE
592              ENDIF
593     
594              AREA_FACE = merge(ONE, AYZ(IJK), CARTESIAN_GRID)
595              I = I_OF(IJK)
596              IJKE = EAST_OF(IJK)
597     
598              DO M= 1, MMAX
599                 EPSA(M) = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
600                 SUM_VXF_SS(M) = ZERO
601                 DO L = 1, MMAX
602                   IF (L .NE. M) THEN
603                      LM = FUNLM(L,M)
604     ! Solids M - All other Solids VolxDrag summation
605                      SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
606                   ENDIF
607                ENDDO
608              ENDDO
609     
610              DO M= 1, MMAX
611     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
612                 NUM_MSol_LSol(M) = ZERO
613                 DEN_MSol_LSol(M)  = ZERO
614                 DO L = 1, MMAX
615                    IF (L .NE. M) THEN
616                       LM = FUNLM(L,M)
617                       SUM_VXF_SS_wt_M = ZERO
618                       DO Lp = 1, MMAX
619                          IF ( (Lp .NE. L) .AND. (Lp .NE. M) ) THEN
620                             LpL = FUNLM(Lp,L)
621     ! Solids L - All other Solids VolxDrag but M summation
622                             SUM_VXF_SS_wt_M =                              &
623                                SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
624                          ENDIF
625                       ENDDO
626                       IF (MOMENTUM_X_EQ(L)) THEN
627                          NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +             &
628                             (VXF_SS(IJK,LM)*EPSA(L) /                      &
629                             ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+    &
630                             SMALL_NUMBER))
631     
632                          DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +             &
633                             VXF_SS(IJK,LM)*(((-AM0(IJK,L))+VXF_GS(IJK,L) + &
634                             SUM_VXF_SS_wt_M)/((-AM0(IJK,L))+VXF_GS(IJK,L)+ &
635                             SUM_VXF_SS(L)+ SMALL_NUMBER ))
636                       ELSE
637                          DEN_MSol_LSol(M) = DEN_MSol_LSol(M)+VXF_SS(IJK,LM)
638                       ENDIF
639                    ENDIF   ! end if (L.ne.M)
640                 ENDDO  ! end do (l=1,mmax)
641     
642     ! Linking velocity correction coefficient to pressure - SOLIDs Phase (Model_A only)
643                 IF (MOMENTUM_X_EQ(M)) THEN
644                    TMPdp = -AM0(IJK,M)+VXF_GS(IJK,M)+DEN_MSol_LSol(M)
645                    IF(abs(TMPdp) > SMALL_NUMBER) THEN
646                       D_E(IJK,M) = P_SCALE*AREA_FACE*                      &
647                          (EPSA(M) + NUM_MSol_LSol(M))/TMPdp
648                    ELSE
649                       D_E(IJK,M) = ZERO
650                    ENDIF
651                 ELSE
652                    D_E(IJK,M) = ZERO
653                 ENDIF
654              ENDDO  ! end do (m=1,mmax)
655           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
656     !$omp end parallel
657     
658           RETURN
659           END SUBROUTINE CALC_D_E_SOLIDS_ONLY
660