  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
#ifdef key_wave_soulsby_bottom_stress
   USE wave,     ONLY : wave_uzk1,wave_vzk1,wave_freq,uorb,fwst
#endif
#ifdef key_wave_qeuler_dynamics
   USE wave,     ONLY : tausurf,wave_tke
#endif
   USE cometdyn
   USE comturb
#if defined key_siggen || defined key_gencoord
   USE toolgeom, ONLY : f_hzu,f_hzw
#endif
#ifdef key_MPI_2D
   USE toolmpi,  ONLY : ex_i_rsh,ex_j_rsh
#endif
#ifdef key_wave_qeuler_diag
   USE wavediag, ONLY : diag_wf_tke,diag_wf_prodtke,diag_wf_dissiptke
   USE wavediag, ONLY : diag_wf_nudtkedz
#endif
#ifdef key_obstructions
   USE obstructions,      ONLY : obst_a3d,obst_t,obst_nbvar
#endif
   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
#ifdef key_wave_soulsby_bottom_stress
   REAL(kind=rsh) :: shw,shcl,shc
#endif
   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
#if defined key_siggen || defined key_gencoord
   REAL(kind=rsh),DIMENSION(kmax)   :: hz,hz_w
#endif
   REAL(kind=rsh) :: sqrtcmu
#ifdef key_wave_qeuler_dynamics
   REAL(kind=rsh) :: var1,var2,var3
#endif
#ifdef key_obstructions
   REAL(kind=rsh) :: obst_akkp1,obst_1mak,obst_1ma,obst_1o1mak
#endif
#ifdef key_wave_qeuler_diag
   REAL(kind=rsh),DIMENSION(kmax) :: ectm
#endif

   !!----------------------------------------------------------------------
   !! * 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,                              &
