  SUBROUTINE turbk(xeuv,uz,vz,pn2,ect)
   
!&E------------------------------------------------------------------------
!&E                 ***  ROUTINE turbk  ***
!&E
!&E ** Purpose : compute the vertical viscosity and diffusion coefficients,
!&E              using the k parameterisation
!&E
!&E ** Description : compute the turbulent kinetic energy (ect),
!&E                  estimate nz and kz from ect and from an anlytical
!&E                  formulation of the mixing length (l)
!&E
!&E ** Called by : stepu, stepv
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : h0,xeuv,nz,z0b,uz,vz,eps,ect,pn2,kz,tosx,tosy
!&E
!&E ** Modified variables : ect,eps,nz
!&E
!&E ** History :
!&E       !  2004-09-15
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-03  (T. Duhaut) corrections
!&E       !  2009-04  (A. Mano) MPI2D
!&E       !  2009-06  (T. Duhaut, F. Dumas) evolution of the scheme
!&E       !  2009-06  (F. Dumas) generalized sigma
!&E       !  2011-10  (V. Garnier) OMP
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2011-12  (S. Petton, F. Dumas) Bottom stress : removing of min, max
!&E       !  2013-08  (A-C. Bennis and F. Dumas) wave-current interaction
!&E       !  2014-07  (F. Ganthy) Introduction of obstruction module
!&E       !  2017-09  (F. Ganthy) Turbulence advection (k-epsilon)
!&E       !  2017-10  (F. Ganthy) Minor change for obstructions (remove call to
!&E       !                       obstruction-induced turbulence coefficient)
!&E
!&E-------------------------------------------------------------------------
!! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d


   USE cometdyn
   USE comturb


   USE toolmpi,  ONLY : ex_i_rsh,ex_j_rsh



   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)    :: xeuv
   REAL(kind=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)    :: uz,vz
   REAL(kind=rsh),DIMENSION(kmax,limin  :limax  ,ljmin  :ljmax  ),INTENT(in)    :: pn2
   REAL(kind=rsh),DIMENSION(kmax,limin  :limax  ,ljmin  :ljmax  ),INTENT(inout) :: ect

!! * Local declarations
   INTEGER        :: k,i,j

   REAL(kind=rsh) :: nzms,nzmi,kfond
   REAL(kind=rsh) :: dtm,hp
   REAL(kind=rsh) :: cd,prod,dest,enplus,tos
   REAL(kind=rsh),DIMENSION(kmax)   :: a,bb,c,y,u0,p0
   REAL(kind=rsh),DIMENSION(0:kmax) :: q0

   REAL(kind=rsh) :: sqrtcmu




!!----------------------------------------------------------------------
!! * Executable part

!$OMP PARALLEL DEFAULT(none) SHARED(turb_2eq_option,icon,         &
!$OMP&  ect,eps,cmu,hm,sigmaect,dsigw,dsigu,halfdt,ectmin,        &
!$OMP&  ig,id,rhoref,                                             &
!$OMP&  ljminp2,ljmaxm1,limin,limax,                              &





!$OMP& sig,                                                       &

!$OMP&  h0,xeuv,nz,kz,z0b,uz,vz,pn2,tosx,tosy)                    &
!$OMP& PRIVATE(i,j,k,gamma,sqrtcmu,dtm,hp,cd,kfond,nzms,nzmi,     &





!$OMP&  a,bb,c,y,prod,dest,enplus,tos,q0,p0,u0)

!------------------------------------------------------------------------------
!  coefficients de dispersion horizontale (pour calcul d advection dispersion)
!------------------------------------------------------------------------------
   
! dans l ideal il faut gerer c1, c2 et c3 et gamma devra disparaitre
   IF (turb_2eq_option == 1) gamma = 1.0_rsh
   IF (turb_2eq_option == 2) gamma = 1.0_rsh
   IF (turb_2eq_option == 3) gamma = 1.0_rsh
   IF (turb_2eq_option == 4) gamma = 1.0_rsh

   sqrtcmu=SQRT(cmu)

   dtm=halfdt

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       hp=h0(i,j)+xeuv(i,j)
       IF (hp > hm) THEN



         nz(:,i,j)=nz(:,i,j)/sigmaect

         k=1
!*****************************


