MODULE TURB

use NUMERICS, only : PROGONKA, IND_STAB_FACT_DB
use NUMERIC_PARAMS, only : vector_length, small_value
use INOUT, only : CHECK_UNIT
use LAKE_DATATYPES, only : ireals, iintegers

contains
SUBROUTINE K_EPSILON(ix, iy, nx, ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, dt, &
& b0, tau_air, tau_wav, tau_i, tau_gr, roughness, fetch)
      
! KT_eq calculates eddy diffusivity (ED) in water coloumn following parameterizations:

! 1) empirical profile of ED;
! 2) Semi-empirical formulations for ED;
! 2) k-epsilon parameterization for ED, including Kolmogorov relation.

use OMP_LIB

use COMPARAMS, only: &
& num_ompthr

use DRIVING_PARAMS , only : &
& nstep_keps, &
& turb_out, &
& Turbpar, &
& path, &
& stabfunc, &
& kepsbc, &
& kwe, &
& M, &
& omp, &
& eos, lindens

use ATMOS, only : &
& uwind, vwind

use ARRAYS_TURB , only : &
& S_integr_positive, &
& S_integr_negative, &
& Gen_integr, &
& Seps_integr_positive, &
& Seps_integr_negative, &
& Geneps_integr, &
& epseps_integr, &
& eps_integr, &
& E_integr, &
& TKE_balance, &
& eps_balance, &
& Seps_integr_positive, &
& Seps_integr_negative, &
& Geneps_integr, &
& epseps_integr, &
& TKE_turb_trans_integr, &
& veg_sw, &
& E_it1, E_it2, &
& E1, E2, &
& eps_it1, eps_it2, &
& eps1, eps2, &
& k2_mid, k2, k2t, &
& E12, eps12, &
& row, row2, &
& k5_mid, k5, &
& L, &
& TF, &
& k3_mid, &
& WU_, WV_, &
& GAMT, GAMU, GAMV, &
& Gen, Gen_seiches, &
& S, &
& F, &
& KT, &
& TKE_turb_trans, &
& Re, &
& C1aup , &
& Ri , Ri_bulk, &
& k_turb_T_flux, &
& H_entrainment, &
& Buoyancy0, signwaveheight, &
& knum, &
& Eseiches

use ARRAYS_GRID, only : &
& ddz, ddz2, ddz05, ddz052, ddz054, &
& dzeta_int, dzeta_05int

use ARRAYS_BATHYM, only : &
& dhw, dhw0, &
& area_int, area_half, l1, h1, &
& bathymwater, vol, botar, ls

use ARRAYS_WATERSTATE, only : &
& Tw1, Tw2, &
& Sal1, Sal2, lamw, lamsal

use ARRAYS, only : &
& um, u2, u1, &
& vm, v2, v1, &
& nstep, &
& time, &
& uv

use PHYS_CONSTANTS, only: &
& row0, roa0, &
& lamw0, &
& cw, &
& g, &
& cw_m_row0, &
& roughness0, niu_wat, &
& alsal

use TURB_CONST

use PHYS_FUNC, only : &
& H13

use WATER_DENSITY, only: &
& water_dens_t_unesco, &
& water_dens_ts, &
& DDENS_DTEMP0, DDENS_DSAL0, &
& SET_DDENS_DTEMP, SET_DDENS_DSAL

use INOUT_PARAMETERS, only : &
& lake_subr_unit_min, &
& lake_subr_unit_max

use NUMERICS, only : STEP

use SEICHES_PARAM, only: SEICHE_ENERGY, C_Deff, cnorm

implicit none

! Input variables

! Reals
real(kind=ireals), intent(in) :: dt
real(kind=ireals), intent(in) :: kor
real(kind=ireals), intent(in) :: h_veg, a_veg, c_veg
real(kind=ireals), intent(in) :: hour
real(kind=ireals), intent(in) :: b0
real(kind=ireals), intent(in) :: tau_air, tau_wav, tau_i, tau_gr
real(kind=ireals), intent(in) :: roughness
real(kind=ireals), intent(in) :: fetch


! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: year, month, day
        
! Local variables     
! Reals
real(kind=ireals) :: C_eps1(M), C_eps2(M), C_eps3(M)
real(kind=ireals) :: CE(M), CEt(M)
real(kind=ireals) :: lam_E(M+1), lam_eps(M+1)

real(kind=ireals) :: a(vector_length)
real(kind=ireals) :: b(vector_length)
real(kind=ireals) :: c(vector_length)
real(kind=ireals) :: d(vector_length)

!real(kind=ireals) :: am(vector_length,2,2)
!real(kind=ireals) :: bm(vector_length,2,2)
!real(kind=ireals) :: cm(vector_length,2,2)
!real(kind=ireals) :: ym(vector_length,2)
!real(kind=ireals) :: dm(vector_length,2)

real(kind=ireals) :: AG(5)

real(kind=ireals) :: dt05
real(kind=ireals) :: wr

real(kind=ireals) :: dhw_keps, dhw0_keps
real(kind=ireals) :: pi
real(kind=ireals) :: ufr
real(kind=ireals) :: dist_surf
real(kind=ireals) :: ACC, ACC2, ACCk
real(kind=ireals) :: dE_it, deps_it
real(kind=ireals) :: xx(1:2), yy(1:2), zz
real(kind=ireals) :: lm
real(kind=ireals) :: ext_lamw
real(kind=ireals) :: month_sec, day_sec, hour_sec
real(kind=ireals) :: al_it
real(kind=ireals) :: maxTKEinput, maxTKEinput_i, maxTKEinput_gr
real(kind=ireals) :: FS, FTKES
real(kind=ireals) :: FB, FTKEB
real(kind=ireals) :: al, alft
real(kind=ireals) :: GAMUN
real(kind=ireals) :: urels, vrels
real(kind=ireals) :: ksurf1(1:2), ksurf2(1:2)
real(kind=ireals) :: Esurf1(1:2), Esurf2(1:2)
real(kind=ireals) :: ceps3
real(kind=ireals) :: vdamp
real(kind=ireals) :: dE_it_max, deps_it_max
real(kind=ireals) :: E_min, eps_min
real(kind=ireals) :: Buoyancy1, dBdz0, z0, bvf2
real(kind=ireals) :: Eseiches_diss

real(kind=ireals) :: tf_integr
     
real(kind=ireals) :: GRADU, GRADV, GRADT, GAMVN
real(kind=ireals) :: TFR
real(kind=ireals) :: COGRTN, COGRUN, COGRVN
real(kind=ireals) :: DTDZH

real(kind=ireals), allocatable :: work(:,:)
real(kind=ireals), allocatable :: rhotemp(:), rhosal(:)

! Integers
integer(kind=iintegers), parameter :: maxiter = 35 !35
integer(kind=iintegers) :: iter, iter_t
integer(kind=iintegers) :: nl
integer(kind=iintegers) :: i, j, k
integer(kind=iintegers) :: keps_coef
! integer(kind=iintegers) :: stabfunc
integer(kind=iintegers) :: time_deriv
integer(kind=iintegers) :: i_entrain, i_entrain_old
integer(kind=iintegers) :: nunit_ = lake_subr_unit_min

integer(kind=iintegers) :: grav_wave_Gill = 0
integer(kind=iintegers) :: seiches_Goudsmit = 0

! Logicals
logical, allocatable :: init_keps(:,:)
logical :: indstab
logical :: ind_bound
logical :: iterat
logical :: firstcall
logical :: cycle_keps
logical :: next_tstep_keps
logical :: smooth
logical :: perdam
logical :: semi_implicit_scheme
logical :: contrgrad

! Characters
character :: tp_name*10
character :: numt*1


data firstcall /.true./

! Externals
real(kind=ireals), external :: DZETA

SAVE

!call TIMEC(0)



firstcall_if : if (firstcall) then

  !$OMP SINGLE

  if (turb_out%par == 1) then
