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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: K_Epsilon_PROP(IER)                                    C
4     !  Purpose: Calculate diffusion coefficeint and sources for K &        C
5     !           Epsilon equations                                          C
6     !                                                                      C
7     !  Author:                                                    Date:    C
8     !  Reviewer:                                                           C
9     !  Modified: S. Benyahia                         Date:May-13-04        C
10     !                                                                      C
11     !                                                                      C
12     !  Literature/Document References: Wilcox, D.C., Turbulence Modeling   C
13     !  for CFD. DCW Industries, Inc. La Canada, Ca. 1994.                  C
14     !                                                                      C
15     !  Variables referenced: None                                          C
16     !  Variables modified: None                                            C
17     !                                                                      C
18     !  Local variables: None                                               C
19     !                                                                      C
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21     !
22           SUBROUTINE K_Epsilon_PROP( IER)
23     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
24     !...Switches: -xf
25     !
26     !-----------------------------------------------
27     !   M o d u l e s
28     !-----------------------------------------------
29           USE param
30           USE param1
31           USE parallel
32           USE physprop
33           USE drag
34           USE run
35           USE output
36           USE geometry
37           USE fldvar
38           USE visc_g
39           USE visc_s
40           USE trace
41           USE indices
42           USE constant
43           Use vshear
44           USE turb
45           USE toleranc
46           USE compar
47           USE TAU_G
48           USE sendrecv
49     
50           USE cutcell
51           USE fun_avg
52           USE functions
53     
54           IMPLICIT NONE
55     !-----------------------------------------------
56     !   G l o b a l   P a r a m e t e r s
57     !-----------------------------------------------
58             DIMENSION D_g(3,3)
59     !-----------------------------------------------
60     !   D u m m y   A r g u m e n t s
61     !-----------------------------------------------
62     !
63     !                       Error index
64           INTEGER          IER
65     !
66     !-----------------------------------------------
67     !-----------------------------------------------
68     !   L o c a l   P a r a m e t e r s
69     !-----------------------------------------------
70           DOUBLE PRECISION, PARAMETER :: F2O3 = 2.D0/3.D0
71     !-----------------------------------------------
72     !   L o c a l   V a r i a b l e s
73     !-----------------------------------------------
74     !
75     !
76     !                      Indices
77           INTEGER          I, J, K, IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
78                           IM, JM, KM
79           INTEGER          IMJPK, IMJMK, IMJKP, IMJKM, IPJKM, IPJMK, IJMKP, &
80                           IJMKM, IJPKM, I1, J1, K1
81     !
82     !                      Solids phase
83           INTEGER          M
84     !
85     !
86     !                      U_g at the north face of the THETA cell-(i, j+1/2, k)
87           DOUBLE PRECISION U_g_N
88     !
89     !                      U_g at the south face of the THETA cell-(i, j-1/2, k)
90           DOUBLE PRECISION U_g_S
91     !
92     !                      U_g at the top face of the THETA cell-(i, j, k+1/2)
93           DOUBLE PRECISION U_g_T
94     !
95     !                      U_g at the bottom face of the THETA cell-(i, j, k-1/2)
96           DOUBLE PRECISION U_g_B
97     !
98     !                      U_g at the center of the THETA cell-(i, j, k)
99     !                      Calculated for Cylindrical coordinates only.
100           DOUBLE PRECISION U_g_C
101     !
102     !                      V_g at the east face of the THETA cell-(i+1/2, j, k)
103           DOUBLE PRECISION V_g_E
104     !
105     !                      V_g at the west face of the THETA cell-(i-1/2, j, k)
106           DOUBLE PRECISION V_g_W
107     !
108     !                      V_g at the top face of the THETA cell-(i, j, k+1/2)
109           DOUBLE PRECISION V_g_T
110     !
111     !                      V_g at the bottom face of the THETA cell-(i, j, k-1/2)
112           DOUBLE PRECISION V_g_B
113     !
114     !                      W_g at the east face of the THETA cell-(i+1/2, j, k)
115           DOUBLE PRECISION W_g_E
116     !
117     !                      W_g at the west face of the THETA cell-(1-1/2, j, k)
118           DOUBLE PRECISION W_g_W
119     !
120     !                      W_g at the north face of the THETA cell-(i, j+1/2, k)
121           DOUBLE PRECISION W_g_N
122     !
123     !                      W_g at the south face of the THETA cell-(i, j-1/2, k)
124           DOUBLE PRECISION W_g_S
125     !
126     !                      W_g at the center of the THETA cell-(i, j, k).
127     !                      Calculated for Cylindrical coordinates only.
128           DOUBLE PRECISION W_g_C
129     !
130     !                      Cell center value of U_g
131           DOUBLE PRECISION UGC
132     !
133     !                      Cell center value of V_g
134           DOUBLE PRECISION VGC
135     !
136     !                      Cell center value of W_g
137           DOUBLE PRECISION WGC
138     !
139     !                      trace_g and eddy viscosity
140           DOUBLE PRECISION Trace_G, Mu_gas_t
141     !
142     !
143           DOUBLE PRECISION C_Eps_Pope, Xsi_Pope, UG(3,3), D_g
144     !
145     !  Production of Turb. Due to shear, Turb Visc, and Constants.
146     !  See Wilcox PP. 89
147           DOUBLE PRECISION Tauij_gDUi_gODxj, C_MU, Sigma_k, Sigma_E, Kappa
148           DOUBLE PRECISION Pos_Tauij_gDUi_gODxj, Neg_Tauij_gDUi_gODxj
149           DOUBLE PRECISION Ceps_1, Ceps_2, C_Eps_3, Check_Log
150           DOUBLE PRECISION Pos_PI_kq_2, Neg_PI_kq_2
151     ! Modif. for Sof Local Var.
152     !
153           INTEGER :: P,Q
154     !-----------------------------------------------
155     
156           IF( .NOT. K_Epsilon) RETURN
157     !
158     ! M should be forced to be equal to one to get some info. from solids phase.
159             M = ONE
160     !
161     ! Add constants. Some of them may not be used now but will be necessary when the Simonin
162     ! model is implemented in the standard MFIX model.
163     ! Most of these constants have the same names and values as the ones defined in Wilcox
164     ! book (turbulence modeling for CFD)
165             C_MU = 9.0D-02
166             Kappa = 0.42D+0
167             Sigma_k = 1.0D0
168             Sigma_E = 1.3D0
169             Ceps_1 = 1.44D0 !should be 1.6 for axisymmetric cases with no Pope Correction.
170             Ceps_2 = 1.92D0
171             C_Eps_3 = 1.2D0 ! for Simonin model
172             C_Eps_Pope = 0.79D0 ! for Pope Correction.
173     
174     !
175     !  ---  Remember to include all the local variables here for parallel
176     !  ---- processing
177     !!!$omp  parallel do private(ijk, L)
178           DO IJK = IJKSTART3, IJKEND3
179              IF (FLUID_AT(IJK)) THEN
180     !
181     !
182     !------------------
183     ! Part copied from calc_mu_g.f to calculate the rate of strain Sii tensor
184     ! of the gas phase.
185     ! Sof@fluent.com
186     !------------------
187     !
188                 I = I_OF(IJK)
189                 J = J_OF(IJK)
190                 K = K_OF(IJK)
191                 IM = IM1(I)
192                 JM = JM1(J)
193                 KM = KM1(K)
194     
195                 IMJK = IM_OF(IJK)
196                 IPJK = IP_OF(IJK)
197                 IJMK = JM_OF(IJK)
198                 IJPK = JP_OF(IJK)
199                 IJKM = KM_OF(IJK)
200                 IJKP = KP_OF(IJK)
201                 IMJPK = IM_OF(IJPK)
202                 IMJMK = IM_OF(IJMK)
203                 IMJKP = IM_OF(IJKP)
204                 IMJKM = IM_OF(IJKM)
205                 IPJKM = IP_OF(IJKM)
206                 IPJMK = IP_OF(IJMK)
207                 IJMKP = JM_OF(IJKP)
208                 IJMKM = JM_OF(IJKM)
209                 IJPKM = JP_OF(IJKM)
210     !
211     !         Find fluid velocity values at faces of the cell
212                 U_G_N = AVG_Y(AVG_X_E(U_G(IMJK),U_G(IJK),I),AVG_X_E(U_G(IMJPK),U_G(&
213                    IJPK),I),J)                       !i, j+1/2, k
214                 U_G_S = AVG_Y(AVG_X_E(U_G(IMJMK),U_G(IJMK),I),AVG_X_E(U_G(IMJK),U_G&
215                    (IJK),I),JM)                      !i, j-1/2, k
216                 U_G_T = AVG_Z(AVG_X_E(U_G(IMJK),U_G(IJK),I),AVG_X_E(U_G(IMJKP),U_G(&
217                    IJKP),I),K)                       !i, j, k+1/2
218                 U_G_B = AVG_Z(AVG_X_E(U_G(IMJKM),U_G(IJKM),I),AVG_X_E(U_G(IMJK),U_G&
219                    (IJK),I),KM)                      !i, j, k-1/2
220                 V_G_E = AVG_X(AVG_Y_N(V_G(IJMK),V_G(IJK)),AVG_Y_N(V_G(IPJMK),V_G(&
221                    IPJK)),I)                         !i+1/2, j, k
222                 V_G_W = AVG_X(AVG_Y_N(V_G(IMJMK),V_G(IMJK)),AVG_Y_N(V_G(IJMK),V_G(&
223                    IJK)),IM)                         !i-1/2, j, k
224                 V_G_T = AVG_Z(AVG_Y_N(V_G(IJMK),V_G(IJK)),AVG_Y_N(V_G(IJMKP),V_G(&
225                    IJKP)),K)                         !i, j, k+1/2
226                 V_G_B = AVG_Z(AVG_Y_N(V_G(IJMKM),V_G(IJKM)),AVG_Y_N(V_G(IJMK),V_G(&
227                    IJK)),KM)                         !i, j, k-1/2
228                 W_G_N = AVG_Y(AVG_Z_T(W_G(IJKM),W_G(IJK)),AVG_Z_T(W_G(IJPKM),W_G(&
229                    IJPK)),J)                         !i, j+1/2, k
230                 W_G_S = AVG_Y(AVG_Z_T(W_G(IJMKM),W_G(IJMK)),AVG_Z_T(W_G(IJKM),W_G(&
231                    IJK)),JM)                         !i, j-1/2, k
232                 W_G_E = AVG_X(AVG_Z_T(W_G(IJKM),W_G(IJK)),AVG_Z_T(W_G(IPJKM),W_G(&
233                    IPJK)),I)                         !i+1/2, j, k
234                 W_G_W = AVG_X(AVG_Z_T(W_G(IMJKM),W_G(IMJK)),AVG_Z_T(W_G(IJKM),W_G(&
235                    IJK)),IM)                         !i-1/2, j, k
236     !
237                 IF (CYLINDRICAL) THEN
238     !                                                !i, j, k
239                    U_G_C = AVG_X_E(U_G(IMJK),U_G(IJK),I)
240     !                                                !i, j, k
241                    W_G_C = AVG_Z_T(W_G(IJKM),W_G(IJK))
242                 ELSE
243                    U_G_C = ZERO
244                    W_G_C = ZERO
245                 ENDIF
246     !
247     !         Calculate velocity components at i, j, k to be used in wall functions
248                 UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
249                 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
250                 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
251     
252     !=======================================================================
253     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
254     !=======================================================================
255     
256                 IF(.NOT.CUT_CELL_AT(IJK)) THEN
257     !
258     !       Velocity derivives computed for the Pope Approximation.
259     !
260                    UG(1,1) = (U_G(IJK)-U_G(IMJK))*ODX(I)
261                    UG(1,2) = (U_G_N - U_G_S)*ODY(J)
262                    UG(1,3) = (U_G_T-U_G_B)*(OX(I)*ODZ(K))-W_G_C*OX(I)
263                    UG(2,1) = (V_G_E-V_G_W)*ODX(I)
264                    UG(2,2) = (V_G(IJK)-V_G(IJMK))*ODY(J)
265                    UG(2,3) = (V_G_T-V_G_B)*(OX(I)*ODZ(K))
266                    UG(3,1) = (W_G_E - W_G_W)*ODX(I)
267                    UG(3,2) = (W_G_N-W_G_S)*ODY(J)
268                    UG(3,3) = (W_G(IJK)-W_G(IJKM))*(OX(I)*ODZ(K)) + U_G_C*OX(I)
269     !
270     !
271     !         Find components of fluid phase strain rate
272     !         tensor, D_g, at center of the cell - (i, j, k)
273                    D_G(1,1) = (U_G(IJK)-U_G(IMJK))*ODX(I)
274                    D_G(1,2) = HALF*((U_G_N - U_G_S)*ODY(J)+(V_G_E-V_G_W)*ODX(I))
275                    D_G(1,3) = HALF*((W_G_E - W_G_W)*ODX(I)+(U_G_T-U_G_B)*(OX(I)*ODZ(K)&
276                       )-W_G_C*OX(I))
277                    D_G(2,1) = D_G(1,2)
278                    D_G(2,2) = (V_G(IJK)-V_G(IJMK))*ODY(J)
279                    D_G(2,3)=HALF*((V_G_T-V_G_B)*(OX(I)*ODZ(K))+(W_G_N-W_G_S)*ODY(J))
280                    D_G(3,1) = D_G(1,3)
281                    D_G(3,2) = D_G(2,3)
282                    D_G(3,3) = (W_G(IJK)-W_G(IJKM))*(OX(I)*ODZ(K)) + U_G_C*OX(I)
283     
284     
285     
286                 ELSE  ! CUT_CELL
287     
288                    CALL CG_CALC_VEL_G_GRAD(IJK,UG)
289     
290                    DO P = 1,3
291                       DO Q = 1,3
292                          D_G(P,Q) = HALF * (UG(P,Q)+UG(Q,P))
293                       ENDDO
294                    ENDDO
295     
296                 ENDIF
297     
298     !=======================================================================
299     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
300     !=======================================================================
301     !!!
302                 Trace_g = D_G(1,1) + D_G(2,2) + D_G(3,3)
303     
304     !
305     ! Sof@fluent.com
306     !
307     ! Pope's correction in 2-D to the Epsilon Equation for a round-Jet
308     ! from Wilcox book, Page 103.
309                 Xsi_Pope = ZERO
310                 DO I1 = 1,3
311                    DO J1 = 1,3
312                       DO K1 = 1,3
313                          Xsi_Pope = Xsi_Pope + (UG(I1,J1) - UG(J1,I1))&
314                                     *(UG(J1,K1) - UG(K1,J1))*(UG(K1,I1) + UG(I1,K1))
315                       END DO
316                    END DO
317                 END DO
318     !
319     
320                 Xsi_Pope = Xsi_Pope/6. * K_Turb_G(IJK)**2/(E_Turb_G(IJK)+Small_number)
321     !
322     !
323     ! This IF statment is to ensure that we are using the turbulent viscosity
324     ! and NOT the effective viscosity.
325     !
326                 IF (MU_GT(IJK) .GE. MU_g(IJK)) THEN
327                    Mu_gas_t = MU_GT(IJK) - MU_g(IJK)
328                 ELSE
329                    Mu_gas_t = ZERO
330                 ENDIF
331     !
332     !        Calculate Tau(i,j)*dUi/dXj (production term in the K Equation
333     
334     !=======================================================================
335     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
336     !=======================================================================
337     
338                 IF(.NOT.CUT_CELL_AT(IJK)) THEN
339     
340                    Tauij_gDUi_gODxj = 2D0*Mu_gas_t*(                             &
341                                       D_G(1,1) * D_G(1,1) +                      &
342                                       D_G(1,2) * (U_G_N - U_G_S)*ODY(J) +        &
343                                       D_G(1,3) * ((U_G_T-U_G_B)*                 &
344                                       (OX(I)*ODZ(K))-W_G_C*OX(I)) +              &
345                                       D_G(2,1) * (V_G_E-V_G_W)*ODX(I) +          &
346                                       D_G(2,2) * D_G(2,2) +                      &
347                                       D_G(2,3) * (V_G_T-V_G_B)*(OX(I)*ODZ(K)) +  &
348                                       D_G(3,1) * (W_G_E - W_G_W)*ODX(I) +        &
349                                       D_G(3,2) * (W_G_N-W_G_S)*ODY(J) +          &
350                                       D_G(3,3) * D_G(3,3)) -                     &
351                                       F2O3 * RO_G(IJK) * K_Turb_G(IJK)*Trace_g   &
352                                      - F2O3 * Mu_gas_t * Trace_g**2
353     
354                 ELSE  ! CUT_CELL       ! This is actually not used because of wall functions in cut cells
355     
356     !               Tauij_gDUi_gODxj = 2D0*Mu_gas_t*(                             &
357     !                                  D_G(1,1) * UG(1,1)  +                      &
358     !                                  D_G(1,2) * UG(1,2)  +                      &
359     !                                  D_G(1,3) * UG(1,3)  +                      &
360     !                                  D_G(2,1) * UG(2,1)  +                      &
361     !                                  D_G(2,2) * UG(2,2)  +                      &
362     !                                  D_G(2,3) * UG(2,3)  +                      &
363     !                                  D_G(3,1) * UG(3,1)  +                      &
364     !                                  D_G(3,2) * UG(3,2)  +                      &
365     !                                  D_G(3,3) * UG(3,3)) -                      &
366     !                                  F2O3 * RO_G(IJK) * K_Turb_G(IJK)*Trace_g   &
367     !                                 - F2O3 * Mu_gas_t * Trace_g**2
368     
369     
370                 ENDIF
371     
372     !=======================================================================
373     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
374     !=======================================================================
375     
376     
377     ! To avoid very small negative numbers
378     !
379     
380                 IF (Tauij_gDUi_gODxj .GE. ZERO) THEN
381                    Pos_Tauij_gDUi_gODxj = Tauij_gDUi_gODxj
382                    Neg_Tauij_gDUi_gODxj = ZERO
383                 ELSE
384                    Pos_Tauij_gDUi_gODxj = ZERO
385                    Neg_Tauij_gDUi_gODxj = Tauij_gDUi_gODxj
386                 ENDIF
387     !
388     !
389     ! Interaction terms in the K-Epsilon equations FOR USE Simonin and Ahmadi models
390     !
391                 IF(SIMONIN) THEN
392                    Pos_PI_kq_2 = F_GS(IJK,1)*K_12(IJK)
393                    Neg_PI_kq_2 = F_GS(IJK,1)*2.0D0* K_Turb_G(IJK)
394                 ELSE IF(AHMADI) THEN
395                    Pos_PI_kq_2 = F_GS(IJK,1)*3.0D0*Theta_m(IJK,M)
396                    Neg_PI_kq_2 = F_GS(IJK,1)*2.0D0* K_Turb_G(IJK)
397                    C_Eps_3 = zero ! no extra terms in epsilon equation for Ahmadi model !
398                 ELSE
399                    Pos_PI_kq_2 = ZERO
400                    Neg_PI_kq_2 = ZERO
401                 ENDIF
402     
403                 IF(K_Turb_G(IJK) > Small_number .AND. E_Turb_G(IJK) > Small_number) THEN
404     !
405     ! Start Adding source terms to the K equation
406     !
407                    K_Turb_G_c (IJK) = (EP_g(IJK) * Pos_Tauij_gDUi_gODxj +   &
408                                        Pos_PI_kq_2 )
409     !
410                    K_Turb_G_p (IJK) =(-EP_g(IJK) * Neg_Tauij_gDUi_gODxj+    &
411                                        EP_g(IJK)*RO_G(IJK)*E_Turb_G(IJK)+   &
412                                        Neg_PI_kq_2)/K_Turb_G(IJK)
413     !
414     ! Implementing wall functions targeted to fluid cells next to walls...
415     ! Setting Source and sink terms in the K equation since the production
416     ! in the K eq. due to shear should include the LOG law of the wall.
417     !
418                    IF(WALL_AT(JP_OF(IJK)).OR.WALL_AT(JM_OF(IJK))) THEN
419     !
420                       Check_Log = LOG(9.81*C_mu**0.25*                          &
421                                   RO_G(IJK)*SQRT(K_Turb_G(IJK))*DY(J)/      &
422                                   2.0D+0/Mu_g(IJK))
423                       IF(Check_Log .LE. ZERO) THEN
424                          K_Turb_G_c (IJK) = zero
425                          K_Turb_G_p (IJK) = zero
426                       ELSE
427                          K_Turb_G_c (IJK) = SQRT(C_mu)*2.D+0/DY(J)*               &
428                                             MAX(ABS(UGC),ABS(WGC))*               &
429                                             EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)     &
430                                             /Check_Log
431     
432                          K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)                 &
433                            *C_mu**0.75*K_Turb_G(IJK)**1.5/DY(J)*2.0D+0/Kappa)     &
434                             )/K_Turb_G(IJK)
435                       ENDIF !for check_log less than zero
436     
437                    ELSE IF(WALL_AT(KP_OF(IJK)).OR.WALL_AT(KM_OF(IJK))) THEN
438     !
439                       Check_Log = LOG(9.81*C_mu**0.25*                            &
440                                   RO_G(IJK)*SQRT(K_Turb_G(IJK))/                  &
441                                   (ODZ(K)*OX(I)*2.0D+0)/Mu_g(IJK))
442                       IF(Check_Log .LE. ZERO) THEN
443                          K_Turb_G_c (IJK) = zero
444                          K_Turb_G_p (IJK) = zero
445                       ELSE
446                          K_Turb_G_c (IJK) = SQRT(C_mu)*(ODZ(K)*OX(I)*2.0D+0)*     &
447                                             MAX(ABS(UGC),ABS(VGC))                &
448                                            *EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)     &
449                                            /Check_Log
450     
451                          K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)                 &
452                                           *C_mu**0.75*K_Turb_G(IJK)**1.5/Kappa*   &
453                                           (ODZ(K)*OX(I)*2.0D+0))                  &
454                                           )/K_Turb_G(IJK)
455                       ENDIF !for check_log less than zero
456     
457                    ENDIF !For walls
458     !
459     !
460     ! For Cylindrical cases, wall_at (IP) is a wall cell, but wall_at (IM) is
461     ! the axis of symmetry where wall functions obviously don't apply.
462     !
463               IF(CYLINDRICAL) THEN
464                 IF (WALL_AT(IP_OF(IJK)))  THEN
465     !!
466                   Check_Log = LOG(9.81*C_mu**0.25*                          &
467                          RO_G(IJK)*SQRT(K_Turb_G(IJK))*DX(I)/               &
468                            2.0D+0/Mu_g(IJK))
469     
470                   IF(Check_Log .LE. ZERO) THEN
471                     K_Turb_G_c (IJK) = zero
472                     K_Turb_G_p (IJK) = zero
473                   ELSE
474                     K_Turb_G_c (IJK) = SQRT(C_mu)*2.D+0/DX(I)*              &
475                                 MAX(ABS(VGC),ABS(WGC))            &
476                                *EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)           &
477                               /Check_Log
478     
479                     K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)                &
480                       *C_mu**0.75*K_Turb_G(IJK)**1.5/DX(I)*2.0D+0/Kappa)    &
481                       )/K_Turb_G(IJK)
482     
483                   ENDIF !for check_log less than zero
484                 ENDIF  ! for wall cells in I direction
485               ELSE IF (WALL_AT(IP_OF(IJK)).OR.WALL_AT(IM_OF(IJK))) THEN
486     !!
487                 Check_Log = LOG(9.81*C_mu**0.25*                            &
488                          RO_G(IJK)*SQRT(K_Turb_G(IJK))*DX(I)/               &
489                            2.0D+0/Mu_g(IJK))
490     
491                 IF(Check_Log .LE. ZERO) THEN
492                   K_Turb_G_c (IJK) = zero
493                   K_Turb_G_p (IJK) = zero
494                 ELSE
495                   K_Turb_G_c (IJK) = SQRT(C_mu)*2.D+0/DX(I)*                &
496                    MAX(ABS(VGC),ABS(WGC))*                        &
497                    EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)                        &
498                                   /Check_Log
499     
500                   K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)                  &
501                       *C_mu**0.75*K_Turb_G(IJK)**1.5/DX(I)*2.0D+0/Kappa)    &
502                     )/K_Turb_G(IJK)
503                 ENDIF !for check_log less than zero
504               ENDIF ! for cylindrical
505     
506     
507     
508              IF(CUT_CELL_AT(IJK)) THEN
509     !
510                 Check_Log = LOG(9.81*C_mu**0.25*                          &
511                                 RO_G(IJK)*SQRT(K_Turb_G(IJK))*DELH_Scalar(IJK)/Mu_g(IJK))
512                 IF(Check_Log .LE. ZERO) THEN
513                   K_Turb_G_c (IJK) = zero
514                   K_Turb_G_p (IJK) = zero
515                 ELSE
516     !              K_Turb_G_c (IJK) = SQRT(C_mu)/DELH_Scalar(IJK)*               &
517     !              MAX(ABS(UGC),ABS(VGC),ABS(WGC))*                        &
518     !              EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)                        &
519     !                            /Check_Log
520     
521     
522                   K_Turb_G_c (IJK) = SQRT(C_mu)/DELH_Scalar(IJK)*               &
523                   DSQRT(UGC**2 + VGC**2 + WGC**2) *                        &
524                   EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)                        &
525                                 /Check_Log
526     
527                   K_Turb_G_p (IJK) = (EP_g(IJK)*RO_G(IJK)                 &
528                     *C_mu**0.75*K_Turb_G(IJK)**1.5)/(DELH_Scalar(IJK)*Kappa)     &
529                      /K_Turb_G(IJK)
530                 ENDIF !for check_log less than zero
531     
532     
533              ENDIF
534     
535     
536     
537     !
538     !            Diffusion coefficient for turbulent kinetic energy (K)
539     !
540                  Dif_K_Turb_G(IJK) = EP_g(IJK)* (MU_G(IJK) + Mu_gas_t /Sigma_k)
541     !
542     ! Add here Dissipation of Turbulence source terms
543     !
544                  E_Turb_G_c (IJK) = (Ceps_1 *EP_g(IJK)*Pos_Tauij_gDUi_gODxj &
545                                     *E_Turb_G(IJK)/K_Turb_G(IJK) +          &
546                                     C_Eps_3*Pos_PI_kq_2*                    &
547                                     E_Turb_G(IJK)/K_Turb_G(IJK))           !&
548     !
549     ! Pope Correction in Xsi_Pope, Add it to E_Turb_G_c to use this option.
550     !
551                                    ! + C_Eps_Pope*RO_G(IJK)*EP_g(IJK)*&
552                                    !  ZMAX(Xsi_Pope)
553     !
554                  E_Turb_G_p (IJK) = -Ceps_1 *EP_g(IJK)*Neg_Tauij_gDUi_gODxj &
555                                      /K_Turb_G(IJK) +                       &
556                                     Ceps_2 * EP_g(IJK) *RO_G(IJK)           &
557                                     *E_Turb_G(IJK)/K_Turb_G(IJK)            &
558                                     + C_Eps_3*(Neg_PI_kq_2)                 &
559                                      /K_Turb_G(IJK)
560     !
561     !
562     !            Diffusion coefficient for Dissipation of turbulent energy (Epsilon)
563     !
564                  Dif_E_Turb_G(IJK) =EP_g(IJK)* (MU_G(IJK) + Mu_gas_t /Sigma_E)
565     
566             ELSE
567               K_Turb_G_c (IJK) = zero
568               K_Turb_G_p (IJK) = zero
569               E_Turb_G_c (IJK) = zero
570               E_Turb_G_p (IJK) = zero
571               Dif_K_Turb_G(IJK) = zero
572               Dif_E_Turb_G(IJK) = zero
573     
574             ENDIF !for K_turb_G and E_Turb_G having very small numbers
575     
576              ENDIF
577           END DO
578     !
579     !
580           RETURN
581           END SUBROUTINE K_Epsilon_PROP
582