! (sig(1)+1.0_rsh)*hp est la distance depuis le fond
!    cd=min((0.4_rsh/log(((sig(1)+1.0_rsh)*hp+z0b(i,j))/z0b(i,j)))**2,0.1_rsh)
!cd=MIN((0.4_rsh/LOG((sig(1)+1.0_rsh)*hp/z0b(i,j)))**2,0.1_rsh)
         cd=(0.4_rsh/LOG((sig(k)+1.0_rsh)*hp/z0b(i,j)+1.0_rsh))**2





         kfond=cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2) &
                  +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2))/(cmu**2)
         nzms=(dsigw(k+1)*nz(k,i,j)+dsigw(k)*nz(k+1,i,j))/(2.0_rsh*dsigu(k+1))
         nzmi=nz(k,i,j)



         c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))

         bb(k)=REAL(dsigw(k),rlg)-c(k)+dtm*nz(k,i,j)/(hp*hp*dsigu(k)) &
               +dtm*dsigw(k)*gamma*eps(k,i,j)/ect(k,i,j)
         a(k)=0.0_rlg


         prod=sigmaect*nz(k,i,j)*0.5_rsh*((vz(k+1,i,j  )-vz(k,i,j  ))**2 &
              +(vz(k+1,i,j-1)-vz(k,i,j-1))**2 &
              +(uz(k+1,i,j  )-uz(k,i,j  ))**2 &
              +(uz(k+1,i-1,j)-uz(k,i-1,j))**2)/(hp*dsigw(k))**2
         enplus=dtm*nz(k,i,j)*kfond/(hp*hp*dsigu(k))

         dest=0.0_rsh
         IF (icon /=0) dest=kz(k,i,j)*pn2(k,i,j)
! dest=MAX(dest , 0.0_rsh)

         y(k)=ect(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest-eps(k,i,j)*(gamma-1.0_rsh))+enplus  




         DO k=2,kmax-1
!*****************************

           nzms=(dsigw(k+1)*nz(k,i,j)+dsigw(k)*nz(k+1,i,j))/(2.0_rsh*dsigu(k+1))


           c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))


           a(k)=c(k-1)
           bb(k)=dsigw(k)-a(k)-c(k)+dtm*dsigw(k)*gamma*eps(k,i,j)/ect(k,i,j)

           prod=sigmaect*nz(k,i,j)*0.5_rsh*((vz(k+1,i,j  )-vz(k,i,j  ))**2 &
                +(vz(k+1,i,j-1)-vz(k,i,j-1))**2 &
                +(uz(k+1,i,j  )-uz(k,i,j  ))**2 &
                +(uz(k+1,i-1,j)-uz(k,i-1,j))**2)/(hp*dsigw(k))**2

           dest=0.0_rsh
           IF (icon /=0) dest=kz(k,i,j)*pn2(k,i,j)
! dest=MAX(dest , 0.0_rsh)

           y(k)=ect(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest-eps(k,i,j)*(gamma-1.0_rsh))


         END DO

!******************
         k=kmax
!******************

         tos=SQRT(tosx(i,j)**2+tosy(i,j)**2)
         a(k)=0.0_rsh
         bb(k)=1.0_rsh


         y(k)=(tos/rhoref)/(cmu**2)



!********************************
!  inversion du systeme 1DV
!********************************

!******************************
!  factorisation LU
!******************************

         q0(1)=-c(1)/bb(1)
         DO k=2,kmax-1
           p0(k)=a(k)*q0(k-1)+bb(k)
           q0(k)=-c(k)/p0(k)
         END DO
         p0(kmax)=a(kmax)*q0(kmax-1)+bb(kmax)

!******************************
!  resolution Lu0=y
!******************************

         u0(1)=y(1)/bb(1)
         DO k=2,kmax
           u0(k)=(y(k)-a(k)*u0(k-1))/p0(k)
         END DO

!******************************
!  resolution Uect=u0
!******************************

         ect(kmax,i,j)=u0(kmax)
         DO k=kmax-1,1,-1
           ect(k,i,j)=q0(k)*ect(k+1,i,j)+u0(k)
         END DO



         ect(:,i,j)=MAX(ect(:,i,j),ectmin)

         nz(:,i,j)=nz(:,i,j)*sigmaect
       END IF ! fin du test hm
     END DO
   END DO
!$OMP END DO

!$OMP END PARALLEL

  END SUBROUTINE turbk
