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