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