File: RELATIVE:/../../../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     !  Purpose: Function for hw                                            C
765     !                                                                      C
766     !  Literature/Document References:                                     C
767     !     See calc_mu_s.f for references on kinetic theory models          C
768     !     Johnson, P. C., and Jackson, R., Frictional-collisional          C
769     !        constitutive relations for granular materials, with           C
770     !        application to plane shearing, Journal of Fluid Mechanics,    C
771     !        1987, 176, pp. 67-93                                          C
772     !     Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating     C
773     !        energy in a collisional grain flow at a flat frictional wall, C
774     !        Physics of Fluids, 1997, 9(10), pp. 2835-2840                 C
775     !                                                                      C
776     !                                                                      C
777     !  Additional Notes:                                                   C
778     !     The JJ granular momentum BC is written as:                       C
779     !     usw/|usw|.tau.n+C+F=0                                            C
780     !       tau = the stress tensor (collisional+frictional)               C
781     !       n = outward normal vector                                      C
782     !       usw = slip velocity with wall                                  C
783     !       C = collisional force                                          C
784     !       F = frictional force                                           C
785     !                                                                      C
786     !    The granular momentum BC is written as the normal vector dot the  C
787     !    stress tensor. Besides the gradient in velocity of phase M, the   C
788     !    stress tensor expression may contain several additional terms     C
789     !    that would need to be accounted for when satisfying the BC. These C
790     !    modifications have NOT been rigorously addressed.                 C
791     !                                                                      C
792     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
793           DOUBLE PRECISION FUNCTION F_HW(g0,EPS,EPG, ep_star_avg, &
794                                          g0EPs_avg,TH,Mu_g_avg,RO_g_avg,ROs_avg, &
795                                          DP_avg,K_12_avg, Tau_12_avg, Tau_1_avg, &
796                                          VREL, VSLIP, M)
797     
798     !-----------------------------------------------
799     ! Modules
800     !-----------------------------------------------
801           USE param
802           USE param1
803           USE constant
804           USE kintheory, only: epm, k_phi, s_star
805           USE physprop
806           USE run
807           USE fldvar
808           USE mpi_utility
809           IMPLICIT NONE
810     !-----------------------------------------------
811     ! Dummy Arguments
812     !-----------------------------------------------
813     ! Radial distribution function of solids phase M with each
814     ! other solids phase
815           DOUBLE PRECISION, INTENT(IN) :: g0(DIMENSION_M)
816     ! Average solids volume fraction of each solids phase
817           DOUBLE PRECISION, INTENT(IN) :: EPS(DIMENSION_M)
818     ! Average solids and gas volume fraction
819           DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
820     ! Sum of eps*G_0
821           DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
822     ! Average theta_m
823           DOUBLE PRECISION, INTENT(INOUT) :: TH (DIMENSION_M)
824     ! Average gas viscosity
825           DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
826     ! Average gas density
827           DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
828     ! Average solids density
829           DOUBLE PRECISION, INTENT(IN) :: ROS_avg(DIMENSION_M)
830     ! Average particle diameter of each solids phase
831           DOUBLE PRECISION, INTENT(IN) :: DP_avg(DIMENSION_M)
832     ! Average cross-correlation K_12 and lagrangian integral time-scale
833           DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
834     ! Magnitude of slip velocity between two phases
835           DOUBLE PRECISION, INTENT(IN) :: VREL
836     ! Slip velocity between wall and particles
837           DOUBLE PRECISION, INTENT(IN) :: VSLIP
838     ! Solids phase index
839           INTEGER, INTENT(IN) :: M
840     !-----------------------------------------------
841     ! Local Variables
842     !-----------------------------------------------
843     ! Solids phase index
844           INTEGER :: LL
845     ! Coefficient of 2nd term
846           DOUBLE PRECISION :: F_2
847     ! Coefficient of 1st term
848           DOUBLE PRECISION :: Mu_s
849     ! Viscosity
850           DOUBLE PRECISION :: Mu
851     ! Bulk viscosity
852           DOUBLE PRECISION :: Mu_b
853     ! Viscosity corrected for interstitial fluid effects
854           DOUBLE PRECISION :: Mu_star
855     ! Solids pressure
856           DOUBLE PRECISION :: Ps
857     ! Reynolds number based on slip velocity
858           DOUBLE PRECISION :: Re_g
859     ! Friction Factor in drag coefficient
860           DOUBLE PRECISION :: C_d
861     ! Drag Coefficient
862           DOUBLE PRECISION :: Beta, DgA
863     ! Constants in Simonin or Ahmadi model
864           DOUBLE PRECISION :: Sigma_c, Tau_2_c, Tau_12_st, Nu_t
865           DOUBLE PRECISION :: Tau_2, zeta_c_2, MU_2_T_Kin, Mu_2_Col
866           DOUBLE PRECISION :: Tmp_Ahmadi_Const
867     ! Variables for Iddir & Arastoopour model
868           DOUBLE PRECISION :: MU_sM_sum, MU_s_MM, MU_s_LM, MU_sM_ip, MU_common_term,&
869                               MU_sM_LM
870           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, D_PL, DPSUMo2
871           DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, Bp_lm
872     ! Variables for Garzo and Dufty model
873           DOUBLE PRECISION :: c_star, zeta0_star, nu_eta_star, press_star, &
874                               gamma_star, eta_k_star, eta_star, eta0
875     ! Variables for GTSH theory
876           DOUBLE PRECISION :: Chi, RdissP, Re_T, Rdiss, GamaAvg, A2_gtshW, &
877                               zeta_star, nu0, nuN, NuK, EDT_s, zeta_avg, etaK, &
878                               Mu_b_avg, mu2_0, mu4_0, mu4_1
879     !-----------------------------------------------
880     ! Functions
881     !-----------------------------------------------
882     ! Variable specularity coefficient
883           DOUBLE PRECISION :: PHIP_JJ
884     !-----------------------------------------------
885     
886     ! This is done here similar to bc_theta to avoid small negative values of
887     ! Theta coming most probably from linear solver
888           IF(TH(M) .LE. ZERO)THEN
889             TH(M) = 1D-8
890             IF (myPE.eq.PE_IO) THEN
891                WRITE(*,*) &
892                  'Warning: Negative granular temp at wall set to 1e-8'
893     !          CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
894             ENDIF
895           ENDIF
896     
897     
898     ! common variables
899           D_PM = DP_avg(M)
900           M_PM = (PI/6.d0)*(D_PM**3)*ROS_avg(M)
901           NU_PM = (EPS(M)*ROS_avg(M))/M_PM
902     
903     ! This is from Wen-Yu correlation, you can put here your own single particle drag
904           Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
905           IF (Re_g .lt. 1000.d0) THEN
906              C_d = (24.d0/(Re_g+SMALL_NUMBER))*(ONE + 0.15d0 * Re_g**0.687d0)
907           ELSE
908              C_d = 0.44d0
909           ENDIF
910           DgA = 0.75d0*C_d*Ro_g_avg*EPG*VREL/(D_PM*EPG**(2.65d0))
911           IF(VREL == ZERO) DgA = LARGE_NUMBER
912           Beta = EPS(M)*DgA !this is equivalent to F_gs(ijk,m)
913     
914     
915     ! Defining viscosity according to each KT for later use in JJ BC
916     ! Also defining solids pressure according to each KT for use if JENKINS
917     ! -------------------------------------------------------------------
918           SELECT CASE (KT_TYPE_ENUM)
919           CASE (LUN_1984)
920              Mu = (5d0*DSQRT(Pi*TH(M))*DP_avg(M)*ROS_avg(M))/96d0
921              Mu_b = (256d0*Mu*EPS(M)*g0EPs_avg)/(5d0*Pi)
922     
923              IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
924                 Mu_star = Mu
925              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
926                 MU_star = ZERO
927              ELSE
928                 Mu_star = ROS_avg(M)*EPS(M)* g0(M)*TH(M)* Mu/ &
929                    (ROS_avg(M)*g0EPs_avg*TH(M) + 2.0d0*DgA/ROS_avg(M)* Mu)
930              ENDIF
931     
932              Mu_s = ((2d0+ALPHA)/3d0)*((Mu_star/(Eta*(2d0-Eta)*&
933                 g0(M)))*(ONE+1.6d0*Eta*g0EPs_avg)*&
934                 (ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
935                 g0EPs_avg)+(0.6d0*Mu_b*Eta))
936     
937     ! defining granular pressure (for Jenkins BC)
938              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
939     
940     
941           CASE (SIMONIN_1996)
942     ! particle relaxation time
943              Tau_12_st = ROS_avg(M)/(DgA+small_number)
944     
945              Sigma_c = (ONE+ C_e)*(3.d0-C_e)/5.d0
946              Tau_2_c = DP_avg(M)/(6.d0*EPS(M)*g0(M)*DSQRT(16.d0*(TH(M)+Small_number)/PI))
947              zeta_c_2= 2.D0/5.D0*(ONE+ C_e)*(3.d0*C_e-ONE)
948              Nu_t =  Tau_12_avg/Tau_12_st
949              Tau_2 = ONE/(2.D0/Tau_12_st+Sigma_c/Tau_2_c)
950              MU_2_T_Kin = (2.0D0/3.0D0*K_12_avg*Nu_t + TH(M) * &
951                   (ONE+ zeta_c_2*EPS(M)*g0(M)))*Tau_2
952              Mu_2_Col = 8.D0/5.D0*EPS(M)*g0(M)*Eta* (MU_2_T_Kin+ &
953                    DP_avg(M)*DSQRT(TH(M)/PI))
954              Mu_s = EPS(M)*ROS_avg(M)*(MU_2_T_Kin + Mu_2_Col)
955     
956     ! defining granular pressure (for Jenkins BC)
957              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
958     
959     
960           CASE (AHMADI_1995)
961     ! particle relaxation time
962              Tau_12_st = ROS_avg(M)/(DgA+small_number)
963     
964              IF(EPS(M) < (ONE-ep_star_avg)) THEN
965                 Tmp_Ahmadi_Const = ONE/&
966                    (ONE+ Tau_1_avg/Tau_12_st * (ONE-EPS(M)/(ONE-ep_star_avg))**3)
967              ELSE
968                 Tmp_Ahmadi_Const = ONE
969              ENDIF
970              Mu_s = Tmp_Ahmadi_Const &
971                 *0.1045D0*(ONE/g0(M)+3.2D0*EPS(M)+12.1824D0*g0(M)*EPS(M)*EPS(M))  &
972                 *DP_avg(M)*ROS_avg(M)* DSQRT(TH(M))
973     
974     ! defining granular pressure (for Jenkins BC)
975              Ps = ROs_avg(M)*EPS(M)*TH(M)*&
976                  ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
977     
978     
979           CASE(IA_2005)
980     ! Use original IA theory if SWITCH_IA is false
981              IF(.NOT. SWITCH_IA) g0EPs_avg = EPS(M)*ROS_avg(M)
982     
983              Mu = (5.d0/96.d0)*D_PM* ROS_avg(M)*DSQRT(PI*TH(M)/M_PM)
984     
985              IF(SWITCH == ZERO .OR. RO_g_avg == ZERO)THEN
986                 Mu_star = Mu
987              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
988                 MU_star = ZERO
989              ELSE
990                 Mu_star = Mu*EPS(M)*g0(M)/ (g0EPs_avg+ 2.0d0*DgA*Mu / &
991                    (ROS_avg(M)**2 *(TH(M)/M_PM)))
992              ENDIF
993     
994              MU_s_MM = (Mu_star/g0(M))*(1.d0+(4.d0/5.d0)*(1.d0+C_E)*g0EPs_avg)**2
995              Mu_sM_sum = ZERO
996     
997              DO LL = 1, SMAX
998                 D_PL = DP_avg(LL)
999                 M_PL = (PI/6.d0)*(D_PL**3.)*ROS_avg(LL)
1000                 MPSUM = M_PM + M_PL
1001                 DPSUMo2 = (D_PM+D_PL)/2.d0
1002                 NU_PL = (EPS(LL)*ROS_avg(LL))/M_PL
1003     
1004                 IF ( LL .eq. M) THEN
1005                    Ap_lm = MPSUM/(2.d0)
1006                    Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
1007                    R1p_lm = ONE/( Ap_lm**1.5 * Dp_lm**3 )
1008                    MU_s_LM = DSQRT(PI)*(DPSUMo2**4 / (48.d0*5.d0))*&
1009                       g0(LL)*(M_PL*M_PM/MPSUM)*(M_PL*M_PM/&
1010                       MPSUM)*((M_PL*M_PM)**1.5)*NU_PM*NU_PL*&
1011                       (1.d0+C_E)*R1p_lm*DSQRT(TH(M))
1012     
1013     ! solids phase 'viscosity' associated with the divergence
1014     ! of solids phase M
1015                    MU_sM_ip = (MU_s_MM + MU_s_LM)
1016     
1017                 ELSE
1018                    Ap_lm = (M_PM*TH(LL)+M_PL*TH(M))/&
1019                       (2.d0)
1020                    Bp_lm = (M_PM*M_PL*(TH(LL)-TH(M) ))/&
1021                       (2.d0*MPSUM)
1022                    Dp_lm = (M_PL*M_PM*(M_PM*TH(M)+M_PL*TH(LL) ))/&
1023                       (2.d0*MPSUM*MPSUM)
1024                    R1p_lm = ( ONE/( Ap_lm**1.5 * Dp_lm**3 ) ) + &
1025                       ( (9.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**4 ) )+&
1026                       ( (30.d0*Bp_lm**4)/( 2.d0*Ap_lm**3.5 * Dp_lm**5 ) )
1027                    MU_common_term = DSQRT(PI)*(DPSUMo2**4 / (48.d0*5.d0))*&
1028                       g0(LL)*(M_PL*M_PM/MPSUM)*(M_PL*M_PM/&
1029                       MPSUM)*((M_PL*M_PM)**1.5)*NU_PM*NU_PL*&
1030                       (1.d0+C_E)*R1p_lm
1031                    MU_sM_LM = MU_common_term*(TH(M)*TH(M)*TH(LL)*TH(LL)*TH(LL) )
1032     
1033     ! solids phase 'viscosity' associated with the divergence
1034     ! of solids phase M
1035                    MU_sM_ip = MU_sM_LM
1036     
1037                 ENDIF
1038                 MU_sM_sum = MU_sM_sum + MU_sM_ip
1039     
1040               ENDDO
1041     
1042     ! Find the term proportional to the gradient in velocity
1043     ! of phase M  (viscosity in the Mth solids phase)
1044               Mu_s = MU_sM_sum
1045     
1046     
1047           CASE(GD_1999)
1048     ! Pressure/Viscosity/Bulk Viscosity
1049     ! Note: k_boltz = M_PM
1050     
1051     ! Find pressure in the Mth solids phase
1052              press_star = 1.d0 + 2.d0*(1.d0+C_E)*EPS(M)*G0(M)
1053     
1054     ! find bulk and shear viscosity
1055              c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
1056                 / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
1057     
1058              zeta0_star = (5.d0/12.d0)*G0(M)*(1.d0 - C_E*C_E) &
1059                 * (1.d0 + (3.d0/32.d0)*c_star)
1060     
1061              nu_eta_star = G0(M)*(1.d0 - 0.25d0*(1.d0-C_E)*(1.d0-C_E)) &
1062                 * (1.d0-(c_star/64.d0))
1063     
1064              gamma_star = (4.d0/5.d0)*(32.d0/PI)*EPS(M)*EPS(M) &
1065                 * G0(M)*(1.d0+C_E) * (1.d0 - (c_star/32.d0))
1066     
1067              eta_k_star = (1.d0 - (2.d0/5.d0)*(1.d0+C_E)*(1.d0-3.d0*C_E) &
1068                 * EPS(M)*G0(M) ) / (nu_eta_star - 0.5d0*zeta0_star)
1069     
1070              eta_star = eta_k_star*(1.d0 + (4.d0/5.d0)*EPS(M)*G0(M) &
1071                 * (1.d0+C_E) ) + (3.d0/5.d0)*gamma_star
1072     
1073              eta0 = 5.0d0*M_PM*DSQRT(TH(M)/PI) / (16.d0*D_PM*D_PM)
1074     
1075     ! added Ro_g = 0 for granular flows (no gas).
1076              IF(SWITCH == ZERO .OR. RO_g_avg == ZERO) THEN
1077                 Mu_star = eta0
1078              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
1079                 Mu_star = ZERO
1080              ELSE
1081                 Mu_star = ROS_avg(M)*EPS(M)*G0(M)*TH(M)*eta0 / &
1082                    ( ROS_avg(M)*EPS(M)*G0(M)*TH(M) + &
1083                    (2.d0*DgA*eta0/ROS_avg(M)) )     ! Note dgA is ~F_gs/ep_s
1084              ENDIF
1085     
1086     !  shear viscosity in Mth solids phase  (add to frictional part)
1087              Mu_s = Mu_star * eta_star
1088     
1089     ! defining granular pressure (for Jenkins BC)
1090              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1091                         ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1092     
1093     
1094           CASE (GTSH_2012)
1095              Chi = g0(M)
1096     
1097              RdissP = one
1098              IF(EPS(M) > small_number) RdissP = &
1099                 one + 3d0*dsqrt(EPS(M)/2d0) + 135d0/64d0*EPS(M)*dlog(EPS(M)) + &
1100                 11.26d0*EPS(M)*(one-5.1d0*EPS(M)+16.57d0*EPS(M)**2-21.77d0*    &
1101                 EPS(M)**3) - EPS(M)*Chi*dlog(epM)
1102     
1103              Re_T = RO_g_avg*D_pm*dsqrt(TH(M)) / Mu_g_avg
1104              Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
1105              Rdiss = RdissP + Re_T * K_phi(EPS(M))
1106              GamaAvg = 3d0*Pi*Mu_g_avg*D_pm*Rdiss
1107              mu2_0 = dsqrt(2d0*pi) * Chi * (one-C_E**2)
1108              mu4_0 = (4.5d0+C_E**2) * mu2_0
1109              mu4_1 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + 2d0*dsqrt(2d0*pi)* &
1110                 Chi*(one+C_E)
1111              A2_gtshW = zero ! for eps = zero
1112     
1113              if((EPS(M)> small_number) .and. (TH(M) > small_number)) then
1114                 zeta_star = 4.5d0*dsqrt(2d0*Pi)*(RO_g_avg/ROs_avg(M))**2*Re_g**2 * &
1115                             S_star(EPS(M),Chi) / (EPS(M)*(one-EPS(M))**2 * Re_T**4)
1116                 A2_gtshW = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
1117                                (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
1118              endif
1119     
1120              eta0 = 0.3125d0/(dsqrt(pi)*D_PM**2)*M_pm*dsqrt(TH(M))
1121              nu0 = (96.d0/5.d0)*(EPS(M)/D_PM)*DSQRT(TH(M)/PI)
1122              nuN = 0.25d0*nu0*Chi*(3d0-C_E)*(one+C_E) * &
1123                 (one+0.7375d0*A2_gtshW)
1124              NuK = nu0*(one+C_E)/3d0*Chi*( one+2.0625d0*(one-C_E)+ &
1125                 ((947d0-579*C_E)/256d0*A2_gtshW) )
1126              EDT_s = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
1127                 (one+0.1875d0*A2_gtshW)*NU_PM*D_PM**2*dsqrt(TH(M))
1128     
1129              if((TH(m) > small_number) .and. (EPS(M) > small_number)) then
1130                 zeta_avg = one/6d0*D_PM*VREL**2*(3d0*pi*Mu_g_avg*D_PM/&
1131                    M_pm)**2 / dsqrt(TH(m)) * S_star(EPS(M), Chi)
1132                 etaK = ROs_avg(M)*EPS(M)*TH(m) / (nuN-0.5d0*( &
1133                    EDT_s-zeta_avg/TH(m) - 2d0*GamaAvg/M_PM)) * &
1134                    (one -0.4d0 * (one+C_E)*(one-3d0*C_E)*EPS(M)*Chi)
1135              else
1136                 etaK = zero
1137              endif
1138     
1139              Mu_b_avg = 25.6d0/pi * EPS(M)**2 * Chi *(one+C_E) * &
1140                 (one - A2_gtshW/16d0)*eta0
1141     
1142              Mu_s = etaK*(one+0.8d0*EPS(M)*Chi*(one+C_E)) + &
1143                 0.6d0*Mu_b_avg
1144     
1145     ! defining granular pressure (for Jenkins BC)
1146              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1147     
1148     
1149           CASE DEFAULT
1150     ! should never hit this
1151              WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1152              WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1153              call mfix_exit(myPE)
1154           END SELECT
1155     
1156     
1157     ! setting the coefficients for JJ BC
1158     ! -------------------------------------------------------------------
1159           SELECT CASE (KT_TYPE_ENUM)
1160     
1161              CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999, GTSH_2012)
1162     ! KT without particle mass in their definition of granular temperature
1163     ! and theta = granular temperature  OR
1164     ! KT with particle mass in their definition of granular temperature
1165     ! and theta = granular temperature/mass
1166     
1167     ! Jenkins BC is implemented for these KT
1168                 IF(JENKINS) THEN
1169                    IF (VSLIP == ZERO) THEN
1170     ! if solids velocity field is initialized to zero, use free slip bc
1171                       F_2 = zero
1172                    ELSE
1173     ! As I understand from soil mechanic papers, the coefficient mu in
1174     ! Jenkins paper is tan_Phi_w. T.  See for example, G.I. Tardos, PT,
1175     ! 92 (1997), 61-74, equation (1). sof
1176     
1177     ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
1178                       F_2 = tan_Phi_w*Ps/VSLIP
1179     
1180                    ENDIF
1181                 ELSE
1182                    IF(.NOT. BC_JJ_M) THEN
1183                       F_2 = (PHIP*DSQRT(3d0*TH(M))*Pi*ROS_avg(M)*EPS(M)*&
1184                          g0(M))/(6d0*(ONE-ep_star_avg))
1185                    ELSE
1186                       F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3d0*TH(M))*Pi*&
1187                          ROS_avg(M)*EPS(M)*g0(M))/(6d0*(ONE-ep_star_avg))
1188                    ENDIF
1189                 ENDIF   ! end if(Jenkins)/else
1190     
1191     
1192             CASE (IA_2005)
1193     ! KTs with particle mass in their definition of granular temperature
1194     ! and theta = granular temperature
1195                 IF(.NOT. BC_JJ_M) THEN
1196                    F_2 = (PHIP*DSQRT(3.d0*TH(M)/M_PM)*PI*ROS_avg(M)*EPS(M)*&
1197                       g0(M))/(6.d0*(ONE-ep_star_avg))
1198                 ELSE
1199                    F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3.d0*TH(M)/M_PM)*PI*&
1200                       ROS_avg(M)*EPS(M)*g0(M))/(6.d0*(ONE-ep_star_avg))
1201                 ENDIF
1202     
1203     
1204              CASE DEFAULT
1205     ! should never hit this
1206                 WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1207                 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1208                 call mfix_exit(myPE)
1209           END SELECT
1210     
1211           F_HW =  F_2/Mu_s
1212     
1213           RETURN
1214           END FUNCTION F_HW
1215     
1216     
1217     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1218     !                                                                      C
1219     !  Subroutine CG_CALC_GRBDRY                                           C
1220     !  Purpose: Cut cell version                                           C
1221     !                                                                      C
1222     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1223     
1224           SUBROUTINE CG_CALC_GRBDRY(IJK, TYPE_OF_CELL, M, L, F_2)
1225     
1226     !-----------------------------------------------
1227     ! Modules
1228     !-----------------------------------------------
1229           USE bc
1230           USE compar
1231           USE constant
1232           USE fldvar
1233           USE fun_avg
1234           USE functions
1235           USE geometry
1236           USE indices
1237           USE param
1238           USE param1
1239           USE physprop
1240           USE rdf
1241           USE run
1242           USE turb
1243           USE visc_s
1244           Use cutcell
1245           use toleranc
1246           IMPLICIT NONE
1247     !-----------------------------------------------
1248     ! Dummy Arguments
1249     !-----------------------------------------------
1250     ! IJK indices
1251           INTEGER, INTENT(IN) :: IJK
1252     
1253           CHARACTER (LEN=*) :: TYPE_OF_CELL
1254     ! Solids phase index
1255           INTEGER, INTENT(IN) :: M
1256     ! Index corresponding to boundary condition
1257           INTEGER, INTENT(IN) ::  L
1258     
1259           DOUBLE PRECISION :: F_2
1260     
1261     !-----------------------------------------------
1262     ! Local Variables
1263     !-----------------------------------------------
1264     ! IJK indices
1265           INTEGER          IJK1, IJK2
1266     ! Other indices
1267           INTEGER          IJK2E, IPJK2, IPJMK2
1268     !
1269           DOUBLE PRECISION :: smallTheta
1270     ! Average scalars
1271           DOUBLE PRECISION EPg_avg, Mu_g_avg, RO_g_avg
1272     ! Average void fraction at packing
1273           DOUBLE PRECISION ep_star_avg
1274     ! Average scalars modified to include all solid phases
1275           DOUBLE PRECISION EPs_avg(DIMENSION_M), DP_avg(DIMENSION_M),&
1276                            TH_avg(DIMENSION_M)
1277     ! Average solids density
1278           DOUBLE PRECISION ROS_AVG(DIMENSION_M)
1279     ! Average Simonin and Ahmadi variables (sof)
1280           DOUBLE PRECISION K_12_avg, Tau_12_avg, Tau_1_avg
1281     ! Average velocities
1282           DOUBLE PRECISION WGC1, WGC2
1283     ! Solids phase index
1284           INTEGER          MM
1285     ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
1286           DOUBLE PRECISION USCM, VSCM,WSCM
1287           DOUBLE PRECISION WSCM1,WSCM2
1288     ! values of U_g, V_g, W_g at appropriate place on boundary wall
1289           DOUBLE PRECISION UGC, VGC, WGC
1290     ! Magnitude of gas-solids relative velocity
1291           DOUBLE PRECISION VREL
1292     ! slip velocity between wall and particles for Jenkins bc (sof)
1293           DOUBLE PRECISION VSLIP
1294     ! radial distribution function at contact
1295           DOUBLE PRECISION g0(DIMENSION_M)
1296     ! Sum of eps*G_0
1297           DOUBLE PRECISION g0EPs_avg
1298     !-----------------------------------------------
1299     
1300     !  Note:  EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).  Hence
1301     !         IJK2 (fluid cell) is used in averages.
1302     !
1303     
1304           smallTheta = (to_SI)**4 * ZERO_EP_S
1305     
1306           SELECT CASE (TYPE_OF_CELL)
1307     
1308              CASE('U_MOMENTUM')
1309     
1310                 IJK2 = EAST_OF(IJK)
1311                 IJK2E = EAST_OF(IJK2)
1312     
1313                 EPg_avg = (VOL(IJK)*EP_g(IJK) + VOL(IJK2)*EP_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1314                 ep_star_avg = (VOL(IJK)*EP_star_array(IJK) + VOL(IJK2)*EP_star_array(IJK2))/(VOL(IJK) + VOL(IJK2))
1315                 Mu_g_avg = (VOL(IJK)*Mu_g(IJK) + VOL(IJK2)*Mu_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1316                 RO_g_avg = (VOL(IJK)*RO_g(IJK) + VOL(IJK2)*RO_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1317     
1318                 K_12_avg = ZERO
1319                 Tau_12_avg = ZERO
1320                 Tau_1_avg = ZERO
1321                 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR.&
1322                    KT_TYPE_ENUM == AHMADI_1995) THEN  ! not converted to CG
1323                    K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
1324                    Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
1325                    Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
1326                 ENDIF
1327     
1328                 g0EPs_avg = ZERO
1329                 DO MM = 1, SMAX
1330                    g0(MM)      = G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM)
1331                    EPs_avg(MM) = (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1332                    DP_avg(MM)  = (VOL(IJK)*D_P(IJK, MM) + VOL(IJK2)*D_P(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1333                    g0EPs_avg   = g0EPs_avg + G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM) &
1334                                * (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1335                    ROS_AVG(MM) = (VOL(IJK)*RO_S(IJK, MM) + VOL(IJK2)*RO_S(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1336     
1337                    TH_avg(MM)  = (VOL(IJK)*Theta_m(IJK,MM) + VOL(IJK2)*Theta_m(IJK2,MM))/(VOL(IJK) + VOL(IJK2))
1338                    IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
1339                 ENDDO
1340     
1341     
1342     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)  ! not converted to CG
1343                 UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
1344                 VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
1345                 WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
1346                              AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
1347                              I_OF(IJK2))
1348                 WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
1349                              AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
1350                              I_OF(IJK1))
1351                 WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
1352                 USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
1353                 VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
1354                 WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
1355                              AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
1356                              I_OF(IJK2))
1357                 WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
1358                              AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
1359                              I_OF(IJK1))
1360                 WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
1361     
1362     ! magnitude of gas-solids relative velocity
1363                 VREL = DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + &
1364                               (WGC - WSCM)**2 )
1365     
1366     ! slip velocity for use in Jenkins bc (sof)
1367                 VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
1368                            + (WSCM-BC_WW_S(L,M))**2 )
1369     
1370                 CALL GET_CG_F2(g0, EPs_avg, EPg_avg, ep_star_avg, &
1371                           g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROS_AVG,&
1372                            DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, &
1373                           VREL, VSLIP, M,F_2)
1374     
1375     
1376              CASE('V_MOMENTUM')
1377     
1378                 IJK2 = NORTH_OF(IJK)
1379     
1380                 EPg_avg = (VOL(IJK)*EP_g(IJK) + VOL(IJK2)*EP_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1381                 ep_star_avg = (VOL(IJK)*EP_star_array(IJK) + VOL(IJK2)*EP_star_array(IJK2))/(VOL(IJK) + VOL(IJK2))
1382                 Mu_g_avg = (VOL(IJK)*Mu_g(IJK) + VOL(IJK2)*Mu_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1383                 RO_g_avg = (VOL(IJK)*RO_g(IJK) + VOL(IJK2)*RO_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1384     
1385                 K_12_avg = ZERO
1386                 Tau_12_avg = ZERO
1387                 Tau_1_avg = ZERO
1388                 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR.&
1389                    KT_TYPE_ENUM == AHMADI_1995) THEN  ! not converted to CG
1390                    Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
1391     ! Simonin only:
1392                    K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
1393                    Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
1394                 ENDIF
1395     
1396                 g0EPs_avg = ZERO
1397                 DO MM = 1, SMAX
1398                    g0(MM)      = G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM)
1399                    EPs_avg(MM) = (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1400                    DP_avg(MM)  = (VOL(IJK)*D_P(IJK, MM) + VOL(IJK2)*D_P(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1401                    g0EPs_avg   = g0EPs_avg + G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM) &
1402                                * (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1403                    ROS_avg(MM) = (VOL(IJK)*RO_S(IJK, MM) + VOL(IJK2)*RO_S(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1404                    TH_avg(MM)  = (VOL(IJK)*Theta_m(IJK,MM) + VOL(IJK2)*Theta_m(IJK2,MM))/(VOL(IJK) + VOL(IJK2))
1405                    IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
1406                 ENDDO
1407     
1408     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)    ! not converted to CG
1409                 UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
1410                 VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
1411                 WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
1412                              AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
1413                              I_OF(IJK2))
1414                 WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
1415                              AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
1416                              I_OF(IJK1))
1417                 WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
1418                 USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
1419                 VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
1420                 WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
1421                              AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
1422                              I_OF(IJK2))
1423                 WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
1424                              AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
1425                              I_OF(IJK1))
1426                 WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
1427     
1428     ! magnitude of gas-solids relative velocity
1429     
1430                 VREL = DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + &
1431                               (WGC - WSCM)**2 )
1432     
1433     ! slip velocity for use in Jenkins bc (sof)
1434                 VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
1435                             + (WSCM-BC_WW_S(L,M))**2 )
1436     
1437     
1438                 CALL GET_CG_F2(g0, EPs_avg, EPg_avg, ep_star_avg, &
1439                           g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROS_AVG,&
1440                           DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, &
1441                           VREL, VSLIP, M,F_2)
1442     
1443     
1444              CASE('W_MOMENTUM')
1445     
1446                 IJK2 = TOP_OF(IJK)
1447     
1448                 EPg_avg = (VOL(IJK)*EP_g(IJK) + VOL(IJK2)*EP_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1449                 ep_star_avg = (VOL(IJK)*EP_star_array(IJK) + VOL(IJK2)*EP_star_array(IJK2))/(VOL(IJK) + VOL(IJK2))
1450                 Mu_g_avg = (VOL(IJK)*Mu_g(IJK) + VOL(IJK2)*Mu_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1451                 RO_g_avg = (VOL(IJK)*RO_g(IJK) + VOL(IJK2)*RO_g(IJK2))/(VOL(IJK) + VOL(IJK2))
1452     
1453                 K_12_avg = ZERO
1454                 Tau_12_avg = ZERO
1455                 Tau_1_avg = ZERO
1456                 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR.&
1457                    KT_TYPE_ENUM == AHMADI_1995) THEN  ! not converted to CG
1458                    Tau_1_avg = AVG_X(Tau_1(IJK2), Tau_1(IJK2E), I_OF(IJK2))
1459     ! Simonon only:
1460                    K_12_avg = AVG_X(K_12(IJK2), K_12(IJK2E), I_OF(IJK2))
1461                    Tau_12_avg = AVG_X(Tau_12(IJK2), Tau_12(IJK2E), I_OF(IJK2))
1462                 ENDIF
1463     
1464                 g0EPs_avg = ZERO
1465                 DO MM = 1, SMAX
1466                    g0(MM)      = G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM)
1467                    EPs_avg(MM) = (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1468                    DP_avg(MM)  = (VOL(IJK)*D_P(IJK, MM) + VOL(IJK2)*D_P(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1469                    ROS_avg(MM) = (VOL(IJK)*RO_S(IJK, MM) + VOL(IJK2)*RO_S(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1470                    g0EPs_avg   = g0EPs_avg + G_0AVG(IJK, IJK, 'X', I_OF(IJK), M, MM) &
1471                                * (VOL(IJK)*EP_s(IJK, MM) + VOL(IJK2)*EP_s(IJK2, MM))/(VOL(IJK) + VOL(IJK2))
1472     
1473                    TH_avg(MM)  = (VOL(IJK)*Theta_m(IJK,MM) + VOL(IJK2)*Theta_m(IJK2,MM))/(VOL(IJK) + VOL(IJK2))
1474                    IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
1475                 ENDDO
1476     
1477     ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)    ! not converted to CG
1478                 UGC  = AVG_Y(U_g(IJK1), U_g(IJK2),J_OF(IJK1))
1479                 VGC  = AVG_X(V_g(IJK1), V_g(IPJMK2),I_OF(IJK1))
1480                 WGC1 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
1481                              AVG_Z_T(W_g(KM_OF(IPJK2)), W_g(IPJK2)),&
1482                              I_OF(IJK2))
1483                 WGC2 = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
1484                              AVG_Z_T(W_g(KM_OF(IPJMK2)), W_g(IPJMK2)),&
1485                              I_OF(IJK1))
1486                 WGC  = AVG_Y(WGC2, WGC1, J_OF(IJK1))
1487                 USCM = AVG_Y(U_s(IJK1,M), U_s(IJK2,M),J_OF(IJK1))
1488                 VSCM = AVG_X(V_s(IJK1,M), V_s(IPJMK2,M),I_OF(IJK1))
1489                 WSCM1= AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M),W_s(IJK2,M)),&
1490                              AVG_Z_T(W_s(KM_OF(IPJK2),M),W_s(IPJK2,M)),&
1491                              I_OF(IJK2))
1492                 WSCM2= AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M),W_s(IJK1,M)),&
1493                              AVG_Z_T(W_s(KM_OF(IPJMK2),M),W_s(IPJMK2,M)),&
1494                              I_OF(IJK1))
1495                 WSCM = AVG_Y(WSCM2, WSCM1, J_OF(IJK1))
1496     
1497     ! magnitude of gas-solids relative velocity
1498                 VREL = DSQRT( (UGC - USCM)**2 + (VGC - VSCM)**2 + &
1499                               (WGC - WSCM)**2 )
1500     
1501     ! slip velocity for use in Jenkins bc (sof)
1502                 VSLIP= DSQRT( (USCM-BC_UW_S(L,M))**2 + (VSCM-BC_VW_S(L,M))**2 &
1503                             + (WSCM-BC_WW_S(L,M))**2 )
1504     
1505     
1506     ! why bother with this since it should not come here...
1507                 CALL GET_CG_F2(g0, EPs_avg, EPg_avg, ep_star_avg, &
1508                           g0EPs_avg, TH_avg, Mu_g_avg, RO_g_avg, ROS_AVG,&
1509                           DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, &
1510                           VREL, VSLIP, M,F_2)
1511     
1512     
1513              CASE DEFAULT
1514                 WRITE(*,*) 'CG_CALC_GRBDRY'
1515                 WRITE(*,*) 'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1516                 WRITE(*,*) 'ACCEPTABLE TYPES ARE:'
1517                 WRITE(*,*) 'U_MOMENTUM'
1518                 WRITE(*,*) 'V_MOMENTUM'
1519                 WRITE(*,*) 'W_MOMENTUM'
1520                 CALL MFIX_EXIT(myPE)
1521               END SELECT
1522     
1523     
1524           RETURN
1525     
1526           END SUBROUTINE CG_CALC_GRBDRY
1527     
1528     
1529     
1530     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1531     !                                                                      C
1532     !  Subroutine: GET_CG_F2                                               C
1533     !  Purpose: Compute F_2 for cut cell version                           C
1534     !                                                                      C
1535     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1536           SUBROUTINE GET_CG_F2(g0,EPS,EPG, ep_star_avg, &
1537                                          g0EPs_avg,TH,Mu_g_avg,RO_g_avg,Ros_avg, &
1538                                          DP_avg,K_12_avg, Tau_12_avg, Tau_1_avg, &
1539                                          VREL, VSLIP, M,F_2)
1540     
1541     !-----------------------------------------------
1542     ! Modules
1543     !-----------------------------------------------
1544           USE param
1545           USE param1
1546           USE constant
1547           USE physprop
1548           USE run
1549           USE fldvar
1550           USE mpi_utility
1551           Use cutcell
1552           IMPLICIT NONE
1553     !-----------------------------------------------
1554     ! Dummy Arguments
1555     !-----------------------------------------------
1556     ! Radial distribution function of solids phase M with each
1557     ! other solids phase
1558           DOUBLE PRECISION, INTENT(IN) :: g0(DIMENSION_M)
1559     ! Average solids volume fraction of each solids phase
1560           DOUBLE PRECISION, INTENT(IN) :: EPS(DIMENSION_M)
1561     ! Average solids and gas volume fraction
1562           DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
1563     ! Sum of eps*G_0
1564           DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
1565     ! Average theta_m
1566           DOUBLE PRECISION, INTENT(INOUT) :: TH (DIMENSION_M)
1567     ! Average gas viscosity
1568           DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
1569     ! Average gas density
1570           DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
1571     !QX: Average solids density
1572           DOUBLE PRECISION, INTENT(IN) :: ROS_avg(DIMENSION_M)
1573     ! Average particle diameter of each solids phase
1574           DOUBLE PRECISION, INTENT(IN) :: DP_avg(DIMENSION_M)
1575     ! Average cross-correlation K_12 and lagrangian integral time-scale
1576           DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
1577     ! Magnitude of slip velocity between two phases
1578           DOUBLE PRECISION, INTENT(IN) :: VREL
1579     ! Slip velocity between wall and particles
1580           DOUBLE PRECISION, INTENT(IN) :: VSLIP
1581     ! Solids phase index
1582           INTEGER, INTENT(IN) :: M
1583     !-----------------------------------------------
1584     ! Local Variables
1585     !-----------------------------------------------
1586     ! Solids pressure
1587           DOUBLE PRECISION :: Ps
1588           DOUBLE PRECISION :: F_2
1589           DOUBLE PRECISION :: D_PM, M_PM
1590     !-----------------------------------------------
1591     ! Functions
1592     !-----------------------------------------------
1593     ! Variable specularity coefficient
1594           DOUBLE PRECISION :: PHIP_JJ
1595     !-----------------------------------------------
1596     
1597     ! This is done here similar to bc_theta to avoid small negative values of
1598     ! Theta coming most probably from linear solver
1599           IF(TH(M) .LE. ZERO)THEN
1600             TH(M) = 1D-8
1601             IF (myPE.eq.PE_IO) THEN
1602                WRITE(*,*) &
1603                   'Warning: Negative granular temp at wall set to 1e-8'
1604     !          CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
1605             ENDIF
1606           ENDIF
1607     
1608     ! common variables
1609           D_PM = DP_avg(M)
1610           M_PM = (PI/6.d0)*(D_PM**3)*ROS_avg(M)
1611     
1612     ! Defining viscosity according to each KT for later use in JJ BC
1613     ! Also defining solids pressure according to each KT for use if JENKINS
1614     ! -------------------------------------------------------------------
1615           SELECT CASE (KT_TYPE_ENUM)
1616           CASE (LUN_1984)
1617     ! defining granular pressure (for Jenkins BC)
1618              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
1619     
1620     
1621           CASE (SIMONIN_1996)
1622     ! defining granular pressure (for Jenkins BC)
1623              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
1624     
1625     
1626           CASE (AHMADI_1995)
1627     ! defining granular pressure (for Jenkins BC)
1628              Ps = ROs_avg(M)*EPS(M)*TH(M)*&
1629                  ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
1630     
1631     
1632           CASE(IA_2005)
1633     ! Use original IA theory if SWITCH_IA is false
1634     ! no ps since no jenkins for this theory
1635     
1636           CASE(GD_1999)
1637     
1638     ! defining granular pressure (for Jenkins BC)
1639              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1640                         ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1641     
1642     
1643           CASE (GTSH_2012)
1644     ! defining granular pressure (for Jenkins BC)
1645              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1646     
1647     
1648           CASE DEFAULT
1649     ! should never hit this
1650              WRITE (*, '(A)') 'CG_CALC_GRBDRY => GET_CG_F2 '
1651              WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1652              call mfix_exit(myPE)
1653           END SELECT
1654     
1655     
1656     ! setting the coefficients for JJ BC
1657     ! -------------------------------------------------------------------
1658           SELECT CASE (KT_TYPE_ENUM)
1659     
1660              CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999, GTSH_2012)
1661     ! KT without particle mass in their definition of granular temperature
1662     ! and theta = granular temperature  OR
1663     ! KT with particle mass in their definition of granular temperature
1664     ! and theta = granular temperature/mass
1665     
1666     ! Jenkins BC is implemented for these KT
1667                 IF(JENKINS) THEN
1668                    IF (VSLIP == ZERO) THEN
1669     ! if solids velocity field is initialized to zero, use free slip bc
1670                       F_2 = zero
1671                    ELSE
1672     ! As I understand from soil mechanic papers, the coefficient mu in
1673     ! Jenkins paper is tan_Phi_w. T.  See for example, G.I. Tardos, PT,
1674     ! 92 (1997), 61-74, equation (1). sof
1675     
1676     ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
1677                       F_2 = tan_Phi_w*Ps/VSLIP
1678     
1679                    ENDIF
1680                 ELSE
1681                    IF(.NOT. BC_JJ_M) THEN
1682                       F_2 = (PHIP*DSQRT(3d0*TH(M))*Pi*ROS_avg(M)*EPS(M)*&
1683                          g0(M))/(6d0*(ONE-ep_star_avg))
1684                    ELSE
1685                       F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3d0*TH(M))*Pi*&
1686                          ROS_avg(M)*EPS(M)*g0(M))/(6d0*(ONE-ep_star_avg))
1687                    ENDIF
1688                 ENDIF   ! end if(Jenkins)/else
1689     
1690              CASE(IA_2005)
1691     ! KTs with particle mass in their definition of granular temperature
1692     ! and theta = granular temperature
1693                 IF(.NOT. BC_JJ_M) THEN
1694                    F_2 = (PHIP*DSQRT(3.d0*TH(M)/M_PM)*PI*ROS_avg(M)*&
1695                       EPS(M)*g0(M))/(6.d0*(ONE-ep_star_avg))
1696                 ELSE
1697                    F_2 = (PHIP_JJ(vslip,th(m))*DSQRT(3.d0*TH(M)/M_PM)*&
1698                       PI*ROS_avg(M)*EPS(M)*g0(M))/(6.d0*(ONE-ep_star_avg))
1699                 ENDIF
1700     
1701     
1702              CASE DEFAULT
1703     ! should never hit this
1704                 WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1705                 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1706                 call mfix_exit(myPE)
1707           END SELECT
1708     
1709     ! all the code used to calculate Mu_s was deleted from this routine
1710     !      F_HW =  F_2/Mu_s  ! Only F_2 is actually needed
1711     
1712     
1713           RETURN
1714           END SUBROUTINE GET_CG_F2
1715     
1716