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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: BC_THETA                                                C
4     !                                                                      C
5     !  Purpose: Calculate boundary conditions for the temperature/         C
6     !     granular energy equation                                         C
7     !                                                                      C
8     !  Author: Kapil Agrawal, Princeton University        Date: 14-MAR-98  C
9     !  Reviewer:                                          Date:            C
10     !                                                                      C
11     !                                                                      C
12     !  Literature/Document References:                                     C
13     !     Johnson, P. C., and Jackson, R., Frictional-collisional          C
14     !        constitutive relations for granular materials, with           C
15     !        application to plane shearing, Journal of Fluid Mechanics,    C
16     !        1987, 176, pp. 67-93                                          C
17     !                                                                      C
18     !  Notes:                                                              C
19     !    If BC_JJ_PS = 3, then set the wall temperature equal to the       C
20     !       adjacent fluid cell temperature (i.e. zero flux)               C
21     !    IF BC_JJ_PS = 1, then set the wall temperature based on           C
22     !       Johnson and Jackson Pseudo-thermal temp BC                     C
23     !    IF BC_JJ_PS !=1, !=3 and >0, then set the wall temperature        C
24     !      based on Johnson and Jackson Pseudo-thermal temp BC without     C
25     !       the generation term                                            C
26     !                                                                      C
27     !                                                                      C
28     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
29     
30           SUBROUTINE BC_THETA(M, A_m, B_m)
31     
32     !-----------------------------------------------
33     ! Modules
34     !-----------------------------------------------
35           USE param
36           USE param1
37           USE parallel
38           USE matrix
39           USE scales
40           USE constant
41           USE toleranc
42           USE run
43           USE physprop
44           USE fldvar
45           USE visc_s
46           USE geometry
47           USE output
48           USE indices
49           USE bc
50           USE compar
51           USE mpi_utility
52           USE fun_avg
53           USE functions
54           IMPLICIT NONE
55     !-----------------------------------------------
56     ! Dummy arguments
57     !-----------------------------------------------
58     ! Solids phase index
59           INTEGER :: M
60     ! Septadiagonal matrix A_m
61           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
62     ! Vector b_m
63           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
64     !-----------------------------------------------
65     ! Local variables
66     !-----------------------------------------------
67     ! Boundary condition
68           INTEGER :: L
69     ! Indices
70           INTEGER :: I,  J, K, I1, I2, J1, J2, K1, K2, IJK,&
71                      IM, JM, KM
72     ! Granular energy coefficient
73           DOUBLE PRECISION :: Gw, Hw, Cw
74     !-----------------------------------------------
75     
76     ! Setup Johnson and Jackson Pseudo-thermal temp B.C.
77           DO L = 1, DIMENSION_BC
78             IF( BC_DEFINED(L) ) THEN
79               IF(BC_TYPE(L) .EQ. 'NO_SLIP_WALL' .OR.&
80                  BC_TYPE(L) .EQ. 'FREE_SLIP_WALL' .OR.&
81                  BC_TYPE(L) .EQ. 'PAR_SLIP_WALL' ) THEN
82                 I1 = BC_I_w(L)
83                 I2 = BC_I_e(L)
84                 J1 = BC_J_s(L)
85                 J2 = BC_J_n(L)
86                 K1 = BC_K_b(L)
87                 K2 = BC_K_t(L)
88     
89                 IF(BC_JJ_PS(L).GT.0) THEN
90                   DO K = K1, K2
91                   DO J = J1, J2
92                   DO I = I1, I2
93     
94                     IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
95                     IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
96                     IJK   = FUNIJK(I, J, K)
97                     IF (FLOW_AT(IJK)) CYCLE !checks for pressure outlets
98                     IM    = Im1(I)
99                     JM    = Jm1(J)
100                     KM    = Km1(K)
101     ! Setting the temperature to zero - recall the center coefficient
102     ! and source vector are negative. The off-diagonal coefficients
103     ! are positive.
104                     A_m(IJK, e, M) =  ZERO
105                     A_m(IJK, w, M) =  ZERO
106                     A_m(IJK, n, M) =  ZERO
107                     A_m(IJK, s, M) =  ZERO
108                     A_m(IJK, t, M) =  ZERO
109                     A_m(IJK, b, M) =  ZERO
110                     A_m(IJK, 0, M) = -ONE
111                     b_m(IJK, M)    =  ZERO
112     
113     ! Checking if the west wall
114                     IF(FLUID_AT(EAST_OF(IJK)))THEN
115                        IF(EP_s(EAST_OF(IJK),M).LE.DIL_EP_s)THEN
116                           A_m(IJK, e, M) = ONE
117                        ELSE
118                           IF (BC_JJ_PS(L).EQ.3) THEN
119     ! Setting the wall temperature equal to the adjacent fluid cell
120     ! temperature (i.e. zero flux)
121                              Gw = 1d0
122                              Hw = 0d0
123                              Cw = 0d0
124                           ELSE
125     ! Setting the wall temperature based on Johnson and Jackson
126     ! Pseudo-thermal temp B.C.
127                              CALL CALC_THETA_BC(IJK,EAST_OF(IJK),'E',&
128                                 M,L,gw,hw,cw)
129     ! Overwriting calculated value of cw
130                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
131                           ENDIF
132     
133                           A_m(IJK, e, M) = -(HALF*hw - oDX_E(I)*gw)
134                           A_m(IJK, 0, M) = -(HALF*hw + oDX_E(I)*gw)
135     
136                           IF (BC_JJ_PS(L) .EQ. 1) THEN
137                             b_m(IJK, M) = -cw
138                           ELSE
139                             b_m(IJK,M) = ZERO
140                           ENDIF
141                        ENDIF
142     
143                     ELSEIF(FLUID_AT(WEST_OF(IJK)))THEN
144                        IF(EP_s(WEST_OF(IJK),M).LE.DIL_EP_s)THEN
145                           A_m(IJK, w, M) = ONE
146                        ELSE
147                          IF (BC_JJ_PS(L).EQ.3) THEN
148                             Gw = 1d0
149                             Hw = 0d0
150                             Cw = 0d0
151                          ELSE
152                             CALL CALC_THETA_BC(IJK,WEST_OF(IJK),'W',&
153                                M,L,gw,hw,cw)
154                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
155                          ENDIF
156     
157                          A_m(IJK, w, M) = -(HALF*hw - oDX_E(IM)*gw)
158                          A_m(IJK, 0, M) = -(HALF*hw + oDX_E(IM)*gw)
159     
160                          IF (BC_JJ_PS(L) .EQ. 1) THEN
161                             b_m(IJK, M) = -cw
162                          ELSE
163                             b_m(IJK,M) = ZERO
164                          ENDIF
165                        ENDIF
166     
167                     ELSEIF(FLUID_AT(NORTH_OF(IJK)))THEN
168                        IF(EP_s(NORTH_OF(IJK),M).LE.DIL_EP_s)THEN
169                           A_m(IJK, n, M) = ONE
170                        ELSE
171                          IF (BC_JJ_PS(L).EQ.3) THEN
172                             Gw = 1d0
173                             Hw = 0d0
174                             Cw = 0d0
175                          ELSE
176                             CALL CALC_THETA_BC(IJK,NORTH_OF(IJK),'N',&
177                                M,L,gw,hw,cw)
178                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
179                          ENDIF
180     
181                          A_m(IJK, n, M) = -(HALF*hw - oDY_N(J)*gw)
182                          A_m(IJK, 0, M) = -(HALF*hw + oDY_N(J)*gw)
183     
184                          IF (BC_JJ_PS(L) .EQ. 1) THEN
185                             b_m(IJK, M) = -cw
186                          ELSE
187                             b_m(IJK,M) = ZERO
188                          ENDIF
189                        ENDIF
190     
191                     ELSEIF(FLUID_AT(SOUTH_OF(IJK)))THEN
192                        IF(EP_s(SOUTH_OF(IJK),M).LE.DIL_EP_s)THEN
193                           A_m(IJK, s, M) = ONE
194                        ELSE
195                          IF (BC_JJ_PS(L).EQ.3) THEN
196                             Gw = 1d0
197                             Hw = 0d0
198                             Cw = 0d0
199                          ELSE
200                             CALL CALC_THETA_BC(IJK,SOUTH_OF(IJK),'S',&
201                                M,L,gw,hw,cw)
202                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
203                          ENDIF
204     
205                          A_m(IJK, s, M) = -(HALF*hw - oDY_N(JM)*gw)
206                          A_m(IJK, 0, M) = -(HALF*hw + oDY_N(JM)*gw)
207     
208                          IF (BC_JJ_PS(L) .EQ. 1) THEN
209                             b_m(IJK, M) = -cw
210                          ELSE
211                             b_m(IJK,M) = ZERO
212                          ENDIF
213                        ENDIF
214     
215                     ELSEIF(FLUID_AT(TOP_OF(IJK)))THEN
216                        IF(EP_s(TOP_OF(IJK),M).LE.DIL_EP_s)THEN
217                           A_m(IJK, t, M) = ONE
218                        ELSE
219                          IF (BC_JJ_PS(L).EQ.3) THEN
220                             Gw = 1d0
221                             Hw = 0d0
222                             Cw = 0d0
223                          ELSE
224                             CALL CALC_THETA_BC(IJK,TOP_OF(IJK),'T',&
225                                M,L,gw,hw,cw)
226                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
227                          ENDIF
228     
229                          A_m(IJK, t, M) = -(HALF*hw - oX(I)*oDZ_T(K)*gw)
230                          A_m(IJK, 0, M) = -(HALF*hw + oX(I)*oDZ_T(K)*gw)
231     
232                          IF (BC_JJ_PS(L) .EQ. 1) THEN
233                             b_m(IJK, M) = -cw
234                          ELSE
235                             b_m(IJK,M) = ZERO
236                          ENDIF
237                        ENDIF
238     
239                     ELSEIF(FLUID_AT(BOTTOM_OF(IJK)))THEN
240                        IF(EP_s(BOTTOM_OF(IJK),M).LE.DIL_EP_s)THEN
241                           A_m(IJK, b , M) = ONE
242                        ELSE
243                          IF (BC_JJ_PS(L).EQ.3) THEN
244                             Gw = 1d0
245                             Hw = 0d0
246                             Cw = 0d0
247                          ELSE
248                             CALL CALC_THETA_BC(IJK,BOTTOM_OF(IJK),'B',&
249                                M,L,gw,hw,cw)
250                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
251                          ENDIF
252     
253                          A_m(IJK, b, M) = -(HALF*hw - oX(I)*oDZ_T(KM)*gw)
254                          A_m(IJK, 0, M) = -(HALF*hw + oX(I)*oDZ_T(KM)*gw)
255     
256                          IF (BC_JJ_PS(L) .EQ. 1) THEN
257                             b_m(IJK, M) = -cw
258                          ELSE
259                             b_m(IJK,M) = ZERO
260                          ENDIF
261                        ENDIF
262                     ENDIF
263     
264                   ENDDO   ! end do over i
265                   ENDDO   ! end do over j
266                   ENDDO   ! end do over k
267     
268                 ENDIF   ! end if (bc_jj_ps>0)
269               ENDIF  ! end if nsw, fsw or psw
270             ENDIF   ! end if (bc_defined)
271           ENDDO   ! end L do loop over dimension_bc
272     
273           RETURN
274           END SUBROUTINE BC_THETA
275     
276     
277     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
278     !                                                                      C
279     !  Subroutine: CALC_THETA_BC                                           C
280     !                                                                      C
281     !  Purpose: Implementation of Johnson & Jackson boundary conditions    C
282     !  for the pseudo-thermal temperature.                                 C
283     !                                                                      C
284     !                                                                      C
285     !  Author: Kapil Agrawal, Princeton University        Date: 14-MAR-98  C
286     !  Reviewer:                                          Date:            C
287     !                                                                      C
288     !                                                                      C
289     !                                                                      C
290     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
291     
292           SUBROUTINE CALC_THETA_BC(IJK1,IJK2,FCELL,M,L,Gw,Hw,Cw)
293     
294     !-----------------------------------------------
295     ! Modules
296     !-----------------------------------------------
297           USE bc
298           USE compar
299           USE constant
300           USE fldvar
301           USE fun_avg
302           USE functions
303           USE geometry
304           USE indices
305           USE mpi_utility
306           USE param
307           USE param1
308           USE physprop
309           USE rdf
310           USE run
311           USE rxns
312           USE toleranc
313           USE turb
314           USE visc_s
315           IMPLICIT NONE
316     !-----------------------------------------------
317     ! Dummy Arguments
318     !-----------------------------------------------
319     ! IJK indices for wall cell (ijk1) and fluid cell (ijk2)
320           INTEGER, INTENT(IN) :: IJK1, IJK2
321     ! The location (e,w,n...) of fluid cell
322           CHARACTER, INTENT(IN) :: FCELL
323     ! Solids phase index
324           INTEGER, INTENT(IN) :: M
325     ! Index corresponding to boundary condition
326           INTEGER, INTENT(IN) ::  L
327     ! Granular energy coefficient
328           DOUBLE PRECISION, INTENT(INOUT) :: Gw, Hw, Cw
329     !-----------------------------------------------
330     ! Local Variables
331     !-----------------------------------------------
332     ! IJK indices for cells
333           INTEGER ::  IJK3
334     ! Solids phase index
335           INTEGER :: MM
336     ! Average scalars
337           DOUBLE PRECISION :: EPg_avg, Mu_g_avg, RO_g_avg, smallTheta
338     ! void fraction at packing
339           DOUBLE PRECISION :: ep_star_avg
340     ! Average scalars modified to include all solid phases
341           DOUBLE PRECISION :: EPs_avg(DIMENSION_M), DP_avg(DIMENSION_M), &
342                               TH_avg(DIMENSION_M)
343     ! Average density of solids phase M
344           DOUBLE PRECISION ROs_avg(DIMENSION_M)
345     ! Average Simonin variables
346           DOUBLE PRECISION :: K_12_avg, Tau_12_avg
347     ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
348           DOUBLE PRECISION :: USCM, VSCM, WSCM
349     ! values of U_g, V_g, W_g at appropriate place on boundary wall
350           DOUBLE PRECISION :: UGC, VGC, WGC
351     ! Magnitude of gas-solids relative velocity
352           DOUBLE PRECISION :: VREL
353     ! Square of slip velocity between wall and particles
354           DOUBLE PRECISION :: VSLIPSQ
355     ! Wall velocity dot outward normal for all solids phases
356           DOUBLE PRECISION :: VWDOTN (DIMENSION_M)
357     ! Gradient in number density at the wall dot with outward normal for
358     ! all solids phases
359           DOUBLE PRECISION :: GNUWDOTN (DIMENSION_M)
360     ! Gradient in granular temperature at the wall dot with outward normal
361     ! for all solids phases
362           DOUBLE PRECISION :: GTWDOTN (DIMENSION_M)
363     ! Error message
364           CHARACTER(LEN=80)     LINE(1)
365     ! Radial distribution function
366           DOUBLE PRECISION :: g0(DIMENSION_M)
367     ! Sum of eps*G_0
368           DOUBLE PRECISION :: g0EPs_avg
369     !-----------------------------------------------
370     ! External functions
371     !-----------------------------------------------
372     ! Variable specularity coefficient
373           DOUBLE PRECISION :: PHIP_JJ
374     !-----------------------------------------------
375     
376     ! Note: EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).
377     !       Hence IJK2 (fluid cell) is used in averages.
378     
379           smallTheta = (to_SI)**4 * ZERO_EP_S
380     
381     ! set location independent quantities
382           EPg_avg = EP_g(IJK2)
383           ep_star_avg = EP_star_array(IJK2)
384           Mu_g_avg = Mu_g(IJK2)
385           RO_g_avg = RO_g(IJK2)
386           g0EPs_avg = ZERO
387     
388     ! added for Simonin model (sof)
389           IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
390              KT_TYPE_ENUM == AHMADI_1995) THEN
391              K_12_avg = K_12(IJK2)
392              Tau_12_avg = Tau_12(IJK2)
393           ELSE
394              K_12_avg = ZERO
395              Tau_12_avg = ZERO
396           ENDIF
397     
398           DO MM = 1, SMAX
399              g0(MM)      = G_0(IJK2, M, MM)
400              EPs_avg(MM) = EP_s(IJK2,MM)
401              DP_avg(MM)  = D_P(IJK2,MM)
402              g0EPs_avg   = g0EPs_avg + G_0(IJK2, M, MM)*EP_s(IJK2,MM)
403              ROs_avg(MM) = RO_S(IJK2,MM)
404           ENDDO
405     
406     
407           IF(FCELL .EQ. 'N')THEN
408              DO MM = 1, SMAX
409                 TH_avg(MM) = AVG_Y(Theta_m(IJK1,MM),Theta_m(IJK2,MM),J_OF(IJK1))
410                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
411     
412     ! added for IA (2005) theory:
413     ! include -1 since normal vector is pointing south (-y)
414     ! velocity at wall (i,j+1/2,k relative to ijk1) dot with the normal
415     ! vector at the wall
416                 VWDOTN(MM) = -1.d0*V_S(IJK1,MM)
417     
418     ! gradient in number density at wall (i,j+1/2,k relative to ijk1) dot
419     ! with the normal vector at the wall. since nu is undefined at ijk1,
420     ! approximate gradient with interior points: ijk2 and i,j+1,k relative
421     ! to ijk2
422                 IJK3 = NORTH_OF(IJK2)
423                 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
424                      (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oDY_N(J_OF(IJK2))
425     
426     ! gradient in granular temperature at wall (i,j+1/2,k relative to ijk1) dot
427     ! with the normal vector at the wall.
428                 GTWDOTN(MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
429                      oDY_N(J_OF(IJK2))
430              ENDDO
431     
432     ! Calculate velocity components at i, j+1/2, k (relative to IJK1)
433              UGC  = AVG_Y(AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
434                          AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
435                          J_OF(IJK1))
436              VGC  = V_g(IJK1)
437              WGC  = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
438                          AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
439                          J_OF(IJK1))
440              USCM = AVG_Y(AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
441                          AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
442                          J_OF(IJK1))
443              VSCM = V_s(IJK1,M)
444              WSCM = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
445                          AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
446                          J_OF(IJK1))
447     
448     ! slip velocity at the wall
449              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
450     
451     
452           ELSEIF(FCELL .EQ. 'S')THEN
453              DO MM = 1, SMAX
454                 TH_avg(MM) = AVG_Y(Theta_m(IJK2, MM),Theta_m(IJK1, MM),J_OF(IJK2))
455                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
456     
457     ! added for IA (2005) theory:
458     ! include 1 since normal vector is pointing north (+y)
459     ! velocity at wall (i,j+1/2,k relative to ijk2) dot with the normal
460     ! vector at the wall.
461                 VWDOTN(MM) = 1.d0*V_S(IJK2,MM)
462     
463     ! gradient in number density at wall (i,j+1/2,k relative to ijk2) dot
464     ! with the normal vector at the wall. since nu is undefined at ijk1,
465     ! approximate gradient with interior points: ijk2 and i,j-1,k relative
466     ! to ijk2
467                 IJK3 = SOUTH_OF(IJK2)
468                 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
469                      (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oDY_N(J_OF(IJK3))
470     
471     ! gradient in granular temperature at wall (i,j+1/2,k relative to ijk2)
472     ! dot with the normal vector at the wall.
473                 GTWDOTN(MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
474                      oDY_N(J_OF(IJK3))
475              ENDDO
476     
477     ! Calculate velocity components at i, j+1/2, k (relative to IJK2)
478              UGC  = AVG_Y(AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
479                          AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
480                          J_OF(IJK2))
481              VGC  = V_g(IJK2)
482              WGC  = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
483                          AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
484                          J_OF(IJK2))
485              USCM = AVG_Y(AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
486                          AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
487                          J_OF(IJK2))
488              VSCM = V_s(IJK2,M)
489              WSCM = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
490                          AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
491                          J_OF(IJK2))
492     
493     ! slip velocity at the wall
494              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
495     
496     
497           ELSEIF(FCELL== 'E')THEN
498               DO MM = 1, SMAX
499                 TH_avg(MM) = AVG_X(Theta_m(IJK1,MM),Theta_m(IJK2,MM),I_OF(IJK1))
500                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
501     
502     ! added for IA (2005) theory:
503     ! include -1 since normal vector is pointing west (-x)
504     ! velocity at wall (i+1/2,j,k relative to ijk1) dot with the normal
505     ! vector at the wall.
506                 VWDOTN(MM) = -1.d0*U_S(IJK1,MM)
507     
508     ! gradient in number density at wall (i+1/2,j,k relative to ijk1) dot
509     ! with the normal vector at the wall. since nu is undefined at ijk1,
510     ! approximate gradient with interior points: ijk2 and i,i+1,k relative
511     ! to ijk2
512                 IJK3 = EAST_OF(IJK2)
513                 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
514                      (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oDX_E(I_OF(IJK2))
515     
516     ! gradient in granular temperature at wall (i+1/2,j,k relative to ijk1)
517     ! dot with the normal vector at the wall.
518                 GTWDOTN(MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
519                      oDX_E(I_OF(IJK2))
520              ENDDO
521     
522     ! Calculate velocity components at i+1/2, j, k (relative to IJK1)
523              UGC  = U_g(IJK1)
524              VGC  = AVG_X(AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
525                          AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
526                          I_OF(IJK1))
527              WGC  = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
528                          AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
529                          I_OF(IJK1))
530              USCM = U_s(IJK1,M)
531              VSCM = AVG_X(AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
532                          AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
533                          I_OF(IJK1))
534              WSCM = AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
535                          AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
536                          I_OF(IJK1))
537     
538     ! slip velocity at the wall
539              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (VSCM-BC_Vw_s(L, M))**2
540     
541     
542           ELSEIF(FCELL== 'W')THEN
543               DO MM = 1, SMAX
544                 TH_avg(MM) = AVG_X(Theta_m(IJK2,MM),Theta_m(IJK1,MM),I_OF(IJK2))
545                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
546     
547     ! added for IA (2005) theory:
548     ! include 1 since normal vector is pointing west (+x)
549     ! velocity at wall (i+1/2,j,k relative to ijk2) dot with the normal
550     ! vector at the wall.
551                 VWDOTN(MM) = 1.d0*U_S(IJK2,MM)
552     
553     ! gradient in number density at wall (i+1/2,j,k relative to ijk2) dot
554     ! with the normal vector at the wall. since nu is undefined at ijk1,
555     ! approximate gradient with interior points: ijk2 and i,i-1,k relative
556     ! to ijk2
557                 IJK3 = WEST_OF(IJK2)
558                 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
559                      (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oDX_E(I_OF(IJK3))
560     
561     ! gradient in granular temperature at wall (i+1/2,j,k relative to ijk2)
562     ! dot with the normal vector at the wall.
563                 GTWDOTN(MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
564                      oDX_E(I_OF(IJK3))
565              ENDDO
566     
567     ! Calculate velocity components at i+1/2, j, k (relative to IJK2)
568              UGC  = U_g(IJK2)
569              VGC  = AVG_X(AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
570                          AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
571                          I_OF(IJK2))
572              WGC  = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
573                          AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
574                          I_OF(IJK2))
575              USCM = U_s(IJK2,M)
576              VSCM = AVG_X(AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
577                          AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
578                          I_OF(IJK2))
579              WSCM = AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
580                          AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
581                          I_OF(IJK2))
582     
583     ! slip velocity at the wall
584              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (VSCM-BC_Vw_s(L, M))**2
585     
586     
587           ELSEIF(FCELL== 'T')THEN
588              DO MM = 1, SMAX
589                 TH_avg(MM) = AVG_Z(Theta_m(IJK1,MM),Theta_m(IJK2,MM),K_OF(IJK1))
590                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
591     
592     ! added for IA (2005) theory:
593     ! include -1 since normal vector is pointing in bottom dir (-z)
594     ! velocity at wall (i,j,k+1/2 relative to ijk1) dot with the normal
595     ! vector at the wall.
596                 VWDOTN(MM) = -1.d0*W_S(IJK1,MM)
597     
598     ! gradient in number density at wall (i,j,k+1/2 relative to ijk1) dot
599     ! with the normal vector at the wall. since nu is undefined at ijk1,
600     ! approximate gradient with interior points: ijk2 and i,j,k+1 relative
601     ! to ijk2
602                 IJK3 = TOP_OF(IJK2)
603                 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
604                      (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oX(I_of(IJK2))*oDZ_T(K_OF(IJK2))
605     
606     ! gradient in granular temperature at wall (i,j,k+1/2 relative to ijk1)
607     ! dot with the normal vector at the wall.
608                 GNUWDOTN(MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
609                      oX(I_of(IJK2))*oDZ_T(K_OF(IJK2))
610              ENDDO
611     
612     ! Calculate velocity components at i, j, k+1/2 (relative to IJK1)
613              UGC  = AVG_Z(AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
614                          AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
615                          K_OF(IJK1))
616              VGC  = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
617                          AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
618                          K_OF(IJK1))
619              WGC  = W_g(IJK1)
620              USCM = AVG_Z(AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
621                          AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
622                          K_OF(IJK1))
623              VSCM = AVG_Z(AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
624                          AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
625                          K_OF(IJK1))
626              WSCM = W_s(IJK1,M)
627     
628     ! slip velocity at the wall
629              VSLIPSQ=(VSCM-BC_Vw_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
630     
631     
632           ELSEIF(FCELL== 'B')THEN
633              DO MM = 1, SMAX
634                 TH_avg(MM) = AVG_Z(Theta_m(IJK2,MM),Theta_m(IJK1,MM),K_OF(IJK2))
635                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
636     
637     ! added for IA (2005) theory:
638     ! include 1 since normal vector is pointing in bottom dir (+z)
639     ! velocity at wall (i,j,k+1/2 relative to ijk2) dot with the normal
640     ! vector at the wall
641                 VWDOTN(MM) = 1.d0*W_S(IJK2,MM)
642     
643     ! gradient in number density at wall (i,j,k+1/2 relative to ijk2) dot
644     ! with the normal vector at the wall. since nu is undefined at ijk1,
645     ! approximate gradient with interior points: ijk2 and i,j,k-1 relative
646     ! to ijk2
647                 IJK3 = BOTTOM_OF(IJK2)
648                 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
649                      (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oX(I_of(IJK3))*oDZ_T(K_OF(IJK3))
650     
651     ! gradient in granular temperature at wall (i,j,k+1/2 relative to ijk2)
652     ! dot with the normal vector at the wall.
653                 GTWDOTN(MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
654                      oX(I_of(IJK3))*oDZ_T(K_OF(IJK3))
655              ENDDO
656     
657     ! Calculate velocity components at i, j, k+1/2 (relative to IJK2)
658              UGC  = AVG_Z(AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
659                          AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
660                          K_OF(IJK2))
661              VGC  = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
662                          AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
663                          K_OF(IJK2))
664              WGC  = W_g(IJK2)
665              USCM = AVG_Z(AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
666                          AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
667                          K_OF(IJK2))
668              VSCM = AVG_Z(AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
669                          AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
670                          K_OF(IJK2))
671              WSCM = W_s(IJK2,M)
672     
673     ! slip velocity at the wall
674              VSLIPSQ=(VSCM-BC_Vw_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
675     
676           ELSE
677              WRITE(LINE,'(A, A)') 'Error: Unknown FCELL'
678              CALL WRITE_ERROR('CALC_THETA_BC', LINE, 1)
679              call exitMPI(myPE)
680           ENDIF
681     
682     ! magnitude of gas-solids relative velocity
683           VREL = DSQRT( (UGC-USCM)**2 + (VGC-VSCM)**2 + (WGC-WSCM)**2 )
684     
685           CALL THETA_Hw_Cw(g0, EPs_avg, EPg_avg, ep_star_avg, &
686                            g0EPs_avg, TH_avg,Mu_g_avg,RO_g_avg, ROs_avg, &
687                            DP_avg, K_12_avg,Tau_12_avg,VREL,VSLIPSQ,VWDOTN,&
688                            GNUWDOTN,GTWDOTN,M,Gw,Hw,Cw,L)
689     
690     ! Output the variable specularity coefficient.  Currently only works
691     ! for one solids phase
692           IF(BC_JJ_M .and. PHIP_OUT_JJ) THEN
693              IF(PHIP_OUT_ITER.eq.1) THEN
694                 ReactionRates(IJK1,1)= PHIP_JJ(dsqrt(VSLIPSQ),TH_avg(1))
695              ENDIF
696           ENDIF
697     
698           RETURN
699           END SUBROUTINE CALC_THETA_BC
700     
701     
702     
703     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
704     !                                                                      C
705     !  Subroutine: SUBROUTINE THETA_HW_CW                                  C
706     !  Purpose: Subroutine for gw, hw and cw                               C
707     !                                                                      C
708     !  Author: Kapil Agrawal, Princeton University         Date: 15-MAR-98 C
709     !  Reviewer:                                           Date:           C
710     !                                                                      C
711     !                                                                      C
712     !  Modified: Sofiane Benyahia, Fluent Inc.             Date: 02-FEB-05 C
713     !  Purpose: Include conductivity defined by Simonin and Ahmadi         C
714     !           Also included Jenkins small frictional limit               C
715     !                                                                      C
716     !  Literature/Document References:                                     C
717     !     See calc_mu_s.f for references on kinetic theory models          C
718     !     Johnson, P. C., and Jackson, R., Frictional-collisional          C
719     !        constitutive relations for granular materials, with           C
720     !        application to plane shearing, Journal of Fluid Mechanics,    C
721     !        1987, 176, pp. 67-93                                          C
722     !     Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating     C
723     !        energy in a collisional grain flow at a flat frictional wall, C
724     !        Physics of Fluids, 1997, 9(10), pp. 2835-2840                 C
725     !                                                                      C
726     !                                                                      C
727     !  Additional Notes:                                                   C
728     !     The JJ granular energy BC is written as: n.q = D-G, where        C
729     !       q = the heat flux, n = outward normal vector                   C
730     !       D = dissipation due to p-w collisions, and                     C
731     !       G = generation due to p-w slip                                 C
732     !     In MFIX this is implemented as k.gradT = G-D, where              C
733     !       k = granular conductivity units (Mass/Length/Time)             C
734     !       T = granular temperature in units (Length/Time)^2              C
735     !     However, the expression for heat flux (q) may contain terms      C
736     !     beyond that of gradient in temperature, such as the dufour       C
737     !     term. A more rigorous implemenation should account for these     C
738     !     additional terms, which is not done here.                        C
739     !                                                                      C
740     !                                                                      C
741     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
742     
743           SUBROUTINE THETA_HW_CW(g0,EPS,EPG, ep_star_avg, &
744                                  g0EPs_avg,TH,Mu_g_avg,RO_g_avg, ROs_avg, DP_avg, &
745                                  K_12_avg,Tau_12_avg,VREL,VSLIPSQ,VWDOTN,&
746                                  GNUWDOTN,GTWDOTN,M,GW,HW,CW,L)
747     
748     !-----------------------------------------------
749     ! Modules
750     !-----------------------------------------------
751           USE param
752           USE param1
753           USE physprop
754           USE run
755           USE constant
756           USE fldvar
757           USE toleranc
758           USE bc
759           USE compar
760           IMPLICIT NONE
761     !-----------------------------------------------
762     ! Dummy Arguments
763     !-----------------------------------------------
764     ! Radial distribution function of solids phase M with each
765     ! other solids phase
766           DOUBLE PRECISION, INTENT(IN) :: g0(DIMENSION_M)
767     ! Average solids volume fraction of each solids phase
768           DOUBLE PRECISION, INTENT(IN) :: EPS(DIMENSION_M)
769     ! Average solids and gas volume fraction
770           DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
771     ! Sum of eps*G_0
772           DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
773     ! Average theta_m
774           DOUBLE PRECISION, INTENT(INOUT) :: TH (DIMENSION_M)
775     ! Average gas viscosity
776           DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
777     ! Average gas density
778           DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
779     ! Average density of solids phase M
780           DOUBLE PRECISION, INTENT(IN) :: ROS_avg(DIMENSION_M)
781     ! Average particle diameter of each solids phase
782           DOUBLE PRECISION, INTENT(IN) :: DP_avg(DIMENSION_M)
783     ! Average cross-correlation K_12 and lagrangian integral time-scale
784           DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg
785     ! Magnitude of slip velocity between two phases
786           DOUBLE PRECISION, INTENT(IN) :: VREL
787     ! Square of slip velocity between wall and particles
788           DOUBLE PRECISION, INTENT(IN) :: VSLIPSQ
789     ! Wall velocity dot outward normal for all solids phases
790           DOUBLE PRECISION, INTENT(IN) :: VWDOTN (DIMENSION_M)
791     ! Gradient in number density at the wall dot with outward
792     ! normal for all solids phases
793           DOUBLE PRECISION, INTENT(IN) :: GNUWDOTN (DIMENSION_M)
794     ! Gradient in temperature at the wall dot with outward
795     ! normal for all solids phases
796           DOUBLE PRECISION, INTENT(IN) :: GTWDOTN (DIMENSION_M)
797     ! Solids phase index
798           INTEGER, INTENT(IN) :: M
799     ! Coefficients in JJ boundary conditions
800     ! Gw = granular conductivity (w particle mass)
801     ! Hw = dissipation due to particle-wall collision
802     ! Cw = generation due to particle-wall slip
803           DOUBLE PRECISION, INTENT(INOUT) :: GW, HW, CW
804     ! Index corresponding to boundary condition
805           INTEGER, INTENT(IN) :: L
806     !-----------------------------------------------
807     ! Local Variables
808     !-----------------------------------------------
809     ! Solids phase index
810           INTEGER :: LL
811     ! Coefficient of 1st term
812           DOUBLE PRECISION :: K_1
813     ! Conductivity
814           DOUBLE PRECISION :: Kgran
815     ! Conductivity corrected for interstitial fluid effects
816           DOUBLE PRECISION :: Kgran_star
817     ! Drag Coefficient
818           DOUBLE PRECISION :: Beta, DgA
819     ! Reynolds number based on slip velocity
820           DOUBLE PRECISION :: Re_g
821     ! Friction Factor in drag coefficient
822           DOUBLE PRECISION :: C_d
823     ! mth solids phase particle diameter, mass and number density
824           DOUBLE PRECISION :: D_PM, M_PM, NU_PM
825     ! Variables for Iddir & Arastoopour model
826           DOUBLE PRECISION :: D_PL, M_PL, NU_PL, MPSUM, DPSUMo2
827           DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R1p_lm, R8p_lm, R9p_lm, &
828                               Bp_lm, R5p_lm, R6p_lm, R7p_lm
829           DOUBLE PRECISION :: K_s_sum, K_s_MM, K_sM_LM
830           DOUBLE PRECISION :: Kvel_s_sum, Kth_s_sum
831           DOUBLE PRECISION :: Knu_s_sum, K_common_term
832     ! Variables for Garzo & Dufty model
833           DOUBLE PRECISION :: c_star, zeta0_star, press_star, eta0, &
834                               kappa0, nu_kappa_star, kappa_k_star, &
835                               kappa_star
836     ! Variables for GTSH theory
837           DOUBLE PRECISION :: Re_T, Chi, vfrac, mu2_0, mu4_0, mu4_1, &
838                               A2_gtshW, zeta_star, nu0, NuK, Kth0, KthK, EDT_s
839     ! Local parameters for Simonin model
840           DOUBLE PRECISION :: Zeta_c, Omega_c, Tau_2_c, Kappa_kin, &
841                               Kappa_Col, Tau_12_st
842     ! Solids pressure and pressure divided by granular temperature
843           DOUBLE PRECISION :: Ps, PsoTheta
844     !-----------------------------------------------
845     ! Function subroutines
846     !-----------------------------------------------
847     ! variable specularity coefficient
848           DOUBLE PRECISION :: PHIP_JJ
849           DOUBLE PRECISION :: S_star
850     !-----------------------------------------------
851     
852           IF(TH(M) .LE. ZERO)THEN
853              TH(M) = 1D-8
854              IF (myPE.eq.PE_IO) THEN
855                 WRITE(*,*)  &
856                    'Warning: Negative granular temp at wall set to 1e-8'
857     !            CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
858              ENDIF
859           ENDIF
860     
861     ! common variables
862           D_PM = DP_avg(M)
863           M_PM = (PI/6.d0)*(D_PM**3.)*ROS_avg(M)
864           NU_PM = (EPS(M)*ROS_avg(M))/M_PM
865     
866     ! This is from Wen-Yu correlation, you can put here your own single particle drag
867           Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
868           IF (Re_g .lt. 1000.d0) THEN
869              C_d = (24.d0/(Re_g+SMALL_NUMBER))*(ONE + 0.15d0 * Re_g**0.687d0)
870           ELSE
871              C_d = 0.44d0
872           ENDIF
873           DgA = 0.75d0*C_d*Ro_g_avg*EPG*VREL/(DP_avg(M)*EPG**(2.65d0))
874           IF(VREL == ZERO) DgA = LARGE_NUMBER
875           Beta = EPS(M)*DgA
876     
877     
878     ! Defining conductivity according to each KT for later use in JJ BC
879     ! Also define solids pressure for use in Jenkins small frictional BC
880     ! -------------------------------------------------------------------
881           SELECT CASE (KT_TYPE_ENUM)
882     
883           CASE (LUN_1984)
884              Kgran = 75d0*ROs_avg(M)*DP_avg(M)*DSQRT(Pi*TH(M))/&
885                 (48*Eta*(41d0-33d0*Eta))
886     
887              IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
888                 Kgran_star = Kgran
889              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
890                 Kgran_star = ZERO
891              ELSE
892                 Kgran_star = ROs_avg(M)*EPS(M)* g0(M)*TH(M)* Kgran/ &
893                    (ROs_avg(M)*g0EPs_avg*TH(M) + &
894                    1.2d0*DgA/ROs_avg(M)* Kgran)
895              ENDIF
896     
897     ! mth solids phase granular conductivity includes particle mass
898              K_1 = Kgran_star/g0(M)*( &
899                 ( ONE + (12d0/5.d0)*Eta*g0EPs_avg ) * &
900                 ( ONE + (12d0/5.d0)*Eta*Eta*(4d0*Eta-3d0) *g0EPs_avg ) + &
901                 (64d0/(25d0*Pi)) * (41d0-33d0*Eta) *(Eta*g0EPs_avg)**2 )
902     
903     ! defining granular pressure (for Jenkins BC)
904              PsoTheta = ROs_avg(M)*EPS(M)*(ONE+4.D0*Eta*g0EPs_avg)
905              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
906     
907     
908           CASE (SIMONIN_1996)
909     ! particle relaxation time
910              Tau_12_st = ROs_avg(M)/(DgA+small_number)
911              Zeta_c  = (ONE+ C_e)*(49.d0-33.d0*C_e)/100.d0
912              Omega_c = 3.d0*(ONE+ C_e)**2 *(2.d0*C_e-ONE)/5.d0
913              Tau_2_c = DP_avg(M)/(6.d0*EPS(M)*g0(M) &
914                 *DSQRT(16.d0*(TH(M)+Small_number)/PI))
915     
916     ! Defining Simonin's Solids Turbulent Kinetic diffusivity: Kappa
917              Kappa_kin = (9.d0/10.d0*K_12_avg*(Tau_12_avg/Tau_12_st) &
918                 + 3.0D0/2.0D0 * TH(M)*(ONE+ Omega_c*EPS(M)*g0(M)))/     &
919                 (9.d0/(5.d0*Tau_12_st) + zeta_c/Tau_2_c)
920     
921              Kappa_Col = 18.d0/5.d0*EPS(M)*g0(M)*Eta* (Kappa_kin+ &
922                 5.d0/9.d0*DP_avg(M)*DSQRT(TH(M)/PI))
923     
924     ! mth solids phase granular conductivity includes particle mass
925              K_1 =  EPS(M)*ROs_avg(M)*(Kappa_kin + Kappa_Col)
926     
927     ! defining granular pressure (for Jenkins BC)
928              PsoTheta = ROs_avg(M)*EPS(M)*(ONE+4.D0*Eta*g0EPs_avg)
929              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
930     
931     
932           CASE (AHMADI_1995)
933     ! mth solids phase granular conductivity includes particle mass
934              K_1 =  0.1306D0*ROs_avg(M)*DP_avg(M)*(ONE+C_e**2)* (  &
935                 ONE/g0(M)+4.8D0*EPS(M)+12.1184D0 *EPS(M)*EPS(M)*g0(M) )*&
936                 DSQRT(TH(M))
937     
938     ! defining granular pressure (for Jenkins BC)
939              PsoTheta = ROs_avg(M)*EPS(M)* &
940                  ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
941              Ps = ROs_avg(M)*EPS(M)*TH(M)*&
942                  ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
943     
944     
945           CASE (IA_2005)
946     ! Use original IA theory if SWITCH_IA is false
947              IF(.NOT. SWITCH_IA) g0EPs_avg = EPS(M)*ROS_avg(M)
948     
949              K_s_sum = ZERO
950              Kvel_s_sum = ZERO
951              Knu_s_sum = ZERO
952              Kth_s_sum = ZERO
953     
954              Kgran = (75.d0/384.d0)*ROs_avg(M)*D_PM*DSQRT(Pi*TH(M)/M_PM)
955     
956              IF(SWITCH == ZERO .OR. RO_g_avg == ZERO)THEN
957                 Kgran_star = Kgran
958              ELSEIF(TH(M)/M_PM .LT. SMALL_NUMBER)THEN
959                 Kgran_star = ZERO
960              ELSE
961                 Kgran_star = Kgran*g0(M)*EPS(M)/ &
962                    (g0EPs_avg+ 1.2d0*DgA*Kgran / (ROS_avg(M)**2 *(TH(M)/M_PM)))
963              ENDIF
964     
965              K_s_MM = (Kgran_star/(M_PM*g0(M)))*&  ! Kth doesn't include the mass.
966                 (1.d0+(3.d0/5.d0)*(1.d0+C_E)*(1.d0+C_E)*g0EPs_avg)**2
967     
968              DO LL = 1, SMAX
969                 D_PL = DP_avg(LL)
970                 M_PL = (PI/6.d0)*(D_PL**3.)*ROS_avg(LL)
971                 MPSUM = M_PM + M_PL
972                 DPSUMo2 = (D_PM+D_PL)/2.d0
973                 NU_PL = (EPS(LL)*ROS_avg(LL))/M_PL
974     
975                 IF ( LL .eq. M) THEN
976                    K_s_sum = K_s_sum + K_s_MM
977     ! Kth_sL_LM is zero when LL=M because it cancels with terms from K_s_LM
978     ! Kvel_s_LM should be zero when LL=M (same solids phase)
979     ! Knu_s_LM should be zero when LL=M (same solids phase)
980                 ELSE
981                    Ap_lm = (M_PM*TH(LL)+M_PL*TH(M))/(2.d0)
982                    Bp_lm = (M_PM*M_PL*(TH(LL)-TH(M) ))/(2.d0*MPSUM)
983                    Dp_lm = (M_PL*M_PM*(M_PM*TH(M)+M_PL*TH(LL) ))/&
984                       (2.d0*MPSUM*MPSUM)
985                    R0p_lm = ( 1.d0/( Ap_lm**1.5 * Dp_lm**2.5 ) )+ &
986                       ( (15.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**2.5 * Dp_lm**3.5 ) )+&
987                       ( (175.d0*(Bp_lm**4))/( 8.d0*Ap_lm**3.5 * Dp_lm**4.5 ) )
988                    R1p_lm = ( 1.d0/( (Ap_lm**1.5)*(Dp_lm**3) ) )+ &
989                       ( (9.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**4 ) )+&
990                       ( (30.d0*Bp_lm**4) /( 2.d0*Ap_lm**3.5 * Dp_lm**5 ) )
991                    R5p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**3 ) )+ &
992                       ( (5.d0*Bp_lm*Bp_lm)/( Ap_lm**3.5 * Dp_lm**4 ) )+&
993                       ( (14.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5 ) )
994                    R6p_lm = ( 1.d0/( Ap_lm**3.5 * Dp_lm**3 ) )+ &
995                       ( (7.d0*Bp_lm*Bp_lm)/( Ap_lm**4.5 * Dp_lm**4 ) )+&
996                       ( (126.d0*Bp_lm**4)/( 5.d0*Ap_lm**5.5 * Dp_lm**5 ) )
997                    R7p_lm = ( 3.d0/( 2.d0*Ap_lm**2.5 * Dp_lm**4 ) )+ &
998                       ( (10.d0*Bp_lm*Bp_lm)/( Ap_lm**3.5 * Dp_lm**5 ) )+&
999                       ( (35.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**6 ) )
1000                    R8p_lm = ( 1.d0/( 2.d0*Ap_lm**1.5 * Dp_lm**4 ) )+ &
1001                       ( (6.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**5 ) )+&
1002                       ( (25.d0*Bp_lm**4)/( Ap_lm**3.5 * Dp_lm**6 ) )
1003                    R9p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**3 ) )+ &
1004                       ( (15.d0*Bp_lm*Bp_lm)/( Ap_lm**3.5 * Dp_lm**4 ) )+&
1005                       ( (70.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5 ) )
1006                    K_common_term = DPSUMo2**3 * M_PL*M_PM/(2.d0*MPSUM)*&
1007                       (1.d0+C_E)*g0(LL) * (M_PM*M_PL)**1.5
1008     
1009     ! solids phase 'conductivity' associated with the
1010     ! gradient in granular temperature of species M
1011                    K_sM_LM = - K_common_term*NU_PM*NU_PL*(&
1012                       ((DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)+&
1013                       ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1014                       (3.d0/2.d0)*R1p_lm)-(&
1015                       ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM/8.d0)*Bp_lm*R6p_lm)+&
1016                       ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1017                       8.d0)*M_PM*R9p_lm)+&
1018                       ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1019                       M_PL*Bp_lm*R7p_lm)+&
1020                       ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1021                       2.d0)*(M_PM/MPSUM)**2 * M_PL*R8p_lm)+&
1022                       ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1023                       R9p_lm)-&
1024                       ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1025                       (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*TH(LL) )*&
1026                       (TH(M)**2 * TH(LL)**3)
1027     
1028     ! These lines were commented because they are not currently used
1029     ! solids phase 'conductivity' associated with the gradient in granular
1030     ! temperature of species L
1031     !               Kth_sL_LM = K_common_term*NU_PM*NU_PL*(&
1032     !                  (-(DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)+&
1033     !                  (-(M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1034     !                  (3.d0/2.d0)*R1p_lm)+(&
1035     !                  ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL/8.d0)*Bp_lm*R6p_lm)+&
1036     !                  ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1037     !                  8.d0)*M_PL*R9p_lm)+&
1038     !                  ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1039     !                  M_PM*Bp_lm*R7p_lm)+&
1040     !                  ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1041     !                  2.d0)*(M_PM*M_PM/(MPSUM*MPSUM))*M_PM*R8p_lm)+&
1042     !                  ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1043     !                  R9p_lm)+&
1044     !                  ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1045     !                  (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*TH(M) )*&
1046     !                  (TH(LL)*TH(LL)*(TH(M)**3.))*(GTWDOTN(LL))
1047     
1048     ! solids phase 'conductivity' associated with the difference in velocities
1049     !               Kvel_s_LM = K_common_term*NU_PM*NU_PL*&
1050     !                  (M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(3.d0*PI/10.d0)*&
1051     !                  R0p_lm*( (TH(M)*TH(LL))**(5./2.) )*&
1052     !                  (VWDOTN(M)-VWDOTN(LL))
1053     
1054     ! solids phase 'conductivity' associated with the difference in the
1055     ! gradient in number densities
1056     !               Knu_s_LM = K_common_term*(((DPSUMo2*DSQRT(PI)/16.d0)*&
1057     !                  Bp_lm*R5p_lm)+((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*&
1058     !                  (DPSUMo2*PI/6.d0)*R1p_lm) )*( (TH(M)*TH(LL))**(3.) )*&
1059     !                  (NU_PM*GNUWDOTN(LL)-NU_PL*GNUWDOTN(M))
1060     !
1061     
1062                    K_s_sum = K_s_sum + K_sM_LM
1063     !               Kth_s_sum = Kth_s_sum + Kth_sL_LM
1064     !               Kvel_s_sum = Kvel_s_sum + Kvel_s_LM
1065     !               Knu_s_sum = Knu_s_sum + Knu_s_LM
1066     
1067                 ENDIF   ! end if( LL .eq. M)/else
1068              ENDDO  ! end do LL = 1, SMAX
1069     
1070     ! mth solids phase granular conductivity DOES NOT include particle mass
1071              K_1 = K_s_sum
1072     
1073     ! Note that the velocity term is not included here because it should
1074     ! become zero when dotted with the outward normal (i.e. no solids
1075     ! flux through the wall; assuming that the solids velocity at the
1076     ! wall in the normal direction is zero)
1077     !         CW = CW + Kvel_s_sum
1078     
1079     ! Note that the gradient in number density at the wall must be
1080     ! approximated with interior points since there is no value of
1081     ! number density associated with the wall location; the ghost
1082     ! cell values are undefined for volume fraction.
1083     !         CW = CW + Knu_s_sum
1084     
1085     ! The gradient in temperature of phase LL at the wall
1086     !         CW = CW + Kth_s_sum
1087     
1088     
1089           CASE (GD_1999)
1090     ! Note: k_boltz = M_PM
1091              D_PM = DP_avg(M)
1092              M_PM = (PI/6.d0)*(D_PM**3.)*ROS_avg(M)
1093              NU_PM = (EPS(M)*ROS_avg(M))/M_PM
1094     
1095     ! Find pressure in the Mth solids phase
1096              press_star = ONE + 2.d0*(1.d0+C_E)*EPS(M)*G0(M)
1097     
1098     ! find conductivity
1099              eta0 = 5.0d0*M_PM*DSQRT(TH(M)/PI) / (16.d0*D_PM*D_PM)
1100     
1101              c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
1102                 / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
1103     
1104              zeta0_star = (5.d0/12.d0)*G0(M)*(1.d0 - C_E*C_E) &
1105                 * (1.d0 + (3.d0/32.d0)*c_star)
1106     
1107              kappa0 = (15.d0/4.d0)*eta0
1108     
1109              nu_kappa_star = (G0(M)/3.d0)*(1.d0+C_E) * ( 1.d0 + &
1110                 (33.d0/16.d0)*(1.d0-C_E) + ((19.d0-3.d0*C_E)/1024.d0)*c_star)
1111     !         nu_mu_star = nu_kappa_star
1112     
1113              kappa_k_star = (2.d0/3.d0)*(1.d0 + 0.5d0*(1.d0+press_star)*c_star + &
1114                 (3.d0/5.d0)*EPS(M)*G0(M)*(1.d0+C_E)*(1.d0+C_E) * &
1115                 (2.d0*C_E - 1.d0 + ( 0.5d0*(1.d0+C_E) - 5.d0/(3*(1.d0+C_E))) * &
1116                 c_star ) ) / (nu_kappa_star - 2.d0*zeta0_star)
1117     
1118              kappa_star = kappa_k_star * (1.d0 + (6.d0/5.d0)*EPS(M)* &
1119                 G0(M)*(1.d0+C_E) ) + (256.d0/25.d0)*(EPS(M)* &
1120                 EPS(M)/PI)*G0(M)*(1.d0+C_E)*(1.d0+(7.d0/32.d0)* &
1121                 c_star)
1122     
1123              IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
1124                 Kgran_star = kappa0
1125              ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
1126                 Kgran_star = ZERO
1127              ELSE
1128                 Kgran_star = ROS_avg(M)*EPS(M)* G0(M)*TH(M)* kappa0/ &
1129                    (ROS_avg(M)*G0(M)*EPS(M)*TH(M) + &
1130                    1.2d0*DgA*kappa0/ROs_avg(M))     ! Note dgA is ~F_gs/ep_s
1131              ENDIF
1132     
1133     ! kgran_star includes particle mass
1134     ! mth solids phase granular conductivity includes particle mass
1135              K_1 = Kgran_star * kappa_star
1136     
1137     ! transport coefficient of the Mth solids phase associated
1138     ! with gradient in volume fraction in heat flux
1139     !         qmu_k_star = 2.d0*( (1.d0+EPS(M)*DG_0DNU(EPS(M)))* &
1140     !            zeta0_star*kappa_k_star + ( (press_star/3.d0) + (2.d0/3.d0)* &
1141     !            EPS(M)*(1.d0+C_E) * (G0(M)+EPS(M)* &
1142     !            DG_0DNU(EPS(M))) )*c_star - (4.d0/5.d0)*EPS(M)* &
1143     !            G0(M)* (1.d0+(EPS(M)/2.d0)*DG_0DNU(EPS(M)))* &
1144     !            (1.d0+C_E) * ( C_E*(1.d0-C_E)+0.25d0*((4.d0/3.d0)+C_E* &
1145     !            (1.d0-C_E))*c_star ) ) / (2.d0*nu_kappa_star-3.d0*zeta0_star)
1146     !         qmu_star = qmu_k_star*(1.d0+(6.d0/5.d0)*EPS(M)*G0(M)*&
1147     !            (1.d0+C_E) )
1148     !         Kphis = (TH(M)*Kgran_star/NU_PM)*qmu_star
1149     
1150     ! defining granular pressure (for Jenkins BC)
1151              PsoTheta = ROs_avg(M)*EPS(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1152                         ! ~ROs_avg(m)*EPS(M)*press_star
1153              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1154                         ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1155     
1156           CASE (GTSH_2012)
1157              D_PM = DP_avg(M)
1158              M_PM = (PI/6.d0)*(D_PM**3.)*ROS_avg(M)
1159              NU_PM = (EPS(M)*ROS_avg(M))/M_PM
1160              Chi = g0(M)
1161              vfrac = EPS(M)
1162     
1163              Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
1164              Re_T = RO_g_avg*D_PM*dsqrt(TH(M)) / Mu_g_avg
1165              mu2_0 = dsqrt(2d0*pi) * Chi * (one-C_E**2)  ! eq. (6.22) GTSH theory
1166              mu4_0 = (4.5d0+C_E**2) * mu2_0              ! eq. (6.23) GTSH theory
1167              mu4_1 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + 2d0*dsqrt(2d0*pi)* &
1168                       Chi*(one+C_E)  ! this is done to avoid /0 in case c_e = 1.0
1169              A2_gtshW = zero ! for eps = zero
1170              if((vfrac> small_number) .and. (TH(M) > small_number)) then
1171                 zeta_star = 4.5d0*dsqrt(2d0*Pi)*(RO_g_avg/ROs_avg(M))**2*Re_g**2 * &
1172                             S_star(vfrac,Chi) / (vfrac*(one-vfrac)**2 * Re_T**4)
1173                 A2_gtshW = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
1174                                (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
1175              endif
1176              eta0 = 0.3125d0/(dsqrt(pi)*D_PM**2)*M_pm*dsqrt(TH(M))
1177              nu0 = (96.d0/5.d0)*(vfrac/D_PM)*DSQRT(TH(M)/PI)
1178              NuK = nu0*(one+C_E)/3d0*Chi*( one+2.0625d0*(one-C_E)+ &
1179                         ((947d0-579*C_E)/256d0*A2_gtshW) )
1180              Kth0 = 3.75d0*eta0/M_pm
1181              EDT_s = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
1182                     (one+0.1875d0*A2_gtshW)*NU_PM*D_PM**2*dsqrt(TH(M))
1183              KthK = zero
1184              if(vfrac > small_number) KthK = 2d0/3d0*Kth0*nu0/(NuK - 2d0*EDT_s) * &
1185                            (one+2d0*A2_gtshW+0.6d0*vfrac*Chi* &
1186                            (one+C_E)**2*(2*C_E-one+A2_gtshW*(one+C_E)))
1187     
1188     ! defining conductivity K_1 at walls equivalent to Kth_s(IJK,M) in calc_mu_s
1189     ! mth solids phase granular conductivity DOES NOT include particle mass
1190              K_1 = KthK*(one+1.2d0*vfrac*Chi*(one+C_E)) + (10.24d0/pi* &
1191                     vfrac**2*Chi*(one+C_E)*(one+0.4375d0*A2_gtshW)*Kth0)
1192     
1193     ! defining granular pressure (for Jenkins BC)
1194              PsoTheta = ROs_avg(M)*EPS(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1195              Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1196     
1197     
1198           CASE DEFAULT
1199     ! should never hit this
1200              WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1201              WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1202              call mfix_exit(myPE)
1203           END SELECT
1204     
1205     
1206     ! setting the coefficients for JJ BC
1207     ! -------------------------------------------------------------------
1208     ! GW = granular conductivity
1209     ! Hw = dissipation due to particle-wall collision
1210     ! Cw = generation due to particle-wall slip
1211           SELECT CASE (KT_TYPE_ENUM)
1212              CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999)
1213     ! KTs without particle mass in their definition of granular temperature
1214     ! and with mass in their conductivity
1215     ! and theta = granular temperature
1216                 GW = K_1
1217     
1218     ! Jenkins BC implemented for these KT models
1219                 IF(JENKINS) THEN
1220                    HW = 3.D0/8.D0*DSQRT(3.D0*TH(M))*((1d0-e_w))*&
1221                         PsoTheta
1222     
1223     ! As I understand from soil mechanic papers, the coefficient mu in
1224     ! Jenkins paper is tan_Phi_w. T.  See for example, G.I. Tardos, PT,
1225     ! 92 (1997), 61-74, equation (1). sof
1226                    CW = tan_Phi_w*tan_Phi_w*(ONE+e_w)*21.d0/16.d0*&
1227                         DSQRT(3.D0*TH(M)) * Ps
1228     
1229                ELSE
1230                    HW = (Pi*DSQRT(3d0)/(4.D0*(ONE-ep_star_avg)))*&
1231                       (1d0-e_w*e_w)*&
1232                       ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))
1233     
1234                    IF(.NOT. BC_JJ_M) THEN
1235                       CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*PHIP*&
1236                          ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1237                    ELSE
1238                       CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*&
1239                          PHIP_JJ(DSQRT(vslipsq),TH(M))*ROs_avg(M)*&
1240                          EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1241                    ENDIF
1242                 ENDIF   ! end if(Jenkins)/else
1243     
1244     
1245              CASE (GTSH_2012)
1246     ! KTs with particle mass in their definition of granular temperature
1247     ! and without mass in their conductivity
1248     ! and theta = granular temperature/mass
1249                 GW = M_PM * K_1
1250     
1251                 HW = (Pi*DSQRT(3d0)/(4.D0*(ONE-ep_star_avg)))*&
1252                    (1d0-e_w*e_w)*&
1253                    ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))
1254     
1255                 IF(.NOT. BC_JJ_M) THEN
1256                    CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*PHIP*&
1257                       ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1258                 ELSE
1259                    CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*&
1260                       PHIP_JJ(DSQRT(vslipsq),TH(M))*ROs_avg(M)*&
1261                       EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1262                 ENDIF
1263     
1264     
1265              CASE (IA_2005)
1266     ! KTs with particle mass in their definition of granular temperature
1267     ! and without mass in their conductivity
1268     ! and theta = granular temperature
1269                 GW = M_PM * K_1
1270     
1271                 HW = (PI*DSQRT(3.d0)/(4.d0*(ONE-ep_star_avg)))*&
1272                    (1.d0-e_w*e_w)*&
1273                    ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M)/M_PM)
1274     
1275                 IF(.NOT. BC_JJ_M) THEN
1276                    CW = (PI*DSQRT(3.d0)/(6.d0*(ONE-ep_star_avg)))*PHIP*&
1277                       ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M)*M_PM)*VSLIPSQ
1278                 ELSE
1279                    CW = (PI*DSQRT(3.d0)/(6.d0*(ONE-ep_star_avg)))*&
1280                       PHIP_JJ(DSQRT(vslipsq),TH(M))*ROs_avg(M)*&
1281                       EPS(M)*g0(M)*DSQRT(TH(M)*M_PM)*VSLIPSQ
1282                 ENDIF
1283     
1284     
1285     
1286              CASE DEFAULT
1287     ! should never hit this
1288                 WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1289                 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1290                 call mfix_exit(myPE)
1291              END SELECT
1292     
1293     
1294           RETURN
1295           END SUBROUTINE THETA_HW_CW
1296     
1297     
1298     
1299     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1300     !                                                                      C
1301     !  Function: PHIP_JJ                                                   C
1302     !  Purpose: Calculate the specularity coefficient for JJ BC            C
1303     !                                                                      C
1304     !  Author: Tingwen Li                                                  C
1305     !                                                                      C
1306     !                                                                      C
1307     !  Literature/Document References:                                     C
1308     !     Li, T., and Benyahia, S., Revisiting Johnson and Jackson         C
1309     !        boundary conditions for granular flows, AIChE Journal, 2012,  C
1310     !        Vol 58 (7), pp. 2058-2068                                     C
1311     !                                                                      C
1312     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1313     
1314     
1315           DOUBLE PRECISION FUNCTION phip_jj(uslip,g_theta)
1316     !-----------------------------------------------
1317     ! Modules
1318     !-----------------------------------------------
1319           use constant, only : e_w, PI, k4phi, phip0
1320           implicit none
1321     !-----------------------------------------------
1322     ! Dummy arguments
1323     !-----------------------------------------------
1324           double precision, INTENT(IN) :: uslip
1325           double precision, INTENT(IN) :: g_theta
1326     !-----------------------------------------------
1327     ! Local variables
1328     !-----------------------------------------------
1329           double precision :: r4phi
1330     !-----------------------------------------------
1331     
1332           !k4phi=7.d0/2.d0*mu4phi*(1.d0+e_w)
1333           r4phi=uslip/dsqrt(3.d0*g_theta)
1334     
1335           IF(r4phi .le. 4.d0*k4phi/7.d0/dsqrt(6.d0*PI)/phip0) THEN
1336              phip_jj=-7.d0*dsqrt(6.d0*PI)*phip0**2/8.d0/k4phi*r4phi+phip0
1337           ELSE
1338              phip_jj=2.d0/7.d0*k4phi/r4phi/dsqrt(6.d0*PI)
1339           ENDIF
1340     
1341           RETURN
1342     
1343           END FUNCTION PHIP_JJ
1344