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

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