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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DRAG_gs                                                 C
4     !  Purpose: Calculate the gas-solids drag coefficient                  C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 20-MAY-96  C
7     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !  Literature/Document References:                                     C
10     !                                                                      C
11     !  Comments:                                                           C
12     !    If changes are made to this file, especially in regard to the     C
13     !    structure of any drag subroutine call, then such changes must     C
14     !    also be made to the analgous routine des_drag_gp in the file      C
15     !    des/drag_fgs.f for consistency!                                   C
16     !                                                                      C
17     !                                                                      C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20     
21           SUBROUTINE DRAG_GS(M, IER)
22     
23     !-----------------------------------------------
24     ! Modules
25     !-----------------------------------------------
26           USE param
27           USE param1
28           USE parallel
29           USE fldvar
30           USE geometry
31           USE indices
32           USE physprop
33           USE run
34           USE constant
35           USE compar
36           USE drag
37           USE sendrecv
38           USE discretelement
39           USE ur_facs
40           USE funits
41           USE mms
42           USE fun_avg
43           USE functions
44           IMPLICIT NONE
45     !-----------------------------------------------
46     ! Dummy arguments
47     !-----------------------------------------------
48     ! Solids phase index
49           INTEGER, INTENT(IN) :: M
50     ! Error index
51           INTEGER, INTENT(INOUT) :: IER
52     !-----------------------------------------------
53     ! Local parameters
54     !-----------------------------------------------
55     
56     !-----------------------------------------------
57     ! Local variables
58     !-----------------------------------------------
59     ! indices
60           INTEGER :: I, IJK, IMJK, IJMK, IJKM
61     ! cell center value of U_sm, V_sm, W_sm
62           DOUBLE PRECISION :: USCM, VSCM, WSCM
63     ! cell center value of x, y and z-particle velocity
64           DOUBLE PRECISION :: USCM_HYS, VSCM_HYS, WSCM_HYS
65     ! cell center value of U_g, V_g, W_g
66           DOUBLE PRECISION :: UGC, VGC, WGC
67     ! magnitude of gas-solids relative velocity
68           DOUBLE PRECISION :: VREL
69     ! gas laminar viscosity redefined here to set viscosity at pressure
70     ! boundaries
71           DOUBLE PRECISION :: Mu
72     ! drag coefficient
73           DOUBLE PRECISION :: DgA
74     ! current value of F_gs (i.e., without underrelaxation)
75           DOUBLE PRECISION F_gstmp
76     ! indices of solids phases (continuous, discrete)
77           INTEGER :: CM, DM, L
78     ! temporary shift of total number of solids phases to account for both
79     ! discrete and continuous solids phases used for the hybrid mdoel
80           INTEGER :: MAXM
81     ! tmp local variable for the particle diameter of solids
82     ! phase M (continuous or discrete)
83           DOUBLE PRECISION :: DP_loc(2*DIM_M)
84     ! tmp local variable for the solids volume fraction of solids
85     ! phase M (continuous or discrete)
86           DOUBLE PRECISION :: EPs_loc(2*DIM_M)
87     ! tmp local variable for the particle density of solids
88     ! phase M (continuous or discrete)
89           DOUBLE PRECISION :: ROs_loc(2*DIM_M)
90     ! correction factors for implementing polydisperse drag model
91     ! proposed by van der Hoef et al. (2005)
92           DOUBLE PRECISION :: F_cor, tmp_sum, tmp_fac
93     ! average particle diameter in polydisperse systems
94           DOUBLE PRECISION :: DPA
95     ! diameter ratio in polydisperse systems
96           DOUBLE PRECISION :: Y_i
97     ! total solids volume fraction
98           DOUBLE PRECISION :: phis
99     ! temporary local variables to use for dummy arguments in subroutines
100     ! void fraction, gas density, gas bulk density, solids volume fraction
101     ! particle diameter, particle density
102           DOUBLE PRECISION :: EPG, ROg, ROPg, EP_SM, DPM, ROs
103     !-----------------------------------------------
104     
105     !!$omp  parallel do default(shared)                                   &
106     !!$omp  private( I,  IJK, IMJK, IJMK, IJKM, DM, MAXM, CM, L,          &
107     !!$omp           UGC, VGC, WGC, USCM, VSCM, WSCM, VREL, USCM_HYS,     &
108     !!$omp           VSCM_HYS, WSCM_HYS, tmp_sum, tmp_fac, Y_i, F_cor,    &
109     !!$omp           EP_SM, EPs_loc, ROs_loc, DP_loc, DPA, DPM, ROs,      &
110     !!$omp           phis, EPg, ROg, ROPg, Mu, DgA, F_gstmp)
111     
112     
113           DO IJK = ijkstart3, ijkend3
114     
115              IF (FLUIDorP_FLOW_AT(IJK)) THEN
116     
117                 I = I_OF(IJK)
118                 IMJK = IM_OF(IJK)
119                 IJMK = JM_OF(IJK)
120                 IJKM = KM_OF(IJK)
121     
122                 MAXM = SMAX
123                 DO CM = 1,SMAX
124                    DP_loc(CM) = D_p(IJK,CM)
125                    EPs_loc(CM) = EP_S(IJK,CM)
126                    ROs_loc(CM) = RO_S(IJK,CM)
127                 ENDDO
128     
129     ! Append DEM phase information to the end for hybrid runs and update
130     ! the number of solids 'phases' to loop over.
131                 IF(DES_CONTINUUM_HYBRID) THEN
132                    MAXM = SMAX + DES_MMAX
133                    DO DM = 1,DES_MMAX
134                       L = SMAX + DM
135                       DP_loc(L) = DES_D_p0(DM)
136                       EPs_loc(L) = DES_ROP_S(IJK,DM)/DES_RO_S(DM)
137                       ROs_loc(L) = DES_RO_S(DM)
138                    ENDDO
139                 ENDIF
140     
141     ! Calculate velocity components at i, j, k
142                 UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
143                 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
144                 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
145     
146                 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
147                 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
148                 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
149     
150     ! magnitude of gas-solids relative velocity
151                 VREL = SQRT((UGC-USCM)**2 + (VGC-VSCM)**2 + (WGC-WSCM)**2)
152     
153     ! Laminar viscosity at a pressure boundary is given the value of the
154     ! fluid cell next to it. This applies just to the calculation of the
155     ! drag, in other routines the value of viscosity at a pressure boundary
156     ! always has a zero value.
157                 IF (P_FLOW_AT(IJK)) THEN
158                    IF( FLUID_AT(EAST_OF(IJK)) ) THEN
159                       Mu = MU_G(EAST_OF(IJK))
160                    ELSEIF ( FLUID_AT(WEST_OF(IJK)) ) THEN
161                       Mu = MU_G(WEST_OF(IJK))
162                    ELSEIF ( FLUID_AT(NORTH_OF(IJK)) ) THEN
163                       Mu = MU_G(NORTH_OF(IJK))
164                    ELSEIF ( FLUID_AT(SOUTH_OF(IJK)) ) THEN
165                       Mu = MU_G(SOUTH_OF(IJK))
166                    ELSEIF ( FLUID_AT(TOP_OF(IJK)) ) THEN
167                       Mu = MU_G(TOP_OF(IJK))
168                    ELSEIF ( FLUID_AT(BOTTOM_OF(IJK)) ) THEN
169                       Mu = MU_G(BOTTOM_OF(IJK))
170                    ENDIF
171                 ELSE
172                    Mu = MU_G(IJK)
173                 ENDIF
174     
175     ! calculate the total solids volume fraction
176                 phis = ZERO
177                 DO L = 1, MAXM
178     ! this is slightly /= one-ep_g due to round-off
179                    phis = phis + EPs_loc(L)
180                 ENDDO
181     
182     ! calculate the average paricle diameter and particle ratio
183                 DPA = ZERO
184                 tmp_sum = ZERO
185                 tmp_fac = ZERO
186                 DO L = 1, MAXM
187                    IF (phis .GT. ZERO) THEN
188                       tmp_fac = EPs_loc(L)/phis
189                       tmp_sum = tmp_sum + tmp_fac/DP_loc(L)
190                     ELSE
191                       tmp_sum = tmp_sum + ONE/DP_loc(L) ! not important, but will avoid NaN's in empty cells
192                     ENDIF
193                 ENDDO
194                 DPA = ONE / tmp_sum
195                 Y_i = DP_loc(M)/DPA
196     
197     ! assign aliases for easy reference
198                 EPg = EP_G(IJK)
199                 ROg = RO_G(IJK)
200                 ROPg = ROP_G(IJK)
201                 EP_SM = EPs_loc(M)
202                 DPM = DP_loc(M)
203                 ROs = ROs_loc(M)
204     
205     
206     ! determine the drag coefficient
207                 IF (EP_SM <= ZERO) THEN
208                    DgA = ZERO
209                 ELSEIF (EPg == ZERO) THEN
210     ! this case will already be caught in most drag subroutines whenever
211     ! RE==0 (for correlations in which RE includes EPg). however, this will
212     ! prevent potential divisions by zero in some models by setting it now.
213                    DgA = ZERO
214     ! Force a ZERO drag coefficient for MMS cases.
215                 ELSEIF (USE_MMS) THEN
216                    DgA = ZERO
217                 ELSE
218                    SELECT CASE(DRAG_TYPE_ENUM)
219     
220                    CASE (SYAM_OBRIEN)
221                       CALL DRAG_SYAM_OBRIEN(DgA,EPG,Mu,ROg,VREL,DPM)
222     
223                    CASE (GIDASPOW)
224                       CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
225     
226                    CASE (GIDASPOW_PCF)
227                       CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
228     
229                    CASE (GIDASPOW_BLEND)
230                       CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
231     
232                    CASE (GIDASPOW_BLEND_PCF)
233                       CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
234     
235                    CASE (WEN_YU)
236                       CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPM)
237     
238                    CASE (WEN_YU_PCF)
239                       CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPA)
240     
241                    CASE (KOCH_HILL)
242                       CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPM,phis)
243     
244                    CASE (KOCH_HILL_PCF)
245                       CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
246     
247                    CASE (BVK)
248                       CALL DRAG_BVK(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
249     
250                    CASE (USER_DRAG)
251                       CALL DRAG_USR(IJK, M, DgA, EPg, Mu, ROg, VREL, DPM,  &
252                          ROs, UGC, VGC, WGC)
253     
254                    CASE (HYS)
255     ! calculate velocity components of each solids phase
256                       USCM_HYS = ZERO
257                       VSCM_HYS = ZERO
258                       WSCM_HYS = ZERO
259                       IF(phis > ZERO) THEN
260                          DO L = 1, MAXM
261                             USCM_HYS = USCM_HYS + EPs_loc(L)*(UGC - &
262                                AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I))
263                             VSCM_HYS = VSCM_HYS + EPs_loc(L)*(VGC - &
264                                AVG_Y_N(V_S(IJMK,L),V_S(IJK,L)))
265                             WSCM_HYS = WSCM_HYS + EPs_loc(L)*(WGC - &
266                                AVG_Z_T(W_S(IJKM,L),W_S(IJK,L)))
267                          ENDDO
268                          USCM_HYS = USCM_HYS/phis
269                          VSCM_HYS = VSCM_HYS/phis
270                          WSCM_HYS = WSCM_HYS/phis
271                       ENDIF
272     ! magnitude of gas-solids relative velocity
273                       VREL = SQRT(USCM_HYS**2 +VSCM_HYS**2 +WSCM_HYS**2)
274     
275                       CALL DRAG_HYS(DgA,EPg,Mu,ROPg,VREL,&
276                            DP_loc(:),DPA,Y_i,EPs_loc(:),phis,M,MAXM,IJK)
277     
278                    CASE DEFAULT
279                       CALL START_LOG
280                       IF(.NOT.DMP_LOG) call open_pe_log(ier)
281                       IF(DMP_LOG) WRITE (*, '(A,A)') &
282                          'Unknown DRAG_TYPE: ', DRAG_TYPE
283                       WRITE (UNIT_LOG, '(A,A)') 'Unknown DRAG_TYPE: ', DRAG_TYPE
284                       CALL END_LOG
285                       CALL mfix_exit(myPE)
286                    END SELECT   ! end selection of drag_type
287                 ENDIF   ! end if/elseif/else (ep_sm <= zero, ep_g==0)
288     
289     ! modify the drag coefficient to account for subgrid domain effects
290                 IF (SUBGRID_TYPE_ENUM .ne. UNDEFINED_SUBGRID_TYPE) THEN
291                    IF (SUBGRID_TYPE_ENUM .EQ. IGCI) THEN
292                       CALL SUBGRID_DRAG_IGCI(DgA,EPg,Mu,ROg,DPM,ROs,IJK)
293                    ELSEIF (SUBGRID_TYPE_ENUM .EQ. MILIOLI) THEN
294                       CALL SUBGRID_DRAG_MILIOLI(DgA,EPg,Mu,ROg,VREL,&
295                            DPM,ROs,IJK)
296                    ENDIF
297                 ENDIF
298     
299     
300     ! Modify drag coefficient to account for possible corrections and
301     ! for differences between Model B and Model A
302                 IF(DRAG_TYPE_ENUM == HYS) THEN
303     ! this drag model is handled differently than the others
304                    IF(Model_B)THEN
305                       F_gstmp = DgA/EPg
306                    ELSE
307                       F_gstmp = DgA
308                    ENDIF
309                 ELSE
310                    IF(DRAG_TYPE_ENUM == GIDASPOW_PCF .OR. &
311                       DRAG_TYPE_ENUM == GIDASPOW_BLEND_PCF .OR. &
312                       DRAG_TYPE_ENUM == WEN_YU_PCF .OR. &
313                       DRAG_TYPE_ENUM == KOCH_HILL_PCF .OR. &
314                       DRAG_TYPE_ENUM == BVK) THEN
315     ! see erratum by Beetstra et al. (2007) : the correction factor differs
316     ! for model A versus model B.
317     ! application of the correction factor for model A is found from
318     ! the correction factor for model B and neglects the Y_i**3 term
319                       IF(Model_B) THEN
320                          IF (M == 1) THEN
321                             F_cor = (EPg*Y_i + phis*Y_i**2)
322                          ELSE
323                             F_cor = (EPg*Y_i + phis*Y_i**2 + &
324                                0.064d0*EPg*Y_i**3)
325                          ENDIF
326                       ELSE
327                          F_cor = Y_i
328                       ENDIF
329                       DgA = ONE/(Y_i*Y_i) * DgA * F_cor
330                    ENDIF
331     
332     ! Calculate the drag coefficient (Model B coeff = Model A coeff/EP_g); all other models, eg., Wen_Yu
333                    IF(Model_B)THEN
334                       F_gstmp = DgA * EP_SM/EPg
335                    ELSE
336                       F_gstmp = DgA * EP_SM                     !3D buoyancy model
337                    ENDIF
338                 ENDIF   !end if/else trim(drag_type=='hys')
339     
340     ! Determine drag force coefficient accounting for any under relaxation
341                 F_GS(IJK,M) = (ONE - UR_F_gs)*F_GS(IJK, M) + &
342                    UR_F_gs*F_gstmp
343     
344                 IF(KT_TYPE_ENUM == GHD_2007) THEN
345                    IF(M==1) THEN
346                       F_gs(IJK,MMAX) = F_gs(IJK,M)
347                    ELSE
348                       F_gs(IJK,MMAX) = F_gs(IJK,MMAX) + F_gs(IJK,M)
349                    ENDIF
350                 ENDIF
351     
352              ELSE   ! .not.(fluidorp_flow_at(ijk)) branch
353     
354                 F_GS(IJK,M) = ZERO
355                 IF(KT_TYPE_ENUM == GHD_2007) F_gs(IJK, MMAX) = ZERO
356     
357              ENDIF   ! end if (fluidorp_flow_at(ijk))
358     
359           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
360     
361           RETURN
362           END SUBROUTINE DRAG_GS
363     
364     
365     
366     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
367     !                                                                      C
368     !  Function(s): C_DsxRe                                                C
369     !  Purpose:                                                            C
370     !     Calculate single sphere drag correlation multiplied by           C
371     !     the Reynolds number or                                           C
372     !     Calculate the single sphere drag correlation                     C
373     !                                                                      C
374     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
375     
376     ! Dalla Valle (1948)
377     !----------------------------------------------------------------->>>
378           DOUBLE PRECISION FUNCTION C_DSXRE_DV(RE)
379           USE param
380           USE param1
381           IMPLICIT NONE
382           DOUBLE PRECISION, INTENT(IN) :: RE ! Reynolds number
383     
384           C_DSXRE_DV = (0.63D0*SQRT(RE) + 4.8D0)**2
385           RETURN
386           END FUNCTION C_DSXRE_DV
387     !-----------------------------------------------------------------<<<
388     
389     ! Turton and Levenspiel (1986)
390     !----------------------------------------------------------------->>>
391           DOUBLE PRECISION FUNCTION C_DSXRE_TL(RE)
392           USE param
393           USE param1
394           IMPLICIT NONE
395           DOUBLE PRECISION, INTENT(IN) :: RE ! Reynolds number
396     
397           C_DSXRE_TL = 24.D0*(1.D0 + 0.173D0*RE**0.657D0) + &
398              0.413D0*RE**2.09D0/(RE**1.09D0 + 16300.D0)
399           RETURN
400           END FUNCTION C_DSXRE_TL
401     !-----------------------------------------------------------------<<<
402     
403     
404     ! Schiller and Naumann (1933)
405     !----------------------------------------------------------------->>>
406           DOUBLE PRECISION FUNCTION C_DS_SN(RE)
407           USE param
408           USE param1
409           IMPLICIT NONE
410           DOUBLE PRECISION, INTENT(IN) :: RE ! Reynolds number
411     
412           C_DS_SN = 24.D0*(1.D0 + 0.15D0*RE**0.687D0)/(RE+SMALL_NUMBER)
413           RETURN
414           END FUNCTION C_DS_SN
415     !-----------------------------------------------------------------<<<
416     
417     
418     
419     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
420     !                                                                      C
421     !  Subroutine: DRAG_SYAM_OBRIEN                                        C
422     !  Purpose: Calculate the gas-solids drag coefficient                  C
423     !                                                                      C
424     !  Literature/Document References:                                     C
425     !     Syamlal M, O'Brien TJ (1988). International Journal of           C
426     !        Multiphase Flow 14: 473-481.                                  C
427     !                                                                      C
428     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
429     
430           SUBROUTINE DRAG_SYAM_OBRIEN(lDgA,EPg,Mug,ROg,VREL,&
431                      DPM)
432     
433     !-----------------------------------------------
434     ! Modules
435     !-----------------------------------------------
436           USE param
437           USE param1
438           USE constant, only : drag_c1, drag_d1, c
439           IMPLICIT NONE
440     !-----------------------------------------------
441     ! Dummy arguments
442     !-----------------------------------------------
443     ! drag coefficient
444           DOUBLE PRECISION, INTENT(OUT) :: ldGA
445     ! gas volume fraction
446           DOUBLE PRECISION, INTENT(IN) :: EPg
447     ! gas laminar viscosity
448           DOUBLE PRECISION, INTENT(IN) :: Mug
449     ! gas density
450           DOUBLE PRECISION, INTENT(IN) :: ROg
451     ! Magnitude of gas-solids relative velocity
452           DOUBLE PRECISION, INTENT(IN) :: VREL
453     ! particle diameter of solids phase M
454           DOUBLE PRECISION, INTENT(IN) :: DPM
455     !-----------------------------------------------
456     ! Local parameters
457     !-----------------------------------------------
458     !     Parameters in the Cluster-effect model
459     !     PARAMETER (a1 = 250.)   !for G_s = 98 kg/m^2.s
460     !     PARAMETER (a1 = 1500.)  !for G_s = 147 kg/m^2.s
461     !     a1 depends upon solids flux.  It has been represented by C(1)
462     !     defined in the data file.
463           DOUBLE PRECISION, PARAMETER :: A2 = 0.005D0
464           DOUBLE PRECISION, PARAMETER :: A3 = 90.0D0
465           DOUBLE PRECISION, PARAMETER :: RE_C = 5.D0
466           DOUBLE PRECISION, PARAMETER :: EP_C = 0.92D0
467     !-----------------------------------------------
468     ! Local variables
469     !-----------------------------------------------
470     ! Variables which are function of EP_g
471           DOUBLE PRECISION :: A, B
472     ! Ratio of settling velocity of a multiparticle system to
473     ! that of a single particle
474           DOUBLE PRECISION :: V_rm
475     ! Reynolds number
476           DOUBLE PRECISION :: RE
477     !-----------------------------------------------
478     ! External functions
479     !-----------------------------------------------
480     ! Single sphere drag coefficient x Re
481           DOUBLE PRECISION, EXTERNAL :: C_DSXRE_DV
482     !-----------------------------------------------
483     
484           IF(Mug > ZERO) THEN
485              RE = DPM*VREL*ROg/Mug
486           ELSE
487              RE = LARGE_NUMBER
488           ENDIF
489     
490     ! Calculate V_rm
491           A = EPg**4.14D0
492           IF (EPg <= 0.85D0) THEN
493              B = drag_c1*EPg**1.28D0
494           ELSE
495             B = EPg**drag_d1
496           ENDIF
497     
498           V_RM=HALF*(A-0.06D0*RE+&
499                SQRT((3.6D-3)*RE*RE+0.12D0*RE*(2.D0*B-A)+A*A) )
500     
501     !------------------Begin cluster correction --------------------------
502     ! uncomment the following four lines ...
503     !       V_RM = V_RM * (ONE + C(1)*&
504     !                      EXP(-A2*(RE-RE_C)**2 - A3*(EPg-EP_C)**2)* &
505     !                      RE*(1. - EPg))
506     !------------------End cluster correction ----------------------------
507     
508           lDgA = 0.75D0*Mug*EPg*C_DSXRE_DV(RE/V_RM)/(V_RM*DPM*DPM)
509     
510           IF (RE == ZERO) lDgA = ZERO
511     
512           RETURN
513           END SUBROUTINE DRAG_SYAM_OBRIEN
514     
515     
516     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
517     !                                                                      C
518     !  Subroutine: DRAG_GIDASPOW                                           C
519     !  Purpose: Calculate the gas-solids drag coefficient                  C
520     !                                                                      C
521     !  Literature/Document References:                                     C
522     !     Ding J, Gidaspow D (1990). AIChE Journal 36: 523-538.            C
523     !                                                                      C
524     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
525     
526           SUBROUTINE DRAG_GIDASPOW(lDgA,EPg,Mug,ROg,ROPg,VREL, DPM)
527     
528     !-----------------------------------------------
529     ! Modules
530     !-----------------------------------------------
531           USE param
532           USE param1
533           IMPLICIT NONE
534     !-----------------------------------------------
535     ! Dummy arguments
536     !-----------------------------------------------
537     ! drag coefficient
538           DOUBLE PRECISION, INTENT(OUT) :: lDgA
539     ! gas volume fraction
540           DOUBLE PRECISION, INTENT(IN) :: EPg
541     ! gas laminar viscosity
542           DOUBLE PRECISION, INTENT(IN) :: Mug
543     ! gas density
544           DOUBLE PRECISION, INTENT(IN) :: ROg
545     ! gas density*EP_g
546           DOUBLE PRECISION, INTENT(IN) :: ROPg
547     ! Magnitude of gas-solids relative velocity
548           DOUBLE PRECISION, INTENT(IN) :: VREL
549     ! particle diameter of solids phase M or
550     ! average particle diameter if PCF
551           DOUBLE PRECISION, INTENT(IN) :: DPM
552     
553     !
554           DOUBLE PRECISION, EXTERNAL :: C_DS_SN
555     
556     !-----------------------------------------------
557     ! Local variables
558     !-----------------------------------------------
559     ! Reynolds number
560           DOUBLE PRECISION :: RE
561     ! Single sphere drag coefficient
562           DOUBLE PRECISION :: C_d
563     !-----------------------------------------------
564     
565     
566     ! Note the presence of gas volume fraction in ROPG
567           RE = merge(DPM*VREL*ROPg/Mug, LARGE_NUMBER, MUg > ZERO)
568     
569     ! Dense phase
570           IF(EPg <= 0.8D0) THEN
571              lDgA = 150D0*(ONE-EPg)*Mug / (EPg*DPM**2) + &
572                     1.75D0*ROg*VREL/DPM
573           ELSE
574     ! Dilute phase - EP_g >= 0.8
575              IF(RE <= 1000D0)THEN
576     ! this could be replaced with the function C_DS_SN
577                 C_d = C_DS_SN(RE)
578              ELSE
579                 C_d = 0.44D0
580              ENDIF
581              lDgA = 0.75D0*C_d*VREL*ROPg*EPg**(-2.65D0) / DPM
582           ENDIF
583     
584           IF (RE == ZERO) lDgA = ZERO
585     
586           RETURN
587           END SUBROUTINE DRAG_GIDASPOW
588     
589     
590     
591     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
592     !                                                                      C
593     !  Subroutine: DRAG_GIDASPOW_BLEND                                     C
594     !  Purpose: Calculate the gas-solids drag coefficient                  C
595     !                                                                      C
596     !  Author: Charles E.A. Finney                        Date: 23-Mar-06  C
597     !  Reviewer: Sreekanth Pannala                        Date: 24-Mar-06  C
598     !                                                                      C
599     !  Literature/Document References:                                     C
600     !     original source unknown:                                         C
601     !     Lathouwers D, Bellan J (2000). Proceedings of the 2000 U.S. DOE  C
602     !        Hydrogen Program Review NREL/CP-570-28890. Available from     C
603     !     http://www.eere.energy.gov/hydrogenandfuelcells/pdfs/28890k.pdf. C
604     !                                                                      C
605     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
606     
607           SUBROUTINE DRAG_GIDASPOW_BLEND(lDgA,EPg,Mug,ROg,ROPg,VREL,&
608                      DPM)
609     
610     !-----------------------------------------------
611     ! Modules
612     !-----------------------------------------------
613           USE param
614           USE param1
615           USE constant, only : PI
616           IMPLICIT NONE
617     !-----------------------------------------------
618     ! Dummy arguments
619     !-----------------------------------------------
620     ! drag coefficient
621           DOUBLE PRECISION, INTENT(OUT) :: lDgA
622     ! gas volume fraction
623           DOUBLE PRECISION, INTENT(IN) :: EPg
624     ! gas laminar viscosity
625           DOUBLE PRECISION, INTENT(IN) :: Mug
626     ! gas density
627           DOUBLE PRECISION, INTENT(IN) :: ROg
628     ! gas density*EP_g
629           DOUBLE PRECISION, INTENT(IN) :: ROPg
630     ! Magnitude of gas-solids relative velocity
631           DOUBLE PRECISION, INTENT(IN) :: VREL
632     ! particle diameter of solids phase M or
633     ! average particle diameter if PCF
634           DOUBLE PRECISION, INTENT(IN) :: DPM
635     !-----------------------------------------------
636     ! Local variables
637     !-----------------------------------------------
638     ! Reynolds number
639           DOUBLE PRECISION :: RE
640     ! Single sphere drag coefficient
641           DOUBLE PRECISION :: C_d
642     ! Gidaspow switch function variables
643           DOUBLE PRECISION :: Ergun
644           DOUBLE PRECISION :: WenYu
645           DOUBLE PRECISION :: PHI_gs
646     !-----------------------------------------------
647     
648           IF(Mug > ZERO) THEN
649     ! Note the presence of gas volume fraction in ROPG
650              RE = DPM*VREL*ROPg/Mug
651           ELSE
652              RE = LARGE_NUMBER
653           ENDIF
654     
655     ! Dense phase - EP_g <= 0.8
656           Ergun = 150D0*(ONE-EPg)*Mug / (EPg*DPM**2) + &
657                    1.75D0*ROg*VREL/DPM
658     ! Dilute phase - EP_g >= 0.8
659           IF(RE <= 1000D0)THEN
660              C_d = (24.D0/(RE+SMALL_NUMBER)) * &
661                (ONE + 0.15D0 * RE**0.687D0)
662           ELSE
663              C_d = 0.44D0
664           ENDIF
665           WenYu = 0.75D0*C_d*VREL*ROPg*EPg**(-2.65D0) / DPM
666     
667     ! Switch function
668           PHI_gs = ATAN(150.D0*1.75D0*(EPg - 0.8D0))/PI + 0.5D0
669     
670     ! Blend the models
671           lDgA = (1.D0-PHI_gs)*Ergun + PHI_gs*WenYu
672           IF (RE == ZERO) lDgA = ZERO
673     
674           RETURN
675           END SUBROUTINE DRAG_GIDASPOW_BLEND
676     
677     
678     
679     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
680     !                                                                      C
681     !  Subroutine: DRAG_WEN_YU                                             C
682     !  Purpose: Calculate the gas-solids drag coefficient                  C
683     !                                                                      C
684     !  Literature/Document References:                                     C
685     !     Wen CY, Yu YH (1966). Chemical Engineering Progress Symposium    C
686     !        Series 62: 100-111.                                           C
687     !                                                                      C
688     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
689     
690           SUBROUTINE DRAG_WEN_YU(lDgA,EPg,Mug,ROPg,VREL,DPM)
691     
692     !-----------------------------------------------
693     ! Modules
694     !-----------------------------------------------
695           USE param
696           USE param1
697           IMPLICIT NONE
698     !-----------------------------------------------
699     ! Dummy arguments
700     !-----------------------------------------------
701     ! drag coefficient
702           DOUBLE PRECISION, INTENT(OUT) :: lDgA
703     ! gas volume fraction
704           DOUBLE PRECISION, INTENT(IN) :: EPg
705     ! gas laminar viscosity
706           DOUBLE PRECISION, INTENT(IN) :: Mug
707     ! gas density*EP_g
708           DOUBLE PRECISION, INTENT(IN) :: ROPg
709     ! Magnitude of gas-solids relative velocity
710           DOUBLE PRECISION, INTENT(IN) :: VREL
711     ! particle diameter of solids phase M or
712     ! average particle diameter if PCF
713           DOUBLE PRECISION, INTENT(IN) :: DPM
714     !-----------------------------------------------
715     ! Local variables
716     !-----------------------------------------------
717     ! Reynolds number
718           DOUBLE PRECISION :: RE
719     ! Single sphere drag coefficient
720           DOUBLE PRECISION :: C_d
721     !-----------------------------------------------
722     
723           IF(Mug > ZERO) THEN
724     ! Note the presence of gas volume fraction in ROPG
725              RE = DPM*VREL*ROPg/Mug
726           ELSE
727              RE = LARGE_NUMBER
728           ENDIF
729     
730           IF(RE <= 1000.0D0)THEN
731              C_d = (24.D0/(RE+SMALL_NUMBER)) * (ONE + 0.15D0*RE**0.687D0)
732           ELSE
733              C_d = 0.44D0
734           ENDIF
735     
736           lDgA = 0.75D0 * C_d * VREL * ROPg * EPg**(-2.65D0) / DPM
737           IF (RE == ZERO) lDgA = ZERO
738     
739           RETURN
740           END SUBROUTINE DRAG_WEN_YU
741     
742     
743     
744     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
745     !                                                                      C
746     !  Subroutine: SUBGRID_DRAG_IGCI                                       C
747     !  Purpose: Calculate subgrid correction to the gas-solids drag        C
748     !           coefficient developed by Wen-Yu                            C
749     !                                                                      C
750     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
751     !                                                                      C
752     !  Revision: 1                                                         C
753     !  Purpose: Minor changes, e.g., fix inconsistenty with analogous      C
754     !     calls in des_drag_gp and with new variable density feature       C
755     !  Author: Janine Carney, June 2013                                    C
756     !                                                                      C
757     !  Literature/Document References:                                     C
758     !     Igci, Y., Pannala, S., Benyahia, S., & Sundaresan S.,            C
759     !        Validation studies on filtered model equations for gas-       C
760     !        particle flows in risers, Industrial & Engineering Chemistry  C
761     !        Research, 2012, 51(4), 2094-2103                              C
762     !                                                                      C
763     !                                                                      C
764     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
765     
766           SUBROUTINE SUBGRID_DRAG_IGCI(lDgA,EPg,Mug,ROg,DPM,ROs,IJK)
767     
768     !-----------------------------------------------
769     ! Modules
770     !-----------------------------------------------
771           USE param
772           USE param1
773           USE run, only : filter_size_ratio, SUBGRID_WALL
774           USE constant, only : GRAVITY
775           USE geometry, only : VOL,AXY,DO_K
776           IMPLICIT NONE
777     !-----------------------------------------------
778     ! Dummy arguments
779     !-----------------------------------------------
780     ! drag coefficient
781           DOUBLE PRECISION, INTENT(INOUT) :: lDgA
782     ! gas volume fraction
783           DOUBLE PRECISION, INTENT(IN) :: EPg
784     ! gas laminar viscosity
785           DOUBLE PRECISION, INTENT(IN) :: Mug
786     ! gas density
787           DOUBLE PRECISION, INTENT(IN) :: ROg
788     ! particle diameter of solids phase M
789           DOUBLE PRECISION, INTENT(IN) :: DPM
790     ! particle density of solids phase M
791           DOUBLE PRECISION, INTENT(IN) :: ROs
792     ! current cell index
793           INTEGER, INTENT(IN) :: IJK
794     !-----------------------------------------------
795     ! Local variables
796     !-----------------------------------------------
797     ! factor to correct the drag for subgrid domain effects
798           DOUBLE PRECISION :: F_Subgrid
799     ! factor to correct the drag for subgrid domain effects arising from
800     ! wall
801           DOUBLE PRECISION :: F_SubGridWall
802     ! particle terminal settling velocity from stokes' formulation
803           DOUBLE PRECISION :: vt
804     ! filter size which is a function of each grid cell volume
805           DOUBLE PRECISION :: filtersize
806     ! inverse Froude number, or dimensionless filter size
807           DOUBLE PRECISION :: Inv_Froude
808     ! total solids volume fraction
809           DOUBLE PRECISION :: EPs
810     ! Variables for Igci model
811           DOUBLE PRECISION :: GG_phip, h_phip, h_phip2, c_function,&
812                               f_filter
813     !-----------------------------------------------
814     
815     ! initialize
816           F_Subgrid = ONE
817           F_SubgridWall = ONE
818     
819     ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
820           vt = GRAVITY*DPM*DPM*(ROs - ROg) / (18.0d0*Mug)
821     ! filter size calculation for each specific gridcell volume
822           IF(DO_K) THEN
823              filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
824           ELSE
825              filtersize = filter_size_ratio * DSQRT(AXY(IJK))
826           ENDIF
827     
828     ! dimensionless inverse of Froude number
829           IF(ABS(vt) > SMALL_NUMBER) THEN
830              Inv_Froude =  filtersize * GRAVITY / vt**2
831           ELSE
832              Inv_Froude =  LARGE_NUMBER
833           ENDIF
834     
835     ! total solids volume fraction
836           EPs = ONE - EPg
837     
838           IF (EPs .LT. 0.0012d0) THEN
839              h_phip = 2.7d0*(EPs**0.234)
840           ELSEIF (EPs .LT. 0.014d0) THEN
841              h_phip = -0.019d0/(EPs**0.455) + 0.963d0
842           ELSEIF (EPs .LT. 0.25d0) THEN
843              h_phip = 0.868d0*EXP((-0.38*EPs)) - &
844                 0.176d0*EXP((-119.2*EPs))
845           ELSEIF (EPs .LT. 0.455d0) THEN
846              h_phip = -4.59d-5*EXP((19.75*EPs)) + &
847                 0.852d0*EXP((-0.268*EPs))
848           ELSEIF (EPs .LE. 0.59d0) THEN
849              h_phip = (EPs - 0.59d0) * (-1501.d0*(EPs**3) + &
850                 2203.d0*(EPs**2) - 1054.d0*EPs + 162.d0)
851           ELSE
852              h_phip=ZERO
853           ENDIF
854     
855           IF (EPs .LT. 0.18d0) THEN
856               GG_phip = (EPs**0.24)*(1.48d0 + EXP(-18.0*EPs))
857           ELSE
858               GG_phip = ONE
859           ENDIF
860     
861     ! a filter function needed in Igci Filtered/subgrid Model [dimensionless]
862           f_filter = (Inv_Froude**1.6) / ((Inv_Froude**1.6)+0.4d0)
863           h_phip2=h_phip*GG_phip
864           c_function=-h_phip2*f_filter
865           F_Subgrid =(ONE + c_function)
866     
867           IF (SUBGRID_WALL) THEN
868              CALL SUBGRID_DRAG_WALL(F_SubgridWall,vt,IJK)
869           ENDIF
870     
871           lDgA = F_SubgridWall*F_Subgrid * lDgA
872     
873           RETURN
874           END SUBROUTINE SUBGRID_DRAG_IGCI
875     
876     
877     
878     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
879     !                                                                      C
880     !  Subroutine: SUBGRID_DRAG_MILIOLI                                    C
881     !  Purpose: Calculate subgrid correction to the gas-solids drag        C
882     !           coefficient developed by Wen-Yu                            C
883     !                                                                      C
884     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
885     !                                                                      C
886     !  Revision: 1                                                         C
887     !  Purpose: Minor changes, e.g., fix inconsistenty with analogous      C
888     !     calls in des_drag_gp and with new variable density feature       C
889     !  Author: Janine Carney, June 2013                                    C
890     !                                                                      C
891     !  Literature/Document References:                                     C
892     !     Milioli, C. C., et al., Filtered two-fluid models of fluidized   C
893     !        gas-particle flows: new constitutive relations, AICHE J,      C
894     !        doi: 10.1002/aic.14130                                        C
895     !                                                                      C
896     !  Comments:                                                           C
897     !     Still needs to be reviewed for accuracy with source material     C
898     !                                                                      C
899     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
900     
901           SUBROUTINE SUBGRID_DRAG_MILIOLI(lDgA,EPg,Mug,ROg,VREL,DPM,ROs,&
902              IJK)
903     
904     !-----------------------------------------------
905     ! Modules
906     !-----------------------------------------------
907           USE param
908           USE param1
909           USE run, only : filter_size_ratio, SUBGRID_WALL
910           USE constant, only : GRAVITY
911           USE geometry, only : VOL,AXY,DO_K
912           IMPLICIT NONE
913     !-----------------------------------------------
914     ! Dummy arguments
915     !-----------------------------------------------
916     ! drag coefficient
917           DOUBLE PRECISION, INTENT(INOUT) :: lDgA
918     ! gas volume fraction
919           DOUBLE PRECISION, INTENT(IN) :: EPg
920     ! gas laminar viscosity
921           DOUBLE PRECISION, INTENT(IN) :: Mug
922     ! gas density
923           DOUBLE PRECISION, INTENT(IN) :: ROg
924     ! Magnitude of gas-solids relative velocity
925           DOUBLE PRECISION, INTENT(IN) :: VREL
926     ! particle diameter of solids phase M
927           DOUBLE PRECISION, INTENT(IN) :: DPM
928     ! particle density of solids phase M
929           DOUBLE PRECISION, INTENT(IN) :: ROs
930     ! current cell index
931           INTEGER, INTENT(IN) :: IJK
932     !-----------------------------------------------
933     ! Local variables
934     !-----------------------------------------------
935     ! factor to correct the drag for subgrid domain effects
936           DOUBLE PRECISION :: F_Subgrid
937     ! factor to correct the drag for subgrid domain effects arising from
938     ! wall
939           DOUBLE PRECISION :: F_SubGridWall
940     ! particle terminal settling velocity from stokes' formulation
941           DOUBLE PRECISION :: vt
942     ! filter size which is a function of each grid cell volume
943           DOUBLE PRECISION :: filtersize
944     ! inverse Froude number, or dimensionless filter size
945           DOUBLE PRECISION :: Inv_Froude
946     ! dimensionless slip velocity = VREL/vt
947           DOUBLE PRECISION :: vslip
948     ! total solids volume fraction
949           DOUBLE PRECISION :: EPs
950     ! Variables for Milioli model
951           DOUBLE PRECISION :: h1, henv, hlin
952     !-----------------------------------------------
953     
954     ! initialize
955           F_Subgrid = ONE
956           F_SubgridWall = ONE
957     
958     ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
959           vt = GRAVITY*DPM*DPM*(ROs - ROg) / (18.0d0*Mug)
960     ! filter size calculation for each specific gridcell volume
961           IF(DO_K) THEN
962              filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
963           ELSE
964              filtersize = filter_size_ratio * DSQRT(AXY(IJK))
965           ENDIF
966     ! dimensionless inverse of Froude number
967           IF(ABS(vt) > SMALL_NUMBER) THEN
968              Inv_Froude =  filtersize * GRAVITY / vt**2
969           ELSE
970              Inv_Froude =  LARGE_NUMBER
971           ENDIF
972     ! total solids volume fractionn
973           EPs = ONE - EPg
974     ! dimensionless slip velocity between gas and solids phase M
975           Vslip = VREL / vt
976     
977           IF (Inv_Froude .LE. 1.028d0) THEN
978              h1 = (1.076d0 + 0.12d0*Vslip - (0.02d0/(Vslip+0.01d0)))*EPs + &
979                 (0.084d0 + 0.09d0*Vslip - (0.01d0/(0.1d0*Vslip+0.01d0)))
980              IF (EPs .LE. 0.53d0) THEN
981                 henv = (6.8d0*(ONE+EPs)*(EPs**0.3)) / &
982                    (10.d0*(EPs**1.5) + 5.d0)
983              ELSEIF (EPs .GT. 0.53d0 .AND. EPs .LE. 0.65d0) THEN
984                 henv = (2.23d0*((0.65d0-EPs)**(0.45))) / &
985                    ((ONE/EPs)-ONE)
986              ELSEIF (EPs .GT. 0.65d0) THEN
987                 henv=ZERO
988              ENDIF
989     
990           ELSEIF (Inv_Froude .GT. 1.028d0 .AND. &
991                   Inv_Froude .LE. 2.056d0) THEN
992              h1 = (1.268d0 - (0.2d0*Vslip) + (0.14d0/(Vslip+0.01d0)))*EPs + &
993                 (0.385d0 + 0.09d0*Vslip - (0.05d0/(0.2d0*Vslip+0.01d0)))
994              IF (EPs .LE. 0.53d0) THEN
995                 henv = (8.6d0*(ONE+EPs)*(EPs**0.2)) / (10.d0*EPs + 6.3d0)
996              ELSEIF (EPs .GT. 0.53d0 .AND. EPs .LE. 0.65d0) THEN
997                 henv = (0.423d0*((0.65d0-EPs)**0.3)) / (ONE-(EPs**0.4))
998              ELSEIF (EPs .GT. 0.65d0) THEN
999                 henv=ZERO
1000              ENDIF
1001     
1002           ELSEIF (Inv_Froude .GT. 2.056d0 .AND. &
1003                   Inv_Froude .LE. 4.112d0) THEN
1004              h1 = ((0.018d0*Vslip + 0.1d0)/(0.14d0*Vslip + 0.01d0))*EPs + &
1005                 (0.9454d0 - (0.09d0/(0.2d0*Vslip + 0.01d0)))
1006              IF (EPs .LE. 0.5d0) THEN
1007                 henv = (7.9d0*(ONE+EPs)*(EPs**0.2)) / &
1008                    (10.d0*(EPs**0.9) + 5.d0)
1009              ELSEIF (EPs .GT. 0.5d0 .AND. EPs .LE. 0.63d0) THEN
1010                 henv = (0.705d0*((0.63d0-EPs)**0.3)) / (ONE-(EPs**0.7))
1011              ELSEIF (EPs .GT. 0.63d0) THEN
1012                 henv=ZERO
1013              ENDIF
1014     
1015           ELSEIF (Inv_Froude .GT. 4.112d0 .AND. &
1016                   Inv_Froude .LE. 8.224d0) THEN
1017              h1 = ((0.05d0*Vslip+0.3d0)/(0.4d0*Vslip+0.06d0))*EPs + &
1018                 (0.9466d0 - (0.05d0/(0.11d0*Vslip+0.01d0)))
1019              IF (EPs .LE. 0.45d0) THEN
1020                 henv = (7.9d0*(ONE+EPs)*(EPs**0.2)) / &
1021                    ((10.d0*(EPs**0.6)) + 3.6d0)
1022              ELSEIF (EPs .GT. 0.45d0 .AND. EPs .LE. 0.57d0) THEN
1023                 henv = (0.78d0*((0.57d0-EPs)**0.2)) / (ONE-(EPs**0.9))
1024              ELSEIF (EPs .GT. 0.57d0) THEN
1025                 henv=ZERO
1026              ENDIF
1027     
1028           ELSEIF (Inv_Froude .GT. 8.224d0 .AND. &
1029                   Inv_Froude .LE. 12.336d0) THEN
1030              h1 = ((1.3d0*Vslip+2.2d0)/(5.2d0*Vslip+0.07d0))*EPs + &
1031                 (0.9363d0-(0.11d0/(0.3d0*Vslip+0.01d0)))
1032              IF (EPs .LE. 0.35d0) THEN
1033                 henv = (7.6d0*(ONE+EPs)*(EPs**0.2)) / &
1034                    ((10.d0*(EPs**0.6)) + 3.3d0)
1035              ELSEIF (EPs .GT. 0.35d0 .AND. EPs .LE. 0.55d0) THEN
1036                 henv = (0.81d0*((0.55d0-EPs)**0.3)) / (ONE-(EPs**0.7))
1037              ELSEIF (EPs .GT. 0.55d0) THEN
1038                 henv=ZERO
1039              ENDIF
1040     
1041           ELSEIF (Inv_Froude .GT. 12.336d0 .AND. &
1042                   Inv_Froude .LE. 16.448d0) THEN
1043              h1 = ((2.6d0*Vslip+4.d0)/(10.d0*Vslip+0.08d0))*EPs + &
1044                 (0.926d0-(0.17d0/(0.5d0*Vslip+0.01d0)))
1045              IF (EPs .LE. 0.25d0) THEN
1046                 henv = (8.4d0*(ONE+EPs)*(EPs**0.2)) / &
1047                    ((10.d0*(EPs**0.5)) + 3.3d0)
1048              ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1049                 henv = (1.01d0*((0.52d0-EPs)**0.03))/(ONE-(EPs**0.9))
1050              ELSEIF (EPs .GT. 0.52d0) THEN
1051                 henv=ZERO
1052              ENDIF
1053     
1054           ELSEIF (Inv_Froude .GT. 16.448d0 .AND. &
1055                   Inv_Froude .LE. 20.56d0) THEN
1056              h1 = ((2.5d0*Vslip+4.d0)/(10.d0*Vslip+0.08d0))*EPs + &
1057                 (0.9261d0-(0.17d0/(0.5d0*Vslip+0.01d0)))
1058              IF (EPs .LE. 0.25d0) THEN
1059                 henv = (8.4d0*(ONE+EPs)*(EPs**0.2)) / &
1060                    ((10.d0*(EPs**0.5)) + 3.3d0)
1061              ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1062                 henv = (1.065d0*((0.52d0-EPs)**0.3))/(ONE-EPs)
1063              ELSEIF (EPs .GT.0.52d0) THEN
1064                 henv=ZERO
1065              ENDIF
1066     
1067           ELSEIF (Inv_Froude .GT. 20.56d0) THEN
1068              h1 = ((1.6d0*Vslip+4.d0)/(7.9d0*Vslip+0.08d0))*EPs + &
1069                 (0.9394d0 - (0.22d0/(0.6d0*Vslip+0.01d0)))
1070              IF (EPs .LE. 0.25d0) THEN
1071                 henv = (9.d0*(ONE+EPs)*(EPs**0.15)) / &
1072                    (10.d0*(EPs**0.45) + 4.2d0)
1073              ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1074                 henv = (0.91d0*((0.52d0-EPs)**0.4))/(ONE-(EPs**0.6))
1075              ELSEIF (EPs .GT. 0.52d0) THEN
1076                 henv=ZERO
1077              ENDIF
1078           ENDIF
1079     
1080           IF (h1 .GT. ZERO) THEN
1081              hlin=h1
1082           ELSE
1083              hlin=ZERO
1084           ENDIF
1085     
1086           IF (Inv_Froude .LT. 1.028d0) THEN
1087     ! for very small filtered size, the drag wont be changed:
1088     ! F_Subgrid = 1.0 - H where H = 0.0
1089              F_Subgrid = ONE
1090           ELSE
1091     ! MIN(henv,hlin) is H in Milioli paper, 2013
1092              F_Subgrid = ONE - MIN(henv,hlin)
1093           ENDIF
1094     
1095     ! Filtered drag = (1 - H)*Microscopic_drag; it is strange Milioli takes EPs
1096     !     F_Subgrid = EPs*(ONE-hmili)
1097     
1098           IF (SUBGRID_WALL) THEN
1099              CALL SUBGRID_DRAG_WALL(F_SubgridWall,vt,IJK)
1100           ENDIF
1101     
1102           lDgA = F_SubgridWall*F_Subgrid * lDgA
1103     
1104     
1105           RETURN
1106           END SUBROUTINE SUBGRID_DRAG_MILIOLI
1107     
1108     
1109     
1110     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1111     !                                                                      C
1112     !  Subroutine: SUBGRID_DRAG_WALL                                       C
1113     !  Purpose: Calculate subgrid correction arising from wall to the      C
1114     !     gas-solids drag coefficient                                      C
1115     !                                                                      C
1116     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
1117     !                                                                      C
1118     !  Revision: 1                                                         C
1119     !  Author: Janine Carney, June 2013                                    C
1120     !                                                                      C
1121     !  Literature/Document References:                                     C
1122     !     Igci, Y., and Sundaresan, S., Verification of filtered two-      C
1123     !        fluid models for gas-particle flows in risers, AICHE J.,      C
1124     !        2011, 57 (10), 2691-2707.                                     C
1125     !                                                                      C
1126     !  Comments: Currently only valid for free-slip wall but no checks     C
1127     !     are made to ensure user has selected free-slip wall when this    C
1128     !     option is invoked                                                C
1129     !                                                                      C
1130     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1131     
1132           SUBROUTINE SUBGRID_DRAG_WALL(lSubgridWall,vt,IJK)
1133     
1134     !-----------------------------------------------
1135     ! Modules
1136     !-----------------------------------------------
1137           USE param
1138           USE param1
1139           USE constant, only : GRAVITY
1140           USE cutcell, only : DWALL
1141           IMPLICIT NONE
1142     !-----------------------------------------------
1143     ! Dummy arguments
1144     !-----------------------------------------------
1145     ! factor to correct the drag for subgrid domain effects arising from
1146     ! wall
1147           DOUBLE PRECISION, INTENT(OUT) :: lSubGridWall
1148     ! particle terminal settling velocity from stokes' formulation
1149           DOUBLE PRECISION, INTENT(IN) :: vt
1150     ! current cell index
1151           INTEGER, INTENT(IN) :: IJK
1152     !-----------------------------------------------
1153     ! Local parameters
1154     !-----------------------------------------------
1155     ! values are only correct for FREE-Slip walls
1156           DOUBLE PRECISION, PARAMETER :: a22=6.0d0, b22=0.295d0
1157     !-----------------------------------------------
1158     ! Local variables
1159     !-----------------------------------------------
1160     ! dimensionless distance to wall
1161           DOUBLE PRECISION :: x_d
1162     !-----------------------------------------------
1163     ! initialize
1164           lSubgridWall = ONE
1165     
1166     ! dimensionless distance to the Wall
1167           x_d = DWALL(IJK) * GRAVITY / vt**2
1168     
1169     ! decrease exponentionally away from the wall
1170     ! more complex model could be implemented with JJ wall model
1171           lSubgridWall = ONE / ( ONE + a22 * (EXP(-b22*x_d)) )
1172     
1173           RETURN
1174           END SUBROUTINE SUBGRID_DRAG_WALL
1175     
1176     
1177     
1178     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1179     !                                                                      C
1180     !  Subroutine: DRAG_KOCH_HILL                                          C
1181     !  Purpose: Calculate the gas-solids drag coefficient                  C
1182     !                                                                      C
1183     !  Author: Clay Sutton (Lehigh University)            Date: 14-Jul-04  C
1184     !                                                                      C
1185     !  Revision: 1                                                         C
1186     !  Author: Sofiane Benyahia                           Date: 21-Jan-05  C
1187     !                                                                      C
1188     !  Literature/Document References:                                     C
1189     !     Benyahia S, Syamlal M, O'Brien TJ (2006). Powder Technology      C
1190     !        162: 166-174.                                                 C
1191     !     Hill RJ, Koch DL, Ladd JC (2001). Journal of Fluid Mechanics     C
1192     !        448: 213-241.                                                 C
1193     !     Hill RJ, Koch DL, Ladd JC (2001). Journal of Fluid Mechanics     C
1194     !        448: 243-278.                                                 C
1195     !                                                                      C
1196     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1197     
1198           SUBROUTINE DRAG_KOCH_HILL(lDgA,EPg,Mug,ROPg,VREL,&
1199                      DPM,DPA,PHIS)
1200     
1201     !-----------------------------------------------
1202     ! Modules
1203     !-----------------------------------------------
1204           USE param
1205           USE param1
1206           IMPLICIT NONE
1207     !-----------------------------------------------
1208     ! Dummy arguments
1209     !-----------------------------------------------
1210     ! drag coefficient
1211           DOUBLE PRECISION, INTENT(OUT) :: lDgA
1212     ! gas volume fraction
1213           DOUBLE PRECISION, INTENT(IN) :: EPg
1214     ! gas laminar viscosity
1215           DOUBLE PRECISION, INTENT(IN) :: Mug
1216     ! gas density*EP_g
1217           DOUBLE PRECISION, INTENT(IN) :: ROPg
1218     ! Magnitude of gas-solids relative velocity
1219           DOUBLE PRECISION, INTENT(IN) :: VREL
1220     ! particle diameter of solids phase M
1221           DOUBLE PRECISION, INTENT(IN) :: DPM
1222     ! average particle diameter if pcf otherwise DPM again
1223           DOUBLE PRECISION, INTENT(IN) :: DPA
1224     ! total solids volume fraction of solids phases
1225           DOUBLE PRECISION, INTENT(IN) :: PHIS
1226     !-----------------------------------------------
1227     ! Local variables
1228     !-----------------------------------------------
1229     ! Reynolds number
1230           DOUBLE PRECISION :: RE
1231     ! transition Reynolds numbers
1232           DOUBLE PRECISION :: Re_Trans_1, Re_Trans_2
1233     ! Stokes Drag Force
1234           DOUBLE PRECISION :: F_STOKES
1235     ! zero Re function for low Reynolds number
1236           DOUBLE PRECISION :: F_0
1237     ! inertial function for low Reynolds number
1238           DOUBLE PRECISION :: F_1
1239     ! zero Re function for high Reynolds number
1240           DOUBLE PRECISION :: F_2
1241     ! inertial function for high Reynolds number
1242           DOUBLE PRECISION :: F_3
1243     ! dimensionless drag force F
1244           DOUBLE PRECISION :: F
1245     ! weighting factor to compute F_0 and F_2
1246           DOUBLE PRECISION :: w
1247     !-----------------------------------------------
1248     
1249     
1250           IF(Mug > ZERO) THEN
1251     ! Note the presence of gas volume fraction in ROPG and factor of 1/2
1252              RE = 0.5D0*DPA*VREL*ROPG/Mug        ! if pcf DPA otherwise DPM
1253           ELSE
1254              RE = LARGE_NUMBER
1255           ENDIF
1256     
1257           F_STOKES = 18.D0*Mug*EPg*EPg/DPM**2    ! use DPM
1258           w = EXP(-10.0D0*(0.4D0-phis)/phis)
1259     
1260           IF(phis > 0.01D0 .AND. phis < 0.4D0) THEN
1261              F_0 = (1.0D0-w) * (1.0D0 + 3.0D0*dsqrt(phis/2.0D0) + &
1262                 135.0D0/64.0D0*phis*LOG(phis) + 17.14D0*phis) / &
1263                 (1.0D0 + 0.681D0*phis - 8.48D0*phis*phis + &
1264                 8.16D0*phis**3) + w*10.0D0*phis/(1.0D0-phis)**3
1265           ELSEIF(phis >= 0.4D0) THEN
1266              F_0 = 10.0D0*phis/(1.0D0-phis)**3
1267           ENDIF
1268     
1269           IF(phis > 0.01D0 .AND. phis <= 0.1D0) THEN
1270             F_1 = dsqrt(2.0D0/phis) / 40.0D0
1271           ELSE IF(phis > 0.1D0) THEN
1272             F_1 = 0.11D0 + 5.1D-04 * exp(11.6D0*phis)
1273           ENDIF
1274     
1275           IF(phis < 0.4D0) THEN
1276             F_2 = (1.0D0-w) * (1.0D0 + 3.0D0*dsqrt(phis/2.0D0) + &
1277                135.0D0/64.0D0*phis*LOG(phis) + 17.89D0*phis) / &
1278                (1.0D0 + 0.681D0*phis - 11.03D0*phis*phis + &
1279                15.41D0*phis**3)+ w*10.0D0*phis/(1.0D0-phis)**3
1280           ELSE
1281              F_2 = 10.0D0*phis/(1.0D0-phis)**3
1282           ENDIF
1283     
1284           IF(phis < 0.0953D0) THEN
1285              F_3 = 0.9351D0*phis + 0.03667D0
1286           ELSE
1287              F_3 = 0.0673D0 + 0.212D0*phis +0.0232D0/(1.0-phis)**5
1288           ENDIF
1289     
1290           Re_Trans_1 = (F_2 - 1.0D0)/(3.0D0/8.0D0 - F_3)
1291           Re_Trans_2 = (F_3 + dsqrt(F_3*F_3 - 4.0D0*F_1 &
1292                *(F_0-F_2))) / (2.0D0*F_1)
1293     
1294           IF(phis <= 0.01D0 .AND. RE <= Re_Trans_1) THEN
1295              F = 1.0D0 + 3.0D0/8.0D0*RE
1296           ELSEIF(phis > 0.01D0 .AND. RE <= Re_Trans_2) THEN
1297              F = F_0 + F_1*RE*RE
1298           ELSEIF(phis <= 0.01D0 .AND. RE > Re_Trans_1 .OR.   &
1299              phis >  0.01D0 .AND. RE > Re_Trans_2) THEN
1300              F = F_2 + F_3*RE
1301           ELSE
1302              F = zero
1303           ENDIF
1304     
1305           lDgA = F * F_STOKES
1306           IF (RE == ZERO) lDgA = ZERO
1307     
1308           RETURN
1309           END SUBROUTINE DRAG_KOCH_HILL
1310     
1311     
1312     
1313     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1314     !                                                                      C
1315     !  Subroutine: DRAG_BVK                                                C
1316     !  Purpose: Calculate the gas-solids drag coefficient                  C
1317     !                                                                      C
1318     !  Literature/Document References:                                     C
1319     !     Beetstra, van der Hoef, Kuipers, Chem. Eng. Science 62           C
1320     !     (Jan 2007)                                                       C
1321     !                                                                      C
1322     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1323     
1324           SUBROUTINE DRAG_BVK(lDgA,EPg,Mug,ROPg,VREL,&
1325                      DPM,DPA,PHIS)
1326     
1327     !-----------------------------------------------
1328     ! Modules
1329     !-----------------------------------------------
1330           USE param
1331           USE param1
1332           IMPLICIT NONE
1333     !-----------------------------------------------
1334     ! Dummy arguments
1335     !-----------------------------------------------
1336     ! drag coefficient
1337           DOUBLE PRECISION, INTENT(OUT) :: lDgA
1338     ! gas volume fraction
1339           DOUBLE PRECISION, INTENT(IN) :: EPg
1340     ! gas laminar viscosity
1341           DOUBLE PRECISION, INTENT(IN) :: Mug
1342     ! gas density*EP_g
1343           DOUBLE PRECISION, INTENT(IN) :: ROPg
1344     ! magnitude of gas-solids relative velocity
1345           DOUBLE PRECISION, INTENT(IN) :: VREL
1346     ! particle diameter of solids phase M or
1347           DOUBLE PRECISION, INTENT(IN) :: DPM
1348     ! average particle diameter
1349           DOUBLE PRECISION, INTENT(IN) :: DPA
1350     ! total solids volume fraction of solids phases
1351           DOUBLE PRECISION, INTENT(IN) :: PHIS
1352     !-----------------------------------------------
1353     ! Local variables
1354     !-----------------------------------------------
1355     ! Reynolds number
1356           DOUBLE PRECISION :: RE
1357     ! Stokes Drag Force
1358           DOUBLE PRECISION :: F_STOKES
1359     ! dimensionless drag force F
1360           DOUBLE PRECISION :: F
1361     !-----------------------------------------------
1362     
1363           IF(Mug > ZERO) THEN
1364     ! Note the presence of gas volume fraction in ROPG
1365              RE = DPA*VREL*ROPg/Mug        ! use DPA
1366           ELSE
1367              RE = LARGE_NUMBER
1368           ENDIF
1369     
1370     ! eq(9) BVK J. fluid. Mech. 528, 2005
1371     ! (this F_Stokes is /= of Koch_Hill by a factor of ep_g)
1372           F_STOKES = 18D0*Mug*EPg/DPM**2   ! use DPM
1373     
1374           F = 10d0*phis/EPg**2 + EPg**2*(ONE+1.5d0*DSQRT(phis))
1375           F = F + 0.413d0*RE/(24.d0*EPg**2) * &
1376                  (ONE/EPg + 3d0*EPg*phis + 8.4d0/RE**0.343) / &
1377                  (ONE+10.d0**(3d0*phis)/RE**(0.5+2.d0*phis))
1378     
1379           lDgA = F*F_STOKES
1380           IF (RE == ZERO) lDgA = ZERO
1381     
1382           RETURN
1383           END SUBROUTINE DRAG_BVK
1384     
1385     
1386     
1387     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1388     !                                                                      C
1389     !  Subroutine: DRAG_HYS                                                C
1390     !  Purpose: Calculate the gas-solids drag coefficient                  C
1391     !                                                                      C
1392     !  Literature/Document References:                                     C
1393     !     Yin, X, Sundaresan, S. (2009). AIChE Journal 55: no 6, 1352-     C
1394     !     1368                                                             C
1395     !                                                                      C
1396     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1397     
1398           SUBROUTINE DRAG_HYS(lDgA,EPg,Mug,ROPg,VREL,&
1399                      DPM,DPA,Y_I,EP_sM,PHIS,M,MAXM,IJK)
1400     
1401     !-----------------------------------------------
1402     ! Modules
1403     !-----------------------------------------------
1404           USE param
1405           USE param1
1406           USE drag, only : beta_ij
1407           USE run, only : LAM_HYS
1408           IMPLICIT NONE
1409     !-----------------------------------------------
1410     ! Dummy arguments
1411     !-----------------------------------------------
1412     ! drag coefficient
1413           DOUBLE PRECISION, INTENT(OUT) :: lDgA
1414     ! gas volume fraction
1415           DOUBLE PRECISION, INTENT(IN) :: EPg
1416     ! gas laminar viscosity
1417           DOUBLE PRECISION, INTENT(IN) :: Mug
1418     ! gas density*EP_g
1419           DOUBLE PRECISION, INTENT(IN) :: ROPg
1420     ! magnitude of gas-solids relative velocity
1421           DOUBLE PRECISION, INTENT(IN) :: VREL
1422     ! local variable for the particle diameter
1423           DOUBLE PRECISION :: DPM(2*DIM_M)
1424     ! average particle diameter
1425           DOUBLE PRECISION, INTENT(IN) :: DPA
1426     ! diameter ratio in polydisperse systems
1427           DOUBLE PRECISION, INTENT(IN) :: Y_i
1428     ! local variable for the solids volume fraction
1429           DOUBLE PRECISION :: EP_sM(2*DIM_M)
1430     ! total solids volume fraction of solids phases
1431           DOUBLE PRECISION, INTENT(IN) :: PHIS
1432     ! current solids phase index and fluid cell index
1433           INTEGER, INTENT(IN) :: M, IJK
1434     ! maximum number of solids phases
1435           INTEGER, INTENT(IN) :: MAXM
1436     !-----------------------------------------------
1437     ! Local variables
1438     !-----------------------------------------------
1439     ! minimum particle diameter in mixture
1440           DOUBLE PRECISION :: DPmin
1441     ! Index for particles of other species
1442           INTEGER :: L
1443     ! Reynolds number
1444           DOUBLE PRECISION :: RE
1445     ! Stokes Drag Force
1446           DOUBLE PRECISION :: F_STOKES
1447     ! dimensionless drag force F
1448           DOUBLE PRECISION :: F
1449     ! Polydisperse correction factor for YS drag relation
1450           DOUBLE PRECISION :: a_YS
1451     ! Lubrication interaction prefactor in YS drag relation
1452           DOUBLE PRECISION :: alpha_YS
1453     ! Friction coefficient for a particle of type i (HYS drag relation)
1454           DOUBLE PRECISION :: beta_i_HYS
1455     ! Friction coefficient for a particle of type j (HYS drag relation)
1456           DOUBLE PRECISION :: beta_j_HYS
1457     ! Stokes drag of a particle of type j
1458           DOUBLE PRECISION :: FSTOKES_j
1459     ! Diameter ratio for particle of type j
1460           DOUBLE PRECISION :: Y_i_J
1461     ! Variable for Beetstra et. al. drag relation
1462           DOUBLE PRECISION :: F_D_BVK
1463     ! Variable for YS drag relation
1464           DOUBLE PRECISION :: F_YS
1465     !-----------------------------------------------
1466     
1467           IF (Mug > ZERO) THEN
1468     ! Note the presence of gas volume fraction in ROPG
1469              RE = DPA*VREL*ROPg/Mug   ! use DPA
1470           ELSE
1471              RE = LARGE_NUMBER
1472           ENDIF
1473     
1474     ! (this F_Stokes is /= of Koch_Hill by a factor of ep_g/ep_sm)
1475     ! (this F_Stokes is /= of BVK by a factor of 1/ep_sm)
1476           F_STOKES = 18D0*Mug*EPg*EP_SM(M)/DPM(M)**2   ! use DPM
1477     
1478     ! Find smallest diameter if number of particle types is greater than 1
1479           Dpmin= DPM(1)
1480           IF (MAXM > 1) THEN
1481              DO L=2,MAXM
1482                 Dpmin = MIN(Dpmin,DPM(L))
1483              ENDDO
1484           ENDIF
1485     
1486           a_YS = 1d0 - 2.66d0*phis + 9.096d0*phis**2 - 11.338d0*phis**3
1487     
1488     ! Calculate the prefactor of the off-diagonal friction coefficient
1489     ! Use default value of lamdba if there are no particle asparities
1490           alpha_YS = 1.313d0*LOG10(DPmin/lam_HYS) - 1.249d0
1491     
1492     
1493     ! Beetstra correction for monodisperse drag
1494           F_D_BVK = ZERO
1495           F = 10.d0*phis/EPg**2 + EPg**2*(ONE+1.5d0*DSQRT(phis))
1496     
1497           IF(RE > ZERO) F_D_BVK = F + 0.413d0*RE/(24.d0*EPg**2)*&
1498              (ONE/EPg + 3.d0*EPg*phis + 8.4d0/RE**0.343d0) / &
1499              (ONE+10**(3.d0*phis)/RE**(0.5d0+2.d0*phis))
1500     
1501     ! YS correction for polydisperse drag
1502           F_YS = 1d0/EPg + (F_D_BVK - 1d0/EPg)*&
1503                            (a_YS*Y_i+(1d0-a_YS)*Y_i**2)
1504           F_YS = F_YS*F_STOKES
1505           beta_i_HYS = F_YS
1506     
1507           DO L= 1,MAXM
1508              IF (L /= M) THEN
1509                 Y_i_J = DPM(L)/DPA
1510                 beta_j_HYS = 1.d0/EPg + (F_D_BVK - 1.d0/EPg) * &
1511                    (a_YS*Y_i_J + (1d0-a_YS)*Y_i_J**2)
1512                 FSTOKES_j = 18.D0*Mug*EP_sM(L)*EPg/&
1513                    DPM(L)**2
1514     
1515                 beta_j_HYS = beta_j_HYS*FSTOKES_j
1516     
1517     ! Calculate off-diagonal friction coefficient
1518                 beta_ij(IJK,M,L) = ZERO
1519     
1520     ! This if statement prevents NaN values from appearing for beta_ij
1521                 IF (EP_sM(M) > ZERO .AND. EP_SM(L) > ZERO) &
1522                    beta_ij(IJK,M,L) = (2.d0*alpha_YS*EP_sM(M)*EP_sM(L))/ &
1523                       (EP_sM(M)/beta_i_HYS + EP_sM(L)/beta_j_HYS)
1524                 F_YS = F_YS + beta_ij(IJK,M,L)
1525     
1526              ENDIF   ! end if (J/=M)
1527           ENDDO   ! end do (J=1,MAXM)
1528     
1529     
1530           lDgA = F_YS
1531     
1532           RETURN
1533           END SUBROUTINE DRAG_HYS
1534