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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_GRBDRY                                             C
4     !  Purpose: Main controller subroutine for calculating coefficients    C
5     !     for momentum boundary conditions using kinetic & frictional      C
6     !     theory                                                           C
7     !                                                                      C
8     !  Author: K. Agrawal & A. Srivastava, Princeton Univ. Date: 19-JAN-98 C
9     !  Reviewer:                                           Date:           C
10     !                                                                      C
11     !                                                                      C
12     !  Literature/Document References:                                     C
13     !                                                                      C
14     !                                                                      C
15     !  Notes:                                                              C
16     !    This routine will not be entered unless granular_energy so we     C
17     !    do not need to check for its condition!                           C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20     
21           SUBROUTINE CALC_GRBDRY(IJK1, IJK2, FCELL, COM, M, L, Gw, Hw, Cw)
22     
23     !-----------------------------------------------
24     ! Modules
25     !-----------------------------------------------
26           USE bc
27           USE compar
28           USE constant
29           USE fldvar
30           USE fun_avg
31           USE functions
32           USE geometry
33           USE indices
34           USE mpi_utility
35           USE param
36           USE param1
37           USE physprop
38           USE rdf
39           USE run
40           USE toleranc
41           USE turb
42           USE visc_s
43           IMPLICIT NONE
44     
45     !-----------------------------------------------
46     ! Dummy Arguments
47     !-----------------------------------------------
48     ! IJK indices for wall cell (ijk1) and fluid cell (ijk2)
49           INTEGER, INTENT(IN) :: IJK1, IJK2
50     ! The location (e,w,n...) of fluid cell
51           CHARACTER, INTENT(IN) :: FCELL
52     ! Velocity component (U, V, W)
53           CHARACTER :: COM
54     ! Solids phase index
55           INTEGER, INTENT(IN) :: M
56     ! Index corresponding to boundary condition
57           INTEGER, INTENT(IN) ::  L
58     ! Wall momentum coefficients:
59     ! 1st term on LHS
60           DOUBLE PRECISION, INTENT(INOUT) :: Gw
61     ! 2nd term on LHS
62           DOUBLE PRECISION, INTENT(INOUT) :: Hw
63     ! all terms appearing on RHS
64           DOUBLE PRECISION, INTENT(INOUT) :: Cw
65     !-----------------------------------------------
66     ! Local Variables
67     !-----------------------------------------------
68     ! Other indices
69           INTEGER :: IJK2E, IPJK2, IPJKM2, IPJKP2, IPJMK2, IPJPK2
70           INTEGER :: IJK2N, IJPK2, IJPKM2, IJPKP2, IMJPK2
71           INTEGER :: IJK2T, IJKP2, IJMKP2, IMJKP2
72     ! Solids phase index
73           INTEGER :: MM
74     !
75           DOUBLE PRECISION :: smallTheta
76     ! Average scalars
77           DOUBLE PRECISION :: EPg_avg, Mu_g_avg, RO_g_avg
78     ! Average void fraction at packing
79           DOUBLE PRECISION :: ep_star_avg
80     ! Average scalars modified to include all solid phases
81           DOUBLE PRECISION :: EPs_avg(DIMENSION_M), DP_avg(DIMENSION_M),&
82                               TH_avg(DIMENSION_M), AVGX1, AVGX2
83           DOUBLE PRECISION ROS_AVG(DIMENSION_M)
84     ! Average Simonin and Ahmadi variables (sof)
85           DOUBLE PRECISION :: K_12_avg, Tau_12_avg, Tau_1_avg
86     ! Average velocities
87     ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
88           DOUBLE PRECISION :: USCM, VSCM,WSCM
89           DOUBLE PRECISION :: USCM1,USCM2,VSCM1,VSCM2,WSCM1,WSCM2
90     ! values of U_g, V_g, W_g at appropriate place on boundary wall
91           DOUBLE PRECISION :: UGC, VGC, WGC
92           DOUBLE PRECISION :: WGC1, WGC2, VGC1, VGC2, UGC1, UGC2
93     ! velocity variables used to standarize dummy argument for different dirs
94           DOUBLE PRECISION :: WVELS
95           DOUBLE PRECISION :: VELS
96     ! del.u
97           DOUBLE PRECISION :: DEL_DOT_U
98     ! S:S
99           DOUBLE PRECISION :: S_DDOT_S
100     ! S_dd (d can be x, y or z)
101           DOUBLE PRECISION :: S_dd
102     ! Magnitude of gas-solids relative velocity
103           DOUBLE PRECISION :: VREL
104     ! slip velocity between wall and particles for Jenkins bc (sof)
105           DOUBLE PRECISION :: VSLIP
106     ! radial distribution function at contact
107           DOUBLE PRECISION :: g0(DIMENSION_M)
108     ! Sum of eps*G_0
109           DOUBLE PRECISION :: g0EPs_avg
110     ! Error message
111           CHARACTER(LEN=80) :: LINE
112     
113     !-----------------------------------------------
114     !  Function subroutines
115     !-----------------------------------------------
116           DOUBLE PRECISION F_HW
117     !-----------------------------------------------
118     
119     ! Note: EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).
120     !       Hence IJK2 (fluid cell) is used in averages.
121     
122           smallTheta = (to_SI)**4 * ZERO_EP_S
123     
124     !-----------------------------------------------
125     ! Calculations for U momentum equation
126           IF (COM .EQ. 'U')THEN
127     
128             IJK2E = EAST_OF(IJK2)
129             IPJK2 = IP_OF(IJK2)
130     
131             EPg_avg = AVG_X(EP_g(IJK2), EP_g(IJK2E), I_OF(IJK2))
132             ep_star_avg = AVG_X(EP_star_array(IJK2), EP_star_array(IJK2E), I_OF(IJK2))
133             Mu_g_avg = AVG_X(Mu_g(IJK2), Mu_g(IJK2E), I_OF(IJK2))
134             RO_g_avg = AVG_X(RO_g(IJK2), RO_g(IJK2E), I_OF(IJK2))
135     
136             K_12_avg = ZERO
137             Tau_12_avg = ZERO
138             Tau_1_avg = ZERO
139             IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
140                KT_TYPE_ENUM == AHMADI_1995) THEN
141                Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
142     ! Simonin only:
143                K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
144                Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
145             ENDIF
146     
147             g0EPs_avg = ZERO
148             DO MM = 1, SMAX
149                g0(MM)      = G_0AVG(IJK2, IJK2E, 'X', I_OF(IJK2), M, MM)
150                EPs_avg(MM) = AVG_X(EP_s(IJK2, MM), EP_s(IJK2E, MM), I_OF(IJK2))
151                DP_avg(MM)  = AVG_X(D_P(IJK2,MM), D_P(IJK2E,MM), I_OF(IJK2))
152                g0EPs_avg   = g0EPs_avg + G_0AVG(IJK2, IJK2E, 'X', I_OF(IJK2), M, MM) &
153                            * AVG_X(EP_s(IJK2, MM), EP_s(IJK2E, MM), I_OF(IJK2))
154                ROs_avg(MM) = AVG_X(RO_S(IJK2, MM), RO_S(IJK2E, MM), I_OF(IJK2))
155     
156                TH_avg(MM) = AVG_X(THETA_M(IJK2,MM), THETA_M(IJK2E,MM), I_OF(IJK2))
157                IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
158             ENDDO
159     
160             WVELS = BC_Uw_s(L,M)
161     
162     
163             IF(FCELL .EQ. 'N')THEN
164               IPJMK2 = JM_OF(IPJK2)
165     
166     ! code modified for some corner cells
167               DO MM = 1, SMAX
168                  AVGX1 = AVG_X(Theta_m(IJK1,MM), Theta_m(IPJMK2,MM), I_OF(IJK1))
169                  AVGX2 = AVG_X(Theta_m(IJK2,MM), Theta_m(IPJK2,MM), I_OF(IJK2))
170                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
171                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
172                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
173                     TH_avg(MM) = smallTheta
174                  ELSE
175                     TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK1))
176                  ENDIF
177               ENDDO
178     
179     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)
180               UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
181               VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
182               WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
183                            AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
184                            I_OF(IJK2))
185               WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
186                            AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
187                            I_OF(IJK1))
188               WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
189               USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
190               VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
191               WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
192                            AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
193                            I_OF(IJK2))
194               WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
195                            AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
196                            I_OF(IJK1))
197               WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
198               VELS = USCM
199     
200             ELSEIF(FCELL .EQ. 'S')THEN
201               IPJPK2= JP_OF(IPJK2)
202     
203     ! code modified for some corner cells
204               DO MM = 1, SMAX
205                  AVGX1 = AVG_X(Theta_m(IJK2,MM),Theta_m(IPJK2,MM),I_OF(IJK2))
206                  AVGX2 = AVG_X(Theta_m(IJK1,MM),Theta_m(IPJPK2,MM),I_OF(IJK1))
207                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
208                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
209                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
210                     TH_avg(MM) = smallTheta
211                  ELSE
212                     TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK2))
213                  ENDIF
214               ENDDO
215     
216     
217     ! Calculate velocity components at i+1/2, j+1/2, k relative to IJK2
218               UGC  = AVG_Y(U_g(IJK2),U_g(IJK1),J_OF(IJK2))
219               VGC  = AVG_X(V_g(IJK2),V_g(IPJK2),I_OF(IJK2))
220               WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
221                            AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
222                            I_OF(IJK2))
223               WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
224                            AVG_Z_T(W_g(KM_OF(IPJPK2)), W_g(IPJPK2)),&
225                            I_OF(IJK1))
226               WGC  = AVG_Y(WGC1, WGC2, J_OF(IJK2))
227               USCM = AVG_Y(U_s(IJK2, M),U_s(IJK1, M),J_OF(IJK2))
228               VSCM = AVG_X(V_s(IJK2, M),V_s(IPJK2, M),I_OF(IJK2))
229               WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
230                            AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
231                            I_OF(IJK2))
232               WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
233                            AVG_Z_T(W_s(KM_OF(IPJPK2),M),W_s(IPJPK2,M)),&
234                            I_OF(IJK1))
235               WSCM = AVG_Y(WSCM1, WSCM2, J_OF(IJK2))
236               VELS = USCM
237     
238             ELSEIF(FCELL .EQ. 'T')THEN
239               IPJKM2= KM_OF(IPJK2)
240     
241     ! code modified for some corner cells
242               DO MM = 1, SMAX
243                  AVGX1 = AVG_X(Theta_m(IJK1,MM),Theta_m(IPJKM2,MM),I_OF(IJK1))
244                  AVGX2 = AVG_X(Theta_m(IJK2,MM),Theta_m(IPJK2,MM),I_OF(IJK2))
245                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
246                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
247                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
248                     TH_avg(MM) = smallTheta
249                  ELSE
250                     TH_avg(MM) = AVG_Z(AVGX1, AVGX2, K_OF(IJK1))
251                  ENDIF
252               ENDDO
253     
254     ! Calculate velocity components at i+1/2,j,k-1/2 relative to IJK2
255               UGC  = AVG_Z(U_g(IJK1), U_g(IJK2), K_OF(IJK1))
256               VGC1 = AVG_X(AVG_Y_N(V_g(JM_OF(IJK2)),V_g(IJK2)),&
257                            AVG_Y_N(V_g(JM_OF(IPJK2)),V_g(IPJK2)),&
258                            I_OF(IJK2))
259               VGC2 = AVG_X(AVG_Y_N(V_g(JM_OF(IJK1)),V_g(IJK1)),&
260                            AVG_Y_N(V_g(JM_OF(IPJKM2)),V_g(IPJKM2)),&
261                            I_OF(IJK1))
262               VGC  = AVG_Z(VGC2,VGC1,K_OF(IJK1))
263               WGC  = AVG_X(W_g(IJK1), W_g(IPJKM2),I_OF(IJK1))
264               USCM = AVG_Z(U_s(IJK1,M), U_s(IJK2,M), K_OF(IJK1))
265               VSCM1= AVG_X(AVG_Y_N(V_s(JM_OF(IJK2),M),V_s(IJK2,M)),&
266                            AVG_Y_N(V_s(JM_OF(IPJK2),M),V_s(IPJK2,M)),&
267                            I_OF(IJK2))
268               VSCM2= AVG_X(AVG_Y_N(V_s(JM_OF(IJK1),M),V_s(IJK1,M)),&
269                            AVG_Y_N(V_s(JM_OF(IPJKM2),M),V_s(IPJKM2,M)),&
270                            I_OF(IJK1))
271               VSCM  = AVG_Z(VSCM2,VSCM1,K_OF(IJK1))
272               WSCM = AVG_X(W_s(IJK1,M), W_s(IPJKM2,M), I_OF(IJK1))
273               VELS = USCM
274     
275             ELSEIF(FCELL .EQ. 'B')THEN
276               IPJKP2= KP_OF(IPJK2)
277     
278     ! code modified for some corner cells
279               DO MM = 1, SMAX
280                  AVGX1 = AVG_X(Theta_m(IJK1,MM), Theta_m(IPJKP2,MM),I_OF(IJK1))
281                  AVGX2 = AVG_X(Theta_m(IJK2,MM), Theta_m(IPJK2,MM),I_OF(IJK2))
282                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
283                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
284                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
285                     TH_avg(MM) = smallTheta
286                  ELSE
287                     TH_avg(MM) = AVG_Z(AVGX1, AVGX2, K_OF(IJK2))
288                  ENDIF
289               ENDDO
290     
291     ! Calculate velocity components at i+1/2,j,k-1/2 relative to IJK1
292               UGC  = AVG_Z(U_g(IJK2), U_g(IJK1), K_OF(IJK2))
293               VGC1 = AVG_X(AVG_Y_N(V_g(JM_OF(IJK2)),V_g(IJK2)),&
294                            AVG_Y_N(V_g(JM_OF(IPJK2)),V_g(IPJK2)),&
295                            I_OF(IJK2))
296               VGC2 = AVG_X(AVG_Y_N(V_g(JM_OF(IJK1)),V_g(IJK1)),&
297                            AVG_Y_N(V_g(JM_OF(IPJKP2)),V_g(IPJKP2)),&
298                            I_OF(IJK1))
299               VGC  = AVG_Z(VGC1, VGC2, K_OF(IJK2))
300               WGC  = AVG_X(W_g(IJK2), W_g(IPJK2),I_OF(IJK2))
301               USCM = AVG_Z(U_s(IJK2, M), U_s(IJK1, M), K_OF(IJK2))
302               VSCM1= AVG_X(AVG_Y_N(V_s(JM_OF(IJK2),M),V_s(IJK2,M)),&
303                            AVG_Y_N(V_s(JM_OF(IPJK2),M),V_s(IPJK2,M)),&
304                            I_OF(IJK2))
305               VSCM2= AVG_X(AVG_Y_N(V_s(JM_OF(IJK1),M),V_s(IJK1,M)),&
306                            AVG_Y_N(V_s(JM_OF(IPJKP2),M),V_s(IPJKP2,M)),&
307                            I_OF(IJK1))
308               VSCM = AVG_Z(VSCM1, VSCM2, K_OF(IJK2))
309               WSCM = AVG_X(W_s(IJK2, M), W_s(IPJK2, M),I_OF(IJK2))
310               VELS = USCM
311     
312             ELSE
313               WRITE(LINE,'(A, A)') 'Error: Unknown FCELL'
314               CALL WRITE_ERROR('CALC_GRBDRY', LINE, 1)
315               CALL exitMPI(myPE)
316             ENDIF
317     
318     !-----------------------------------------------
319     ! Calculations for V momentum equation
320           ELSEIF (COM .EQ. 'V')THEN
321     
322             IJK2N = NORTH_OF(IJK2)
323             IJPK2 = JP_OF(IJK2)
324     
325             EPg_avg = AVG_Y(EP_g(IJK2), EP_g(IJK2N), J_OF(IJK2))
326             ep_star_avg = AVG_Y(EP_star_array(IJK2), EP_star_array(IJK2N), J_OF(IJK2))
327             Mu_g_avg = AVG_Y(Mu_g(IJK2), Mu_g(IJK2N), J_OF(IJK2))
328             RO_g_avg = AVG_Y(RO_g(IJK2), RO_g(IJK2N), J_OF(IJK2))
329     
330             K_12_avg = ZERO
331             Tau_12_avg = ZERO
332             Tau_1_avg = ZERO
333             IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
334                KT_TYPE_ENUM == AHMADI_1995) THEN
335                Tau_1_avg = AVG_Y(Tau_1(IJK2), Tau_1(IJK2N), J_OF(IJK2))
336     ! Simonin only:
337                K_12_avg = AVG_Y(K_12(IJK2), K_12(IJK2N), J_OF(IJK2))
338                Tau_12_avg = AVG_Y(Tau_12(IJK2), Tau_12(IJK2N), J_OF(IJK2))
339             ENDIF
340     
341     
342             g0EPs_avg = ZERO
343             DO MM = 1, SMAX
344                g0(MM)      = G_0AVG(IJK2, IJK2N, 'Y', J_OF(IJK2), M, MM)
345                EPs_avg(MM) = AVG_Y(EP_s(IJK2, MM), EP_s(IJK2N, MM), J_OF(IJK2))
346                DP_avg(MM)  = AVG_Y(D_p(IJK2,MM), D_p(IJK2N,MM), J_OF(IJK2))
347                g0EPs_avg   = g0EPs_avg + G_0AVG(IJK2, IJK2N, 'Y', J_OF(IJK2), M, MM) &
348                             * AVG_Y(EP_s(IJK2, MM), EP_s(IJK2N, MM), J_OF(IJK2))
349                ROS_avg(MM) = AVG_Y(RO_S(IJK2, MM), RO_S(IJK2N, MM), J_OF(IJK2))
350                TH_avg(MM) = AVG_Y(THETA_M(IJK2,MM), THETA_M(IJK2N,MM), J_OF(IJK2))
351                IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
352             ENDDO
353     
354             WVELS = BC_Vw_s(L, M)
355     
356             IF(FCELL .EQ. 'T')THEN
357               IJPKM2 = KM_OF(IJPK2)
358     
359               DO MM = 1, SMAX
360                  AVGX1 = AVG_Z(Theta_m(IJK1,MM), Theta_m(IJK2,MM), K_OF(IJK1))
361                  AVGX2 = AVG_Z(Theta_m(IJPKM2,MM), Theta_m(IJPK2,MM), K_OF(IJPKM2))
362                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
363                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
364                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
365                     TH_avg(MM) = smallTheta
366                  ELSE
367                     TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK1))
368                  ENDIF
369               ENDDO
370     
371     ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK1)
372               UGC1 = AVG_X_E(&
373                              AVG_Y(U_g(IM_OF(IJK1)), U_g(IM_OF(IJPKM2)),&
374                              J_OF(IM_OF(IJK1))),&
375                              AVG_Y(U_g(IJK1), U_g(IJPKM2),J_OF(IJK1)),&
376                              I_OF(IJK1))
377               UGC2 = AVG_X_E(&
378                              AVG_Y(U_g(IM_OF(IJK2)), U_g(IM_OF(IJPK2)),&
379                              J_OF(IM_OF(IJK2))),&
380                              AVG_Y(U_g(IJK2), U_g(IJPK2),J_OF(IJK2)),&
381                              I_OF(IJK2))
382               UGC  = AVG_Z(UGC1, UGC2, K_OF(IJK1))
383               VGC  = AVG_Z(V_g(IJK1), V_g(IJK2),K_OF(IJK1))
384               WGC  = AVG_Y(W_g(IJK1), W_g(IJPKM2), J_OF(IJK1))
385               USCM1= AVG_X_E(&
386                            AVG_Y(U_s(IM_OF(IJK1),M),U_s(IM_OF(IJPKM2),M),&
387                            J_OF(IM_OF(IJK1))),&
388                            AVG_Y(U_s(IJK1,M), U_s(IJPKM2,M),J_OF(IJK1)),&
389                            I_OF(IJK1))
390               USCM2= AVG_X_E(&
391                            AVG_Y(U_s(IM_OF(IJK2),M),U_s(IM_OF(IJPK2),M),&
392                            J_OF(IM_OF(IJK2))),&
393                            AVG_Y(U_s(IJK2,M), U_s(IJPK2,M),J_OF(IJK2)),&
394                            I_OF(IJK2))
395               USCM = AVG_Z(USCM1, USCM2, K_OF(IJK1))
396               VSCM = AVG_Z(V_s(IJK1,M), V_s(IJK2,M),K_OF(IJK1))
397               WSCM = AVG_Y(W_s(IJK1,M), W_s(IJPKM2,M), J_OF(IJK1))
398               VELS = VSCM
399     
400             ELSEIF(FCELL .EQ. 'B')THEN
401               IJPKP2 = KP_OF(IJPK2)
402     
403               DO MM = 1, SMAX
404                  AVGX1 = AVG_Z(Theta_m(IJK2,MM), Theta_m(IJK1,MM), K_OF(IJK2))
405                  AVGX2 = AVG_Z(Theta_m(IJPK2,MM), Theta_m(IJPKP2,MM), K_OF(IJPK2))
406                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
407                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
408                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
409                     TH_avg(MM) = smallTheta
410                  ELSE
411                     TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK2))
412                  ENDIF
413               ENDDO
414     
415     ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK2)
416               UGC1 = AVG_X_E(&
417                              AVG_Y(U_g(IM_OF(IJK1)), U_g(IM_OF(IJPKP2)),&
418                              J_OF(IM_OF(IJK1))),&
419                              AVG_Y(U_g(IJK1), U_g(IJPKP2),J_OF(IJK1)),&
420                              I_OF(IJK1))
421               UGC2 = AVG_X_E(&
422                              AVG_Y(U_g(IM_OF(IJK2)), U_g(IM_OF(IJPK2)),&
423                              J_OF(IM_OF(IJK2))),&
424                              AVG_Y(U_g(IJK2), U_g(IJPK2),J_OF(IJK2)),&
425                              I_OF(IJK2))
426               UGC  = AVG_Z(UGC2, UGC1, K_OF(IJK2))
427               VGC  = AVG_Z(V_g(IJK2), V_g(IJK1),K_OF(IJK2))
428               WGC  = AVG_Y(W_g(IJK2), W_g(IJPK2), J_OF(IJK2))
429               USCM1= AVG_X_E(&
430                            AVG_Y(U_s(IM_OF(IJK1),M),U_s(IM_OF(IJPKP2),M),&
431                            J_OF(IM_OF(IJK1))),&
432                            AVG_Y(U_s(IJK1,M), U_s(IJPKP2,M),J_OF(IJK1)),&
433                            I_OF(IJK1))
434               USCM2= AVG_X_E(&
435                            AVG_Y(U_s(IM_OF(IJK2),M),U_s(IM_OF(IJPK2),M),&
436                            J_OF(IM_OF(IJK2))),&
437                            AVG_Y(U_s(IJK2,M), U_s(IJPK2,M),J_OF(IJK2)),&
438                            I_OF(IJK2))
439               USCM = AVG_Z(USCM2, USCM1, K_OF(IJK2))
440               VSCM = AVG_Z(V_s(IJK2,M), V_s(IJK1,M),K_OF(IJK2))
441               WSCM = AVG_Y(W_s(IJK2,M), W_s(IJPK2,M), J_OF(IJK2))
442               VELS = VSCM
443     
444             ELSEIF(FCELL .EQ. 'E')THEN
445               IMJPK2= IM_OF(IJPK2)
446     
447               DO MM = 1, SMAX
448                  AVGX1 = AVG_X(Theta_m(IJK1,MM),Theta_m(IJK2,MM),I_OF(IJK1))
449                  AVGX2 = AVG_X(Theta_m(IMJPK2,MM),Theta_m(IJPK2,MM),I_OF(IMJPK2))
450                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
451                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
452                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
453                     TH_avg(MM) = smallTheta
454                  ELSE
455                     TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK1))
456                  ENDIF
457               ENDDO
458     
459     ! Calculate velocity components at i+1/2,j+1/2,k relative to IJK1
460               UGC  = AVG_Y(U_g(IJK1), U_g(IMJPK2), J_OF(IJK1))
461               VGC  = AVG_X(V_g(IJK1), V_g(IJK2), I_OF(IJK1))
462               WGC1 = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
463                            AVG_Z_T(W_g(KM_OF(IMJPK2)), W_g(IMJPK2)),&
464                            J_OF(IJK1))
465               WGC2 = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
466                            AVG_Z_T(W_g(KM_OF(IJPK2)), W_g(IJPK2)),&
467                            J_OF(IJK2))
468               WGC  = AVG_X(WGC1, WGC2, I_OF(IJK1))
469               USCM = AVG_Y(U_s(IJK1,M), U_s(IMJPK2,M), J_OF(IJK1))
470               VSCM = AVG_X(V_s(IJK1, M), V_s(IJK2, M), I_OF(IJK1))
471               WSCM1= AVG_Y(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
472                            AVG_Z_T(W_s(KM_OF(IMJPK2),M), W_s(IMJPK2,M)),&
473                            J_OF(IJK1))
474               WSCM2 = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
475                            AVG_Z_T(W_s(KM_OF(IJPK2),M), W_s(IJPK2,M)),&
476                            J_OF(IJK2))
477               WSCM  = AVG_X(WSCM1, WSCM2, I_OF(IJK1))
478               VELS = VSCM
479     
480             ELSEIF(FCELL .EQ. 'W')THEN
481               IPJPK2= IP_OF(IJPK2)
482     
483               DO MM = 1, SMAX
484                  AVGX1 = AVG_X(Theta_m(IJK2,MM),Theta_m(IJK1,MM),I_OF(IJK2))
485                  AVGX2 = AVG_X(Theta_m(IJPK2,MM),Theta_m(IPJPK2,MM),I_OF(IJPK2))
486                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
487                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
488                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
489                     TH_avg(MM) = smallTheta
490                  ELSE
491                     TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK2))
492                  ENDIF
493               ENDDO
494     
495     ! Calculate velocity components at i+1/2,j+1/2,k relative to IJK2
496               UGC  = AVG_Y(U_g(IJK2), U_g(IJPK2), J_OF(IJK2))
497               VGC  = AVG_X(V_g(IJK2), V_g(IJK1), I_OF(IJK2))
498               WGC1 = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
499                            AVG_Z_T(W_g(KM_OF(IPJPK2)), W_g(IPJPK2)),&
500                            J_OF(IJK1))
501               WGC2 = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
502                            AVG_Z_T(W_g(KM_OF(IJPK2)), W_g(IJPK2)),&
503                            J_OF(IJK2))
504               WGC  = AVG_X(WGC2, WGC1, I_OF(IJK2))
505               USCM = AVG_Y(U_s(IJK2,M), U_s(IJPK2,M), J_OF(IJK2))
506               VSCM = AVG_X(V_s(IJK2, M), V_s(IJK1, M), I_OF(IJK2))
507               WSCM1= AVG_Y(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
508                            AVG_Z_T(W_s(KM_OF(IPJPK2),M), W_s(IPJPK2,M)),&
509                            J_OF(IJK1))
510               WSCM2 = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
511                            AVG_Z_T(W_s(KM_OF(IJPK2),M), W_s(IJPK2,M)),&
512                            J_OF(IJK2))
513               WSCM  = AVG_X(WSCM2, WSCM1, I_OF(IJK2))
514               VELS = VSCM
515     
516             ELSE
517               WRITE(LINE,'(A, A)') 'Error: Unknown FCELL'
518               CALL WRITE_ERROR('CALC_GRBDRY', LINE, 1)
519               CALL exitMPI(myPE)
520             ENDIF
521     
522     
523     !-----------------------------------------------
524     ! Calculations for W momentum equation
525           ELSEIF (COM .EQ. 'W')THEN
526             IJK2T = TOP_OF(IJK2)
527             IJKP2 = KP_OF(IJK2)
528     
529             EPg_avg = AVG_Z(EP_g(IJK2), EP_g(IJK2T), K_OF(IJK2))
530             Mu_g_avg = AVG_Z(Mu_g(IJK2), Mu_g(IJK2T), K_OF(IJK2))
531             RO_g_avg = AVG_Z(RO_g(IJK2), RO_g(IJK2T), K_OF(IJK2))
532             ep_star_avg = AVG_Z(EP_star_array(IJK2), EP_star_array(IJK2T), K_OF(IJK2))
533     
534             K_12_avg = ZERO
535             Tau_12_avg = ZERO
536             Tau_1_avg = ZERO
537             IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
538                KT_TYPE_ENUM == AHMADI_1995) THEN
539                Tau_1_avg = AVG_Z(Tau_1(IJK2), Tau_1(IJK2T), K_OF(IJK2))
540     ! Simonin only:
541                K_12_avg = AVG_Z(K_12(IJK2), K_12(IJK2T), K_OF(IJK2))
542                Tau_12_avg = AVG_Z(Tau_12(IJK2), Tau_12(IJK2T), K_OF(IJK2))
543             ENDIF
544     
545             g0EPs_avg = ZERO
546             DO MM = 1, SMAX
547                g0(MM)      = G_0AVG(IJK2, IJK2T, 'Z', K_OF(IJK2), M, MM)
548                EPs_avg(MM) = AVG_Z(EP_s(IJK2,MM), EP_s(IJK2T,MM), K_OF(IJK2))
549                DP_avg(MM)  = AVG_Z(D_p(IJK2,MM), D_p(IJK2T,MM), K_OF(IJK2))
550                g0EPs_avg   = g0EPs_avg + G_0AVG(IJK2, IJK2T, 'Z', K_OF(IJK2), M, MM) &
551                            * AVG_Z(EP_s(IJK2, MM), EP_s(IJK2T, MM), K_OF(IJK2))
552                ROs_avg(MM) = AVG_Z(RO_S(IJK2,MM), RO_S(IJK2T,MM), K_OF(IJK2))
553                TH_avg(MM) = AVG_Z(THETA_M(IJK2,MM), THETA_M(IJK2T,MM), K_OF(IJK2))
554                IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
555             ENDDO
556     
557             WVELS = BC_Ww_s(L,M)
558     
559             IF(FCELL .EQ. 'N')THEN
560               IJMKP2 = JM_OF(IJKP2)
561     
562               DO MM = 1, SMAX
563                  AVGX1 = AVG_Z(Theta_m(IJK1,MM), Theta_m(IJMKP2,MM), K_OF(IJK1))
564                  AVGX2 = AVG_Z(Theta_m(IJK2,MM), Theta_m(IJKP2,MM), K_OF(IJK2))
565                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
566                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
567                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
568                    TH_avg(MM) = smallTheta
569                  ELSE
570                    TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK1))
571                  ENDIF
572               ENDDO
573     
574     ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK1)
575               UGC1 = AVG_X_E(&
576                              AVG_Z(U_g(IM_OF(IJK1)), U_g(IM_OF(IJMKP2)),&
577                              K_OF(IM_OF(IJK1)) ),&
578                              AVG_Z(U_g(IJK1), U_g(IJMKP2), K_OF(IJK1)),&
579                              I_OF(IJK1))
580               UGC2 = AVG_X_E(&
581                              AVG_Z(U_g(IM_OF(IJK2)), U_g(IM_OF(IJKP2)),&
582                              K_OF(IM_OF(IJK2))),&
583                              AVG_Z(U_g(IJK2), U_g(IJKP2), K_OF(IJK2)),&
584                              I_OF(IJK2))
585               UGC  = AVG_Y(UGC1, UGC2, J_OF(IJK1))
586               VGC  = AVG_Z(V_g(IJK1), V_g(IJMKP2),K_OF(IJK1))
587               WGC  = AVG_Y(W_g(IJK1), W_g(IJK2), J_OF(IJK1))
588               USCM1= AVG_X_E(&
589                            AVG_Z(U_s(IM_OF(IJK1),M),U_s(IM_OF(IJMKP2),M),&
590                            K_OF(IM_OF(IJK1))),&
591                            AVG_Z(U_s(IJK1,M), U_s(IJMKP2,M),K_OF(IJK1)),&
592                            I_OF(IJK1))
593               USCM2= AVG_X_E(&
594                            AVG_Z(U_s(IM_OF(IJK2),M),U_s(IM_OF(IJKP2),M),&
595                            K_OF(IM_OF(IJK2))),&
596                            AVG_Z(U_s(IJK2,M), U_s(IJKP2,M),K_OF(IJK2)),&
597                            I_OF(IJK2))
598               USCM = AVG_Y(USCM1, USCM2, J_OF(IJK1))
599               VSCM = AVG_Z(V_s(IJK1,M), V_s(IJMKP2,M),K_OF(IJK1))
600               WSCM = AVG_Y(W_s(IJK1,M), W_s(IJK2,M), J_OF(IJK1))
601               VELS = WSCM
602     
603             ELSEIF(FCELL .EQ. 'S')THEN
604               IJPKP2 = JP_OF(IJKP2)
605     
606               DO MM = 1, SMAX
607                  AVGX1 = AVG_Z(Theta_m(IJK2,MM), Theta_m(IJKP2,MM), K_OF(IJK2))
608                  AVGX2 = AVG_Z(Theta_m(IJK1,MM), Theta_m(IJPKP2,MM), K_OF(IJK1))
609                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
610                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
611                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
612                    TH_avg(MM) = smallTheta
613                  ELSE
614                    TH_avg(MM) = AVG_Y(AVGX1, AVGX2, J_OF(IJK2))
615                  ENDIF
616                ENDDO
617     
618     ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK2)
619               UGC1 = AVG_X_E(&
620                              AVG_Z(U_g(IM_OF(IJK1)), U_g(IM_OF(IJPKP2)),&
621                              K_OF(IM_OF(IJK1))),&
622                              AVG_Z(U_g(IJK1), U_g(IJPKP2), K_OF(IJK1)),&
623                              I_OF(IJK1))
624               UGC2 = AVG_X_E(&
625                              AVG_Z(U_g(IM_OF(IJK2)), U_g(IM_OF(IJKP2)),&
626                              K_OF(IM_OF(IJK2))),&
627                              AVG_Z(U_g(IJK2), U_g(IJKP2), K_OF(IJK2)),&
628                              I_OF(IJK2))
629               UGC  = AVG_Y(UGC2, UGC1, J_OF(IJK2))
630               VGC  = AVG_Z(V_g(IJK2), V_g(IJKP2),K_OF(IJK2))
631               WGC  = AVG_Y(W_g(IJK2), W_g(IJK1), J_OF(IJK2))
632               USCM1= AVG_X_E(&
633                            AVG_Z(U_s(IM_OF(IJK1),M),U_s(IM_OF(IJPKP2),M),&
634                            K_OF(IM_OF(IJK1))),&
635                            AVG_Z(U_s(IJK1,M), U_s(IJPKP2,M),K_OF(IJK1)),&
636                            I_OF(IJK1))
637               USCM2= AVG_X_E(&
638                            AVG_Z(U_s(IM_OF(IJK2),M),U_s(IM_OF(IJKP2),M),&
639                            K_OF(IM_OF(IJK2))),&
640                            AVG_Z(U_s(IJK2,M), U_s(IJKP2,M),K_OF(IJK2)),&
641                            I_OF(IJK2))
642               USCM = AVG_Y(USCM2, USCM1, J_OF(IJK2))
643               VSCM = AVG_Z(V_s(IJK2,M), V_s(IJKP2,M),K_OF(IJK2))
644               WSCM = AVG_Y(W_s(IJK2,M), W_s(IJK1,M), J_OF(IJK2))
645               VELS = WSCM
646     
647             ELSEIF(FCELL .EQ. 'E')THEN
648               IMJKP2 = IM_OF(IJKP2)
649     
650               DO MM = 1, SMAX
651                  AVGX1 = AVG_X(Theta_m(IJK1,MM),Theta_m(IJK2,MM),I_OF(IJK1))
652                  AVGX2 = AVG_X(Theta_m(IMJKP2,MM),Theta_m(IJKP2,MM),I_OF(IMJKP2))
653                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
654                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
655                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
656                    TH_avg(MM) = smallTheta
657                  ELSE
658                    TH_avg(MM) = AVG_Z(AVGX1, AVGX2, K_OF(IJK1))
659                  ENDIF
660               ENDDO
661     
662     ! Calculate velocity components at i+1/2,j,k+1/2 relative to IJK1
663               UGC  = AVG_Z(U_g(IJK1), U_g(IMJKP2), K_OF(IJK1))
664               VGC1 = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK1)),V_g(IJK1)),&
665                            AVG_Y_N(V_g(JM_OF(IMJKP2)),V_g(IMJKP2)),&
666                            K_OF(IJK1))
667               VGC2 = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK2)),V_g(IJK2)),&
668                            AVG_Y_N(V_g(JM_OF(IJKP2)),V_g(IJKP2)),&
669                            K_OF(IJK2))
670               VGC  = AVG_X(VGC1,VGC2,I_OF(IJK1))
671               WGC  = AVG_X(W_g(IJK1), W_g(IJK2),I_OF(IJK1))
672               USCM = AVG_Z(U_s(IJK1,M), U_s(IMJKP2,M), K_OF(IJK1))
673               VSCM1= AVG_Z(AVG_Y_N(V_s(JM_OF(IJK1),M),V_s(IJK1,M)),&
674                            AVG_Y_N(V_s(JM_OF(IMJKP2),M),V_s(IMJKP2,M)),&
675                            K_OF(IJK1))
676               VSCM2= AVG_Z(AVG_Y_N(V_s(JM_OF(IJK2),M),V_s(IJK2,M)),&
677                            AVG_Y_N(V_s(JM_OF(IJKP2),M),V_s(IJKP2,M)),&
678                            K_OF(IJK2))
679               VSCM  = AVG_X(VSCM1,VSCM2,I_OF(IJK1))
680               WSCM = AVG_X(W_s(IJK1,M), W_s(IJK2,M), I_OF(IJK1))
681               VELS = WSCM
682     
683             ELSEIF(FCELL .EQ. 'W')THEN
684               IPJKP2= IP_OF(IJKP2)
685     
686               DO MM = 1, SMAX
687                  AVGX1 = AVG_X(Theta_m(IJK2,MM),Theta_m(IJK1,MM),I_OF(IJK2))
688                  AVGX2 = AVG_X(Theta_m(IJKP2,MM),Theta_m(IPJKP2,MM),I_OF(IJKP2))
689                  IF(AVGX1 < ZERO .AND. AVGX2 > ZERO) AVGX1 = AVGX2
690                  IF(AVGX2 < ZERO .AND. AVGX1 > ZERO) AVGX2 = AVGX1
691                  IF(AVGX1 < ZERO .AND. AVGX2 < ZERO) THEN
692                    TH_avg(MM) = smallTheta
693                  ELSE
694                    TH_avg(MM) = AVG_Z(AVGX1, AVGX2, K_OF(IJK2))
695                  ENDIF
696               ENDDO
697     
698     ! Calculate velocity components at i+1/2,j,k+1/2 relative to IJK2
699               UGC  = AVG_Z(U_g(IJK2), U_g(IJKP2), K_OF(IJK2))
700               VGC1 = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK1)),V_g(IJK1)),&
701                            AVG_Y_N(V_g(JM_OF(IPJKP2)),V_g(IPJKP2)),&
702                            K_OF(IJK1))
703               VGC2 = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK2)),V_g(IJK2)),&
704                            AVG_Y_N(V_g(JM_OF(IJKP2)),V_g(IJKP2)),&
705                            K_OF(IJK2))
706               VGC  = AVG_X(VGC2,VGC1,I_OF(IJK2))
707               WGC  = AVG_X(W_g(IJK2), W_g(IJK1),I_OF(IJK2))
708               USCM = AVG_Z(U_s(IJK2,M), U_s(IJKP2,M), K_OF(IJK2))
709               VSCM1= AVG_Z(AVG_Y_N(V_s(JM_OF(IJK1),M),V_s(IJK1,M)),&
710                            AVG_Y_N(V_s(JM_OF(IPJKP2),M),V_s(IPJKP2,M)),&
711                            K_OF(IJK1))
712               VSCM2= AVG_Z(AVG_Y_N(V_s(JM_OF(IJK2),M),V_s(IJK2,M)),&
713                            AVG_Y_N(V_s(JM_OF(IJKP2),M),V_s(IJKP2,M)),&
714                            K_OF(IJK2))
715               VSCM  = AVG_X(VSCM2,VSCM1,I_OF(IJK2))
716               WSCM = AVG_X(W_s(IJK2,M), W_s(IJK1,M), I_OF(IJK2))
717               VELS = WSCM
718     
719             ELSE
720               WRITE(LINE,'(A, A)') 'Error: Unknown FCELL'
721               CALL WRITE_ERROR('CALC_GRBDRY', LINE, 1)
722               CALL exitMPI(myPE)
723             ENDIF
724     
725     
726           ELSE
727              WRITE(LINE,'(A, A)') 'Error: Unknown COM'
728              CALL WRITE_ERROR('CALC_GRBDRY', LINE, 1)
729              CALL exitMPI(myPE)
730           ENDIF
731     
732     ! magnitude of gas-solids relative velocity
733           VREL =&
734              DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + (WGC - WSCM)**2 )
735     
736     ! slip velocity for use in Jenkins bc (sof)
737           VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
738              + (WSCM-BC_WW_S(L,M))**2 )
739     
740           IF (FRICTION .AND. (ONE-EP_G(IJK2))>EPS_F_MIN) THEN
741             CALL CALC_S_DDOT_S(IJK1, IJK2, FCELL, COM, M, DEL_DOT_U,&
742                S_DDOT_S, S_dd)
743     
744             CALL CALC_Gw_Hw_Cw(g0(M), EPs_avg(M), EPg_avg, ep_star_avg, &
745                g0EPs_avg, TH_avg(M), Mu_g_avg, RO_g_avg, ROs_avg(M), &
746                DP_avg(M), K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP,&
747                DEL_DOT_U, S_DDOT_S, S_dd, VELS, WVELS, M, gw, hw, cw)
748           ELSE
749              GW = 1D0
750              Hw = F_Hw(g0, EPs_avg, EPg_avg, ep_star_avg, &
751                 g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROs_avg, &
752                 DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP, M)
753              CW = HW*WVELS
754           ENDIF
755     
756     
757           RETURN
758           END SUBROUTINE CALC_GRBDRY
759     
760     
761     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
762     !                                                                      C
763     !  Function: F_HW                                                      C
764     !                                                                      C
765     !  Purpose: Function for hw                                            C
766     !                                                                      C
767     !  Author: K. Agrawal & A. Srivastava, Princeton Univ. Date: 24-JAN-98 C
768     !  Reviewer:                                           Date:           C
769     !                                                                      C
770     !  Modified: Sofiane Benyahia, Fluent Inc.             Date: 02-FEB-05 C
771     !  Purpose: Include conductivity defined by Simonin and Ahmadi         C
772     !           Also included Jenkins small frictional limit               C
773     !                                                                      C
774     !  Literature/Document References:                                     C
775     !     See calc_mu_s.f for references on kinetic theory models          C
776     !     Johnson, P. C., and Jackson, R., Frictional-collisional          C
777     !        constitutive relations for granular materials, with           C
778     !        application to plane shearing, Journal of Fluid Mechanics,    C
779     !        1987, 176, pp. 67-93                                          C
780     !     Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating     C
781     !        energy in a collisional grain flow at a flat frictional wall, C
782     !        Physics of Fluids, 1997, 9(10), pp. 2835-2840                 C
783     !                                                                      C
784     !                                                                      C
785     !  Additional Notes:                                                   C
786     !     The JJ granular momentum BC is written as:                       C
787     !     usw/|usw|.tau.n+C+F=0                                            C
788     !       tau = the stress tensor (collisional+frictional)               C
789     !       n = outward normal vector                                      C
790     !       usw = slip velocity with wall                                  C
791     !       C = collisional force                                          C
792     !       F = frictional force                                           C
793     !                                                                      C
794     !    The granular momentum BC is written as the normal vector dot the  C
795     !    stress tensor. Besides the gradient in velocity of phase M, the   C
796     !    stress tensor expression may contain several additional terms     C
797     !    that would need to be accounted for when satisfying the BC. These C
798     !    modifications have NOT been rigorously addressed.                 C
799     !                                                                      C
800     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
801           DOUBLE PRECISION FUNCTION F_HW(g0,EPS,EPG, ep_star_avg, &
802                                          g0EPs_avg,TH,Mu_g_avg,RO_g_avg,ROs_avg, &
803                                          DP_avg,K_12_avg, Tau_12_avg, Tau_1_avg, &
804                                          VREL, VSLIP, M)
805     
806     !-----------------------------------------------
807     ! Modules
808     !-----------------------------------------------
809           USE param
810           USE param1
811           USE constant
812           USE physprop
813           USE run
814           USE fldvar
815           USE mpi_utility
816           IMPLICIT NONE
817     !-----------------------------------------------
818     ! Dummy Arguments
819     !-----------------------------------------------
820     ! Radial distribution function of solids phase M with each
821     ! other solids phase
822           DOUBLE PRECISION, INTENT(IN) :: g0(DIMENSION_M)
823     ! Average solids volume fraction of each solids phase
824           DOUBLE PRECISION, INTENT(IN) :: EPS(DIMENSION_M)
825     ! Average solids and gas volume fraction
826           DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
827     ! Sum of eps*G_0
828           DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
829     ! Average theta_m
830           DOUBLE PRECISION, INTENT(INOUT) :: TH (DIMENSION_M)
831     ! Average gas viscosity
832           DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
833     ! Average gas density
834           DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
835     ! Average solids density
836           DOUBLE PRECISION, INTENT(IN) :: ROS_avg(DIMENSION_M)
837     ! Average particle diameter of each solids phase
838           DOUBLE PRECISION, INTENT(IN) :: DP_avg(DIMENSION_M)
839     ! Average cross-correlation K_12 and lagrangian integral time-scale
840           DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
841     ! Magnitude of slip velocity between two phases
842           DOUBLE PRECISION, INTENT(IN) :: VREL
843     ! Slip velocity between wall and particles
844           DOUBLE PRECISION, INTENT(IN) :: VSLIP
845     ! Solids phase index
846           INTEGER, INTENT(IN) :: M
847     !-----------------------------------------------
848     ! Local Variables
849     !-----------------------------------------------
850     ! Solids phase index
851           INTEGER :: LL
852     ! Coefficient of 2nd term
853           DOUBLE PRECISION :: F_2
854     ! Coefficient of 1st term
855           DOUBLE PRECISION :: Mu_s
856     ! Viscosity
857           DOUBLE PRECISION :: Mu
858     ! Bulk viscosity
859           DOUBLE PRECISION :: Mu_b
860     ! Viscosity corrected for interstitial fluid effects
861           DOUBLE PRECISION :: Mu_star
862     ! Solids pressure
863           DOUBLE PRECISION :: Ps
864     ! Reynolds number based on slip velocity
865           DOUBLE PRECISION :: Re_g
866     ! Friction Factor in drag coefficient
867           DOUBLE PRECISION :: C_d
868     ! Drag Coefficient
869           DOUBLE PRECISION :: Beta, DgA
870     ! Constants in Simonin or Ahmadi model
871           DOUBLE PRECISION :: Sigma_c, Tau_2_c, Tau_12_st, Nu_t
872           DOUBLE PRECISION :: Tau_2, zeta_c_2, MU_2_T_Kin, Mu_2_Col
873           DOUBLE PRECISION :: Tmp_Ahmadi_Const
874     ! Variables for Iddir & Arastoopour model
875           DOUBLE PRECISION :: MU_sM_sum, MU_s_MM, MU_s_LM, MU_sM_ip, MU_common_term,&
876                               MU_sM_LM
877           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, D_PL, DPSUMo2
878           DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, Bp_lm
879     ! Variables for Garzo and Dufty model
880           DOUBLE PRECISION :: c_star, zeta0_star, nu_eta_star, press_star, &
881                               gamma_star, eta_k_star, eta_star, eta0
882     ! Variables for GTSH theory
883           DOUBLE PRECISION :: Chi, RdissP, Re_T, Rdiss, GamaAvg, A2_gtshW, &
884                               zeta_star, nu0, nuN, NuK, EDT_s, zeta_avg, etaK, &
885                               Mu_b_avg, mu2_0, mu4_0, mu4_1
886     !-----------------------------------------------
887     ! Functions
888     !-----------------------------------------------
889     ! Variable specularity coefficient
890           DOUBLE PRECISION :: PHIP_JJ
891           DOUBLE PRECISION :: S_star
892           DOUBLE PRECISION :: K_phi
893     !-----------------------------------------------
894     
895     ! This is done here similar to bc_theta to avoid small negative values of
896     ! Theta coming most probably from linear solver
897           IF(TH(M) .LE. ZERO)THEN
898             TH(M) = 1D-8
899             IF (myPE.eq.PE_IO) THEN
900                WRITE(*,*) &
901                  'Warning: Negative granular temp at wall set to 1e-8'
902     !          CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
903             ENDIF
904           ENDIF
905     
906     
907     ! common variables
908           D_PM = DP_avg(M)
909           M_PM = (PI/6.d0)*(D_PM**3)*ROS_avg(M)
910           NU_PM = (EPS(M)*ROS_avg(M))/M_PM
911     
912     ! This is from Wen-Yu correlation, you can put here your own single particle drag
913           Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
914           IF (Re_g .lt. 1000.d0) THEN
915              C_d = (24.d0/(Re_g+SMALL_NUMBER))*(ONE + 0.15d0 * Re_g**0.687d0)
916           ELSE
917              C_d = 0.44d0
918           ENDIF
919           DgA = 0.75d0*C_d*Ro_g_avg*EPG*VREL/(D_PM*EPG**(2.65d0))
920           IF(VREL == ZERO) DgA = LARGE_NUMBER
921           Beta = EPS(M)*DgA !this is equivalent to F_gs(ijk,m)
922     
923     
924     ! Defining viscosity according to each KT for later use in JJ BC
925     ! Also defining solids pressure according to each KT for use if JENKINS
926     ! -------------------------------------------------------------------
927           SELECT CASE (KT_TYPE_ENUM)
928           CASE (LUN_1984)
929              Mu = (5d0*DSQRT(Pi*TH(M))*DP_avg(M)*ROS_avg(M))/96d0
930              Mu_b = (256d0*Mu*EPS(M)*g0EPs_avg)/(5d0*Pi)
931     
932              IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
933                 Mu_star = Mu
934              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
935                 MU_star = ZERO
936              ELSE
937                 Mu_star = ROS_avg(M)*EPS(M)* g0(M)*TH(M)* Mu/ &
938                    (ROS_avg(M)*g0EPs_avg*TH(M) + 2.0d0*DgA/ROS_avg(M)* Mu)
939              ENDIF
940     
941              Mu_s = ((2d0+ALPHA)/3d0)*((Mu_star/(Eta*(2d0-Eta)*&
942                 g0(M)))*(ONE+1.6d0*Eta*g0EPs_avg)*&
943                 (ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
944                 g0EPs_avg)+(0.6d0*Mu_b*Eta))
945     
946     ! defining granular pressure (for Jenkins BC)
947              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
948     
949     
950           CASE (SIMONIN_1996)
951     ! particle relaxation time
952              Tau_12_st = ROS_avg(M)/(DgA+small_number)
953     
954              Sigma_c = (ONE+ C_e)*(3.d0-C_e)/5.d0
955              Tau_2_c = DP_avg(M)/(6.d0*EPS(M)*g0(M)*DSQRT(16.d0*(TH(M)+Small_number)/PI))
956              zeta_c_2= 2.D0/5.D0*(ONE+ C_e)*(3.d0*C_e-ONE)
957              Nu_t =  Tau_12_avg/Tau_12_st
958              Tau_2 = ONE/(2.D0/Tau_12_st+Sigma_c/Tau_2_c)
959              MU_2_T_Kin = (2.0D0/3.0D0*K_12_avg*Nu_t + TH(M) * &
960                   (ONE+ zeta_c_2*EPS(M)*g0(M)))*Tau_2
961              Mu_2_Col = 8.D0/5.D0*EPS(M)*g0(M)*Eta* (MU_2_T_Kin+ &
962                    DP_avg(M)*DSQRT(TH(M)/PI))
963              Mu_s = EPS(M)*ROS_avg(M)*(MU_2_T_Kin + Mu_2_Col)
964     
965     ! defining granular pressure (for Jenkins BC)
966              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
967     
968     
969           CASE (AHMADI_1995)
970     ! particle relaxation time
971              Tau_12_st = ROS_avg(M)/(DgA+small_number)
972     
973              IF(EPS(M) < (ONE-ep_star_avg)) THEN
974                 Tmp_Ahmadi_Const = ONE/&
975                    (ONE+ Tau_1_avg/Tau_12_st * (ONE-EPS(M)/(ONE-ep_star_avg))**3)
976              ELSE
977                 Tmp_Ahmadi_Const = ONE
978              ENDIF
979              Mu_s = Tmp_Ahmadi_Const &
980                 *0.1045D0*(ONE/g0(M)+3.2D0*EPS(M)+12.1824D0*g0(M)*EPS(M)*EPS(M))  &
981                 *DP_avg(M)*ROS_avg(M)* DSQRT(TH(M))
982     
983     ! defining granular pressure (for Jenkins BC)
984              Ps = ROs_avg(M)*EPS(M)*TH(M)*&
985                  ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
986     
987     
988           CASE(IA_2005)
989     ! Use original IA theory if SWITCH_IA is false
990              IF(.NOT. SWITCH_IA) g0EPs_avg = EPS(M)*ROS_avg(M)
991     
992              Mu = (5.d0/96.d0)*D_PM* ROS_avg(M)*DSQRT(PI*TH(M)/M_PM)
993     
994              IF(SWITCH == ZERO .OR. RO_g_avg == ZERO)THEN
995                 Mu_star = Mu
996              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
997                 MU_star = ZERO
998              ELSE
999                 Mu_star = Mu*EPS(M)*g0(M)/ (g0EPs_avg+ 2.0d0*DgA*Mu / &
1000                    (ROS_avg(M)**2 *(TH(M)/M_PM)))
1001              ENDIF
1002     
1003              MU_s_MM = (Mu_star/g0(M))*(1.d0+(4.d0/5.d0)*(1.d0+C_E)*g0EPs_avg)**2
1004              Mu_sM_sum = ZERO
1005     
1006              DO LL = 1, SMAX
1007                 D_PL = DP_avg(LL)
1008                 M_PL = (PI/6.d0)*(D_PL**3.)*ROS_avg(LL)
1009                 MPSUM = M_PM + M_PL
1010                 DPSUMo2 = (D_PM+D_PL)/2.d0
1011                 NU_PL = (EPS(LL)*ROS_avg(LL))/M_PL
1012     
1013                 IF ( LL .eq. M) THEN
1014                    Ap_lm = MPSUM/(2.d0)
1015                    Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
1016                    R1p_lm = ONE/( Ap_lm**1.5 * Dp_lm**3 )
1017                    MU_s_LM = DSQRT(PI)*(DPSUMo2**4 / (48.d0*5.d0))*&
1018                       g0(LL)*(M_PL*M_PM/MPSUM)*(M_PL*M_PM/&
1019                       MPSUM)*((M_PL*M_PM)**1.5)*NU_PM*NU_PL*&
1020                       (1.d0+C_E)*R1p_lm*DSQRT(TH(M))
1021     
1022     ! solids phase 'viscosity' associated with the divergence
1023     ! of solids phase M
1024                    MU_sM_ip = (MU_s_MM + MU_s_LM)
1025     
1026                 ELSE
1027                    Ap_lm = (M_PM*TH(LL)+M_PL*TH(M))/&
1028                       (2.d0)
1029                    Bp_lm = (M_PM*M_PL*(TH(LL)-TH(M) ))/&
1030                       (2.d0*MPSUM)
1031                    Dp_lm = (M_PL*M_PM*(M_PM*TH(M)+M_PL*TH(LL) ))/&
1032                       (2.d0*MPSUM*MPSUM)
1033                    R1p_lm = ( ONE/( Ap_lm**1.5 * Dp_lm**3 ) ) + &
1034                       ( (9.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**4 ) )+&
1035                       ( (30.d0*Bp_lm**4)/( 2.d0*Ap_lm**3.5 * Dp_lm**5 ) )
1036                    MU_common_term = DSQRT(PI)*(DPSUMo2**4 / (48.d0*5.d0))*&
1037                       g0(LL)*(M_PL*M_PM/MPSUM)*(M_PL*M_PM/&
1038                       MPSUM)*((M_PL*M_PM)**1.5)*NU_PM*NU_PL*&
1039                       (1.d0+C_E)*R1p_lm
1040                    MU_sM_LM = MU_common_term*(TH(M)*TH(M)*TH(LL)*TH(LL)*TH(LL) )
1041     
1042     ! solids phase 'viscosity' associated with the divergence
1043     ! of solids phase M
1044                    MU_sM_ip = MU_sM_LM
1045     
1046                 ENDIF
1047                 MU_sM_sum = MU_sM_sum + MU_sM_ip
1048     
1049               ENDDO
1050     
1051     ! Find the term proportional to the gradient in velocity
1052     ! of phase M  (viscosity in the Mth solids phase)
1053               Mu_s = MU_sM_sum
1054     
1055     
1056           CASE(GD_1999)
1057     ! Pressure/Viscosity/Bulk Viscosity
1058     ! Note: k_boltz = M_PM
1059     
1060     ! Find pressure in the Mth solids phase
1061              press_star = 1.d0 + 2.d0*(1.d0+C_E)*EPS(M)*G0(M)
1062     
1063     ! find bulk and shear viscosity
1064              c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
1065                 / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
1066     
1067              zeta0_star = (5.d0/12.d0)*G0(M)*(1.d0 - C_E*C_E) &
1068                 * (1.d0 + (3.d0/32.d0)*c_star)
1069     
1070              nu_eta_star = G0(M)*(1.d0 - 0.25d0*(1.d0-C_E)*(1.d0-C_E)) &
1071                 * (1.d0-(c_star/64.d0))
1072     
1073              gamma_star = (4.d0/5.d0)*(32.d0/PI)*EPS(M)*EPS(M) &
1074                 * G0(M)*(1.d0+C_E) * (1.d0 - (c_star/32.d0))
1075     
1076              eta_k_star = (1.d0 - (2.d0/5.d0)*(1.d0+C_E)*(1.d0-3.d0*C_E) &
1077                 * EPS(M)*G0(M) ) / (nu_eta_star - 0.5d0*zeta0_star)
1078     
1079              eta_star = eta_k_star*(1.d0 + (4.d0/5.d0)*EPS(M)*G0(M) &
1080                 * (1.d0+C_E) ) + (3.d0/5.d0)*gamma_star
1081     
1082              eta0 = 5.0d0*M_PM*DSQRT(TH(M)/PI) / (16.d0*D_PM*D_PM)
1083     
1084     ! added Ro_g = 0 for granular flows (no gas).
1085              IF(SWITCH == ZERO .OR. RO_g_avg == ZERO) THEN
1086                 Mu_star = eta0
1087              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
1088                 Mu_star = ZERO
1089              ELSE
1090                 Mu_star = ROS_avg(M)*EPS(M)*G0(M)*TH(M)*eta0 / &
1091                    ( ROS_avg(M)*EPS(M)*G0(M)*TH(M) + &
1092                    (2.d0*DgA*eta0/ROS_avg(M)) )     ! Note dgA is ~F_gs/ep_s
1093              ENDIF
1094     
1095     !  shear viscosity in Mth solids phase  (add to frictional part)
1096              Mu_s = Mu_star * eta_star
1097     
1098     ! defining granular pressure (for Jenkins BC)
1099              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1100                         ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1101     
1102     
1103           CASE (GTSH_2012)
1104              Chi = g0(M)
1105     
1106              RdissP = one
1107              IF(EPS(M) > small_number) RdissP = &
1108                 one + 3d0*dsqrt(EPS(M)/2d0) + 135d0/64d0*EPS(M)*dlog(EPS(M)) + &
1109                 11.26d0*EPS(M)*(one-5.1d0*EPS(M)+16.57d0*EPS(M)**2-21.77d0*    &
1110                 EPS(M)**3) - EPS(M)*Chi*dlog(epM)
1111     
1112              Re_T = RO_g_avg*D_pm*dsqrt(TH(M)) / Mu_g_avg
1113              Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
1114              Rdiss = RdissP + Re_T * K_phi(EPS(M))
1115              GamaAvg = 3d0*Pi*Mu_g_avg*D_pm*Rdiss
1116              mu2_0 = dsqrt(2d0*pi) * Chi * (one-C_E**2)
1117              mu4_0 = (4.5d0+C_E**2) * mu2_0
1118              mu4_1 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + 2d0*dsqrt(2d0*pi)* &
1119                 Chi*(one+C_E)
1120              A2_gtshW = zero ! for eps = zero
1121     
1122              if((EPS(M)> small_number) .and. (TH(M) > small_number)) then
1123                 zeta_star = 4.5d0*dsqrt(2d0*Pi)*(RO_g_avg/ROs_avg(M))**2*Re_g**2 * &
1124                             S_star(EPS(M),Chi) / (EPS(M)*(one-EPS(M))**2 * Re_T**4)
1125                 A2_gtshW = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
1126                                (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
1127              endif
1128     
1129              eta0 = 0.3125d0/(dsqrt(pi)*D_PM**2)*M_pm*dsqrt(TH(M))
1130              nu0 = (96.d0/5.d0)*(EPS(M)/D_PM)*DSQRT(TH(M)/PI)
1131              nuN = 0.25d0*nu0*Chi*(3d0-C_E)*(one+C_E) * &
1132                 (one+0.7375d0*A2_gtshW)
1133              NuK = nu0*(one+C_E)/3d0*Chi*( one+2.0625d0*(one-C_E)+ &
1134                 ((947d0-579*C_E)/256d0*A2_gtshW) )
1135              EDT_s = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
1136                 (one+0.1875d0*A2_gtshW)*NU_PM*D_PM**2*dsqrt(TH(M))
1137     
1138              if((TH(m) > small_number) .and. (EPS(M) > small_number)) then
1139                 zeta_avg = one/6d0*D_PM*VREL**2*(3d0*pi*Mu_g_avg*D_PM/&
1140                    M_pm)**2 / dsqrt(TH(m)) * S_star(EPS(M), Chi)
1141                 etaK = ROs_avg(M)*EPS(M)*TH(m) / (nuN-0.5d0*( &
1142                    EDT_s-zeta_avg/TH(m) - 2d0*GamaAvg/M_PM)) * &
1143                    (one -0.4d0 * (one+C_E)*(one-3d0*C_E)*EPS(M)*Chi)
1144              else
1145                 etaK = zero
1146              endif
1147     
1148              Mu_b_avg = 25.6d0/pi * EPS(M)**2 * Chi *(one+C_E) * &
1149                 (one - A2_gtshW/16d0)*eta0
1150     
1151              Mu_s = etaK*(one+0.8d0*EPS(M)*Chi*(one+C_E)) + &
1152                 0.6d0*Mu_b_avg
1153     
1154     ! defining granular pressure (for Jenkins BC)
1155              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1156     
1157     
1158           CASE DEFAULT
1159     ! should never hit this
1160              WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1161              WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1162              call mfix_exit(myPE)
1163           END SELECT
1164     
1165     
1166     ! setting the coefficients for JJ BC
1167     ! -------------------------------------------------------------------
1168           SELECT CASE (KT_TYPE_ENUM)
1169     
1170              CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999, GTSH_2012)
1171     ! KT without particle mass in their definition of granular temperature
1172     ! and theta = granular temperature  OR
1173     ! KT with particle mass in their definition of granular temperature
1174     ! and theta = granular temperature/mass
1175     
1176     ! Jenkins BC is implemented for these KT
1177                 IF(JENKINS) THEN
1178                    IF (VSLIP == ZERO) THEN
1179     ! if solids velocity field is initialized to zero, use free slip bc
1180                       F_2 = zero
1181                    ELSE
1182     ! As I understand from soil mechanic papers, the coefficient mu in
1183     ! Jenkins paper is tan_Phi_w. T.  See for example, G.I. Tardos, PT,
1184     ! 92 (1997), 61-74, equation (1). sof
1185     
1186     ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
1187                       F_2 = tan_Phi_w*Ps/VSLIP
1188     
1189                    ENDIF
1190                 ELSE
1191                    IF(.NOT. BC_JJ_M) THEN
1192                       F_2 = (PHIP*DSQRT(3d0*TH(M))*Pi*ROS_avg(M)*EPS(M)*&
1193                          g0(M))/(6d0*(ONE-ep_star_avg))
1194                    ELSE
1195                       F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3d0*TH(M))*Pi*&
1196                          ROS_avg(M)*EPS(M)*g0(M))/(6d0*(ONE-ep_star_avg))
1197                    ENDIF
1198                 ENDIF   ! end if(Jenkins)/else
1199     
1200     
1201             CASE (IA_2005)
1202     ! KTs with particle mass in their definition of granular temperature
1203     ! and theta = granular temperature
1204                 IF(.NOT. BC_JJ_M) THEN
1205                    F_2 = (PHIP*DSQRT(3.d0*TH(M)/M_PM)*PI*ROS_avg(M)*EPS(M)*&
1206                       g0(M))/(6.d0*(ONE-ep_star_avg))
1207                 ELSE
1208                    F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3.d0*TH(M)/M_PM)*PI*&
1209                       ROS_avg(M)*EPS(M)*g0(M))/(6.d0*(ONE-ep_star_avg))
1210                 ENDIF
1211     
1212     
1213              CASE DEFAULT
1214     ! should never hit this
1215                 WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1216                 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1217                 call mfix_exit(myPE)
1218           END SELECT
1219     
1220           F_HW =  F_2/Mu_s
1221     
1222           RETURN
1223           END FUNCTION F_HW
1224     
1225     
1226     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1227     !                                                                      C
1228     !  Subroutine CG_CALC_GRBDRY                                           C
1229     !  Purpose: Cut cell version                                           C
1230     !                                                                      C
1231     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1232     
1233           SUBROUTINE CG_CALC_GRBDRY(IJK, TYPE_OF_CELL, M, L, F_2)
1234     
1235     !-----------------------------------------------
1236     ! Modules
1237     !-----------------------------------------------
1238           USE bc
1239           USE compar
1240           USE constant
1241           USE fldvar
1242           USE fun_avg
1243           USE functions
1244           USE geometry
1245           USE indices
1246           USE param
1247           USE param1
1248           USE physprop
1249           USE rdf
1250           USE run
1251           USE turb
1252           USE visc_s
1253           Use cutcell
1254           use toleranc
1255           IMPLICIT NONE
1256     !-----------------------------------------------
1257     ! Dummy Arguments
1258     !-----------------------------------------------
1259     ! IJK indices
1260           INTEGER, INTENT(IN) :: IJK
1261     
1262           CHARACTER (LEN=*) :: TYPE_OF_CELL
1263     ! Solids phase index
1264           INTEGER, INTENT(IN) :: M
1265     ! Index corresponding to boundary condition
1266           INTEGER, INTENT(IN) ::  L
1267     
1268           DOUBLE PRECISION :: F_2
1269     
1270     !-----------------------------------------------
1271     ! Local Variables
1272     !-----------------------------------------------
1273     ! IJK indices
1274           INTEGER          IJK1, IJK2
1275     ! Other indices
1276           INTEGER          IJK2E, IPJK2, IPJMK2
1277     !
1278           DOUBLE PRECISION :: smallTheta
1279     ! Average scalars
1280           DOUBLE PRECISION EPg_avg, Mu_g_avg, RO_g_avg
1281     ! Average void fraction at packing
1282           DOUBLE PRECISION ep_star_avg
1283     ! Average scalars modified to include all solid phases
1284           DOUBLE PRECISION EPs_avg(DIMENSION_M), DP_avg(DIMENSION_M),&
1285                            TH_avg(DIMENSION_M)
1286     ! Average solids density
1287           DOUBLE PRECISION ROS_AVG(DIMENSION_M)
1288     ! Average Simonin and Ahmadi variables (sof)
1289           DOUBLE PRECISION K_12_avg, Tau_12_avg, Tau_1_avg
1290     ! Average velocities
1291           DOUBLE PRECISION WGC1, WGC2
1292     ! Solids phase index
1293           INTEGER          MM
1294     ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
1295           DOUBLE PRECISION USCM, VSCM,WSCM
1296           DOUBLE PRECISION WSCM1,WSCM2
1297     ! values of U_g, V_g, W_g at appropriate place on boundary wall
1298           DOUBLE PRECISION UGC, VGC, WGC
1299     ! Magnitude of gas-solids relative velocity
1300           DOUBLE PRECISION VREL
1301     ! slip velocity between wall and particles for Jenkins bc (sof)
1302           DOUBLE PRECISION VSLIP
1303     ! radial distribution function at contact
1304           DOUBLE PRECISION g0(DIMENSION_M)
1305     ! Sum of eps*G_0
1306           DOUBLE PRECISION g0EPs_avg
1307     !-----------------------------------------------
1308     
1309     !  Note:  EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).  Hence
1310     !         IJK2 (fluid cell) is used in averages.
1311     !
1312     
1313           smallTheta = (to_SI)**4 * ZERO_EP_S
1314     
1315           SELECT CASE (TYPE_OF_CELL)
1316     
1317              CASE('U_MOMENTUM')
1318     
1319                 IJK2 = EAST_OF(IJK)
1320                 IJK2E = EAST_OF(IJK2)
1321     
1322                 EPg_avg = (VOL(IJK)*EP_g(IJK) + VOL(IJK2)*EP_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1323                 ep_star_avg = (VOL(IJK)*EP_star_array(IJK) + VOL(IJK2)*EP_star_array(IJK2))/(VOL(IJK) + VOL(IJK2))
1324                 Mu_g_avg = (VOL(IJK)*Mu_g(IJK) + VOL(IJK2)*Mu_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1325                 RO_g_avg = (VOL(IJK)*RO_g(IJK) + VOL(IJK2)*RO_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1326     
1327                 K_12_avg = ZERO
1328                 Tau_12_avg = ZERO
1329                 Tau_1_avg = ZERO
1330                 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR.&
1331                    KT_TYPE_ENUM == AHMADI_1995) THEN  ! not converted to CG
1332                    K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
1333                    Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
1334                    Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
1335                 ENDIF
1336     
1337                 g0EPs_avg = ZERO
1338                 DO MM = 1, SMAX
1339                    g0(MM)      = G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM)
1340                    EPs_avg(MM) = (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1341                    DP_avg(MM)  = (VOL(IJK)*D_P(IJK, MM) + VOL(IJK2)*D_P(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1342                    g0EPs_avg   = g0EPs_avg + G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM) &
1343                                * (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1344                    ROS_AVG(MM) = (VOL(IJK)*RO_S(IJK, MM) + VOL(IJK2)*RO_S(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1345     
1346                    TH_avg(MM)  = (VOL(IJK)*Theta_m(IJK,MM) + VOL(IJK2)*Theta_m(IJK2,MM))/(VOL(IJK) + VOL(IJK2))
1347                    IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
1348                 ENDDO
1349     
1350     
1351     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)  ! not converted to CG
1352                 UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
1353                 VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
1354                 WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
1355                              AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
1356                              I_OF(IJK2))
1357                 WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
1358                              AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
1359                              I_OF(IJK1))
1360                 WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
1361                 USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
1362                 VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
1363                 WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
1364                              AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
1365                              I_OF(IJK2))
1366                 WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
1367                              AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
1368                              I_OF(IJK1))
1369                 WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
1370     
1371     ! magnitude of gas-solids relative velocity
1372                 VREL = DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + &
1373                               (WGC - WSCM)**2 )
1374     
1375     ! slip velocity for use in Jenkins bc (sof)
1376                 VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
1377                            + (WSCM-BC_WW_S(L,M))**2 )
1378     
1379                 CALL GET_CG_F2(g0, EPs_avg, EPg_avg, ep_star_avg, &
1380                           g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROS_AVG,&
1381                            DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, &
1382                           VREL, VSLIP, M,F_2)
1383     
1384     
1385              CASE('V_MOMENTUM')
1386     
1387                 IJK2 = NORTH_OF(IJK)
1388     
1389                 EPg_avg = (VOL(IJK)*EP_g(IJK) + VOL(IJK2)*EP_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1390                 ep_star_avg = (VOL(IJK)*EP_star_array(IJK) + VOL(IJK2)*EP_star_array(IJK2))/(VOL(IJK) + VOL(IJK2))
1391                 Mu_g_avg = (VOL(IJK)*Mu_g(IJK) + VOL(IJK2)*Mu_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1392                 RO_g_avg = (VOL(IJK)*RO_g(IJK) + VOL(IJK2)*RO_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1393     
1394                 K_12_avg = ZERO
1395                 Tau_12_avg = ZERO
1396                 Tau_1_avg = ZERO
1397                 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR.&
1398                    KT_TYPE_ENUM == AHMADI_1995) THEN  ! not converted to CG
1399                    Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
1400     ! Simonin only:
1401                    K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
1402                    Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
1403                 ENDIF
1404     
1405                 g0EPs_avg = ZERO
1406                 DO MM = 1, SMAX
1407                    g0(MM)      = G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM)
1408                    EPs_avg(MM) = (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1409                    DP_avg(MM)  = (VOL(IJK)*D_P(IJK, MM) + VOL(IJK2)*D_P(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1410                    g0EPs_avg   = g0EPs_avg + G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM) &
1411                                * (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1412                    ROS_avg(MM) = (VOL(IJK)*RO_S(IJK, MM) + VOL(IJK2)*RO_S(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1413                    TH_avg(MM)  = (VOL(IJK)*Theta_m(IJK,MM) + VOL(IJK2)*Theta_m(IJK2,MM))/(VOL(IJK) + VOL(IJK2))
1414                    IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
1415                 ENDDO
1416     
1417     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)    ! not converted to CG
1418                 UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
1419                 VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
1420                 WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
1421                              AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
1422                              I_OF(IJK2))
1423                 WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
1424                              AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
1425                              I_OF(IJK1))
1426                 WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
1427                 USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
1428                 VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
1429                 WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
1430                              AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
1431                              I_OF(IJK2))
1432                 WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
1433                              AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
1434                              I_OF(IJK1))
1435                 WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
1436     
1437     ! magnitude of gas-solids relative velocity
1438     
1439                 VREL = DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + &
1440                               (WGC - WSCM)**2 )
1441     
1442     ! slip velocity for use in Jenkins bc (sof)
1443                 VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
1444                             + (WSCM-BC_WW_S(L,M))**2 )
1445     
1446     
1447                 CALL GET_CG_F2(g0, EPs_avg, EPg_avg, ep_star_avg, &
1448                           g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROS_AVG,&
1449                           DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, &
1450                           VREL, VSLIP, M,F_2)
1451     
1452     
1453              CASE('W_MOMENTUM')
1454     
1455                 IJK2 = TOP_OF(IJK)
1456     
1457                 EPg_avg = (VOL(IJK)*EP_g(IJK) + VOL(IJK2)*EP_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1458                 ep_star_avg = (VOL(IJK)*EP_star_array(IJK) + VOL(IJK2)*EP_star_array(IJK2))/(VOL(IJK) + VOL(IJK2))
1459                 Mu_g_avg = (VOL(IJK)*Mu_g(IJK) + VOL(IJK2)*Mu_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1460                 RO_g_avg = (VOL(IJK)*RO_g(IJK) + VOL(IJK2)*RO_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1461     
1462                 K_12_avg = ZERO
1463                 Tau_12_avg = ZERO
1464                 Tau_1_avg = ZERO
1465                 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR.&
1466                    KT_TYPE_ENUM == AHMADI_1995) THEN  ! not converted to CG
1467                    Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
1468     ! Simonon only:
1469                    K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
1470                    Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
1471                 ENDIF
1472     
1473                 g0EPs_avg = ZERO
1474                 DO MM = 1, SMAX
1475                    g0(MM)      = G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM)
1476                    EPs_avg(MM) = (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1477                    DP_avg(MM)  = (VOL(IJK)*D_P(IJK, MM) + VOL(IJK2)*D_P(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1478                    ROS_avg(MM) = (VOL(IJK)*RO_S(IJK, MM) + VOL(IJK2)*RO_S(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1479                    g0EPs_avg   = g0EPs_avg + G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM) &
1480                                * (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1481     
1482                    TH_avg(MM)  = (VOL(IJK)*Theta_m(IJK,MM) + VOL(IJK2)*Theta_m(IJK2,MM))/(VOL(IJK) + VOL(IJK2))
1483                    IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
1484                 ENDDO
1485     
1486     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)    ! not converted to CG
1487                 UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
1488                 VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
1489                 WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
1490                              AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
1491                              I_OF(IJK2))
1492                 WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
1493                              AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
1494                              I_OF(IJK1))
1495                 WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
1496                 USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
1497                 VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
1498                 WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
1499                              AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
1500                              I_OF(IJK2))
1501                 WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
1502                              AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
1503                              I_OF(IJK1))
1504                 WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
1505     
1506     ! magnitude of gas-solids relative velocity
1507                 VREL = DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + &
1508                               (WGC - WSCM)**2 )
1509     
1510     ! slip velocity for use in Jenkins bc (sof)
1511                 VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
1512                             + (WSCM-BC_WW_S(L,M))**2 )
1513     
1514     
1515     ! why bother with this since it should not come here...
1516                 CALL GET_CG_F2(g0, EPs_avg, EPg_avg, ep_star_avg, &
1517                           g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROS_AVG,&
1518                           DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, &
1519                           VREL, VSLIP, M,F_2)
1520     
1521     
1522              CASE DEFAULT
1523                 WRITE(*,*) 'CG_CALC_GRBDRY'
1524                 WRITE(*,*) 'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1525                 WRITE(*,*) 'ACCEPTABLE TYPES ARE:'
1526                 WRITE(*,*) 'U_MOMENTUM'
1527                 WRITE(*,*) 'V_MOMENTUM'
1528                 WRITE(*,*) 'W_MOMENTUM'
1529                 CALL MFIX_EXIT(myPE)
1530               END SELECT
1531     
1532     
1533           RETURN
1534     
1535           END SUBROUTINE CG_CALC_GRBDRY
1536     
1537     
1538     
1539     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1540     !                                                                      C
1541     !  Subroutine: GET_CG_F2                                               C
1542     !  Purpose: Compute F_2 for cut cell version                           C
1543     !                                                                      C
1544     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1545           SUBROUTINE GET_CG_F2(g0,EPS,EPG, ep_star_avg, &
1546                                          g0EPs_avg,TH,Mu_g_avg,RO_g_avg,Ros_avg, &
1547                                          DP_avg,K_12_avg, Tau_12_avg, Tau_1_avg, &
1548                                          VREL, VSLIP, M,F_2)
1549     
1550     !-----------------------------------------------
1551     ! Modules
1552     !-----------------------------------------------
1553           USE param
1554           USE param1
1555           USE constant
1556           USE physprop
1557           USE run
1558           USE fldvar
1559           USE mpi_utility
1560           Use cutcell
1561           IMPLICIT NONE
1562     !-----------------------------------------------
1563     ! Dummy Arguments
1564     !-----------------------------------------------
1565     ! Radial distribution function of solids phase M with each
1566     ! other solids phase
1567           DOUBLE PRECISION, INTENT(IN) :: g0(DIMENSION_M)
1568     ! Average solids volume fraction of each solids phase
1569           DOUBLE PRECISION, INTENT(IN) :: EPS(DIMENSION_M)
1570     ! Average solids and gas volume fraction
1571           DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
1572     ! Sum of eps*G_0
1573           DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
1574     ! Average theta_m
1575           DOUBLE PRECISION, INTENT(INOUT) :: TH (DIMENSION_M)
1576     ! Average gas viscosity
1577           DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
1578     ! Average gas density
1579           DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
1580     !QX: Average solids density
1581           DOUBLE PRECISION, INTENT(IN) :: ROS_avg(DIMENSION_M)
1582     ! Average particle diameter of each solids phase
1583           DOUBLE PRECISION, INTENT(IN) :: DP_avg(DIMENSION_M)
1584     ! Average cross-correlation K_12 and lagrangian integral time-scale
1585           DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
1586     ! Magnitude of slip velocity between two phases
1587           DOUBLE PRECISION, INTENT(IN) :: VREL
1588     ! Slip velocity between wall and particles
1589           DOUBLE PRECISION, INTENT(IN) :: VSLIP
1590     ! Solids phase index
1591           INTEGER, INTENT(IN) :: M
1592     !-----------------------------------------------
1593     ! Local Variables
1594     !-----------------------------------------------
1595     ! Solids pressure
1596           DOUBLE PRECISION :: Ps
1597           DOUBLE PRECISION :: F_2
1598           DOUBLE PRECISION :: D_PM, M_PM
1599     !-----------------------------------------------
1600     ! Functions
1601     !-----------------------------------------------
1602     ! Variable specularity coefficient
1603           DOUBLE PRECISION :: PHIP_JJ
1604     !-----------------------------------------------
1605     
1606     ! This is done here similar to bc_theta to avoid small negative values of
1607     ! Theta coming most probably from linear solver
1608           IF(TH(M) .LE. ZERO)THEN
1609             TH(M) = 1D-8
1610             IF (myPE.eq.PE_IO) THEN
1611                WRITE(*,*) &
1612                   'Warning: Negative granular temp at wall set to 1e-8'
1613     !          CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
1614             ENDIF
1615           ENDIF
1616     
1617     ! common variables
1618           D_PM = DP_avg(M)
1619           M_PM = (PI/6.d0)*(D_PM**3)*ROS_avg(M)
1620     
1621     ! Defining viscosity according to each KT for later use in JJ BC
1622     ! Also defining solids pressure according to each KT for use if JENKINS
1623     ! -------------------------------------------------------------------
1624           SELECT CASE (KT_TYPE_ENUM)
1625           CASE (LUN_1984)
1626     ! defining granular pressure (for Jenkins BC)
1627              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
1628     
1629     
1630           CASE (SIMONIN_1996)
1631     ! defining granular pressure (for Jenkins BC)
1632              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
1633     
1634     
1635           CASE (AHMADI_1995)
1636     ! defining granular pressure (for Jenkins BC)
1637              Ps = ROs_avg(M)*EPS(M)*TH(M)*&
1638                  ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
1639     
1640     
1641           CASE(IA_2005)
1642     ! Use original IA theory if SWITCH_IA is false
1643     ! no ps since no jenkins for this theory
1644     
1645           CASE(GD_1999)
1646     
1647     ! defining granular pressure (for Jenkins BC)
1648              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1649                         ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1650     
1651     
1652           CASE (GTSH_2012)
1653     ! defining granular pressure (for Jenkins BC)
1654              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1655     
1656     
1657           CASE DEFAULT
1658     ! should never hit this
1659              WRITE (*, '(A)') 'CG_CALC_GRBDRY => GET_CG_F2 '
1660              WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1661              call mfix_exit(myPE)
1662           END SELECT
1663     
1664     
1665     ! setting the coefficients for JJ BC
1666     ! -------------------------------------------------------------------
1667           SELECT CASE (KT_TYPE_ENUM)
1668     
1669              CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999, GTSH_2012)
1670     ! KT without particle mass in their definition of granular temperature
1671     ! and theta = granular temperature  OR
1672     ! KT with particle mass in their definition of granular temperature
1673     ! and theta = granular temperature/mass
1674     
1675     ! Jenkins BC is implemented for these KT
1676                 IF(JENKINS) THEN
1677                    IF (VSLIP == ZERO) THEN
1678     ! if solids velocity field is initialized to zero, use free slip bc
1679                       F_2 = zero
1680                    ELSE
1681     ! As I understand from soil mechanic papers, the coefficient mu in
1682     ! Jenkins paper is tan_Phi_w. T.  See for example, G.I. Tardos, PT,
1683     ! 92 (1997), 61-74, equation (1). sof
1684     
1685     ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
1686                       F_2 = tan_Phi_w*Ps/VSLIP
1687     
1688                    ENDIF
1689                 ELSE
1690                    IF(.NOT. BC_JJ_M) THEN
1691                       F_2 = (PHIP*DSQRT(3d0*TH(M))*Pi*ROS_avg(M)*EPS(M)*&
1692                          g0(M))/(6d0*(ONE-ep_star_avg))
1693                    ELSE
1694                       F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3d0*TH(M))*Pi*&
1695                          ROS_avg(M)*EPS(M)*g0(M))/(6d0*(ONE-ep_star_avg))
1696                    ENDIF
1697                 ENDIF   ! end if(Jenkins)/else
1698     
1699              CASE(IA_2005)
1700     ! KTs with particle mass in their definition of granular temperature
1701     ! and theta = granular temperature
1702                 IF(.NOT. BC_JJ_M) THEN
1703                    F_2 = (PHIP*DSQRT(3.d0*TH(M)/M_PM)*PI*ROS_avg(M)*&
1704                       EPS(M)*g0(M))/(6.d0*(ONE-ep_star_avg))
1705                 ELSE
1706                    F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3.d0*TH(M)/M_PM)*&
1707                       PI*ROS_avg(M)*EPS(M)*g0(M))/(6.d0*(ONE-ep_star_avg))
1708                 ENDIF
1709     
1710     
1711              CASE DEFAULT
1712     ! should never hit this
1713                 WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1714                 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1715                 call mfix_exit(myPE)
1716           END SELECT
1717     
1718     ! all the code used to calculate Mu_s was deleted from this routine
1719     !      F_HW =  F_2/Mu_s  ! Only F_2 is actually needed
1720     
1721     
1722           RETURN
1723           END SUBROUTINE GET_CG_F2
1724     
1725