#ifdef key_obstructions
!$OMP&  obst_a3d,obst_t,                                          &
#endif
#ifdef key_wave_soulsby_bottom_stress
!$OMP& wave_uzk1,wave_vzk1,uorb,wave_freq,fwst,                   &
#endif
#ifdef key_wave_qeuler_dynamics
!$OMP& tausurf,wave_tke,                                          &
#endif
#ifdef key_wave_qeuler_diag
!$OMP& diag_wf_tke,diag_wf_prodtke,diag_wf_dissiptke,diag_wf_nudtkedz, &
#endif
#if defined key_siggen || defined key_gencoord
#else
!$OMP& sig,                                                       &
#endif
!$OMP&  h0,xeuv,nz,kz,z0b,uz,vz,pn2,tosx,tosy)                    &
!$OMP& PRIVATE(i,j,k,gamma,sqrtcmu,dtm,hp,cd,kfond,nzms,nzmi,     &
#if defined key_siggen || defined key_gencoord
!$OMP& hz,hz_w,                                                   &
#endif
#ifdef key_wave_qeuler_dynamics
!$OMP& var1,var2,var3,                                            &
#endif
#ifdef key_wave_soulsby_bottom_stress
!$OMP& shcl,shc,shw,                                              &
#endif
#ifdef key_wave_qeuler_diag
!$OMP& ectm,                                                      &
#endif
#ifdef key_obstructions
!$OMP& obst_akkp1,obst_1mak,obst_1ma,obst_1o1mak,                 & 
#endif
!$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

#if defined key_siggen || defined key_gencoord
         ! estimation au centre de la maille
         DO k=1,kmax
           hz(k)=f_hzu(h0(i,j),xeuv(i,j),k,i,j)
         END DO
         DO k=1,kmax-1
           hz_w(k)=f_hzw(h0(i,j),xeuv(i,j),k,i,j)
         END DO
#endif

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

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

#if defined key_siggen || defined key_gencoord
         ! on suppose que le premier point de u est au milieu de la premiere maille.
         !    cd=min((0.4_rsh/log((0.5_rsh*dsigu(k)*hz(k)+z0b(i,j))/z0b(i,j)))**2,0.1_rsh)
         !cd=MIN((0.4_rsh/LOG(0.5_rsh*dsigu(k)*hz(k)/z0b(i,j)))**2,0.1_rsh)
         cd=(0.4_rsh/LOG(0.5_rsh*dsigu(k)*hz(k)/z0b(i,j)+1.0_rsh))**2
         ! toute l equation d ect est multipliee par dsigw(k). Attention au Hz_w qui reste qui est de k aussi.
#else
         ! (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
#endif
#ifdef key_no_bottom_stress
         cd=0.0_rsh
#endif
#ifdef key_wave_soulsby_bottom_stress
         ! bottom stress linked to quasi-eulerian current (=cd*U)
         ! hauteur depuis le fond
#if defined key_siggen || defined key_gencoord
         shcl=(0.4_rsh/log(0.5_rsh*dsigu(k)*hz(k)/z0b(i,j)+1.0_rsh))**2*sqrt(wave_uzk1(i,j)**2+wave_vzk1(i,j)**2)
#else
         shcl=(0.4_rsh/log(hp*(sig(k)+1.0_rsh)/z0b(i,j)+1.0_rsh))**2*sqrt(wave_uzk1(i,j)**2+wave_vzk1(i,j)**2)
#endif
         ! norm of the bottom stress linked to quasi-eulerian current (=cd*U2)
         shc=shcl*sqrt(wave_uzk1(i,j)**2+wave_vzk1(i,j)**2)
         shw=0.0_rsh
         IF (uorb(i,j) /= 0._rsh) THEN
           ! bottom stress linked to waves only
           shw=0.5_rsh*fwst(i,j)*(uorb(i,j))**2
         END IF
!
         cd=1.0e-10
         IF ( (shw+shc)*wave_uzk1(i,j)*wave_vzk1(i,j) /= 0.0_rsh ) &
!        IF ( shw+shc /= 0.0_rsh ) &
           cd=shcl*(1.0_rsh+1.2_rsh*((shw/(shw+shc))**3.2_rsh)) &
           /sqrt(0.5_rsh*(wave_uzk1(i,j)**2+wave_uzk1(i-1,j)**2+wave_vzk1(i,j)**2+wave_vzk1(i,j-1)**2))
!VALtoACB pourquoi cette division par variable moyennees ?
#endif
#ifdef key_wave_yu_slinn_bottom_stress
           cd=max(2.0_rsh/pi*cflin*uorb(i,j) &
           /sqrt(0.5_rsh*(wave_uzk1(i,j)**2+wave_uzk1(i-1,j)**2+wave_vzk1(i,j)**2+wave_vzk1(i,j-1)**2)),1.0e-10)
#endif

         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)
#ifdef key_obstructions
         obst_akkp1 = (obst_a3d(obst_nbvar+2,k,i,j)*dsigw(k+1)+obst_a3d(obst_nbvar+2,k+1,i,j)*dsigw(k))/(2.0_rsh*dsigu(k+1))
         IF (obst_a3d(obst_nbvar+2,k,i,j).EQ.0.0_rsh)THEN
           obst_1ma    = 1.0_rsh
         ELSE
           obst_1ma    = (1.0_rsh-obst_a3d(obst_nbvar+2,k,i,j))
         ENDIF
         IF (obst_akkp1.EQ.0.0_rsh)THEN
           obst_1mak   = 1.0_rsh
           obst_1o1mak = 1.0_rsh
         ELSE
           obst_1mak   = (1.0_rsh-obst_akkp1)
           obst_1o1mak = (1.0_rsh/obst_1mak)
         ENDIF
         nzms = nzms * obst_1mak
         nzmi = nzmi * obst_1ma
#endif

#if defined key_siggen || defined key_gencoord
         c(k)=-dtm*nzms/(hz_w(k)*hz(k+1)*dsigu(k+1))
#ifdef key_obstructions
         c(k) = c(k) * obst_1o1mak
#endif
         bb(k)=dsigw(k)-c(k)+dtm*nz(k,i,j)/(hz_w(k)*hz(k)*dsigu(k)) &
               +dtm*dsigw(k)*gamma*eps(k,i,j)/ect(k,i,j)
         a(k)=0.0_rlg
#else
         c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))
#ifdef key_obstructions
         c(k) = c(k) * obst_1o1mak
#endif
         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
#endif
#if defined key_siggen || defined key_gencoord
         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)/(hz_w(k)*dsigw(k))**2
         enplus=dtm*nz(k,i,j)*kfond/(hz_w(k)*hz(k)*dsigu(k))
#else
         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))
#endif
         dest=0.0_rsh
         IF (icon /=0) dest=kz(k,i,j)*pn2(k,i,j)
           ! dest=MAX(dest , 0.0_rsh)
