!
! Copyright (C) 1991-2004  ; All Rights Reserved ; Colorado State University
! Colorado State University Research Foundation ; ATMET, LLC
! 
! This file is free software; you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software 
! Foundation; either version 2 of the License, or (at your option) any later version.
! 
! This software is distributed in the hope that it will be useful, but WITHOUT ANY 
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 
! PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this 
! code; if not, write to the Free Software Foundation, Inc., 
! 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!======================================================================================

Subroutine tkescl (m1,m2,m3,m4,ia,iz,ja,jz &
          ,tkep,tket,epsp &
          ,vt3da,vt3dc,vt3dh,vt3di,vt3dj &
          ,scr1,scr2,rtgt &
          ,vt3dd,vt3de,dxt &
          ,ustar,patch_area,dn0)     
!_STC  +-------------------------------------------------------------------+
!_STC  \  Authors: S. Trini Castelli, E. Ferrero                           \
!_STC  \  First version:        1998                                       \
!_STC  \  Last updated version: 2000                                       \
!_STC  \  Version for RAMS4.3 : 2001                                       \
!_STC  +-------------------------------------------------------------------+
!_STC  \  This routine calculates the mixing coefficients based on the     \
!_STC  \    prognostic TKE and diagnostic length scale, according to a     \
!_STC  \    standard E-l closure model                                     \
!_STC  +-------------------------------------------------------------------+

use mem_grid
use mem_scratch
use ke_coms
use rconstants

implicit none

integer :: m1,m2,m3,m4,ia,iz,ja,jz
real, dimension(m1,m2,m3) :: tkep,tket,epsp  &
          ,vt3da,vt3dc,vt3dh,vt3di,vt3dj,scr1,scr2  &
          ,vt3dd,vt3de,dn0
real, dimension(m2,m3) ::  rtgt,dxt
real, dimension(m2,m3,m4) ::  ustar,patch_area

real :: coef2
real :: dudz_1,dudz_2,dvdz_1,dvdz_2,dvelmin_1,dvelmin_2 &
        ,psizil_1q,psizil_2q,psizil_q
real :: al0,al0_zil,scl,sumtkz,sumtk,tkep_k2,dzloc,dpsi2dz,scl_max
real, external :: ssumvect
integer :: k,i,j,np

coef2= 1./(c_eps**(2./3.))

!  - C_EPS,C_MI,ALF_TKE,ALF_EPS,ALF_THT assigned in rconstants.h
!  - aymptotic mixing length option assigned in rconstants.h
!  - Constant value (Ying, 1992) for the asymptotic mixing length assigned 
!  in rconstants.h

      do j=ja,jz
         do i=ia,iz

!_STC....................................
! Calculation of the asymptotic lenght scale: OPTIONS !
! MELLOR-YAMADA 1974: pre-calculation
!_STC....................................

            if(iopzl.eq.2) then

               do k=2,m1-1
                  vctr30(k) = sqrt(2.0 * tkep(k,i,j))
                  vctr31(k)=.75*vctr30(k)/sqrt(max(1.e-20,vt3dj(k,i,j)))
                  vctr1(k)=(zt(k)-zm(1))*rtgt(i,j)
                  dzloc=(zm(k)-zm(k-1))*rtgt(i,j)
                  vctr33(k)=vctr30(k)*dzloc
                  vctr32(k)=vctr33(k)*vctr1(k)
               enddo

               sumtkz = ssumvect(m1-2,vctr32(2),1)
               sumtk  = ssumvect(m1-2,vctr33(2),1)
               al0=.1*sumtkz/sumtk

            endif


            do k=2,m1-1


!_STC....................................
!_STC-----> Ying's value
!_STC....................................

               if(iopzl == 1) then    ! YING's length scale
                  scl=vonk*vctr1(k)/(1.+vonk*vctr1(k)/al0_const) 

!_STC....................................
!_STC-----> Mellor-Yamada 1974 formulation
!_STC....................................

               elseif(iopzl == 2) then
                  scl=min(vonk*vctr1(k)/(1.+vonk*vctr1(k)/al0),vctr31(k))

!_STC................................................
!_STC-----> Zilitinkevich-Laikhtman 1965 formulation
!_STC................................................
!_5Mar02     trasformo la formulazione di Zilitinkevich da funzione di
!_5Mar02     psi a funzione di psi**2

               elseif(iopzl == 3) then ! opzione Zilitinkevich