!    write (numt,'(i1)') Turbpar%par
!    open(2114,file=path(1:len_trim(path))//'results/'// &
!    & 'err_progon.dat',  status = 'unknown')
!    open(2115,file=path(1:len_trim(path))// &
!    & 'results/'//'E-eps.dat', &
!    & status = 'unknown')
!    open(2116,file=path(1:len_trim(path))// &
!    & 'results/'//'E-eps1.dat', &
!    & status = 'unknown')
!    open (2117,file=path(1:len_trim(path))//'results/' &
!    & //'turb'//numt//'.dat',      status = 'unknown')
!    open (2118,file=path(1:len_trim(path))//'results/' &
!    & //'temp'//numt//'.dat',      status = 'unknown')
!    open (2119,file=path(1:len_trim(path))//'results/' &
!    & //'uv'//numt//'.dat',        status = 'unknown')
!    open (2120,file=path(1:len_trim(path))//'results/' &
!    & //'KT'//numt//'.dat',        status = 'unknown') 
!    open (2121,file=path(1:len_trim(path))//'results/' &
!    & //'dE1'//numt//'.dat',       status = 'unknown') 
!    open (2122,file=path(1:len_trim(path))//'results/' &
!    & //'dE2'//numt//'.dat',       status = 'unknown') 
    call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,nunit_)
    open (nunit_,file=path(1:len_trim(path))//'results/debug/err_keps_iter.dat', &
    & status = 'unknown')
    write(unit=nunit_,fmt=*) 'Errors of iterative process in k-epsilon solver'
!    open (unit=1, file='test.dat', status='unknown')
  endif
  
  allocate (init_keps(1:nx, 1:ny) )
  init_keps(:,:) = .false.
       
  month_sec   = 30.*24.*60.*60.
  day_sec     = 24*60.*60.
  hour_sec    = 60.*60.

  pi          = 4.*atan(1.e0_ireals)

  AL    = g/row0
  ALFT  = 1.d0 
  
! Parameters of numerical scheme
  semi_implicit_scheme = .true. ! No iterations in k-epsilon, semi-implicit scheme
  al_it   = 0.2 !0.8
  ACC     = 1.d-20 !1.d-20
  ACC2    = 1.d-20 !1.d-20
  ACCk    = 1.d-20 !1.d-20
  knum    = 0.

  z0 = 1.d-2
  
  eps_min = 1.d-12 ! The range for dissipation rate values 10**(-11) - 10**(-6) m**2/s**3 (Wuest and Lorke, 2003)
  E_min   = sqrt(eps_min*1.E-5*lamw0*(cw_m_row0*CEt0)**(-1) ) ! This expression ensures that
                                                         ! eddy diffusivity in the decaying
                                                         ! turbulence (E -> Emin, eps -> eps_min)
                                                         ! is much less than the molecular diffusivity

  dE_it_max = 1.d-11 !1.d-9
  deps_it_max = 1.d-15 !1.d-12

  smooth = .false.
  perdam = .false.
  vdamp  = 1.d-3 !1.d-3
 
  contrgrad = .false.
!  AG(1)   = 0._ireals
!  AG(2)   = 0._ireals
!  AG(3)   = 0._ireals
! AG(4)=1. !0.
! AG(5)=1. !0.

  keps_coef = 1

! keps_coef = 1 - standard empirical coefficients in epsilon-equation
! keps_coef = 2 - the coefficient by (Aupoix et al., 1989) is used
  
! stabfunc = 1  - stability functions CE and CEt are constants
! stabfunc = 2  - stability functions CE and CEt are dependent 
! on shear and buoyancy (Canuto et al., 2001)
! stabfunc = 3  - stability functions CE and CEt are dependent
! only on buoyancy (Galperin et al., 1988)

       
! FRICTION: VEGETATION

  if (h_veg>0) then
    print*, 'Vegetation friction parameterization &
    & is not operational currently: STOP'
    STOP
  endif

  veg_sw = 0.
  do i = 1, M+1
    if (h1-dzeta_int(i)*h1 < h_veg) veg_sw(i) = 1.
  enddo

  !$OMP END SINGLE

endif firstcall_if


dt05 = 0.5d0*dt

!$OMP SINGLE
allocate (work(1:M+1,1:2))
allocate (rhotemp(1:M+1), rhosal(1:M+1))
!$OMP END SINGLE

  !$OMP MASTER
  
! Implementation of the Krank-Nikolson numerical scheme
! for k-epsilon parameterization using simple iterations to
! solve the set of nonlinear finite-difference equations

!  k_turb_T_flux used here from the previous time step

  ! Mixed layer defined from the gradient of heat turbulent flux
  !i_entrain = M
  !do i = M-1, 1, -1
  !  if (abs(k_turb_T_flux(i+1)-k_turb_T_flux(i))/abs(k_turb_T_flux(i+1) + 1.d-20) &
  !  & > 0.2 ) then ! 0.2 - quite arbitrary value 
  !    i_entrain = i
  !    exit
  !  endif
  !enddo

  ! Seiche energy evolution, TKE production due to seiche dissipation
  if (seiches_Goudsmit == 1) then
    call SEICHE_ENERGY &
    & (M,bathymwater,ls,tau_air-tau_wav,sqrt(uwind*uwind+vwind*vwind),vol,dt,Eseiches,Eseiches_diss)
    !A part of seiche energy dissipation feeding TKE (Goudsmit et al. 2002, eg. (24))
    do i = 1, M
      bvf2 = max(AL*(row(i+1) - row(i))/(h1*ddz(i)),0._ireals)
      Gen_seiches(i) = (1. - 10.*sqrt(C_Deff))*bvf2/(row0*cnorm*botar)*&
      & (area_int(i) - area_int(i+1))/area_half(i)*Eseiches_diss
    enddo
  else
    Gen_seiches(1:M) = 0.
  endif

  ! Mixed layer defined from TKE gradient
  i_entrain = minloc(E1(1:M),1)
  !do i = 1, M
  !  if (E1(i) < 1.E-6) then
  !    i_entrain = i
  !    exit
  !  endif
  !enddo

  !do i = i_entrain, 1, -1
  !  if (E1(i)/E1(i_entrain) > 5.E+0_ireals) then
  !    i_entrain = i
  !    exit
  !  endif
  !enddo

  H_entrainment = dzeta_05int(i_entrain) * h1
  !$OMP END MASTER
  !$OMP FLUSH (i_entrain,H_entrainment,work)

  ! Bulk Richardson number for mixed layer
  Ri_bulk = g/row0*(row(i_entrain)-row(1))*(dzeta_int(i_entrain)*h1)/ &
  & ( (u1(i_entrain) - u1(1))**2 + (v1(i_entrain) - v1(1))**2 + small_value)


  ! Significant wave height
  if (kepsbc%par == 2 .or. kepsbc%par == 3) then
    signwaveheight = H13(tau_air,fetch)
  endif

  dhw_keps  = dhw  / real(nstep_keps%par)
  dhw0_keps = dhw0 / real(nstep_keps%par)

  !$OMP DO
  ! This time interpolation ensures correspondence between mean energy
  ! viscous dissipation in Crank-Nickolson scheme for momentum equation and TKE shear production
  ! when semi-implicit scheme is used (semi_implicit_scheme = .true.)
  do i = 1, M+1
    Um(i) = 0.5d0 * (u2(i) + u1(i))
    Vm(i) = 0.5d0 * (v2(i) + v1(i))
  enddo

! The TKE injection from waves breaking is limited
! by the wind kinetic energy input to waves development
  maxTKEinput = tau_wav/row0*sqrt(uwind*uwind + vwind*vwind)
!  maxTKEinput_i = tau_i/row0*sqrt(Um(1)*Um(1)+Vm(1)*Vm(1))
!  maxTKEinput_gr = tau_gr/row0*sqrt(Um(M+1)*Um(M+1)+Vm(M+1)*Vm(M+1))
      
!  DTDZH=(Tw1(2)-Tw1(1))/(h1*ddz(1))
  
  next_tstep_keps = .false.
      
  keps_mode: do while (.not.next_tstep_keps)


  if (.not.init_keps(ix,iy)) then
! Initialization mode of k-epsilon parameterization: time derivatives are neglected  
    time_deriv = 0
  else
! Evolution mode of k-epsilon parameterization: full equations are solved
    time_deriv = 1
  endif
  
  !$OMP DO
  do i = 1, M+1
    E_it1(i)   = E1(i)
    eps_it1(i) = eps1(i)
    k2_mid(i)  = 0._ireals
    k2(i)      = 0._ireals
    k2t(i)     = 0._ireals
  enddo


  iter_t = 0

  cycle_keps = .true.
  
  if (semi_implicit_scheme) then
    iterat = .false.
  else
    iterat = .true.
  endif  

  ! Setting density derivatives on temperature and salinity
  do i = 1, M+1
    rhotemp(i) = SET_DDENS_DTEMP(eos%par,lindens%par,Tw1(i),Sal1(i))
    rhosal(i) = SET_DDENS_DSAL(eos%par,lindens%par,Tw1(i),Sal1(i))
  enddo

!call TIMEC(1)

  iterat_keps: do while (cycle_keps)

! The cycle implements iterations to
! solve the set of nonlinear finite-difference equations
! of k-epsilon parameterization
   

   !$OMP DO
   do i = 1, M+1
     E12(i)   = 0.5d0*(E1(i)   + E_it1(i))
     eps12(i) = 0.5d0*(eps1(i) + eps_it1(i))
   enddo

 
!    E12   =  E_it1
!    eps12 =  eps_it1

! Calculating the stability functions
    !$OMP DO
    do i = 1, M
      if (stabfunc%par == 1) then
         lam_eps(i) = lam_eps0
         lam_E  (i) = lam_E0
         CE     (i) = CE0
         CEt    (i) = CEt0
      else
        work(i,1) = E12(i)*E12(i)/((eps12(i) + ACC2)*(eps12(i) + ACC2))* & 
        &        (rhotemp(i) * & ! Assuming Pr = Sc
        &        ( Tw1(i+1) - Tw1(i) ) + &
        &         alsal*rhosal(i) * &
        &        ( Sal1(i+1) - Sal1(i)) ) / &
        &        (h1*ddz(i)) *AL ! ~ bouyancy production of TKE
        work(i,2) = E12(i)*E12(i)/((eps12(i) + ACC2)*(eps12(i) + ACC2))* &
        & ( (U1(i+1) - U1(i))*(U1(i+1) - U1(i)) + &
        &   (V1(i+1) - V1(i))*(V1(i+1) - V1(i)) ) / &
        &   (h1*h1*ddz(i)*ddz(i)) ! ~ shear production of TKE
        if (stabfunc%par == 2) then
          CE(i)  = CE_CANUTO (work(i,2), work(i,1))
          CEt(i) = CEt_CANUTO(work(i,2), work(i,1))
        elseif (stabfunc%par == 3) then
          CE(i)  = sqrt(2.d0)*CL_K_KL_MODEL*SMOMENT_GALPERIN(work(i,1))
          CEt(i) = sqrt(2.d0)*CL_K_KL_MODEL*SHEAT_GALPERIN(work(i,1))
        endif
        lam_eps(i) = CE(i)/sigmaeps
        lam_E  (i) = CE(i)/sigmaE
      endif
    enddo

! call TIMEC(2)

! Calculating eddy diffusivities for kinetic energy and its dissipation
! in the middle between half levels
    !$OMP DO
    do i = 2, M
!      E_mid=0.5d0*(E12(i-1)+E12(i))
!      eps_mid=0.5d0*(eps12(i-1)+eps12(i))
      work(i,1) = INTERPOLATE_TO_INTEGER_POINT(E12(i-1),E12(i),ddz(i-1),ddz(i)) ! interpolated TKE
      work(i,2) = INTERPOLATE_TO_INTEGER_POINT(eps12(i-1),eps12(i),ddz(i-1),ddz(i)) ! interpolated dissipation
      xx(1) = INTERPOLATE_TO_INTEGER_POINT(lam_eps(i-1),lam_eps(i),ddz(i-1),ddz(i))
      xx(2) = INTERPOLATE_TO_INTEGER_POINT(lam_E(i-1),lam_E(i),ddz(i-1),ddz(i))
      if (iter_t == 0) then
        k2_mid(i) = work(i,1)*work(i,1)/(work(i,2) + ACCk)
      else
        k2_mid(i) = k2_mid(i)*al_it + work(i,1)*work(i,1)/(work(i,2) + ACCk)*(1-al_it)
      endif     
      k5_mid(i) = niu_wat + max(k2_mid(i)*xx(1), min_visc/sigmaeps)
      k3_mid(i) = niu_wat + max(k2_mid(i)*xx(2), min_visc/sigmaE)
    enddo
! call TIMEC(3)
    !$OMP DO
    do i = 1, M
      if (iter_t == 0) then
        k2(i)  = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))
        k2t(i) = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_diff/CEt(i))
      else
        k2(i)  = k2(i) *al_it + &
        & max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))*(1-al_it)
        k2t(i) = k2t(i)*al_it + &
        & max(E12(i)*E12(i)/(eps12(i) + ACCk), min_diff/CEt(i))*(1-al_it)
      endif
      k5(i) = niu_wat + k2(i)*lam_eps(i)
      l(i) = E12(i)**(1.5d0)/(eps12(i) + ACC)
      TF(i) = sqrt(E12(i))/(l(i) + ACC)
    enddo
