#ifdef key_turbadv
   SUBROUTINE turbpsi(xeadv,xeuv,uz,vz,pn2,psi,ect)
#else
  SUBROUTINE turbpsi(xeuv,uz,vz,pn2,psi,ect)
#endif

   !&E----------------------------------------------------------------------------------------------------
   !&E                 ***  ROUTINE turbpsi  ***
   !&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,z0b,nz,uz,vz,eps,ect,pn2,psi,tosx,tosy
   !&E
   !&E ** Modified variables : psi,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, V. Garnier) MPI2D
   !&E       !  2009-06    (T. Duhaut) Choice of BC (surface and bottom): Flux or Dirichlet (default)
   !&E       !  2009-06    (F. Dumas) generalized sigma coordinates
   !&E       !  2009-12    (T. Duhaut, F. Dumas) evolution
   !&E       !  2011-01    (A.-C. Bennis) correction dtm instead of dt
   !&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       !  2012-02  (A.C Bennis, F. Dumas) Mixing length due to surface waves and bottom rugosity
   !&E                   Warner et al.  (2005) Ocean Modelling : Performance of four turbulence closure models...
   !&E       !  2012-02  (F. Dumas) Initialization of nz, kz at kmax (+ other corrections)
   !&E       !  2013-08  (A-C. Bennis and F. Dumas) wave-current coupling
   !&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_freq,uorb,wave_uzk1,wave_vzk1,fwst
#endif
#ifdef key_wave_qeuler_dynamics
   USE wave,     ONLY : peps,tausurf
#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_psi,diag_wf_prodpsi,diag_wf_dissippsi
   USE wavediag, ONLY : diag_wf_nudpsidz
#endif
#ifdef key_obstructions
   USE obstructions,      ONLY : obst_a3d,obst_t,obst_tau,obst_nbvar
#endif
#ifdef key_turbadv
   USE turbadvobc3dapply
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)    :: xeuv
#ifdef key_turbadv
   REAL(kind=rsh),DIMENSION(limin:limax,ljmin:ljmax),INTENT(in)            :: xeadv
#endif
   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) :: psi,ect

   !! * Local declarations
#ifdef key_wave_qeuler_diag
   REAL(kind=rsh),DIMENSION(kmax) :: psim
#endif
   INTEGER        :: k,i,j,kkmax
#ifdef key_wave_soulsby_bottom_stress
   REAL(kind=rsh) :: shw,shcl,shc
#endif
   REAL(kind=rsh) :: nzms,nzmi,kfond
   REAL(kind=rsh) :: sqrtcmu,dtm,hp
   REAL(kind=rsh) :: cd,prod,dest,enplus
   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) :: kZb,kZs,FpsiS,db,ds
#ifdef key_obstructions
   REAL(KIND=rsh) :: taupsi_free,obst_akkp1,obst_1mak,obst_1ma,obst_1o1mak
#endif

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

