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