! call TIMEC(4)
    !$OMP DO
    do i = 2, M-1
      WU_(i) = (E_it1(i+1)-E_it1(i-1))*(U2(i)-U2(i-1))/ &
      & (h1*h1*ddz052(i)*ddz(i-1) )
      WV_(i) = (E_it1(i+1)-E_it1(i-1))*(V2(i)-V2(i-1))/ &
      & (h1*h1*ddz052(i)*ddz(i-1) )
    enddo
! call TIMEC(5) 
 
    if (contrgrad) then
      !$OMP SINGLE ! the loop is not parallelized because the option contrgrad is normally not used      
      do i = 2, M
        GRADT = (Tw1(i+1)-Tw1(i-1))/(h1*ddz2(i-1))
        GRADU = (U2(i)-U2(i-1))/(h1*ddz(i-1) )
        GRADV = (V2(i)-V2(i-1))/(h1*ddz(i-1) )
        GAMUN = - CONUV*(WU_(i+1)-WU_(i-1))/(h1*ddz2(i-1))
        GAMVN = - CONUV*(WV_(i+1)-WV_(i-1))/(h1*ddz2(i-1))
        COGRTN = DTDZH * 0.1
        TFR=E_it1(i)/(l(i)*l(i)+ACC) + acc
        COGRUN = GAMUN/TFR
        COGRVN = GAMVN/TFR 
        if(gradT < 0.) then
          GAMT(i) = - AG(3)*cogrtn
          GAMU(i) = - AG(1)*cogrun
          GAMV(i) = - AG(2)*cogrvn
        else
          GAMT(i) = - AG(3)*(AL*GRADT**2+CON0*TFR*COGRTN) / &
          &            (AL*GRADT+CON0*TFR)
          GAMU(i) = - AG(1)*(AL*(CON2*GRADT+(CON2-1.)*GAMT(i))*GRADU + &
          &            CON1*TFR*COGRUN)/(AL*GRADT+CON1*TFR)
          GAMV(i) = - AG(2)*(AL*(CON2*GRADT+(CON2-1.)*GAMT(i))*GRADV + &
          &            CON1*TFR*COGRVN)/(AL*GRADT+CON1*TFR)
        endif
      enddo
      GAMT(1)=0.
      GAMU(1)=0.
      GAMV(1)=0.
      !$OMP END SINGLE
    else
      !$OMP DO
      do i = 1, M+1
        GAMT(i) = 0.
        GAMU(i) = 0.
        GAMV(i) = 0.
      enddo
    endif

!	call TIMEC(6)

!    zz = OMP_GET_WTIME()
    !$OMP DO
    do i = 1, M
      Gen(i) = ((Um(i+1) - Um(i))*(Um(i+1) - Um(i) + h1*ddz(i)*GAMU(i)) + &
      &        (Vm(i+1) - Vm(i))*(Vm(i+1) - Vm(i) + h1*ddz(i)*GAMV(i))) / &
      &        (h1**2*ddz(i)**2)*CE(i)
!      S(i) = - ( 0.5d0*(row(i+1) + row2(i+1) - row(i) - row2(i)) / &
!      &       (h1*ddz(i)) + GAMT(i))*ALFT*AL*CEt(i)
      S(i) = - ( 0.5d0*( rhotemp(i) * & ! Assuming Pr = Sc
      &        (Tw1(i+1) + Tw2(i+1) - Tw1(i) - Tw2(i)) + &
      &                  alsal*rhosal(i) * &
      &        (Sal1(i+1) + Sal2(i+1) - Sal1(i) - Sal2(i)) ) / &
      &        (h1*ddz(i)) + GAMT(i) )*ALFT*AL*CEt(i)
      F(i) = Gen(i) + S(i)
      ! Adding shear production of turbulence by gravitational waves in stable stratification
      Gen(i) = Gen(i) - grav_wave_Gill*grav_wave_Gill_const*STEP(-S(i))*S(i)*CE(i)/CEt(i) 
    enddo
!    print*, OMP_GET_WTIME() - zz, OMP_GET_THREAD_NUM()

!   Sensitivity test    
!    S(1) = min(S(1), 0._ireals)
    
    if (iter_t == 0) then
      Buoyancy1 = S(1)*k2t(1)
      dBdz0 = (Buoyancy1 - Buoyancy0)/(0.5*h1*ddz(1))
    endif
      
! Solution of the equation for TKE (turbulent kinetic energy)

! Boundary condition at the top
    !$OMP SECTIONS
    !$OMP SECTION
    if (l1 == 0) then
      if (kepsbc%par < 4) then
        FTKES = TKE_FLUX_SHEAR(tau_air,kwe%par,maxTKEinput,kepsbc%par)
      elseif (kepsbc%par == 4) then ! The new boundary condition, 
                                    ! implemented currently only for open water case
        FTKES = TKE_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment)
      endif
    elseif (l1 > 0 .and. l1 <= L0) then
      FTKES =        b0*TKE_FLUX_SHEAR(tau_i,0._ireals,maxTKEinput,1) + &
      &    (1.d0-b0)*TKE_FLUX_SHEAR(tau_air,kwe%par,maxTKEinput,kepsbc%par)
    elseif (l1 > L0) then
      FTKES = TKE_FLUX_SHEAR(tau_i,0._ireals,maxTKEinput,1)
    endif
    xx(1) = 0.5*(area_int(2)/area_half(1)*k3_mid(2)*dt / &
    &  (h1**2*ddz(1)*ddz05(1)) + &
    &  time_deriv*dzeta_05int(1)*dhw_keps/(h1*ddz05(1) ) - &
    &  time_deriv*dhw0_keps/(h1*ddz05(1) ) )
    yy(1) = dt/(h1*ddz(1))
    b(1) = - xx(1)
    c(1) = - (xx(1) + time_deriv*1.d0) - (abs(S(1)) - S(1))/(TF(1) + ACC)*dt05 - TF(1)*dt
    d(1) = - time_deriv*E1(1) - xx(1)*(E1(2) - E1(1)) - &
    & area_int(1)/area_half(1)*yy(1)*FTKES - &
    & k2t(1)*(S(1) + abs(S(1)))*dt05 - dt*(k2(1)*Gen(1) + Gen_seiches(1))!+eps_it1(i)*dt
    

! Boundary condition at the bottom
    !$OMP SECTION       
    FTKEB = - TKE_FLUX_SHEAR(tau_gr,0._ireals,maxTKEinput,1)
    xx(2) = 0.5*(area_int(M)/area_half(M)*k3_mid(M)*dt / &
    & (h1**2*ddz(M)*ddz05(M-1) ) - &
    & time_deriv*dzeta_05int(M)*dhw_keps/(h1*ddz05(M-1) ) + &
    & time_deriv*dhw0_keps/(h1*ddz05(M-1) ) )
    yy(2) = dt/(h1*ddz(M))
    a(M) = - xx(2)
    c(M) = - (xx(2) + time_deriv*1.d0) - (abs(S(M)) - S(M))/(TF(M) + ACC)*dt05 - &
    & TF(M)*dt
    d(M) = - time_deriv*E1(M) - xx(2)*(E1(M-1) - E1(M)) + &
    & area_int(M+1)/area_half(M)*yy(2)*FTKEB - &
    & k2t(M)*(S(M) + abs(S(M)))*dt05 - dt*(k2(M)*Gen(M) + Gen_seiches(M)) !+eps_it1(i)*dt

    !$OMP END SECTIONS

!	call TIMEC(4)

! The coefficients of linear algebraic equations in the internal points of the mesh
    !$OMP DO
    do i = 2, M-1
      a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
      & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
      & area_int(i)/area_half(i)*k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )

      b(i) = -time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) + &
      & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
      & area_int(i+1)/area_half(i)*k3_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )

      c(i) = a(i) + b(i) - time_deriv*1.d0-(abs(S(i))-S(i))/(TF(i)+ACC)*dt05 - TF(i)*dt

      d(i) = -time_deriv*E1(i)-k2t(i)*(S(i)+abs(S(i)))*dt05 - dt*(k2(i)*Gen(i) + Gen_seiches(i))
      d(i) = d(i)-dt*(h1**2*ddz(i)*ddz2(i))**(-1)* &
      & (area_int(i+1)/area_half(i)*k3_mid(i+1)*(E1(i+1)-E1(i)))   
      d(i) = d(i)+dt*(h1**2*ddz(i)*ddz2(i-1))**(-1)* &
      & (area_int(i)/area_half(i)*k3_mid(i)*(E1(i)-E1(i-1)))
      d(i) = d(i)-time_deriv*dzeta_05int(i)*dhw_keps* &
      & (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1)) 
      d(i) = d(i)+time_deriv*dhw0_keps* &
      & (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1)) 
    enddo
    ind_bound = .false.
    !$OMP MASTER    
    call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
    !$OMP END MASTER
    !$OMP FLUSH(indstab) 
    if (indstab .eqv. .false.) then
      !$OMP DO
      do i = 2, M-1
        a(i) = - k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
        b(i) = - k3_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
      enddo
    endif
     
    !$OMP MASTER
    call PROGONKA (a,b,c,d,E_it2,1,M)
    E_it2(M+1) = 0.
                
    call CHECK_MIN_VALUE(E_it2, M, E_min) 

    dE_it = maxval(abs(E_it1-E_it2))
    !$OMP END MASTER
    !$OMP FLUSH (E_it2,dE_it)