!$OMP PARALLEL DEFAULT(none) SHARED(psi,ect,eps,h0,xeuv,uz,vz,nz,kz,pn2,    &
!$OMP& limin,limax,ljminp2,ljmaxm1,ig,id,dsigu,dsigw,                       &
!$OMP& z0b,z0s,turb_2eq_option,cmu,hm,halfdt,sigmapsi,psimin,icon,          &
#ifdef key_obstructions
!$OMP& obst_a3d,obst_t,obst_tau,                                            &
#endif
#ifdef key_wave_soulsby_bottom_stress
!$OMP& wave_freq,uorb,wave_uzk1,wave_vzk1,fwst,                             &
#endif
#ifdef key_wave_qeuler_dynamics
!$OMP& peps,tausurf,                                                        &
#endif
#ifdef key_wave_qeuler_diag
!$OMP& diag_wf_psi,diag_wf_prodpsi,diag_wf_dissippsi,diag_wf_nudpsidz,      &
#endif
#if defined key_siggen || defined key_gencoord
#else
!$OMP& sig,sigw,                                                            &
#endif
!$OMP& c1turb,c2turb,c3mturb,c3pturb)                                       &
!$OMP& PRIVATE(i,j,k,gamma,l_bc_turb,kkmax,sqrtcmu,dtm,hp,enplus,           &
!$OMP& kZb,kZs,cd,kfond,db,ds,fwall,nzms,nzmi,dest,prod,FpsiS,              &
#if defined key_siggen || defined key_gencoord
!$OMP& hz,hz_w,                                                             &
#endif
#ifdef key_wave_soulsby_bottom_stress
!$OMP& shw,shcl,shc,                                                        &
#endif
#ifdef key_obstructions
!$OMP& taupsi_free,obst_akkp1,obst_1mak,obst_1ma,obst_1o1mak,               &
#endif
!$OMP& a,bb,c,y,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 ! k-kl
   IF ( turb_2eq_option == 2) gamma = 1.0_rsh ! k-epsilon
   IF ( turb_2eq_option == 3) gamma = 1.0_rsh ! k-omega
   IF ( turb_2eq_option == 4) gamma = 1.0_rsh ! generic length scale

   l_bc_turb=.true.       ! Von Neumann boundary condition (we prescribe a flux in the middle of the surface or bottom cell)
   ! l_bc_turb=.false.    ! Dirichlet boundary condition (we prescribe the value at the interface (bottom or surface)

   IF (l_bc_turb) THEN
     kkmax=kmax-1
   ELSE
     kkmax=kmax
   END IF

   sqrtcmu =SQRT(cmu)
   fwall=1.0_rsh      ! private variable, so must be initialized

   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

         ! What a mess !!!! grrrr ! cval
         ! for the surface : Dirichlet BC uses kZb nd Neumann BC uses FpsiS
         ! for the bottom : Dirichlet and Neumann BC use the same variable kZb

         ! kappa*z at the bottom 
         IF ( l_bc_turb) THEN          ! Neumann bc (flux) at u position
#if defined key_siggen || defined key_gencoord
           kZb = 0.4_rsh*(0.5_rsh*dsigu(1)*hz(1)+z0b(i,j))
#else
           kZb = 0.4_rsh*((sig(1)+1.0_rsh)*hp+z0b(i,j))
#endif
         ELSE                          ! Dirichlet (value) at the interface (bottom)
           kZb = 0.4_rsh*z0b(i,j)
         END IF

         ! kappa*z at the surface
         kZs = 0.4_rsh*z0s(i,j)

         ! nu/sigmapsi
         nz(:,i,j)=nz(:,i,j)/sigmapsi

         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 quiest 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
         ! hauteur d eau 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
         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))
#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 = epsilon at the bottom (for Dirichlet BC at least and for Neumann BC d epsilon /dz)....
         ! No comment for the choice of the name.... ! cval
         IF ( turb_2eq_option == 1 ) THEN   ! k-kl
           IF ( l_bc_turb) THEN
             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))*0.4_rsh/(cmu**2)
           ELSE
             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))*kZb/(cmu**2)
           END IF
#if defined key_siggen || defined key_gencoord
           db=hz(k)*dsigu(k)
#else
           db=hp*dsigu(k)
#endif
           ds=hp-db
           fwall = 1.0_rsh+1.33_rsh*( psi(k,i,j)*(hp/(db*ds))  &
                                /(ect(k,i,j)*0.4_rsh) )**2

         ELSE IF ( turb_2eq_option == 2 ) THEN   ! k-epsilon
           IF ( l_bc_turb) THEN
             kfond =-0.4_rsh*( (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)))**(1.5_rsh) )/kZb**2
           ELSE
             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)))**(1.5_rsh) )/kZb
           END IF
         ELSE IF ( turb_2eq_option == 3 ) THEN   ! k-omega
           IF ( l_bc_turb) THEN
             kfond =-0.4_rsh*SQRT( 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)) )/kZb**2/(cmu**2)
           ELSE
             kfond = SQRT( 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)) )/kZb/(cmu**2)
           END IF
         ELSE IF ( turb_2eq_option == 4 ) THEN    ! generic length scale
           IF ( l_bc_turb) THEN
             kfond =-0.4_rsh*(2.0_rsh/3.0_rsh)*( 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)) ) &
                                              /kZb**(5.0_rsh/3.0_rsh)
           ELSE
             kfond = SQRT( 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)) )/kZb**(2.0_rsh/3.0_rsh)
           END IF
         END IF

         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))
#else
         c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))
#endif
         a(k)=0.0_rsh

#ifdef key_obstructions
         c(k) = c(k) * obst_1o1mak
#endif

#if defined key_siggen || defined key_gencoord
         prod=c1turb*sigmapsi*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=c1turb*sigmapsi*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) THEN
           IF (pn2(k,i,j) < 0.0_rsh) dest=c3pturb*kz(k,i,j)*pn2(k,i,j)
           IF (pn2(k,i,j) > 0.0_rsh) dest=c3mturb*kz(k,i,j)*pn2(k,i,j)
         END IF

