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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CALC_d_n                                                !
4     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
5     !                                                                      !
6     !  Purpose: calculate coefficients linking velocity correction to      !
7     !           pressure correction -- North                               !
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_N(A_M, VXF_GS, VXF_SS, D_N, IER)
16     
17     ! Global Variables:
18     !---------------------------------------------------------------------//
19     ! Number of solids phases.
20           use physprop, only: MMAX
21     ! Flag: Solve the Y-Momentum Equations
22           use run, only: MOMENTUM_Y_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(OUT) :: D_N(DIMENSION_3, 0:DIMENSION_M)
52     ! Integer error flag.
53           INTEGER, INTENT(OUT) :: 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_Y_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_Y_MOMENTUM = any(MOMENTUM_Y_EQ(1:MMAX))
82     
83           IF(MOMENTUM_Y_EQ(0))THEN
84              IF(ANY_SOLIDS_Y_MOMENTUM) THEN
85                 CALL CALC_D_N_GAS_AND_SOLIDS(AM0, VxF_GS, VxF_SS, D_N)
86              ELSE
87                 CALL CALC_D_N_GAS_ONLY(AM0, VxF_GS, D_N)
88              ENDIF
89     
90           ELSEIF(ANY_SOLIDS_Y_MOMENTUM) THEN
91              CALL CALC_D_N_SOLIDS_ONLY(AM0, VxF_GS, VxF_SS, D_N)
92           ENDIF
93     
94           deallocate(AM0)
95     
96           RETURN
97           END SUBROUTINE CALC_D_N
98     
99     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
100     !                                                                      !
101     !  Subroutine: CALC_D_N_GAS_AND_SOLIDS                                 !
102     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
103     !                                                                      !
104     !  Purpose: calculate coefficients linking velocity correction to      !
105     !           pressure correction -- North                               !
106     !                                                                      !
107     !  CASE: The gas phase Y-momentum equation is solved and at least one  !
108     !    solids phase (m>0) y-momentum equation is solved.                 !
109     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
110           SUBROUTINE CALC_D_N_GAS_AND_SOLIDS(AM0, VXF_GS, VXF_SS, D_N)
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 Y-Momentum Equations
119           use run, only: MOMENTUM_Y_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 XZ face.
127           use geometry, only: AXZ
128     ! Function to average across Y face.
129           use fun_avg, only: AVG_Y
130     ! Fluid grid loop bounds.
131           use compar, only: IJKStart3, IJKEnd3
132     ! Function to convert to phases to single array, IJK of cell to north
133           use functions, only: FUNLM, NORTH_OF
134     ! Flags: Impermeable surface and mass flow at north face
135           use functions, only: IP_AT_N, MFLOW_AT_N
136     ! Indices: J index of cell
137           use indices, only: J_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_N(DIMENSION_3, 0:DIMENSION_M)
161     
162     ! Local variables:
163     !---------------------------------------------------------------------//
164     ! Usual Indices
165           INTEGER :: LM, M, L, LPL, LP
166           INTEGER :: J, IJK, IJKN
167     ! Average volume fraction at momentum cell centers
168           DOUBLE PRECISION :: EPGA
169     ! Average solid volume fraction at momentum cell centers
170           DOUBLE PRECISION :: EPSA(DIMENSION_M)
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,j,ijkn,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_n,mmax,cartesian_grid,axz,j_of,ep_g,vxf_gs,vxf_ss,momentum_y_eq,am0,model_b,p_scale)
206     !$omp do
207           DO IJK = IJKSTART3, IJKEND3
208              IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN
209                 D_N(IJK,0:MMAX) = ZERO
210                 CYCLE
211              ENDIF
212     
213              AREA_FACE = merge(ONE, AXZ(IJK), CARTESIAN_GRID)
214              J = J_OF(IJK)
215              IJKN = NORTH_OF(IJK)
216              EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
217     
218              SUM_VXF_GS = ZERO
219              DO M= 1, MMAX
220                 EPSA(M) = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
221     ! Gas - All Solids VolxDrag summation
222                 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
223                 SUM_VXF_SS(M) = ZERO
224                 DO L = 1, MMAX
225                    IF (L .NE. M) THEN
226                       LM = FUNLM(L,M)
227     ! Solids M - All other Solids VolxDrag summation
228                       SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
229                    ENDIF
230                 ENDDO
231              ENDDO
232     
233     ! Calculating DEN_MGas and NUM_MGas
234              NUM_MGas = ZERO
235              DEN_MGas  = ZERO
236              DO M= 1, MMAX
237                 IF(MOMENTUM_Y_EQ(M)) THEN
238                    NUM_MGas = NUM_MGas + (EPSA(M)*VXF_GS(IJK,M) /          &
239                       ((-AM0(IJK,M)) + VXF_GS(IJK,M) + SUM_VXF_SS(M) +     &
240                       SMALL_NUMBER))
241                    DEN_MGas = DEN_MGas + (VXF_GS(IJK,M)*                   &
242                       ((-AM0(IJK,M)) + SUM_VXF_SS(M)) /                    &
243                       ((-AM0(IJK,M)) + VXF_GS(IJK,M) + SUM_VXF_SS(M) +     &
244                       SMALL_NUMBER))
245                 ELSE
246                    DEN_MGas = DEN_MGas + VXF_GS(IJK,M)
247                 ENDIF
248              ENDDO
249     
250     ! Model B
251     ! -------------------------------------------------------------------
252              IF (MODEL_B) THEN
253     ! Linking velocity correction coefficient to pressure - GAS Phase
254                 TMPdp = -AM0(IJK,0) + DEN_MGas
255                 IF(abs(TMPdp) > SMALL_NUMBER ) THEN
256                    D_N(IJK,0) = P_SCALE*AREA_FACE/TMPdp
257                 ELSE
258                    D_N(IJK,0) = ZERO
259                 ENDIF
260     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
261                 DO M = 1, MMAX
262                    IF(MOMENTUM_Y_EQ(M)) THEN
263                       TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M)
264                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
265                          D_N(IJK,M) = D_N(IJK,0)*VXF_GS(IJK,M) / TMPdp
266                       ELSE
267                          D_N(IJK,M) = ZERO
268                       ENDIF
269                    ELSE
270                       D_N(IJK,M) = ZERO
271                    ENDIF
272                 ENDDO
273     
274     ! Model A
275     ! -------------------------------------------------------------------
276              ELSE
277     
278     ! Linking velocity correction coefficient to pressure - GAS Phase
279                 TMPdp = -AM0(IJK,0) + DEN_MGas
280                 IF(abs(TMPdp) > SMALL_NUMBER) THEN
281                    D_N(IJK,0) = P_SCALE*AREA_FACE*(EPGA+NUM_MGas) / TMPdp
282                 ELSE
283                    D_N(IJK,0) = ZERO
284                 ENDIF
285     
286                 DO M= 1, MMAX
287     ! calculating NUM_MSol_LGas and DEN_MSol_LGas
288                    NUM_MSol_LGas = VXF_GS(IJK,M)*EPGA/                     &
289                       ((-AM0(IJK,0))+SUM_VXF_GS+SMALL_NUMBER)
290                    DEN_MSol_LGas = VXF_GS(IJK,M)*(                         &
291                       ((-AM0(IJK,0))+(SUM_VXF_GS - VXF_GS(IJK,M))) /       &
292                       ((-AM0(IJK,0))+SUM_VXF_GS + SMALL_NUMBER) )
293     
294     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
295                    NUM_MSol_LSol(M) = ZERO
296                    DEN_MSol_LSol(M)  = ZERO
297                    DO L = 1, MMAX
298                       IF (L .NE. M) THEN
299                          LM = FUNLM(L,M)
300                          SUM_VXF_SS_wt_M = ZERO
301                          DO Lp = 1, MMAX
302                             IF((Lp /= L) .AND. (Lp /= M)) THEN
303                                LpL = FUNLM(Lp,L)
304     ! Solids L - All other Solids VolxDrag but M summation
305                                SUM_VXF_SS_wt_M = SUM_VXF_SS_wt_M +         &
306                                   VXF_SS(IJK,LpL)
307                             ENDIF
308                          ENDDO
309     
310                          IF(MOMENTUM_Y_EQ(L)) THEN
311                             NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +          &
312                                (VXF_SS(IJK,LM)*EPSA(L) /                   &
313                                ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
314                                SMALL_NUMBER ))
315     
316                             DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +          &
317                                VXF_SS(IJK,LM)*(((-AM0(IJK,L)) +            &
318                                VXF_GS(IJK,L) + SUM_VXF_SS_wt_M ) /         &
319                                ((-AM0(IJK,L)) + VXF_GS(IJK,L)+             &
320                                SUM_VXF_SS(L) + SMALL_NUMBER) )
321                          ELSE
322                             DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +          &
323                                VXF_SS(IJK,LM)
324                          ENDIF
325                       ENDIF  ! end if (l.ne.m)
326                    ENDDO  ! end do (l=1,mmax)
327     
328     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
329                    IF(MOMENTUM_Y_EQ(M)) THEN
330                       TMPdp = -AM0(IJK,M) + DEN_MSol_LGas + DEN_MSol_LSol(M)
331     
332                       IF(abs(TMPdp) > SMALL_NUMBER) THEN
333                          D_N(IJK,M) = P_SCALE*AREA_FACE*(EPSA(M) +         &
334                             NUM_MSol_LSol(M) + NUM_MSol_LGas) / TMPdp
335                       ELSE
336                          D_N(IJK,M) = ZERO
337                       ENDIF
338                    ELSE
339                       D_N(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_N_GAS_AND_SOLIDS
349     
350     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
351     !                                                                      !
352     !  Subroutine: CALC_D_N_GAS_ONLY                                       !
353     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
354     !                                                                      !
355     !  Purpose: Calculate coefficients linking velocity correction to      !
356     !  pressure correction -- North                                        !
357     !                                                                      !
358     !  CASE: Only the gas phase (M=0) Y-momentnum equation is solved; the  !
359     !  solids Y-momentum equations are NOT solved. This routine is where a !
360     !  coupled DEM simulation should be directed for proper evaluation of  !
361     !  pressure correction terms.                                          !
362     !                                                                      !
363     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
364           SUBROUTINE CALC_D_N_GAS_ONLY(AM0, VXF_GS, D_N)
365     
366     ! Global Variables:
367     !---------------------------------------------------------------------//
368     ! Volume fractions of gas and solids phases.
369          use fldvar, only: EP_G
370     ! Number of solids phases.
371          use physprop, only: MMAX
372     ! Flag: Use MODEL_B
373          use run, only: MODEL_B
374     ! Pressure scale factor
375          use scales, only: P_SCALE
376     ! Flag: Cartesian grid simulation
377           use cutcell, only: CARTESIAN_GRID
378     ! Area of cell's XZ face.
379           use geometry, only: AXZ
380     ! Volume of V-momentum cell.
381           use geometry, only: VOL_V
382     ! Function to average across Y face.
383           use fun_avg, only: AVG_Y
384     ! Fluid grid loop bounds.
385           use compar, only: IJKStart3, IJKEnd3
386     ! Flags: Impermeable surface and mass flow at north face, IJK of cell to north
387           use functions, only: IP_AT_N, MFLOW_AT_N, NORTH_OF
388     ! Indices: J index of cell
389           use indices, only: J_OF
390     ! Flag and variables for QMOM implementation.
391           use qmom_kinetic_equation, only: QMOMK, QMOMK_NN
392           use qmom_kinetic_equation, only: QMOMK_F_GS, QMOMK_F_GS
393     
394     ! Global Parameters:
395     !---------------------------------------------------------------------//
396     ! Size of IJK arrays and size of solids phase arrays.
397           use param, only: DIMENSION_3, DIMENSION_M
398     ! Size of MxM upper triangular matrix.
399           use param1, only: DIMENSION_LM
400     ! Double precision parameters.
401           use param1, only: ZERO, SMALL_NUMBER, ONE
402     
403           IMPLICIT NONE
404     
405     ! Dummy arguments
406     !---------------------------------------------------------------------//
407     ! Temporary variable to store A matrix for local manipulation
408           DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
409     ! Volume x average at momentum cell centers
410           DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
411     ! Coefficients for pressure correction equation. These coefficients are
412     ! initialized to zero in the subroutine time_march before the time loop
413           DOUBLE PRECISION, INTENT(INOUT) :: d_n(DIMENSION_3, 0:DIMENSION_M)
414     
415     ! Local variables:
416     !---------------------------------------------------------------------//
417     ! Usual Indices
418           INTEGER :: M
419           INTEGER :: INN
420           INTEGER :: J, IJK, IJKN
421     ! Average volume fraction at momentum cell centers
422           DOUBLE PRECISION :: EPGA
423     ! local alias for face area
424           DOUBLE PRECISION :: AREA_FACE
425     ! sum of All Solid-Gas VolxDrag
426           DOUBLE PRECISION :: SUM_VXF_GS
427     ! Temp variable for double precision values.
428           DOUBLE PRECISION :: TMPdp
429     
430     !......................................................................!
431     
432     !$omp parallel default(none) &
433     !$omp          private(ijk,area_face,j,ijkn,epga,sum_vxf_gs,tmpdp) &
434     !$omp          shared(ijkstart3,ijkend3,d_n,mmax,cartesian_grid,axz,j_of,ep_g,vxf_gs,am0,model_b,p_scale,qmomk,vol_v,qmomk_f_gs)
435     !$omp do
436           DO IJK = IJKSTART3, IJKEND3
437     
438              IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN
439                 D_N(IJK,0) = ZERO
440                 CYCLE
441              ENDIF
442     
443              AREA_FACE = merge(ONE, AXZ(IJK), CARTESIAN_GRID)
444              J = J_OF(IJK)
445              IJKN = NORTH_OF(IJK)
446              EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
447     
448              SUM_VXF_GS = ZERO
449              IF (.NOT. QMOMK) THEN
450                 DO M= 1, MMAX
451     ! Gas - All Solids VolxDrag summation
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_V(IJK)*AVG_Y(          &
458                          QMOMK_F_GS(INN,IJK,M),QMOMK_F_GS(INN,IJKN,M),J)
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_N(IJK,0) = P_SCALE*AREA_FACE/TMPdp
467                 ELSE
468                    D_N(IJK,0) = P_SCALE*AREA_FACE*EPGA/TMPdp
469                 ENDIF   !end if/else branch Model_B/Model_A
470              ELSE
471                 D_N(IJK,0) = ZERO
472              ENDIF
473     
474           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
475     !$omp end parallel
476     
477           RETURN
478           END SUBROUTINE CALC_D_N_GAS_ONLY
479     
480     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
481     !                                                                      !
482     !  Subroutine: CALC_D_N_SOLIDS_ONLY                                    !
483     !  Author: M. Syamlal                                 Date: 21-JUN-96  !
484     !                                                                      !
485     !  Purpose: calculate coefficients linking velocity correction to      !
486     !           pressure correction -- North                               !
487     !                                                                      !
488     !  CASE: At least one solids phase Y-momentum equation is solved but   !
489     !  the gas phase (m=0) Y-momentum equation is NOT solved.              !
490     !                                                                      !
491     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
492           SUBROUTINE CALC_D_N_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_N)
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 Y-Momentum Equations
501          use run, only: MOMENTUM_Y_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 XZ face.
509           use geometry, only: AXZ
510     ! Function to average across Y face.
511           use fun_avg, only: AVG_Y
512     ! Fluid grid loop bounds.
513           use compar, only: IJKStart3, IJKEnd3
514     ! Function to convert to phases to single array, IJK of cell to north
515           use functions, only: FUNLM, NORTH_OF
516     ! Flags: Impermeable surface and mass flow at north face
517           use functions, only: IP_AT_N, MFLOW_AT_N
518     ! Indices: J index of cell
519           use indices, only: J_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_N(DIMENSION_3, 0:DIMENSION_M)
543     
544     ! Local variables:
545     !---------------------------------------------------------------------//
546     ! Usual Indices
547           INTEGER :: LM, M, L, LPL, LP
548           INTEGER :: J, IJK, IJKN
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,j,ijkn,epsa,lm,sum_vxf_ss,tmpdp,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
572     !$omp          shared(ijkstart3,ijkend3,d_n,mmax,cartesian_grid,axz,j_of,ep_g,vxf_gs,vxf_ss,momentum_y_eq,am0,model_b,p_scale)
573     !$omp do
574           DO IJK = IJKSTART3, IJKEND3
575              IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK) .OR. MODEL_B) THEN
576                 D_N(IJK,1:MMAX) = ZERO
577                 CYCLE
578              ENDIF
579     
580              AREA_FACE = merge(ONE, AXZ(IJK), CARTESIAN_GRID)
581              J = J_OF(IJK)
582              IJKN = NORTH_OF(IJK)
583     
584              DO M= 1, MMAX
585                 EPSA(M) = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
586                 SUM_VXF_SS(M) = ZERO
587                 DO L = 1, MMAX
588                    IF (L .NE. M) THEN
589                       LM = FUNLM(L,M)
590     ! Solids M - All other Solids VolxDrag summation
591                       SUM_VXF_SS(M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
592                    ENDIF
593                 ENDDO
594              ENDDO
595     
596              DO M= 1, MMAX
597     ! calculating NUM_MSol_LSol and DEN_MSol_LSol
598                 NUM_MSol_LSol(M) = ZERO
599                 DEN_MSol_LSol(M)  = ZERO
600                 DO L = 1, MMAX
601                    IF (L .NE. M) THEN
602                       LM = FUNLM(L,M)
603                       SUM_VXF_SS_wt_M = ZERO
604                       DO Lp = 1, MMAX
605                          IF((Lp /= L) .AND. (Lp /= M)) THEN
606                             LpL = FUNLM(Lp,L)
607     ! Solids L - All other Solids VolxDrag but M summation
608                             SUM_VXF_SS_wt_M=SUM_VXF_SS_wt_M+VXF_SS(IJK,LpL)
609                          ENDIF
610                       ENDDO
611                       IF(MOMENTUM_Y_EQ(L)) THEN
612                          NUM_MSol_LSol(M) = NUM_MSol_LSol(M) +             &
613                             (VXF_SS(IJK,LM)*EPSA(L)/((-AM0(IJK,L)) +       &
614                             VXF_GS(IJK,L)+SUM_VXF_SS(L)+SMALL_NUMBER ))
615                          DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +             &
616                             VXF_SS(IJK,LM)*(((-AM0(IJK,L)) + VXF_GS(IJK,L)+&
617                             SUM_VXF_SS_wt_M )/((-AM0(IJK,L))+VXF_GS(IJK,L)+&
618                             SUM_VXF_SS(L) + SMALL_NUMBER))
619                       ELSE
620                          DEN_MSol_LSol(M) = DEN_MSol_LSol(M)+VXF_SS(IJK,LM)
621                       ENDIF
622                    ENDIF   ! end if (l.ne.m)
623                 ENDDO  ! end do (l=1,mmax)
624     
625     ! Linking velocity correction coefficient to pressure - SOLIDs Phase
626                 IF (MOMENTUM_Y_EQ(M)) THEN
627                    TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M) + DEN_MSol_LSol(M)
628                    IF(abs(TMPdp) > SMALL_NUMBER) THEN
629                       D_N(IJK,M) = P_SCALE*AREA_FACE* (EPSA(M) +           &
630                          NUM_MSol_LSol(M))/TMPdp
631                    ELSE
632                       D_N(IJK,M) = ZERO
633                    ENDIF
634                 ELSE
635                    D_N(IJK,M) = ZERO
636                 ENDIF
637              ENDDO  ! end do (m=1,mmax)
638           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
639     !$omp end parallel
640     
641           RETURN
642           END SUBROUTINE CALC_D_N_SOLIDS_ONLY
643