! C1 is given following Satyanarayana et al., 1999; Aupoix et al., 1989

    if (keps_coef == 1) then
      !$OMP DO PRIVATE (yy,zz,ceps3)
      do i = 1, M
        yy(1) = 0.5*(1. + sign(1._ireals,S(i)))
        zz    = 0.5*(1. - sign(1._ireals,S(i)))
        ceps3 = ceps3_stable*zz + ceps3_unstable*yy(1)
!        ceps3 = 1.14d0
!         ceps3 < -0.4 should be used for stable stratification (Burchard, 2002)
!         ceps3 = 1. ! following (Kochergin and Sklyar, 1992)
!         ceps3 = 1.14d0 ! Baum and Capony (1992)
!         Other values for ceps3:
!         0.<ceps3<0.29 for stable and ceps3 = 1.44 for unstable conditions (Rodi, 1987);
!         ceps3 >= 1 (Kochergin and Sklyar, 1992)
        C_eps1(i) = ceps1
        C_eps2(i) = ceps2
        C_eps3(i) = ceps3
      enddo
    elseif (keps_coef == 2) then
      !$OMP DO
      do i = 1, M 
        Re(i) = ACC + (2*E_it1(i)/3.)**2/(eps_it1(i)*niu_wat + ACC)
        C1aup(i) = Co/(1.d0 + 0.69d0*(2.d0 - Co)/sqrt(Re(i)))
        C_eps1(i) = C1aup(i)
        C_eps2(i) = C1aup(i)
        C_eps3(i) = C1aup(i)
      enddo  
    endif

!	call TIMEC(5)

! The solution of the equation for dissipation rate

! dist_surf is the distance from the surface, 
! where the b.c. for epsilon are formulated

!    dist_surf = 0._ireals !0.5d0*ddz(1)*h1

! Boundary condition at the top
    !$OMP SECTIONS
    !$OMP SECTION
    if (l1 == 0) then
      if (kepsbc%par < 4) then
!        FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
!        & (roughness0 + dist_surf)**(-2)
!        FS = DISSIPATION_FLUX_SHEAR(tau_air,roughness0)

!       Extrapolation of the turbulent diffusivity to the boundary level, assuming
!       logarithmic profile
        ksurf1(1) = max(niu_wat, CE(1)*k2(1)-kar*sqrt(tau_air/row0)*0.5*ddz(1)*h1)
        
        ! ksurf1(1) = kar*sqrt(tau_air/row0)*z0 !roughness0
        ! ksurf1(1) = CE(1)*k2(1)
        Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)
        FS = DISSIPATION_FLUX_SHEAR(tau_air,ksurf1(1),Esurf1(1),z0, &
        & kwe%par,fetch,maxTKEinput,kepsbc%par,signwaveheight)
      elseif (kepsbc%par == 4) then ! The new boundary condition, 
                                    ! implemented currently only for open water case
        Esurf1(1) = E_it1(1)
        FS = DISSIPATION_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment,dBdz0,Esurf1(1))
      endif
    elseif (l1>0 .and. l1<=L0) then
!      FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
!      & (0.0._ireals + dist_surf)**(-2) !0.01 m - roughness of ice 
!      FS =     b0*DISSIPATION_FLUX_SHEAR(tau_i,1.d-3) + & ! 1.d-3 roughness of ice
!      & (1.d0-b0)*DISSIPATION_FLUX_SHEAR(tau_air,roughness0)

!     Extrapolation of the turbulent diffusivity to the boundary level, assuming
!     logarithmic profile
      ksurf1(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_air/row0)*0.5*ddz(1)*h1)
      ksurf2(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_i/row0)*0.5*ddz(1)*h1)
      
      ! ksurf1(1) = kar*sqrt(tau_air/row0)*z0 !*roughness0
      ! ksurf2(1) = kar*sqrt(tau_i/row0)*z0 !1.d-3 is the roughness of ice bottom
!      ksurf1(1) = CE(1)*k2(1)
!      ksurf2(1) = CE(1)*k2(1)
      Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)
      Esurf2(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf2(1)
      FS =     b0*DISSIPATION_FLUX_SHEAR(tau_i,ksurf2(1),Esurf2(1),z0,0._ireals,fetch,maxTKEinput,1,signwaveheight) + &
      & (1.d0-b0)*DISSIPATION_FLUX_SHEAR(tau_air,ksurf1(1),Esurf1(1),z0, &
      & kwe%par,fetch,maxTKEinput,kepsbc%par,signwaveheight)
    elseif (l1 > L0) then
!      FS = DISSIPATION_FLUX_SHEAR(tau_i,1.d-3)

!     Extrapolation of the turbulent diffusivity to the boundary level, assuming
!     logarithmic profile
      ksurf2(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_i/row0)*0.5*ddz(1)*h1)
      
      ! ksurf2(1) = kar*sqrt(tau_i/row0)*z0 !1.d-3 is the roughness of ice bottom
      ! ksurf2(1) = CE(1)*k2(1)
      Esurf2(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf2(1)
      FS = DISSIPATION_FLUX_SHEAR(tau_i,ksurf2(1),Esurf2(1),z0,0._ireals,fetch,maxTKEinput,1,signwaveheight)
    endif 
    xx(1) = 0.5*(area_int(2)/area_half(1)*k5_mid(2)*dt / &
    & (h1**2*ddz(1)*ddz05(1)) + &
    & time_deriv*dzeta_05int(1)*dhw_keps/(h1*ddz05(1) ) - &
    & time_deriv*dhw0_keps/(h1*ddz05(1) ) )
    yy(1) = dt/(h1*ddz(1))
    b(1) = - xx(1)
    c(1) = - (xx(1) + time_deriv*1.d0) - C_eps2(1)*TF(1)*dt + &
    & C_eps3(1)*(- abs(S(1)) + S(1))/(TF(1) + ACC)*dt05
    d(1) = - time_deriv*eps1(1) - xx(1)*(eps1(2) - eps1(1)) - &
    & area_int(1)/area_half(1)*yy(1)*FS - C_eps3(1) * &
    & (abs(S(1)) + S(1))*TF(1)*k2t(1)*dt05 - &
    & C_eps1(1)*TF(1)*dt*(k2(1)*Gen(1) + Gen_seiches(1)) !+eps_it1(i)*dt
 
!   Boundary condition at the bottom
    !$OMP SECTION
!   dist_surf = 0._ireals ! ddz(M)/2*h1
    
!   FB = - CE0**(0.75d0)*k5(M)* &
!   & E_it1(M)**(1.5d0)/kar*(0.0._ireals + dist_surf)**(-2) 
!   FB = - DISSIPATION_FLUX_SHEAR(tau_gr,1.d-3) ! 1.d-3 m - rougnness of bottom

!   Extrapolation of the turbulent diffusivity to the boundary level, assuming
!   logarithmic profile
    ksurf2(2) = max(niu_wat, CE(M)*k2(M) - kar*sqrt(tau_gr/row0)*0.5*ddz(M)*h1)
 
    ! ksurf2(2) = kar*sqrt(tau_gr/row0)*z0 !1.d-3 is the roughness of bottom
    ! ksurf2(2) = CE(M)*k2(M)
    Esurf2(2) = E_it1(M) - 0.5*FTKEB*sigmaE*h1*ddz(M)/ksurf2(2)
    FB = - DISSIPATION_FLUX_SHEAR(tau_gr,ksurf2(2),Esurf2(2),z0,0._ireals,fetch,maxTKEinput,1,signwaveheight)
    xx(2) = 0.5*(area_int(M)/area_half(M)*k5_mid(M)*dt / &
    & (h1**2*ddz(M)*ddz05(M-1) ) - &
    & time_deriv*dzeta_05int(M)*dhw_keps/(h1*ddz05(M-1) ) + &
    & time_deriv*dhw0_keps/(h1*ddz05(M-1) ))
    yy(2) = dt/(h1*ddz(M))
    a(M) = - xx(2)
    c(M) = - (xx(2) + time_deriv*1.d0) - C_eps2(M)*TF(M)*dt + &
    & C_eps3(M)*( - abs(S(M)) + S(M))/(TF(M) + ACC)*dt05
    d(M) = - time_deriv*eps1(M) - xx(2)*(eps1(M-1) - eps1(M)) + &
    & area_int(M+1)/area_half(M)*yy(2)*FB - C_eps3(M) * &
    & (abs(S(M)) + S(M))*TF(M)*k2t(M)*dt05 - &
    & C_eps1(M)*TF(M)*dt*(k2(M)*Gen(M) + Gen_seiches(M))
    !$OMP END SECTIONS

!	call TIMEC(6)

   !$OMP DO
    do i = 2, M-1
      a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
      & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
      & area_int(i)/area_half(i)*k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )

      b(i) = -time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) + &
      & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
      & area_int(i+1)/area_half(i)*k5_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )

      c(i) = a(i) + b(i) - time_deriv*1.d0 - C_eps2(i)*TF(i)*dt + &
      & C_eps3(i)*(-abs(S(i))+S(i))/(TF(i)+ACC)*dt05

      d(i) = - time_deriv*eps1(i) - C_eps3(i)*(abs(S(i)) + S(i))*TF(i)*k2t(i)*dt05 - &
      & C_eps1(i)*TF(i)*dt*(k2(i)*Gen(i) + Gen_seiches(i))
      d(i) = d(i)-dt*(h1**2*ddz(i)*ddz2(i))**(-1)* &
      & (area_int(i+1)/area_half(i)*k5_mid(i+1)*(eps1(i+1)-eps1(i)))
      d(i) = d(i)+dt*(h1**2*ddz(i)*ddz2(i-1))**(-1)* &
      & (area_int(i)/area_half(i)*k5_mid(i)*(eps1(i)-eps1(i-1)))
      d(i) = d(i)-time_deriv*dzeta_05int(i)*dhw_keps* &
      & (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
      d(i) = d(i)+time_deriv*dhw0_keps* &
      & (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
    enddo
    ind_bound = .false.
    !$OMP MASTER
    call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
    !$OMP END MASTER
    !$OMP FLUSH (indstab)
    if (indstab .eqv. .false.) then
      !$OMP DO
      do i = 2, M-1
        a(i) = - k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1))
        b(i) = - k5_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i))
      enddo
    endif
      
    !$OMP MASTER
    call PROGONKA (a,b,c,d,eps_it2,1,M)
    eps_it2(M+1) = 0.
 
    call CHECK_MIN_VALUE(eps_it2, M, eps_min)

    deps_it = maxval(abs(eps_it1-eps_it2))
    !$OMP END MASTER
    !$OMP FLUSH(eps_it2,deps_it)