#ifdef key_wave_qeuler_diag
         diag_wf_prodpsi(k,i,j)=prod
#endif

         IF (l_bc_turb) THEN
#if defined key_siggen || defined key_gencoord
           enplus=dtm*nzmi*kfond/hz_w(k)
#else
           enplus=dtm*nzmi*kfond/hp
#endif
           bb(k)=dsigw(k)-c(k)+dtm*dsigw(k)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)
           IF (psi(k,i,j)<1.e-08) psi(k,i,j)=1e-08
           IF (ect(k,i,j)<1.e-08) ect(k,i,j)=1e-08
           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j) &
                -enplus+dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh) 
         ELSE
#if defined key_siggen || defined key_gencoord
           bb(k)=dsigw(k)-c(k)+dtm*nzmi/(REAL(hz_w(k)*hz(k),rsh)*dsigu(k)*0.5_rsh) &
                 +dtm*dsigw(k)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)
           enplus=dtm*nzmi*kfond/(hz_w(k)*hz(k)*dsigu(k)*0.5_rlg)
#else
           bb(k)=dsigw(k)-c(k)+dtm*nzmi/(hp*hp*dsigu(k)*0.5_rsh)        &
                 +dtm*dsigw(k)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)
           enplus=dtm*nzmi*kfond/(hp*hp*dsigu(k)*0.5_rsh)
#endif
           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)                &
                +enplus+dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh) 
         END IF

#ifdef key_obstructions
         taupsi_free = 1.0_rsh / ((1.0_rsh/c2turb)*(ect(k,i,j)/psi(k,i,j)))
         y(k) = y(k) + obst_t(k,i,j)*MAX(obst_tau(k,i,j),taupsi_free)*dtm*dsigw(k)
#endif

#ifdef key_wave_qeuler_dynamics
         y(k)=y(k)+peps(k,i,j)*dtm*dsigw(k)
#endif
#ifdef key_wave_qeuler_diag
         diag_wf_dissippsi(k,i,j)=-c2turb*eps(k,i,j)*fwall/ect(k,i,j)*psi(k,i,j)*(gamma-1.0_rsh)
#endif

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

           IF ( turb_2eq_option == 1 ) THEN
#if defined key_siggen || defined key_gencoord
             db=db+hz(k)*dsigu(k)
#else
             db=db+hp*dsigu(k)
#endif
             ds=hp-db
             fwall = 1.0_rsh+1.33_rsh*( psi(k,i,j)*(hp/(db*ds))  &
                                   /(ect(k,i,j)*0.4_rsh) )**2
           END IF

           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)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)
#if defined key_siggen || defined key_gencoord
           prod=c1turb*sigmapsi*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=c1turb*sigmapsi*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) THEN
             IF (pn2(k,i,j) < 0.0_rsh) dest=c3pturb*kz(k,i,j)*pn2(k,i,j)
             IF (pn2(k,i,j) > 0.0_rsh) dest=c3mturb*kz(k,i,j)*pn2(k,i,j)
           END IF
#ifdef key_wave_qeuler_diag
           diag_wf_prodpsi(k,i,j)=prod
#endif
           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)     &
                       +dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh)  

#ifdef key_obstructions
           taupsi_free = 1.0_rsh / ((1.0_rsh/c2turb)*(ect(k,i,j)/psi(k,i,j)))
           y(k) = y(k) + obst_t(k,i,j)*MAX(obst_tau(k,i,j),taupsi_free)*dtm*dsigw(k)
#endif

#ifdef key_wave_qeuler_dynamics
           y(k)=y(k)+peps(k,i,j)*dtm*dsigw(k)
#endif
#ifdef key_wave_qeuler_diag
           diag_wf_dissippsi(k,i,j)=-c2turb*eps(k,i,j)*fwall/ect(k,i,j)*psi(k,i,j)*(gamma-1.0_rsh)
#endif

         END DO

         !******************
         k=kkmax
         !******************

         IF ( l_bc_turb) THEN
#if defined key_siggen || defined key_gencoord
           prod=c1turb*sigmapsi*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=c1turb*sigmapsi*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) THEN
             IF (pn2(k,i,j) < 0.0_rsh) dest=c3pturb*kz(k,i,j)*pn2(k,i,j)
             IF (pn2(k,i,j) > 0.0_rsh) dest=c3mturb*kz(k,i,j)*pn2(k,i,j)
           END IF
