File: RELATIVE:/../../../mfix.git/model/k_epsilon_prop.f

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