!	call TIMEC(7)

!    write(*,*) 'iter_t = ', iter_t

    if (iterat) then
      if ((dE_it > dE_it_max .or. deps_it > deps_it_max) .and. iter_t < maxiter) then
        !$OMP DO
        do i = 1, M+1
          E_it1(i)   = E_it1(i)   * al_it + E_it2(i)   * (1-al_it)
          eps_it1(i) = eps_it1(i) * al_it + eps_it2(i) * (1-al_it)
        enddo
        !$OMP MASTER
        iter = iter + 1
        iter_t = iter_t + 1
        if (iter == 8) iter = 0
        !$OMP END MASTER
        !$OMP FLUSH (iter,iter_t)
      else
        if (dE_it < 1.E-6 .and. deps_it < 1.E-7) then
          !$OMP DO
          do i = 1, M+1
            E2(i) = E_it2(i)
            eps2(i) = eps_it2(i)
          enddo
          iter = 0
          iter_t = 0
          nl = 2
          cycle_keps = .false.
        else
         !$OMP DO
          do i = 1, M+1
            eps_it1(i) = eps1(i)
            E_it1(i) = E1(i)
          enddo
          iterat = .false.
          iter_t = 0
          nl = 1
        endif
      endif
    else
      !$OMP DO
      do i = 1, M+1
        E2(i) = E_it2(i)
        eps2(i) = eps_it2(i)
      enddo
!      E2   = E_it1   * al_it + E_it2   * (1-al_it)
!      eps2 = eps_it1 * al_it + eps_it2 * (1-al_it)
      iter = 0
      iter_t = 0
      cycle_keps = .false.
    endif
!    if (dE_it<1.E-10.or.deps_it<1.E-13) nl=3
              
  enddo iterat_keps

!  call TIMEC(8)

  
  if (.not.init_keps(ix,iy)) then
! K-epsilon solver is implemented in initialization mode, i.e.
! to obtain initial profiles of E1 and eps1. 
! Time derivatives in k-epslilon equations are neglected in this case.
    !$OMP DO
    do i = 1, M+1
      E1(i) = E2(i)
      eps1(i) = eps2(i)
    enddo
    init_keps(ix,iy) = .true.
  else
! The case when k-epsilon solver is implemented in evolution mode, i.e.
! E2 and epsilon2 are obtained at the next time step.
    next_tstep_keps = .true.
  endif

  enddo keps_mode

    
! Smoothing of the k and epsilon profiles
  !$OMP MASTER
  if (smooth) then
    call VSMOP_LAKE(E2,E2,lbound(E2,1),ubound(E2,1),1,M,vdamp,perdam)
    call CHECK_MIN_VALUE(E2, M, E_min)
    call VSMOP_LAKE(eps2,eps2,lbound(eps2,1),ubound(eps2,1),1,M,vdamp,perdam)
    call CHECK_MIN_VALUE(eps2, M, eps_min)
  endif
  !$OMP END MASTER
  !$OMP FLUSH (E2,eps2)

  !$OMP DO
  do i = 1, M+1
    E12(i) = 0.5d0*(E1(i) + E2(i))
    eps12(i) = 0.5d0*(eps1(i) + eps2(i))
  enddo

  if (stabfunc%par /= 1) then
    do i = 1, M
        work(i,1) = E2(i)*E2(i)/((eps2(i) + ACC2)*(eps2(i) + ACC2))* & 
        &        (rhotemp(i) * & ! Assuming Pr = Sc
        &        ( Tw2(i+1) - Tw2(i) ) + &
        &         alsal*rhosal(i) * &
        &        ( Sal2(i+1) - Sal2(i)) ) / &
        &        (h1*ddz(i)) *AL ! ~ bouyancy production of TKE
        work(i,2) = E2(i)*E2(i)/((eps2(i) + ACC2)*(eps2(i) + ACC2))* &
        & ( (U2(i+1) - U2(i))*(U2(i+1) - U2(i)) + &
        &   (V2(i+1) - V2(i))*(V2(i+1) - V2(i)) ) / &
        &   (h1*h1*ddz(i)*ddz(i)) ! ~ shear production of TKE
        if (stabfunc%par == 2) then
          CEt(i) = CEt_CANUTO(work(i,2), work(i,1))
        elseif (stabfunc%par == 3) then
          CEt(i) = sqrt(2.d0)*CL_K_KL_MODEL*SHEAT_GALPERIN(work(i,1))
        endif
    enddo
  endif

  !$OMP DO
  do i = 1, M
    KT(i) = lam_T*max(CEt(i)*E2(i)**2/(eps2(i) + ACCk),min_diff)*cw_m_row0 !E2, eps2
  enddo

! Diagnostic calculation
  !Turbulent transport of TKE
  !$OMP DO
  do i = 2, M-1
    TKE_turb_trans(i) = 1.d0/(area_half(i)*h1*h1*ddz(i))* &
    & (area_int(i+1)*k3_mid(i+1)*0.5*(E2(i+1) + E1(i+1) - E2(i  ) - E1(i) ) / ddz05(i)  - &
    & area_int (i)  *k3_mid(i)  *0.5*(E2(i)   + E1(i)   - E2(i-1) - E1(i-1)) / ddz05(i-1) )
  enddo

  TKE_turb_trans(1) = 1.d0/(area_half(1)*h1*h1*ddz(1)) * &
  & (area_int(2)*k3_mid(2)*0.5*(E2(2) + E1(2) - E2(1) - E1(1)) / ddz05(1) + &
  &  area_int(1)*FTKES*h1 )
  TKE_turb_trans(M) = 1.d0/(area_half(M)*h1*h1*ddz(M))* &
  & ( - area_int(M+1)*FTKEB*h1  - &
  & area_int (M)  *k3_mid(M)*0.5 *(E2(M) + E1(M) - E2(M-1) - E1(M-1)) / ddz05(M-1) )

! Diagnostic calculations
  !$OMP DO
  do i = 1, M
    S  (i) = S  (i)*k2t(i)
    Gen(i) = Gen(i)*k2(i)
  enddo

  
  S_integr_positive = 0._ireals
  S_integr_negative = 0._ireals
  Gen_integr = 0._ireals
  TKE_turb_trans_integr = 0._ireals
  
  Seps_integr_positive = 0._ireals
  Seps_integr_negative = 0._ireals
  Geneps_integr = 0._ireals
  epseps_integr = 0._ireals
  
  eps_integr = 0._ireals
  E_integr = 0._ireals
  TF_integr = 0._ireals

  !$OMP SINGLE
  
  do i = 1, i_entrain !M
    xx(1) = ddz(i)*h1
    yy(1) = 0.5*(1. + sign(1._ireals,S(i)))
    zz = 0.5*(1. - sign(1._ireals,S(i)))

    S_integr_positive = S_integr_positive + yy(1)*S(i)*xx(1)
    Seps_integr_positive = Seps_integr_positive + &
    & yy(1)*S(i)*C_eps3(i)*eps2(i)/E2(i)*xx(1)

    S_integr_negative = S_integr_negative + zz*S(i)*xx(1)
    Seps_integr_negative = Seps_integr_negative + &
    & zz*S(i)*C_eps3(i)*eps2(i)/E2(i)*xx(1)

    Gen_integr = Gen_integr + Gen(i)*xx(1)

    TKE_turb_trans_integr = TKE_turb_trans_integr + TKE_turb_trans(i)*xx(1)

    Geneps_integr = Geneps_integr + Gen(i)*C_eps1(i)*eps2(i)/E2(i)*xx(1)
    epseps_integr = epseps_integr + eps2(i)*C_eps2(i)*eps2(i)/E2(i)*xx(1)
    eps_integr = eps_integr + TF(i)*E2(i)*xx(1)
    E_integr = E_integr + E2(i)*xx(1)
    TF_integr = TF_integr + TF(i)*xx(1)
  enddo
  TKE_balance = Gen_integr + S_integr_positive + S_integr_negative - eps_integr
  eps_balance = Geneps_integr + Seps_integr_positive + Seps_integr_negative - epseps_integr
  
  ! Averaging
  S_integr_positive = S_integr_positive/H_entrainment
  S_integr_negative = S_integr_negative/H_entrainment
  Gen_integr = Gen_integr/H_entrainment
  TKE_turb_trans_integr = TKE_turb_trans_integr/H_entrainment

  Seps_integr_positive = Seps_integr_positive/H_entrainment 
  Seps_integr_negative = Seps_integr_negative/H_entrainment 
  Geneps_integr = Geneps_integr/H_entrainment 
  epseps_integr = epseps_integr/H_entrainment 
  
  eps_integr = eps_integr/H_entrainment 
  E_integr = E_integr/H_entrainment 
  TF_integr = TF_integr/H_entrainment 

  TKE_balance = TKE_balance/H_entrainment 
  eps_balance = eps_balance/H_entrainment   

  ! Debug output
  if (firstcall) then
    i_entrain_old = i_entrain
  else