!... calcolo della derivata in z e del valore della lunghezza l0


                  dudz_1=vt3dd(k,i,j)
                  dudz_2=vt3dd(k+1,i,j)
                  dvdz_1=vt3de(k,i,j)
                  dvdz_2=vt3de(k+1,i,j)
   
!.... valore minimo per le derivate in z di u e v
   
                  dvelmin_1=0.001/dzm(k)
                  dvelmin_2=0.001/dzm(k+1)
   
                  if(abs(dudz_1)<dvelmin_1 &
                            .and. dudz_1>=0.) dudz_1=dvelmin_1
                  if(abs(dudz_1)<dvelmin_1 &
                            .and. dudz_1<0.) dudz_1=-dvelmin_1
                  if(abs(dvdz_1)<dvelmin_1 &
                            .and. dvdz_1>=0.) dvdz_1=dvelmin_1
                  if(abs(dvdz_1)<dvelmin_1 &
                            .and. dvdz_1<0.) dvdz_1=-dvelmin_1
                  if(abs(dudz_2)<dvelmin_2 &
                            .and. dudz_2>=0.) dudz_2=dvelmin_2
                  if(abs(dudz_2)<dvelmin_2 &
                            .and. dudz_2<0.) dudz_2=-dvelmin_2
                  if(abs(dvdz_2)<dvelmin_2 &
                            .and. dvdz_2>=0.) dvdz_2=dvelmin_2
                  if(abs(dvdz_2)<dvelmin_2 &
                            .and. dvdz_2<0.) dvdz_2=-dvelmin_2
   
                  psizil_1q=dudz_1*dudz_1 + dvdz_1*dvdz_1  &
                                - alf_tht*vt3dj(k,i,j)
                  psizil_2q=dudz_2*dudz_2 + dvdz_2*dvdz_2  &
                                - alf_tht*vt3dj(k+1,i,j)

                  psizil_q=0.5*(psizil_1q+psizil_2q)

! 'special' unlucky case control !
                  if(psizil_2q == psizil_1q) then
                   psizil_1q=2*dvelmin_1*dvelmin_1 - alf_tht*vt3dj(k,i,j)
                   psizil_2q=2*dvelmin_2*dvelmin_2 - alf_tht*vt3dj(k+1,i,j)
                  endif
   
                  dpsi2dz=(psizil_2q-psizil_1q)*dzm(k)/rtgt(i,j)

                  if(dpsi2dz.eq.0.) stop 'dpsi2dz = 0. chissa perche!'

                  if(dpsi2dz.ne.0.) then
 
                  al0_zil = -4.* vonk *(psizil_q/dpsi2dz)    
                  if(al0_zil.lt.0.) al0_zil=abs(al0_zil)
               endif
   
               scl=vonk*vctr1(k)/(1. + vonk*vctr1(k)/al0_zil ) ! Zily

!_8Mar02   maximum limit for SCL: copied from Deardoff scheme

               scl_max = 1.5*(rtgt(i,j)*(zm(k) - zm(k-1)))**(1./3.) &
                   /dxt(i,j) ** (2./3.)
               if (scl .gt. scl_max) scl = scl_max

   
            endif     !fine variante Zilitinkevic
!----------------------------------------------------------------------------
!_STC................................................
!_STC	T.K.E. equation terms
!_STC................................................
!
!_STC................................................
!_STC......vertical diffusion coefficient assignment
!_STC................................................
           scr1(k,i,j)=c_mi*scl*sqrt(tkep(k,i,j))    ! Km

!_STC................................................
!_STC   ......production term
!_STC................................................

           vctr25(k)=scr1(k,i,j)*vt3dh(k,i,j)-0.333333*tkep(k,i,j)  &
                *(vt3da(k,i,j)+vt3dc(k,i,j)+scr2(k,i,j))

!_STC................................................
!_STC	......buoyancy term
!_STC................................................

           vctr27(k)=-(scr1(k,i,j)*alf_tht)*vt3dj(k,i,j)


!_STC................................................
!_STC	......dissipation term (diagnostic)
!_STC................................................

           vctr28(k)=c_eps*tkep(k,i,j)*sqrt(max(1.e-20,tkep(k,i,j)))/scl
             
           epsp(k,i,j)= vctr28(k)

