File: RELATIVE:/../../../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     
205     !$omp          private(ijk,area_face,k,ijkt,epga,sum_vxf_gs,epsa,lm,sum_vxf_ss,den_mgas,num_mgas,         &
206     !$omp                  tmpdp,num_msol_lgas,den_msol_lgas,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
207     !$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)
208     !$omp do
209           DO IJK = IJKSTART3, IJKEND3
210     
211              IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN   !impermeable
212                 DO M= 0, MMAX
213                    D_T(IJK,M) = ZERO
214                 ENDDO
215                 CYCLE
216              ENDIF
217     
218              AREA_FACE = merge(ONE, AXY(IJK), CARTESIAN_GRID)
219              K = K_OF(IJK)
220              IJKT = TOP_OF(IJK)
221              EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
222     
223              SUM_VXF_GS = ZERO
224              DO M= 1, MMAX
225                 EPSA(M) = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
226     ! Gas - All Solids VolxDrag summation
227                 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
228                 SUM_VXF_SS(M) = ZERO
229                 DO L = 1, MMAX
230                    IF (L .NE. M) THEN
231                       LM = FUNLM(L,M)
232     ! Solid M - All other Solids VolxDrag summation
233                       SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
234                    ENDIF
235                 ENDDO
236              ENDDO
237     
238     ! calculating DEN_MGas and NUM_MGas
239              DEN_MGas  = ZERO
240              NUM_MGas = ZERO
241              DO M= 1, MMAX
242                 IF(MOMENTUM_Z_EQ(M)) THEN
243                    NUM_MGas = NUM_MGas + (EPSA(M)*VXF_GS(IJK,M)/           &
244                       ((-AM0(IJK,M))+VXF_GS(IJK,M)+SUM_VXF_SS(M) +         &
245                       SMALL_NUMBER))
246                    DEN_MGas = DEN_MGas + (VXF_GS(IJK,M)*                   &
247                      ((-AM0(IJK,M))+SUM_VXF_SS(M))/                        &
248                      ((-AM0(IJK,M))+VXF_GS(IJK,M)+SUM_VXF_SS(M)+           &
249                      SMALL_NUMBER))
250                 ELSE
251                    DEN_MGas = DEN_MGas + VXF_GS(IJK,M)
252                 ENDIF
253              ENDDO
254     
255     ! Model B
256     ! -------------------------------------------------------------------
257              IF (MODEL_B) THEN
258     ! Linking velocity correction coefficient to pressure - GAS Phase
259                 TMPdp = -AM0(IJK,0) + DEN_MGas
260                 IF(abs(TMPdp) > SMALL_NUMBER) THEN
261                    D_T(IJK,0) = P_SCALE*AREA_FACE/TMPdp
262                 ELSE
263                    D_T(IJK,0) = ZERO
264                 ENDIF
265     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
266                 DO M = 1, MMAX
267                    IF(MOMENTUM_Z_EQ(M)) THEN
268                       TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M)
269                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
270                          D_T(IJK,M) = D_T(IJK,0)*VXF_GS(IJK,M)/TMPdp
271                       ELSE
272                          D_T(IJK,M) = ZERO
273                       ENDIF
274                    ELSE
275                       D_T(IJK,M) = ZERO
276                    ENDIF
277                 ENDDO
278     
279     ! Model A
280     ! -------------------------------------------------------------------
281              ELSE
282     
283     ! Linking velocity correction coefficient to pressure - GAS Phase
284                 TMPdp = -AM0(IJK,0)+DEN_MGas
285                 IF(abs(TMPdp) > SMALL_NUMBER) THEN
286                    D_T(IJK,0) = P_SCALE*AREA_FACE*(EPGA+NUM_MGas)/TMPdp
287                 ELSE
288                    D_T(IJK,0) = ZERO
289                 ENDIF
290     
291                 DO M= 1, MMAX
292     ! calculating NUM_MSol_LGas and DEN_MSol_LGas
293                    NUM_MSol_LGas = VXF_GS(IJK,M)*EPGA/                     &
294                       ((-AM0(IJK,0))+SUM_VXF_GS+SMALL_NUMBER)
295                    DEN_MSol_LGas = VXF_GS(IJK,M)*(                         &
296                       ((-AM0(IJK,0)) + (SUM_VXF_GS - VXF_GS(IJK,M)))/      &
297                       ((-AM0(IJK,0)) + SUM_VXF_GS+SMALL_NUMBER))
298     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
299                    NUM_MSol_LSol(M) = ZERO
300                    DEN_MSol_LSol(M)  = ZERO
301                    DO L = 1, MMAX
302                       IF (L .NE. M) THEN
303                          LM = FUNLM(L,M)
304                          SUM_VXF_SS_wt_M = ZERO
305                          DO Lp = 1, MMAX
306                             IF( (Lp /= L) .AND. (Lp /= M) ) THEN
307                                LpL = FUNLM(Lp,L)
308     ! Solids L - All other Solids VolxDrag but M summation
309                                SUM_VXF_SS_wt_M =                           &
310                                   SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
311                             ENDIF
312                          ENDDO
313     
314                          IF(MOMENTUM_Z_EQ(L)) THEN
315                             NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +          &
316                                (VXF_SS(IJK,LM)*EPSA(L)/((-AM0(IJK,L)) +    &
317                                VXF_GS(IJK,L)+SUM_VXF_SS(L)+ SMALL_NUMBER))
318     
319                             DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +          &
320                                VXF_SS(IJK,LM)*(((-AM0(IJK,L)) +            &
321                                VXF_GS(IJK,L)+SUM_VXF_SS_wt_M )/            &
322                                ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
323                                SMALL_NUMBER))
324                          ELSE
325                             DEN_MSol_LSol(M) =                             &
326                                DEN_MSol_LSol(M) + VXF_SS(IJK,LM)
327                          ENDIF
328                       ENDIF  ! end if (l.ne.m)
329                    ENDDO  ! end do (l=1,mmax)
330     
331     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
332                    IF(MOMENTUM_Z_EQ(M)) THEN
333                       TMPdp = -AM0(IJK,M) + DEN_MSol_LGas + DEN_MSol_LSol(M)
334                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
335                          D_T(IJK,M) = P_SCALE*AREA_FACE*(EPSA(M) +         &
336                             NUM_MSol_LSol(M) + NUM_MSol_LGas)/TMPdp
337                       ELSE
338                          D_T(IJK,M) = ZERO
339                       ENDIF
340                    ELSE
341                       D_T(IJK,M) = ZERO
342                    ENDIF
343                 ENDDO   ! end do (m=1,mmax)
344     
345              ENDIF    !end if/else branch Model_B/Model_A
346           ENDDO  ! end do loop (ijk=ijkstart3,ijkend3)
347     !$omp end parallel
348     
349           RETURN
350           END SUBROUTINE CALC_D_T_GAS_AND_SOLIDS
351     
352     
353     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
354     !                                                                      !
355     !  Subroutine: CALC_D_T_GAS_ONLY                                       !
356     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
357     !                                                                      !
358     !  Purpose: Calculate coefficients linking velocity correction to      !
359     !           pressure correction -- Top                                 !
360     !                                                                      !
361     !  CASE: Only the gas phase (M=0) Z-momentnum equation is solved; the  !
362     !  solids Z-momentum equations are NOT solved. This routine is where a !
363     !  coupled DEM simulation should be directed for proper evaluation of  !
364     !  pressure correction terms.                                          !
365     !                                                                      !
366     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
367           SUBROUTINE CALC_D_T_GAS_ONLY(AM0, VXF_GS, D_T)
368     
369     ! Global Variables:
370     !---------------------------------------------------------------------//
371     ! Volume fractions of gas and solids phases.
372          use fldvar, only: EP_G
373     ! Number of solids phases.
374          use physprop, only: MMAX
375     ! Flag: Use MODEL_B
376          use run, only: MODEL_B
377     ! Pressure scale factor
378          use scales, only: P_SCALE
379     ! Flag: Cartesian grid simulation
380           use cutcell, only: CARTESIAN_GRID
381     ! Area of cell's XY face.
382           use geometry, only: AXY
383     ! Volume of W-momentum cell.
384           use geometry, only: VOL_W
385     ! Function to average across Z face.
386           use fun_avg, only: AVG_Z
387     ! Fluid grid loop bounds.
388           use compar, only: IJKStart3, IJKEnd3
389     ! Flags: Impermeable surface and mass flow at top face, IJK of cell to top
390           use functions, only: IP_AT_T, MFLOW_AT_T, TOP_OF
391     ! Indices: K index of cell
392           use indices, only: K_OF
393     ! Flag and variables for QMOM implementation.
394           use qmom_kinetic_equation, only: QMOMK, QMOMK_NN
395           use qmom_kinetic_equation, only: QMOMK_F_GS, QMOMK_F_GS
396     
397     ! Global Parameters:
398     !---------------------------------------------------------------------//
399     ! Size of IJK arrays and size of solids phase arrays.
400           use param, only: DIMENSION_3, DIMENSION_M
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_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,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