!    if (i_entrain /= i_entrain_old) then
!      write (1,'(i10,f12.2,6e15.5)') nstep, H_entrainment, eps_integr, TF_integr, FS, &
!      & Seps_integr_positive, Seps_integr_negative, epseps_integr
!      i_entrain_old = i_entrain
!    endif
  endif
  ! End debug output
  
!  print*, 'S+ = ', S_integr_positive, 'S- = ', S_integr_negative
!  print*, 'eps_integr = ', eps_integr

!  do i = 2, M-1
!    work(i,1) = dhw_keps/(dt*h1)*dzeta_05int(i)*(E2(i+1)-E2(i-1))/ &
!    & (0.5d0*ddz(i-1)+ddz(i)+0.5d0*ddz(i+1))/Buoyancy0
!    work(i,2) = dhw0_keps/(dt*h1)*(E2(i+1)-E2(i-1))/ &
!    & (0.5d0*ddz(i-1)+ddz(i)+0.5d0*ddz(i+1))/Buoyancy0
!  enddo

! Output
  if (turb_out%par == 1) then
!  do i=1,M+1
!    write (2115,7) time/hour_sec+12.-3*24, dzeta_05int(i)*h1, &
!    & E2(i)*1.E9,eps2(i)*1.E9, KT(i),tau_air*1.E3,Tw1(i)
!  enddo
!  write (2116,8) time/hour_sec+12.-3*24, E2(1)*1.E9, &
!  & eps2(1)*1.E9, KT(1)
    if ((dE_it > 1.d-7 .or. deps_it > 1.d-9).and.(.not.semi_implicit_scheme)) then
      write (nunit_, *) nstep, nl, dE_it, deps_it
    endif
  endif

7 format (f10.4,f8.2,2f12.1,3f8.2)
8 format (f10.4,2f12.1,f8.2)      

  deallocate (work)
  deallocate (rhotemp, rhosal)
                  
  !$OMP END SINGLE
    
  !$OMP DO
  do i = 1, M+1     
    E1(i) = E2(i)
    eps1(i) = eps2(i)
!    write(*,*) 'hihi', i, OMP_GET_THREAD_NUM()
  enddo
  
! write(*,*) 'K_EPSILON finished', nstep  
 !STOP 
     
firstcall = .false.

END SUBROUTINE K_EPSILON


SUBROUTINE CHECK_MIN_VALUE(E,N,threshold)
implicit none

! Input variables
! Integers
integer(kind=iintegers), intent(in) :: N
! Reals
real(kind=ireals), intent(in) :: threshold

! Input/output variables
real(kind=ireals), intent(inout) :: E(1:N)

integer(kind=iintegers) :: i ! Loop index

do i = 1, N 
  E(i) = max(E(i),threshold)
enddo
END SUBROUTINE CHECK_MIN_VALUE


real(kind=ireals) FUNCTION CE_CANUTO(lambda_M, lambda_N)

! The function CE_CANUTO calculates stability function for momentum
! according to Canuto et al., 2001

use TURB_CONST, only: &
& CEcoef01, &
& CEcoef02, &
& CEcoef03, &

& CEcoef11, &
& CEcoef12, &
& CEcoef13, &
& CEcoef14, &
& CEcoef15, &
& CEcoef16

implicit none

! Upper limit for stability function for momentum (Galperin et al., 1988)
real(kind=ireals), parameter :: upper_limit = 0.46d0 
real(kind=ireals), parameter :: small_number = 1.d-20 
real(kind=ireals), parameter :: lambda_N_low_limit = - 5.9d0

real(kind=ireals), intent(in) :: lambda_M 
real(kind=ireals), intent(in) :: lambda_N

real(kind=ireals) :: lambda_N1

! Putting lower limit for lambda_N, that is defined from the condition
! d(CE_CANUTO)/d(lambda_N)<0

lambda_N1 = max(lambda_N_low_limit,lambda_N)

CE_CANUTO = &
& (CEcoef01 + CEcoef02*lambda_N1 - CEcoef03*lambda_M)/ &
& (CEcoef11 + CEcoef12*lambda_N1 + CEcoef13*lambda_M + &
&  CEcoef14*lambda_N1*lambda_N1 + CEcoef15*lambda_N1*lambda_M - & 
&  CEcoef16*lambda_M*lambda_M)

CE_CANUTO = max(small_number,CE_CANUTO)
CE_CANUTO = min(upper_limit,CE_CANUTO)

END FUNCTION CE_CANUTO



real(kind=ireals) FUNCTION CEt_CANUTO(lambda_M, lambda_N)

! The function CEt_CANUTO calculates stability function for scalars
! according to Canuto et al., 2001

use TURB_CONST , only: &
& CEtcoef01, &
& CEtcoef02, &
& CEtcoef03, &

& CEcoef11, &
& CEcoef12, &
& CEcoef13, &
& CEcoef14, &
& CEcoef15, &
& CEcoef16

implicit none
! Upper limit for stability function for scalars (Galperin et al., 1988)
real(kind=ireals), parameter:: upper_limit = 0.61d0
real(kind=ireals), parameter:: small_number = 1.d-20 
real(kind=ireals), parameter:: lambda_N_low_limit = - 5.9d0   

real(kind=ireals), intent(in):: lambda_M
real(kind=ireals), intent(in):: lambda_N
 
real(kind=ireals) lambda_N1
 
! Putting lower limit for lambda_N, that is defined from the condition
! d(CEt_CANUTO)/d(lambda_N)<0 
lambda_N1 = max(lambda_N,lambda_N_low_limit)

CEt_CANUTO = &
& (CEtcoef01 + CEtcoef02*lambda_N1 + CEtcoef03*lambda_M) / &
& (CEcoef11 + CEcoef12*lambda_N1 + CEcoef13*lambda_M + &
& CEcoef14*lambda_N1**2 + CEcoef15*lambda_N1*lambda_M - &
& CEcoef16*lambda_M**2)

CEt_CANUTO = max(small_number,CEt_CANUTO)
CEt_CANUTO = min(upper_limit,CEt_CANUTO)
 
END FUNCTION CEt_CANUTO
      
      
real(kind=ireals) FUNCTION SMOMENT_GALPERIN(lambda_N)
 
! The function SMOMENT_GALPERIN calculates the equilibrium stability function
! for momentum according to Galperin et al., 1988
 
use TURB_CONST!, only: &
!& A1, &
!& B1, &
!& C1, &
!& A2, &
!& B2, &
!& CL_K_KL_MODEL 

implicit none

! Upper limit for stability function (Galperin et al., 1988) 
! real(kind=ireals), parameter:: upper_limit = 0.46d0 ! This value is in terms of CE
real(kind=ireals), parameter:: small_number = 1.d-2
! Below this limit for lambda_N stability function has negative values    
!real(kind=ireals), parameter:: lambda_N_lower_limit = - 1.8d0 
real(kind=ireals), parameter:: GH_MAX = (A2*(12.d0*A1 + B1 + 3.d0*B2))**(-1) ! Galperin et al., 1988
real(kind=ireals), parameter:: GH_MIN = - 0.28d0 ! Galperin et al., 1988
 
real(kind=ireals), intent(in):: lambda_N

real(kind=ireals) :: GH
 
GH = - lambda_N ! According to Mellor and Yamada, 1982 GH has the opposite sign to lambda_N
! The scaling of lambda_N according to Mellor and Yamada, 1982 
GH = CL_K_KL_MODEL**2*0.5d0*GH
! Implementing limitations Galperin et al., 1988
GH = min(GH,GH_MAX)
GH = max(GH,GH_MIN)
 
SMOMENT_GALPERIN = 1.d0 - 3.d0 * C1 - 6.d0 * A1 / B1 - &
& 3.d0 * A2 * GH * ( (B2 - 3.d0*A2) * (1.d0 - 6.d0*A1 / B1) - &
& 3.d0 * C1 * (B2 + 6.d0 * A1) )
SMOMENT_GALPERIN = SMOMENT_GALPERIN / &
& ( (1.d0 - 3.d0 * A2 * GH * (6.d0 * A1 + B2) ) * &
&   (1.d0 - 9.d0 * A1 * A2 * GH) )
SMOMENT_GALPERIN = SMOMENT_GALPERIN * A1
 
SMOMENT_GALPERIN = max(small_number,SMOMENT_GALPERIN)
!SMOMENT_GALPERIN = min(upper_limit, SMOMENT_GALPERIN)

END FUNCTION SMOMENT_GALPERIN
      
      
real(kind=ireals) FUNCTION SHEAT_GALPERIN(lambda_N)
 
! The function SHEAT_GALPERIN calculates the equilibrium stability function
! for heat according to Galperin et al., 1988
 
use TURB_CONST!, only: &
!& A1, &
!& B1, &
!& C1, &
!& A2, &
!& B2, &
!& CL_K_KL_MODEL

implicit none

! Upper limit for stability function (Galperin et al., 1988)
! real(kind=ireals), parameter:: upper_limit = 0.61d0 ! This value is in terms of CEt
real(kind=ireals), parameter:: small_number = 1.d-2 
! Below this limit for lambda_N stability function has negative values
! real(kind=ireals), parameter:: lambda_N_lower_limit = - 1.8d0 
real(kind=ireals), parameter:: GH_MAX = (A2*(12.d0*A1 + B1 + 3.d0*B2))**(-1) ! Galperin et al., 1988
real(kind=ireals), parameter:: GH_MIN = - 0.28d0 ! Galperin et al., 1988

real(kind=ireals), intent(in):: lambda_N
 
real(kind=ireals) :: GH
 