#ifdef key_wave_qeuler_diag
           diag_wf_prodpsi(k,i,j)=prod
#endif

           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
           ELSE
             obst_1mak   = (1.0_rsh-obst_akkp1)
           ENDIF
           nzms = nzms * obst_1mak
#endif

           IF ( turb_2eq_option == 1 ) THEN
#if defined key_siggen || defined key_gencoord
             db=db+hz(k)*dsigu(k)
#else
             db=db+hp*dsigu(k)
#endif
             ds=hp-db
             fwall = 1.0_rsh+1.33_rsh*( psi(k,i,j)*(hp/(db*ds))  &
                                   /(ect(k,i,j)*0.4_rsh) )**2

             FpsiS = -ect(kmax,i,j)*0.4_rsh
           ELSE IF ( turb_2eq_option == 2 ) THEN
#if defined key_siggen || defined key_gencoord
             FpsiS = (cmu**3)*(ect(kmax,i,j)**1.5_rsh)/0.4_rsh/((REAL(hz(kmax),rsh)*dsigu(kmax)*0.5_rsh+z0s(i,j))**2)
#ifdef key_wave_qeuler_dynamics
             FpsiS = FpsiS &
                     +(cmu**3)*1.5_rsh*(1.0_rsh*(cmu**3)*100.0_rsh*(sqrt(tausurf(i,j))**3))  &
                     /(0.09_rsh*(0.4_rsh**2)*((REAL(hz(kmax),rsh)*dsigu(kmax)*0.5_rsh+z0s(i,j))**2))
#endif
#else
!             FpsiS = (cmu**3)*(ect(kmax,i,j)**1.5_rsh)/0.4_rsh/((hp*sigw(kmax-1)*0.5_rsh)**2)
             FpsiS = (cmu**3)*(ect(kmax,i,j)**1.5_rsh)/0.4_rsh/((-1.0_rsh*hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**2)
#ifdef key_wave_qeuler_dynamics
! ref: see eq. (3.129) in Burchard (2002)
             FpsiS = FpsiS &
                     +(cmu**3)*1.5_rsh*(1.0_rsh*(cmu**3)*100.0_rsh*(sqrt(tausurf(i,j))**3))  &
                     /(0.09_rsh*(0.4_rsh**2)*((-1.0_rsh*hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**2))
#endif
#endif
           ELSE IF ( turb_2eq_option == 3 ) THEN
!cval pourquoi n y a t il pas waveqeuler la ?
#if defined key_siggen || defined key_gencoord
             FpsiS = SQRT(ect(kmax,i,j))/0.4_rsh/((REAL(hz(kmax),rsh)*dsigu(kmax)*0.5_rsh+z0s(i,j))**2)/cmu
#else
             FpsiS = SQRT(ect(kmax,i,j))/0.4_rsh/((-hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**2)/cmu
#endif
           ELSE IF ( turb_2eq_option == 4 ) THEN
!cval pourquoi n y a t il pas waveqeuler la ?
#if defined key_siggen || defined key_gencoord
             FpsiS = (2.0_rsh/3.0_rsh)*(cmu**2)*ect(kmax,i,j)/(0.4_rsh**(2.0_rsh/3.0_rsh)) &
                    /((REAL(hz(kmax),rsh)*dsigu(kmax)*0.5_rsh+z0s(i,j))**(5.0_rsh/3.0_rsh))
#else
             FpsiS = (2.0_rsh/3.0_rsh)*(cmu**2)*ect(kmax,i,j)/(0.4_rsh**(2.0_rsh/3.0_rsh)) &
                    /((-hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**(5.0_rsh/3.0_rsh))
#endif
           END IF

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

#if defined key_siggen || defined key_gencoord
           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)         &
                +dtm*nzms*FpsiS/hz_w(k)                                                    &
                +dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh)
#else
           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)         &
                +dtm*nzms*FpsiS/hp                                                         &
                +dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh) 
#endif

#ifdef key_obstructions
           taupsi_free = 1.0_rsh / ((1.0_rsh/c2turb)*(ect(k,i,j)/eps(k,i,j)))
           y(k) = y(k) + obst_t(k,i,j)*MAX(obst_tau(k,i,j),taupsi_free)*dtm*dsigw(k)
#endif