!_STC................................................
!_STC	T.K.E. equation
!_STC................................................
           tket(k,i,j)=tket(k,i,j)+vctr25(k)+vctr27(k)-vctr28(k)

!_STC....................................................................
!_STC   Assignment of diffusion coefficients to the correspondent arrays
!_STC	Km(oriz) = Km(vert) : isotropy
!_STC....................................................................
!
           scr1(k,i,j)=scr1(k,i,j)*dn0(k,i,j)   ! account for density
           scr2(k,i,j)=scr1(k,i,j)   ! used in DIFFVEL as horiz. diff. coef.
           vt3dh(k,i,j)=scr1(k,i,j)*alf_tht   ! HEAT diffusion coefficient
           vt3di(k,i,j)=scr1(k,i,j)*alf_tke   ! TKE diffusion coefficient Ke

        enddo
        
! Bottom boundary condition        

      tkep_k2=0.
      do np=1,npatch
         tkep_k2=tkep_k2+(ustar(i,j,np)**2)*coef2  &
                         *patch_area(i,j,np)
      enddo
      tket(2,i,j)=(tkep_k2-tkep(2,i,j))/dtlt
   
     enddo
  enddo

return
END SUBROUTINE tkescl

!##############################################################################
Subroutine tkeeps (m1,m2,m3,m4,ia,iz,ja,jz  &
          ,tkep,tket,epsp,epst &
          ,vt3da,vt3dc,vt3dh,vt3di,vt3dj,scr1,scr2 &
          ,rtgt,ustar,patch_area,dn0)