! lambda_N1 = - max(lambda_N,lambda_N_lower_limit)
! The scaling of lambda_N according to Mellor and Yamada, 1982
GH = - lambda_N ! According to Mellor and Yamada, 1982 GH has the opposite sign to lambda_N
! The scaling of lambda_N according to Mellor and Yamada, 1982 
GH = CL_K_KL_MODEL**2*0.5d0*GH
! Implementing limitations Galperin et al., 1988
GH = min(GH,GH_MAX)
GH = max(GH,GH_MIN)
 
SHEAT_GALPERIN = 1.d0 - 6.d0 * A1 / B1
SHEAT_GALPERIN = SHEAT_GALPERIN / &
& (1.d0 - 3.d0 * A2 * GH * (6.d0 * A1 + B2) )
SHEAT_GALPERIN = SHEAT_GALPERIN * A2
 
SHEAT_GALPERIN = max(small_number,SHEAT_GALPERIN)
! SHEAT_GALPERIN = min(upper_limit, SHEAT_GALPERIN)
 
END FUNCTION SHEAT_GALPERIN


SUBROUTINE VSMOP_LAKE(f,fs,nmin,nmax,k0,k1,gama,perdam)

! only the perturbation is smoothed (if perdam is .true.)
!
! 4th-order vertical smoothing. 2nd order for grid-points adjacent
! to the boundary. no smoothing for boundary points.
!
! the second order coeficient is calculated from the fourth order one,
! imposing similar damping for 2 grid-length waves
!
! fourth order damping: response function
! r(lx,ly=)1-16*gama*sin(pi*dx/lx)**4

implicit none

! Input variables
! Integers
integer(kind=iintegers), intent(in) :: nmin
integer(kind=iintegers), intent(in) :: nmax
integer(kind=iintegers), intent(in) :: k0
integer(kind=iintegers), intent(in) :: k1

! Reals
real(kind=ireals), intent(in) :: fs(nmin:nmax)
real(kind=ireals), intent(in) :: gama

! Logicals
logical, intent(in) :: perdam

! Input/output variables
! Reals
real(kind=ireals), intent(inout) :: f(nmin:nmax)

! Local variables
! Reals
real(kind=ireals) :: work(nmin:nmax)
real(kind=ireals) :: um6g
real(kind=ireals) :: gasc
real(kind=ireals) :: gasc05
real(kind=ireals) :: umgasc

! Integers
integer(kind=iintegers) :: i ! Loop index

um6g = 1. - 6.*gama
gasc = 8.*gama
gasc05 = 0.5*gasc
umgasc = 1. - gasc

if (perdam) then
  do i = k0, k1
    f(i) = f(i) - fs(i)
  enddo
endif

do i = k0+2, k1-2
  work(i) = um6g*f(i) - gama*(f(i+2) + f(i-2) - 4.*(f(i+1) + f(i-1)) )
enddo

do i = k0+1, k1-1, k1-k0-2
  work(i) = umgasc*f(i) + gasc05*(f(i-1) + f(i+1))
enddo

do i = k0+1, k1-1
  f(i) = work(i)
enddo

if (perdam) then
  do i = k0, k1
    f(i) = f(i) + fs(i)
  enddo
endif

RETURN
END SUBROUTINE VSMOP_LAKE


FUNCTION INTERPOLATE_TO_INTEGER_POINT(fim05,fip05,ddzim1,ddzi)

! Interpolates linearly the grid function 
! from two half-integer points to integer point between them

implicit none

real(kind=ireals) :: INTERPOLATE_TO_INTEGER_POINT

! Input variables
real(kind=ireals), intent(in) :: fim05  ! The function at (i-0.5)-th level
real(kind=ireals), intent(in) :: fip05  ! The function at (i+0.5)-th level
real(kind=ireals), intent(in) :: ddzim1 ! The thickness of (i-1)-th layer
real(kind=ireals), intent(in) :: ddzi   ! The thickness of i-th layer

INTERPOLATE_TO_INTEGER_POINT = fim05 + (fip05 - fim05)*ddzim1/(ddzim1 + ddzi)

END FUNCTION INTERPOLATE_TO_INTEGER_POINT


FUNCTION TKE_FLUX_BUOY(Buoyancy_flux,roughness_length,H_entrainment)

! The function TKE_FLUX_BUOY calculates the flux of turbulent kinetic energy
! at the surface in the case of free convection, i.e. no wind stress but buoyancy flux
! given at the surface (Stepanenko and Lykosov, 2009).

use TURB_CONST

real(kind=ireals) :: TKE_FLUX_BUOY

! Input variables
real(kind=ireals), intent(in) :: Buoyancy_flux
real(kind=ireals), intent(in) :: roughness_length
real(kind=ireals), intent(in) :: H_entrainment

! Local variables
real(kind=ireals), save :: const

logical, save :: firstcall = .true.

if (firstcall) then
  const = 2.d0/3.d0*CE0*kar*kar/(sigmaE*CL*CL)
endif

!TKE_FLUX_BUOY = - const*Buoyancy_flux*roughness_length * &
!&               (1.d0-2.d0*roughness_length/H_entrainment)

TKE_FLUX_BUOY = 0._ireals

if (firstcall) firstcall = .false.
END FUNCTION TKE_FLUX_BUOY


FUNCTION DISSIPATION_FLUX_BUOY(Buoyancy_flux,roughness_length,H_entrainment, dBdz0, TKE0)

! The function DISSIPATION_FLUX_BUOY calculates the flux of turbulent kinetic energy
! dissipation rate at the surface in the case of free convection, 
! i.e. no wind stress but buoyancy flux given at the surface (Stepanenko and Lykosov, 2009).

use TURB_CONST

real(kind=ireals) :: DISSIPATION_FLUX_BUOY

! Input variables
real(kind=ireals), intent(in) :: Buoyancy_flux
real(kind=ireals), intent(in) :: roughness_length
real(kind=ireals), intent(in) :: H_entrainment
real(kind=ireals), intent(in) :: dBdz0 ! The vertical gradient of Buoyancy flux at the water surface
real(kind=ireals), intent(in) :: TKE0 ! The turbulent kinetic energy at the water surface

! Local variables
real(kind=ireals), save :: const
real(kind=ireals) :: wstar

logical, save :: firstcall = .true.

if (firstcall) then
  const = CE0*(kar/CL)**(4.d0/3.d0)/sigmaeps
endif

! Flux obtained neglecting diffusion term in TKE equation
! DISSIPATION_FLUX_BUOY = const*(Buoyancy_flux*roughness_length)**(4.d0/3.d0) / &
! &                       H_entrainment

! Flux obtained including diffusion term in TKE equation
!wstar = (Buoyancy_flux*H_entrainment)**(1.d0/3.d0)
!DISSIPATION_FLUX_BUOY = CE0*C_wstar*C_wstar*Buoyancy_flux*wstar/sigmaeps

DISSIPATION_FLUX_BUOY = - CE0*TKE0*TKE0*dBdz0/((Buoyancy_flux+small_value)*sigmaeps)

if (firstcall) firstcall = .false.
END FUNCTION DISSIPATION_FLUX_BUOY


FUNCTION TKE_FLUX_SHEAR(momentum_flux,kwe,maxTKEinput,nbc)

! The function TKE_FLUX_SHEAR calculates the flux of turbulent kinetic energy
! at the surface in the case of shear lorarithmic flow, i.e. no buoyancy flux but wind stress
! given at the surface (Burchard and Petersen, 1999).

use PHYS_CONSTANTS, only : &
& row0

use TURB_CONST, alpham_=>alpham!, only : &
!& kar, sigmaE, &
!& CE0, kar_tilde

implicit none

real(kind=ireals) :: TKE_FLUX_SHEAR

! Input variables
real(kind=ireals), intent(in) :: momentum_flux
real(kind=ireals), intent(in) :: kwe
real(kind=ireals), intent(in) :: maxTKEinput

integer(kind=iintegers), intent(in) :: nbc

!Local variables
real(kind=ireals), parameter :: alpham = (1.5*CE0**0.5*sigmaE*kar_tilde**(-2))**(0.5)
real(kind=ireals), parameter :: const0 = (1.5*sigmaE)**0.5*CE0**0.25
real(kind=ireals), parameter :: const = 2./3.*CE0**(-0.5)*kar*const0*alpham/sigmaE

select case (nbc)
  case(1)
    TKE_FLUX_SHEAR = 0._ireals ! for logarithmic profile (shear only)  
  case(2)
    ! TKE_FLUX_SHEAR = min(kwe*(momentum_flux/row0)**1.5d0,maxTKEinput) ! for breaking waves case 
                                                                      ! Craig and Banner (1994)
    TKE_FLUX_SHEAR = kwe*(momentum_flux/row0)**1.5d0
  case(3)
    ! Generalized b.c. for sheared flow with wave breaking (Burchard, 2002)
    TKE_FLUX_SHEAR = min(const*kwe*(momentum_flux/row0)**1.5d0,maxTKEinput)
end select

END FUNCTION TKE_FLUX_SHEAR


FUNCTION DISSIPATION_FLUX_SHEAR(momentum_flux,kturb,Esurf_in,z0_in, &
& kwe,fetch,maxTKEinput,nbc,signwaveheight)

! The function DISSIPATION_FLUX_SHEAR calculates the flux of turbulent kinetic energy
! dissipation rate at the surface in the case of shear logarithmic flow, 
! i.e. no buoyancy flux but wind stress given at the surface (Burchard, 2000).

! Note, that turbulent kinetic energy (TKE) and turbulent diffusivity have to be taken
! at the surface. The TKE may be taken at the first level below the surface, assuming
! d(TKE)/dz = 0, that is satisfied for logarithmic surface layer. 
! The appropriate extrapolation of turbulent diffusivity (kturb) from the first level below
! the surface to the surface itself is using the linear profile of kturb:

! kturb ~ karman_constant*friction_velocity*z

! that is exact one for logarithmic profile.

use PHYS_CONSTANTS, only : &
& row0
use TURB_CONST !, only : sigmaeps, CE0, kar, sigmaE

use PHYS_FUNC, only : &
& H13

implicit none

real(kind=ireals) :: DISSIPATION_FLUX_SHEAR