#ifdef key_wave_qeuler_diag
         diag_wf_prodtke(k,i,j)=prod
         diag_wf_dissiptke(k,i,j) = -eps(k,i,j)*(gamma-1.0_rsh)
#endif
         y(k)=ect(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest-eps(k,i,j)*(gamma-1.0_rsh))+enplus  
#ifdef key_wave_qeuler_dynamics
         y(k)=y(k)+wave_tke(k,i,j)*dtm*dsigw(k)
#endif

#ifdef key_obstructions
         y(k) = y(k) + obst_t(k,i,j)*dtm*dsigw(k)
#endif

         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))
#ifdef key_obstructions
           obst_akkp1 = (obst_a3d(obst_nbvar+2,k,i,j)*dsigw(k+1)+obst_a3d(obst_nbvar+2,k+1,i,j)*dsigw(k))/(2.0_rsh*dsigu(k+1))
           IF (obst_akkp1.EQ.0.0_rsh)THEN
             obst_1mak   = 1.0_rsh
             obst_1o1mak = 1.0_rsh
           ELSE
             obst_1mak   = (1.0_rsh-obst_akkp1)
             obst_1o1mak = (1.0_rsh/obst_1mak)
           ENDIF
           nzms = nzms * obst_1mak
#endif
#if defined key_siggen || defined key_gencoord
           c(k)=-dtm*nzms/(hz_w(k)*hz(k+1)*dsigu(k+1))
#else
           c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))
#endif
#ifdef key_obstructions
           c(k) = c(k) * obst_1o1mak
#endif
           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)
#if defined key_siggen || defined key_gencoord
           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)/(hz_w(k)*dsigw(k))**2
#else
           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
#endif
           dest=0.0_rsh
           IF (icon /=0) dest=kz(k,i,j)*pn2(k,i,j)
             ! dest=MAX(dest , 0.0_rsh)
#ifdef key_wave_qeuler_diag
           diag_wf_prodtke(k,i,j)=prod
           diag_wf_dissiptke(k,i,j) = -eps(k,i,j)*(gamma-1.0_rsh)
#endif
           y(k)=ect(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest-eps(k,i,j)*(gamma-1.0_rsh))
#ifdef key_wave_qeuler_dynamics
           y(k)=y(k)+wave_tke(k,i,j)*dtm*dsigw(k)
#endif
#ifdef key_obstructions
           y(k) = y(k) + obst_t(k,i,j)*dtm*dsigw(k)
#endif
         END DO

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

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

#ifdef key_wave_qeuler_dynamics
! VALtoACB reference pour cela ???
         var1=0.09_rsh
         var2=1.5_rsh
         var3=100.0_rsh
         y(k) = tausurf(i,j)/sqrt(var1)*  &
                  (1.0_rsh+sqrt(var2)*(var1**0.25_rsh)*var3)**(0.66_rsh)
#else
         y(k)=(tos/rhoref)/(cmu**2)
#endif
#ifdef key_obstructions
         y(k) = y(k) + obst_t(k,i,j)*dtm*dsigw(k)
#endif 

         !********************************
         !  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
         !******************************
#ifdef key_wave_qeuler_diag
         ectm(1:kmax)=ect(1:kmax,i,j)
#endif
         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

#ifdef key_wave_qeuler_diag
         diag_wf_tke(1:kmax,i,j)=ect(1:kmax,i,j)
         k=1
         diag_wf_nudtkedz(k,i,j)=(                  c(k)*ect(k+1,i,j) &
                                 -(a(k)+c(k))*ect(k,i,j))/(dtm*dsigw(k))
         diag_wf_dissiptke(k,i,j) =diag_wf_dissiptke(k,i,j)-gamma*eps(k,i,j)/ectm(k)*ect(k,i,j)
         DO k=2,kmax-1
            diag_wf_nudtkedz(k,i,j)=(a(k)*ect(k-1,i,j)+c(k)*ect(k+1,i,j) &
                                  -(a(k)+c(k))*ect(k,i,j))/(dtm*dsigw(k))
            diag_wf_dissiptke(k,i,j)=diag_wf_dissiptke(k,i,j)-gamma*eps(k,i,j)/ectm(k)*ect(k,i,j)
         END DO
#endif

         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