!_STC  +-------------------------------------------------------------------+
!_STC  \  Authors: S. Trini Castelli, E. Ferrero                           \
!_STC  \  First version:        1998                                       \
!_STC  \  Last updated version: 2000                                       \
!_STC  \  Version for RAMS4.3 : 2001                                       \
!_STC  +-------------------------------------------------------------------+
!_STC  \  This routine calculates the mixing coefficients based on the     \
!_STC  \  prognostic TKE and prognostic TKE dissipation, according         \
!_STC  \  for instance, to the RODI (1980)                                 \
!_STC  \  (see also DETERING-ETLING's paper (BLM, 1985, n. 33) )           \
!_STC  +-------------------------------------------------------------------+

use mem_grid
use mem_scratch
use ke_coms
use rconstants

implicit none

integer :: m1,m2,m3,m4,ia,iz,ja,jz
real, dimension(m1,m2,m3) :: tkep,tket,epsp,epst  &
          ,vt3da,vt3dc,vt3dh,vt3di,vt3dj,scr1,scr2,dn0
real, dimension(m2,m3) :: rtgt
real, dimension(m2,m3,m4) ::  ustar,patch_area

integer :: k,i,j,np
real :: epsp_k2,tkep_k2
real :: coef1,coef_km_sqr

coef_km_sqr = sqrt(coef_km)
!
!_STC................................................................
!     Empirical constants and coefficients for the E-eps closure
!           C_MI,C1_EPS,C2_EPS,ALF_TKE,ALF_EPS,ALF_THT
!     assigned in ke_coms
!_STC................................................................


do j=ja,jz
   do i=ia,iz

      coef1=1./(vonk*(zt(2)-zm(1))*rtgt(i,j))

      do k=2,m1-1         !! loop extended to level 2

!_STC................................................
!_STC	......vertical diffusion coefficient
!_STC................................................
!
         scr1(k,i,j)=coef_km*(tkep(k,i,j)**2.)/epsp(k,i,j)
!
!_STC................................................
!_STC   wind shear production term: P_s
!_STC................................................
!
         vctr25(k)=scr1(k,i,j)*vt3dh(k,i,j)-0.333333*tkep(k,i,j)  &
                  *(vt3da(k,i,j)+vt3dc(k,i,j)+scr2(k,i,j))

!_STC................................................
!_STC  buoyancy term
!_STC    VT3DJ(K,I,J) = Brunt-Vaisala frequency from routine BRUVAIS
!_STC................................................

         vctr27(k)=-(scr1(k,i,j)*alf_tht)*vt3dj(k,i,j)

!_STC................................................
!_STC  dissipation term (from prognostic equation at a former time step)
!_STC................................................

         vctr28(k)= epsp(k,i,j)

!_STC................................................
!_STC     equation for the TKE 
!_STC................................................

         tket(k,i,j)=tket(k,i,j)+vctr25(k)+vctr27(k)-vctr28(k)
!_STC.....................................................................
!_STC    equation for the TKE dissipation: RODI
!_STC
!_STC in the "vertical layers" the lateral component of the buoyancy 
!_STC is orthogonal to gravity, so that there is not latera contribution
!_STC from the buoyancy and the correction with c3_eps vanishes. (see Rodi)
!_STC.....................................................................

         epst(k,i,j)= epst(k,i,j) &
           + c1_eps*epsp(k,i,j)*(vctr25(k)+vctr27(k))/tkep(k,i,j) &
           - c2_eps*epsp(k,i,j)*epsp(k,i,j)/tkep(k,i,j)
!
!
!_STC.....................................................................
!_STC  Array assignment of the diffusion coefficients for the
!_STC  scalar diffusion. In SCR2 is assigned the diffusion coefficient
!_STC  for the dissipation tkeps.
!_STC.....................................................................

         scr1(k,i,j)=scr1(k,i,j)*dn0(k,i,j)   ! account for density
         scr2(k,i,j)=scr1(k,i,j)            ! used in DIFFVEL as horizontal
                                            !      diffusion coefficient
         vt3dh(k,i,j)=scr1(k,i,j)*alf_tht   ! HEAT diffusion coefficient
         vt3di(k,i,j)=scr1(k,i,j)*alf_tke   ! TKE diffusion coefficient Ke

if(i==25.and.j==20.and.ngrid==ngrids.and.k<10) print '(a,i3,12e14.5)','tkeeps:',k  &
,tkep(k,i,j),epsp(k,i,j),scr1(k,i,j),tket(k,i,j),vctr25(k),vctr27(k),-vctr28(k)

      enddo
      
! bottom boundary conditions
      epsp_k2=0.
      tkep_k2=0.
      do np=1,m4
         tkep_k2=tkep_k2+(ustar(i,j,np)**2)*coef_km_sqr  &
             *patch_area(i,j,np)
!altern       tkep_k2=tkep_k2+(ustar(i,j,np)**2.)*patch_area(i,j,np)/C_ust_tk
         epsp_k2=epsp_k2+(ustar(i,j,np)**3)*coef1  &
                *patch_area(i,j,np)
!alternative  epsp_k2=c_eps*(tkep_k2**1.5)*coef1
      enddo
      tket(2,i,j)=(tkep_k2-tkep(2,i,j))/dtlt
      epst(2,i,j)=(epsp_k2-epsp(2,i,j))/dtlt

   enddo
enddo

if(ngrid==ngrids)print*,'-------------------------------',time,coef_km

return
END SUBROUTINE tkeeps

!##############################################################################
Subroutine mxtked (m1,m2,m3,ia,iz,ja,jz  &
   ,tkep,tket,vt3da,vt3dc,vt3dh,vt3dj,scr1,scr2,dxt,rtgt)
!  +-------------------------------------------------------------------+
!  \  this routine calculates the mixing coefficients based on the     \
!  \    prognostic tke and diagnostic length scale, according to the   \
!  \    deardorff scheme.                                              \
!  +-------------------------------------------------------------------+

use mem_grid
use mem_scratch
use rconstants

implicit none

integer :: m1,m2,m3,ia,iz,ja,jz
real, dimension(m1,m2,m3) :: tkep,tket,vt3da,vt3dc,vt3dh,vt3dj,scr1,scr2
real, dimension(m2,m3) :: dxt,rtgt
integer :: i,j,k
real :: sqrttkep,tket2,c1,sclu,scl

do j = ja,jz
   do i = ia,iz
      tket2 = tket(2,i,j)
      c1 = 3.375 * rtgt(i,j) / dxt(i,j) ** 2  ! 3.375 = 1.5 ** 3
      do k = 2,m1-1
         sqrttkep = sqrt(tkep(k,i,j))
         sclu = (c1 * (zm(k) - zm(k-1))) ** .333333
         if (vt3dj(k,i,j) .lt. 1.e-10) then
            scl = sclu                                   ! unstable case
         else
            scl = .76 * sqrt(tkep(k,i,j) / vt3dj(k,i,j))   !stable case
            if (scl .gt. sclu) scl = sclu
         endif

! now, vt3dh is strain rate and scr2 is dw/dz

         scr1(k,i,j) = .1 * scl * sqrttkep
         vctr25(k) = scr1(k,i,j) * vt3dh(k,i,j)  &
            - .333333 * tkep(k,i,j)  &
            * (vt3da(k,i,j) + vt3dc(k,i,j) + scr2(k,i,j))
         vt3dh(k,i,j) = scr1(k,i,j) * (1. + 2. * scl / sclu)
         scr2(k,i,j) = scr1(k,i,j)

! now, vt3dh is vkh and scr2 is hkm

         vctr26(k) = .19 + .51 * scl / sclu
         vctr27(k) = - vt3dh(k,i,j) * vt3dj(k,i,j)
         vctr28(k) = vctr26(k) * tkep(k,i,j) * sqrttkep / scl
         tket(k,i,j) = tket(k,i,j) + vctr25(k) + vctr27(k)  &
            - vctr28(k)
      enddo
      tket(2,i,j) = tket2 + vctr25(3) + vctr27(2) - vctr28(2)

   enddo
enddo

return
END SUBROUTINE mxtked

!##############################################################################
Subroutine tkemy (m1,m2,m3,ia,iz,ja,jz,jd  &
   ,tkep,tket,vt3dh,vt3di,vt3dj,scr1,rtgt  &
   ,theta,dn0,up,vp,wp,sflux_u,sflux_v,sflux_w,sflux_t,tkep2)

use mem_grid
use mem_scratch
use rconstants

implicit none

integer :: m1,m2,m3,ia,iz,ja,jz,jd
real, dimension(m1,m2,m3) :: tkep,tket,vt3dh,vt3di,vt3dj,scr1,theta,dn0  &
   ,up,vp,wp
real, dimension(m2,m3) :: rtgt,sflux_u,sflux_v,sflux_w,sflux_t
real, dimension(m1) :: tkep2
integer :: i,j,k
real :: a1,a2,b1,b2,c1,aux1,aux2,rf1,rf2,rf3,rf4,wght1,wght3,sumtkz,sumtk  &
   ,al0,tket2,ri,rf,shr,smr,tker,qq,ssmf,shf,sh0,ssm,aux,gm,gh,sm1,sm2  &
   ,sh1,sh2,dzloc
real, external :: ssumvect

data a1,a2,b1,b2,c1/0.92,0.74,16.6,10.1,0.08/
data aux1,aux2/0.758964199,2.58286747/
data rf1,rf2,rf3,rf4/1.,0.191232309,0.223117196,0.234067819/

!        7 - mellor and yamada (after andre et al, 1978)

!Saleeby(2014):Check older code for possible experimental values

wght3=1.
wght1=1.0-wght3

do j=ja,jz
   do i=ia,iz
     do k=2,m1-1
         tkep2(k) = 2.0 * tkep(k,i,j)
         vctr30(k) = sqrt(tkep2(k))
!Upper limit for length scale "L" in stable conditions (Andre et al 1978)
         vctr31(k)=.75*vctr30(k)/sqrt(max(1.e-20,vt3dj(k,i,j)))
         vctr1(k)=(zt(k)-zm(1))*rtgt(i,j)
         dzloc=(zm(k)-zm(k-1))*rtgt(i,j)
         vctr33(k)=vctr30(k)*dzloc
         vctr32(k)=vctr33(k)*vctr1(k)
      enddo

      sumtkz = ssumvect(m1-2,vctr32(2),1)
      sumtk  = ssumvect(m1-2,vctr33(2),1)
      al0=.1*sumtkz/sumtk
      tket2=tket(2,i,j)

      do k=2,m1-1
!Turbulent length scale "L" after Mellor-Yamada (1982)
         vctr9(k)=min(vonk*vctr1(k)/(1.+vonk*vctr1(k)/al0)  &
              ,vctr31(k))

! --- for growing turbulence use helfand and labraga's modified sm and sh (sh0)
!Brunt-Vaisala contribution over vertical strain
         ri=min(vt3dj(k,i,j)/max(vt3di(k,i,j),1.e-11),0.190)
         rf=min(0.6588*(ri+0.1776-sqrt(ri*(ri-0.3221)  &
              +0.03156)),0.16)
         shr=aux2*(rf-rf2)/(rf-rf1)
         smr=aux1*(rf-rf4)/(rf-rf3)*shr
         tker=max(16.6*vctr9(k)*vctr9(k)*(smr*vt3di(k,i,j)  &
            -shr*vt3dj(k,i,j)),2.*tkmin)
         if(tker.gt.tkep2(k) )then
            qq=sqrt(tkep2(k)/tker)
            ssmf=qq*smr
            shf=qq*shr
            sh0=shf
            ssm=ssmf
         else

! --- for decaying turbulence use mellor and yamada's sm and sh (sh0)
            aux=vctr9(k)*vctr9(k)/tkep2(k) ! L^2 / 2e
            gm=aux*vt3di(k,i,j) ! L^2 / 2e  * vertical strain
            gh=-aux*vt3dj(k,i,j) ! L^2 / 2e  * Brunt-vaisala contribution

            sm1=0.6992-9.33948672*gh
            sm2=1.-(36.7188-187.4408515*gh+88.83949824*gm)*gh  &
                 +5.0784*gm
            ssm=sm1/max(sm2,1e-10) !Sm nondimensional eddy diffusivity

            sh1=0.74-4.0848*ssm*gm
            sh2=1.-30.5916*gh
            sh0=sh1/max(sh2,1e-10) !Sh nondimensional eddy diffusivity
         endif

!Vertical eddy diffusivity for momentum VKM
         scr1(k,i,j)=vctr9(k)*vctr30(k)*ssm
!Vertical eddy diffusivity for heat (scalars) VKH
         vt3dh(k,i,j)=vctr9(k)*vctr30(k)*sh0
!Ps + Pb or (shear production term + buoyancy production term in
!prognostic TKE equation
!(VKM*[(du/dz)^2 + (dv/dz)^2])  - (VKH*(-g/theta * dtheta/dz))
         vctr5(k)=scr1(k,i,j)*vt3di(k,i,j)  &
            -vt3dh(k,i,j)*vt3dj(k,i,j)
         tket(k,i,j)=tket(k,i,j) + 0.5 * vctr5(k)  &
            -tkep2(k) * sqrt(max(1.e-20,tkep2(k)))/(vctr9(k)*16.6)
!Density weighting for VKM and VKH
         scr1(k,i,j)=scr1(k,i,j)*dn0(k,i,j)
         vt3dh(k,i,j)=vt3dh(k,i,j)*dn0(k,i,j)
!Vertical eddy diffusivity for TKE
         vt3di(k,i,j)=0.2*vctr9(k)*vctr30(k)*dn0(k,i,j)
!print*,'VHM= ',scr1(k,i,j),'  VKH= ',vt3dh(k,i,j)
      enddo
      if(nstbot.eq.1)then

         tket(2,i,j) =tket2  &
            + 0.5 * wght3 * vctr5(2+1)   &
            + wght1 *  &
             (-(sflux_u(i,j)*(up(2,i,j)+up(2,i-1,j))  &
               +sflux_v(i,j)*(vp(2,i,j)+vp(2,i,j-jd))  &
               +sflux_w(i,j)*wp(2,i,j))/vctr1(2)  &
               +sflux_t(i,j)*g/theta(2,i,j))/dn0(2,i,j) &
            
            - tkep2(2)*sqrt(max(1.e-20,tkep2(2)))/(vctr9(2)*16.6)

      endif
      
   enddo
enddo

return
END SUBROUTINE tkemy

!##############################################################################
Subroutine tkeinit (n1,n2,n3)

use mem_grid
use mem_turb
use ke_coms
use rconstants

implicit none

integer :: n1,n2,n3
integer :: i,j,k
real :: epsmin

!        limit the values to the minimum

if( allocated(turb_g(ngrid)%tkep) ) then
   do j = 1,n3
      do i = 1,n2
         do k = 1,n1
            turb_g(ngrid)%tkep(k,i,j) = max(tkmin,turb_g(ngrid)%tkep(k,i,j))
         enddo
      enddo
   enddo
endif
if( allocated(turb_g(ngrid)%epsp) ) then
   epsmin=(c_eps*(tkmin**1.5))/(vonk*zt(2))
   do j = 1,n3
      do i = 1,n2
         do k = 1,n1
            turb_g(ngrid)%epsp(k,i,j) = max(epsmin,turb_g(ngrid)%epsp(k,i,j))
         enddo
      enddo
   enddo
endif

return
END SUBROUTINE tkeinit