! Input variables
real(kind=ireals), intent(in) :: momentum_flux
real(kind=ireals), intent(in) :: kturb
real(kind=ireals), intent(in) :: Esurf_in
real(kind=ireals), intent(in) :: z0_in
real(kind=ireals), intent(in) :: kwe
real(kind=ireals), intent(in) :: fetch
real(kind=ireals), intent(in) :: maxTKEinput
real(kind=ireals), intent(in) :: signwaveheight

integer(kind=iintegers), intent(in) :: nbc ! switch for b.c.s

! Local variables
real(kind=ireals), parameter :: small_number = 1.d-20
real(kind=ireals), save :: const
real(kind=ireals), save :: const1, const2, const3, a
real(kind=ireals) :: const4, const5
real(kind=ireals) :: z0s, xx, Esurf

logical, save :: firstcall = .true.

if (firstcall) then
  const = 1.d0/sigmaeps
  const1 = CE0**(0.75d0)/(sigmaeps*kar)
  const2 = const1/kar
  const3 = 1.5*sigmaE*CE0**(0.75d0)/CE0 ! must be CE instead of CE0 in denominator
  a = 1. ! Switch for shear (0 or 1)
endif

z0s = 1.d-2 !signwaveheight*0.85 ! 0.85 - empirical constant

!DISSIPATION_FLUX_SHEAR = (momentum_flux/row0)*(momentum_flux/row0) * &
!& const/z0_in

select case (nbc)
  case(1)
    ! From logarithmic profile (sheared flow without waves breaking)
    Esurf = Esurf_in !momentum_flux/row0/CE0**0.5 ! TKE at the surface
    DISSIPATION_FLUX_SHEAR = const1*kturb*Esurf**(1.5d0) / &
    & (z0_in*z0_in + small_number) ! roughness_length -> z0s?
  case(2)
    ! Dissipation flux in flow with breaking waves (Burchard, 2001)
    const4 = (1.5*sigmaE)**0.5*CE0**0.25*kwe
    ! xx = min(kwe*(momentum_flux/row0)**1.5,maxTKEinput)
    Esurf = Esurf_in !momentum_flux/row0/(CE0**0.5)*const4**(2./3.) ! TKE at the surface
    xx = kwe*(momentum_flux/row0)**1.5
    DISSIPATION_FLUX_SHEAR = const2*(momentum_flux/row0)**0.5*kar*z0s*const4**(1./3.)* &
    & (const3*xx + kar*Esurf**1.5) / (z0s*z0s + small_number)
  case(3)
    ! Generalized b.c. for sheared flow with wave breaking (Burchard, 2002)
    const4 = (1.5*sigmaE)**0.5*CE0**0.25*kwe
    const5 = sigmaeps**(-1)*(a + const4)**(1./3.)*((a + const4)/(z0s + small_number) + alpham*const4)
    DISSIPATION_FLUX_SHEAR = const5*(momentum_flux/row0)**2
end select

if (firstcall) firstcall = .false.
END FUNCTION DISSIPATION_FLUX_SHEAR


SUBROUTINE ED_TEMP_HOSTETLER(wind,frict_vel,zref,phi180,levels,row,ed_temp,N_levels)

! Subroutine ED_TEMP_HOSTETLER calculates the eddy diffusivity for temperature
! in a lake according to Technical Description of the Community Land Model (2004), p. 160

use PHYS_CONSTANTS, only : &
& kappa, & ! Karman constant, n/d
& g ! Acceleration due to gravity, m/s**2

implicit none

! Input variables
real(kind=ireals), intent(in) :: wind ! Module of the wind speed, m/s
real(kind=ireals), intent(in) :: frict_vel ! Friction velocity, m/s
real(kind=ireals), intent(in) :: zref ! The height in the atmosphere, where the wind is measured (computed), m
real(kind=ireals), intent(in) :: phi180 ! The latitude, degrees

integer(kind=iintegers), intent(in) :: N_levels ! The number of levels in water column
real(kind=ireals), intent(in) :: levels(1:N_levels+1) ! The level depths, m
real(kind=ireals), intent(in) :: row(1:N_levels+1) ! Water density, kg/m**3

! Output variables
real(kind=ireals), intent(out) :: ed_temp(1:N_levels) ! The eddy diffusivity for temperature, m**2/s

! Local variables
real(kind=ireals), parameter :: w_star_const = 1.2d-3
real(kind=ireals), parameter :: k_star_const = 6.6, k_star_power_const = -1.84
real(kind=ireals), parameter :: PrandtL0 = 1.
real(kind=ireals), parameter :: Ri_const1 = 37., Ri_const2 = 20.
real(kind=ireals), parameter :: Brunt_Vaisala2_const = 40.
real(kind=ireals), parameter :: z0_water = 1.d-3 ! Water surface roughness, m
real(kind=ireals), parameter :: wind2_min = 1. ! Minimal wind speed at 2 meters
real(kind=ireals), parameter :: pi = 3.14159265 ! pi constant
real(kind=ireals), parameter :: small_value = 1.e-20

real(kind=ireals) :: phi
real(kind=ireals) :: wind2
real(kind=ireals) :: k_star, w_star
real(kind=ireals) :: Brunt_Vaisala2, Richardson
real(kind=ireals) :: level

integer(kind=iintegers) :: i

phi = phi180/180.*pi
wind2 = max(frict_vel/kappa*log(2./z0_water),wind2_min)
w_star = w_star_const*wind2
k_star = k_star_const*wind2**(k_star_power_const)*sqrt(abs(sin(phi)))

do i = 1, N_levels
  level = 0.5*(levels(i) + levels(i+1))
  Brunt_Vaisala2 = 2.*g/(row(i)+row(i+1))*(row(i+1)-row(i))/(levels(i+1)-levels(i))
  Richardson = &
  & (- 1. + sqrt(1. + &
  & max(Brunt_Vaisala2_const*Brunt_Vaisala2*kappa**2*level**2/ &
  & (w_star**2*exp( - 2.*k_star*level) + small_value), -1.d0) &
  & ))/Ri_const2
  ed_temp(i) = kappa*w_star*level*exp(-k_star*level)/ &
  & (PrandtL0*(1.+Ri_const1*Richardson**2) )
enddo

END SUBROUTINE ED_TEMP_HOSTETLER


SUBROUTINE ED_TEMP_HOSTETLER2(wind,frict_vel,zref,phi180,levels,row,ed_temp,N_levels)

! Subroutine ED_TEMP_HOSTETLER calculates the eddy diffusivity for temperature
! in a shallow lake (pond), where bottom and top boundary layers intersect
! in the middle of a water column. Hendersson-Sellers formulations are applied
! for bottom and top boundary layers, excluding Ekman terms.

use PHYS_CONSTANTS, only : &
& kappa, & ! Karman constant, n/d
& g ! Acceleration due to gravity, m/s**2

implicit none

! Input variables
real(kind=ireals), intent(in) :: wind ! Module of the wind speed, m/s
real(kind=ireals), intent(in) :: frict_vel ! Friction velocity, m/s
real(kind=ireals), intent(in) :: zref ! The height in the atmosphere, where the wind is measured (computed), m
real(kind=ireals), intent(in) :: phi180 ! The latitude, degrees

integer(kind=iintegers), intent(in) :: N_levels ! The number of levels in water column
real(kind=ireals), intent(in) :: levels(1:N_levels+1) ! The level depths, m
real(kind=ireals), intent(in) :: row(1:N_levels+1) ! Water density, kg/m**3

! Output variables
real(kind=ireals), intent(out) :: ed_temp(1:N_levels) ! The eddy diffusivity for temperature, m**2/s

! Local variables
real(kind=ireals), parameter :: w_star_const = 1.2d-3
real(kind=ireals), parameter :: k_star_const = 6.6, k_star_power_const = -1.84
real(kind=ireals), parameter :: PrandtL0 = 1.
real(kind=ireals), parameter :: Ri_const1 = 37., Ri_const2 = 20.
real(kind=ireals), parameter :: Brunt_Vaisala2_const = 40.
real(kind=ireals), parameter :: z0_water = 1.d-3 ! Water surface roughness, m
real(kind=ireals), parameter :: wind2_min = 1. ! Minimal wind speed at 2 meters
real(kind=ireals), parameter :: pi = 3.14159265 ! pi constant
real(kind=ireals), parameter :: small_value = 1.e-20

real(kind=ireals) :: phi
real(kind=ireals) :: wind2
real(kind=ireals) :: k_star, w_star
real(kind=ireals) :: Brunt_Vaisala2, Richardson
real(kind=ireals) :: level

integer(kind=iintegers) :: i

phi = phi180/180.*pi
wind2 = max(frict_vel/kappa*log(2./z0_water),wind2_min)
w_star = w_star_const*wind2
!k_star = k_star_const*wind2**(k_star_power_const)*sqrt(abs(sin(phi)))

do i = 1, N_levels
  level = 0.5*(levels(i) + levels(i+1))
  Brunt_Vaisala2 = 2.*g/(row(i)+row(i+1))*(row(i+1)-row(i))/(levels(i+1)-levels(i))
  if ( level <= 0.5*levels(N_levels+1) ) then
    Richardson = &
    & (- 1. + sqrt(1. + &
    & max(Brunt_Vaisala2_const*Brunt_Vaisala2*kappa**2*level**2/ &
    & (w_star**2 + small_value), -1._ireals) ))/Ri_const2
    ed_temp(i) = kappa*w_star*level/(PrandtL0*(1.+Ri_const1*Richardson**2) )
  else
    Richardson = &
    & (- 1. + sqrt(1. + &
    & max(Brunt_Vaisala2_const*Brunt_Vaisala2*kappa**2*(levels(N_levels+1) - level)**2/ &
    & (w_star**2 + small_value), -1._ireals) ))/Ri_const2
    ed_temp(i) = kappa*w_star*(levels(N_levels+1) - level)/(PrandtL0*(1.+Ri_const1*Richardson**2) )
  endif 
enddo

END SUBROUTINE ED_TEMP_HOSTETLER2


END MODULE TURB
