File: N:\mfix\model\calc_grbdry.f

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