#ifdef key_wave_qeuler_dynamics
           y(k)=y(k)+peps(k,i,j)*dtm*dsigw(k)
#endif
#ifdef key_wave_qeuler_diag
           diag_wf_dissippsi(k,i,j)=-c2turb*eps(k,i,j)*fwall/ect(k,i,j)*psi(k,i,j)*(gamma-1.0_rsh)
#endif

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

           a(k)=0.0_rsh
           bb(k)=1.0_rsh
           IF (turb_2eq_option==1) THEN
             y(k) = ect(k,i,j)*kZs
           ELSE IF (turb_2eq_option==2) THEN
#ifdef key_wave_qeuler_dynamics
              if (kZs.ne.0.0) then
!VALtoACB reference pour cela ?
              y(k) = (((sqrt(tausurf(i,j)))**3)/kZs)*(1.0_rsh+sqrt(1.5_rsh) &
                              *(0.09_rsh**0.25_rsh)*100.0_rsh)
              else
              y(k)=1.0e-10
              end if
#else
             y(k) = (cmu**3)*(ect(k,i,j)**1.5_rsh)/kZs
#endif

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

           ELSE IF (turb_2eq_option==3) THEN
             y(k) = SQRT(ect(k,i,j))/(kZs*cmu)
           ELSE IF (turb_2eq_option==4) THEN
             y(k) = (cmu**2)*ect(k,i,j)/(kZs**(2.0_rsh/3.0_rsh))
           END IF

         ELSE
           ! same code than above since kkmax=kmax if we use dirichlet boundary condition (l_bc_turb=.false.)
           a(k)=0.0_rsh
           bb(k)=1.0_rsh
           IF (turb_2eq_option==1) THEN
             y(k) = ect(k,i,j)*kZs
           ELSE IF (turb_2eq_option==2) THEN
             y(k) = (cmu**3)*(ect(k,i,j)**1.5_rsh)/kZs
           ELSE IF (turb_2eq_option==3) THEN
             y(k) = SQRT(ect(k,i,j))/(kZs*cmu)
           ELSE IF (turb_2eq_option==4) THEN
             y(k) = (cmu**2)*ect(k,i,j)/(kZs**(2.0_rsh/3.0_rsh))
           END IF
         END IF

         !********************************
         !  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 Upsi=u0
         !******************************
         psi(kmax,i,j)=u0(kmax)
         DO k=kmax-1,1,-1
           psi(k,i,j)=q0(k)*psi(k+1,i,j)+u0(k)
         END DO

#ifdef key_wave_qeuler_diag
           diag_wf_psi(1:kmax,i,j)=psi(1:kmax,i,j)
           k=1
           diag_wf_nudpsidz(k,i,j)=(                  c(k)*psi(k+1,i,j) &
                                  -(a(k)+c(k))*psi(k,i,j))/(dtm*dsigw(k))
           diag_wf_dissippsi(k,i,j)=diag_wf_dissippsi(k,i,j)-c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)*psi(k,i,j)
           do k=2,kmax-1
             diag_wf_nudpsidz(k,i,j)=(a(k)*psi(k-1,i,j)+c(k)*psi(k+1,i,j) &
                                    -(a(k)+c(k))*psi(k,i,j))/(dtm*dsigw(k))
             diag_wf_dissippsi(k,i,j)=diag_wf_dissippsi(k,i,j)-c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)*psi(k,i,j)
           END DO
#endif

         psi(:,i,j)=MAX(psi(:,i,j),psimin)
         nz(:,i,j)=nz(:,i,j)*sigmapsi       ! since we have divided by sigmapsi at the entrance of the routine
       END IF
     END DO
   END DO
!$OMP END DO

!$OMP END PARALLEL

#ifdef key_turbadv
   !-----------------------------------------------------------------------------
   ! TURBULENCE ADVECTION
   CALL turbadv_uzvzwz(xeadv,xeuv,uz,vz)
#ifdef key_turbadv_wcompact
     !schema compact pour advection de ect & psi sur la verticale
     CALL turbadv_wcompact(xeadv,xeuv,ect,psi)
#else
     !schema centre par defaut pour advection de ect & psi sur la verticale
     CALL turbadv_wcen(xeadv,xeuv,ect,psi)
#endif
   CALL turbadv_obc3d_apply(uzsw,vzsw,ect,psi)
#endif
   !-----------------------------------------------------------------------------

  END SUBROUTINE turbpsi
