  SUBROUTINE turbkl(xeuv,uz,vz,boz,pn2,ect)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE turbkl_gasp  ***
!&E
!&E ** Purpose : compute the vertical viscosity and diffusion coefficients,
!&E              using the k-l parameterisation by Gaspard
!&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,boz,ect,lk,z0b,le,uz,vz,nz,kz,tosx,tosy
!&E
!&E ** Modified variables : le, lk, kz, nz, ect
!&E
!&E ** Reference : Gaspar et al, 1990
!&E
!&E ** History :
!&E       !  2004-09-15  Original code
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-01  (P. Marsaleix, F. Dumas) square estimate of shear before spatial average
!&E       !  2009-03  (T. Duhaut) corrections
!&E       !  2009-04  (A. Mano, V. Garnier) MPI2D + removing of twisted open boundary option
!&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       !  2012-02  (F. Dumas) Initialization of nz, kz at kmax
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters,ONLY : rsh,rlg,imin,imax,ljmin,ljmax,ljminm1,kmax,ljminp2,ljmaxm1,    &
                         limin,limax,liminm1,limaxp1,ljminm1,jmin,jmax,limaxp2,ljmaxp2, &
                         liminm2,ljminm2,ljmaxp1,liminp2,limaxm1
   USE parameters,ONLY : l_testcase,num_testcase,epsilon
   USE parameters,ONLY : ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north
   USE comvars2d, ONLY : ig,id,jh,jb,nbouc,h0,hm,z0b,z0s,l_obc_cycl_x,l_obc_cycl_y
   USE comvars2d, ONLY : halfdt
   USE comvars3d, ONLY : dsigu,sig,sigw,nz,kz,dsigw,rhoref
   USE comturb,   ONLY : cmu,emin
   USE cometdyn,  ONLY : tosx,tosy
   USE toolgeom,  ONLY : f_zw
   USE obccombine, ONLY : l_obc_south, l_obc_north, l_obc_east,l_obc_west


   USE toolmpi,   ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_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,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in)    :: boz
   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,n,kk
   REAL(kind=rsh) :: nzms,kfond,immer
   REAL(kind=rsh) :: sqrtcmu,eps0,dtm,hp,ckl,zce,epote,xrod,aux,ri1,ri2
   REAL(kind=rsh) :: cd,epsplus,prod,dest,enplus,tos,xlme,xlpe,Prt,lk

   REAL(kind=rlg),DIMENSION(kmax)   :: a,bb,c,y,u0,p0

   REAL(kind=rlg),DIMENSION(0:kmax) :: q0
   REAL(kind=rlg),DIMENSION(kmax) :: le
   REAL(kind=rsh) :: zlim

!!----------------------------------------------------------------------
!! * Executable part
 
!------------------------------------------------------------------------------
!  coefficients de dispersion horizontale (pour calcul d advection dispersion)
!------------------------------------------------------------------------------
 
!$OMP PARALLEL default(none) shared(ljmin,ljmax,ljminm1,     &
!$OMP& ljminp2,ljmaxm1,limin,limax,liminm1,limaxm1,          &
!$OMP& limaxp2,ljmaxp2,liminm2,ljminm2,liminp2,              &
!$OMP& ig,id,nbouc,h0,hm,z0b,z0s,halfdt,emin,                &
!$OMP& dsigu,nz,kz,dsigw,rhoref,tosx,tosy,                   &
!$OMP& l_obc_south, l_obc_north, l_obc_east,l_obc_west,      &

!$OMP& limaxp1,ljmaxp1,                                      &


!$OMP& sig,sigw,                                             &

!$OMP& xeuv,uz,vz,boz,pn2,ect,l_obc_cycl_x,l_obc_cycl_y)     &

!$OMP& PRIVATE(k,i,j,n,kk,nzms,kfond,cmu,immer,              &
!$OMP& sqrtcmu,eps0,dtm,hp,ckl,zce,epote,xrod,aux,ri1,ri2,   &
!$OMP& cd,epsplus,prod,dest,enplus,tos,xlme,xlpe,Prt,lk,     &


!$OMP& a,bb,c,y,u0,p0,q0,le,zlim)

   cmu=0.09_rsh
   sqrtcmu=SQRT(cmu)
   eps0=0.166_rsh
   dtm=halfdt
   ckl =0.1_rsh
   zce=1.0_rsh/0.7_rsh
   ri1=0.0_rsh
   ri2=0.5_rsh

!$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


!calcul des longueurs de melange
!*******************************

         DO k=1,kmax-1
           xrod=(boz(k,i,j)+boz(k+1,i,j))*0.5_rsh
           xlpe=0.0_rsh
           epote=0.0_rsh
           n=k

           DO kk=k+1,kmax

             aux=epote+dsigu(kk)*hp*(boz(kk,i,j)-xrod)
             IF (aux >= ect(k,i,j)) GOTO 40
             n=kk
             epote=epote+dsigu(kk)*hp*(boz(kk,i,j)-xrod)
             xlpe=xlpe+dsigu(kk)*hp

           END DO

40         IF (n==k) THEN

             xlpe=SQRT(2.0_rsh*(sig(k+1)-sigw(k))*hp/(boz(k+1,i,j)-xrod)*ect(k,i,j))

           ELSE
             IF (n.NE.kmax) THEN
               xlpe=xlpe+(ect(k,i,j)-epote)/(boz(n+1,i,j)-xrod)
             END IF
           END IF

           xlme=0.0_rsh
           epote=0.0_rsh
           n=k

           DO kk=k,1,-1

             aux=epote+dsigu(kk)*hp*(xrod-boz(kk,i,j))
             IF (aux >= ect(k,i,j)) GOTO 50
             n=kk
             epote=epote+dsigu(kk)*hp*(xrod-boz(kk,i,j))
             xlme=xlme+dsigu(kk)*hp

           END DO

50         IF (kk==k) THEN

             xlme=SQRT(2.0_rsh*(sigw(k)-sig(k))*hp/(xrod-boz(k,i,j))*ect(k,i,j))

           ELSE
             IF (n.NE.1) THEN
               xlme=xlme+(ect(k,i,j)-epote)/(xrod-boz(n-1,i,j))
             END IF
           END IF

           immer = - ( f_zw(h0(i,j),xeuv(i,j),k,i,j) - xeuv(i,j) )
           IF ( immer <= z0s(i,j) )  xlpe = MAX(xlpe,z0s(i,j))

           immer = h0(i,j) + f_zw(h0(i,j),xeuv(i,j),k,i,j)
           IF ( immer <= z0b(i,j) ) xlme = MAX(xlme,z0b(i,j))

           le(k)=MAX(REAL(SQRT(xlme*xlpe),rlg),epsilon)
           lk=MIN(xlme,xlpe)
           nz(k,i,j)=ckl*lk*SQRT(ect(k,i,j))
           kz(k,i,j)=ckl*lk*SQRT(ect(k,i,j))

         END DO

! Estimate of the viscosity and diffusivity at the surface  ! What can we do otherwise ????
! We can also use the intermediate variable with no vertical gradient of uz and vz
! or use a logarithmic gradient of velocity ?
         nz(kmax,i,j)=nz(kmax-1,i,j)
         kz(kmax,i,j)=kz(kmax-1,i,j)

         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/MAX(LOG((sig(1)+1.0_rsh)*hp/z0b(i,j)),epsilon))**2,0.1_rsh)
         cd=(0.4_rsh/LOG((sig(k)+1.0_rsh)*hp/z0b(i,j)+1.0_rsh))**2


         epsplus=dtm*dsigw(k)*SQRT(ect(k,i,j))/le(k)/zce
! kfond=(g1*sqrt(uinf**2+vinf**2)+g2*cd*(uinf**2+vinf**2))/sqrtcmu
         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))/sqrtcmu

! c est une interpolation lineaire des nz aux point u sur la verticale.
! correct si on considere que le point u est au milieu des deux points w.
         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))
         bb(k)=REAL(dsigw(k),rlg)-c(k)+dtm*nz(k,i,j)/(hp*hp*dsigu(k))+1.3_rlg*epsplus
         a(k)=0.0_rlg

!    ri=(boz(k+1,i,j)-boz(k,i,j))/(hp*dsigw(k))/ &
!       max(((vsup-vinf)**2+(usup-uinf)**2)/(hp*dsigw(k))**2,10.0e-8)
         Prt=1.0_rsh
!    if(ri >= ri1 .and. ri < ri2) then
!     Prt= (1.0_rsh*(ri2-ri)+5.0_rsh*(ri-ri1))/(ri2-ri1)
!    end if
!    if(ri >= ri2 ) then
!    Prt=5.0_rsh
!   end if
         kz(k,i,j)=kz(k,i,j)/Prt

         prod=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=kz(k,i,j)*pn2(k,i,j)
         y(k)=ect(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)+enplus+epsplus*ect(k,i,j)*0.3_rlg

         DO k=2,kmax-1

!*****************************

! interpolation lineaire a voir dans le detail.
           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)
           epsplus=dtm*dsigw(k)*SQRT(ect(k,i,j))/le(k)/zce
           bb(k)=REAL(dsigw(k),rlg)-a(k)-c(k)+1.3_rlg*epsplus
!     usup=0.5_rsh*(uz(k+1,i,j)+uz(k+1,i-1,j))  ! cval attention faire carre avant moyenne
!     vsup=0.5_rsh*(vz(k+1,i,j)+vz(k+1,i,j-1))  ! cval attention faire carre avant moyenne
!     ri=(boz(k+1,i,j)-boz(k,i,j))/(hp*dsigw(k))/ &
!       max(((vsup-vinf)**2+(usup-uinf)**2)/(hp*dsigw(k))**2,10.0e-8)
           Prt=1.0_rsh
!     if(ri >=ri1 .and. ri < ri2) then
!      Prt= (1.0_rsh*(ri2-ri)+5.0_rsh*(ri-ri1))/(ri2-ri1)
!     end if
!     if(ri >=ri2 ) then
!      Prt=5.0_rsh
!     end if

           kz(k,i,j)=kz(k,i,j)/Prt

           prod=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=kz(k,i,j)*pn2(k,i,j)
           y(k)=ect(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)+epsplus*ect(k,i,j)*0.3_rlg

         END DO

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

         tos=SQRT(tosx(i,j)**2+tosy(i,j)**2)
         a(k)=0.0_rlg
         bb(k)=1.0_rlg
         y(k)=(tos/rhoref)/sqrtcmu

!********************************
!  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

         DO k=1,kmax
           ect(k,i,j)=MAX(ect(k,i,j),emin)
         END DO

       END IF

     END DO
   END DO
!$OMP END DO

!ect3(:,:,:)=ect(:,:,:)

!$OMP SINGLE
   IF ( l_testcase .AND. num_testcase==22 ) THEN
     zlim=50.0_rsh

     DO k=kmax-1,1,-1
       IF ( ect(k,imax/2,2) > 0.00001_rsh ) THEN
         zlim=50.0_rsh*(dsigu(1)*0.5_rsh +(k-1)*dsigu(1))
       ELSE
         IF ( ect(k,imax/2,2)-ect(k+1,imax/2,2) < 0.0_rsh ) THEN
           WRITE(150,*) nbouc*15.0_rsh/3600.0_rsh,(50.0_rsh-zlim)
           GOTO 45
         ELSE
           zlim=dsigu(1)*0.5_rsh +(k-1)*dsigu(1)
         END IF
       END IF
     END DO

45   IF (MOD(nbouc,54) ==0) THEN
       DO k=kmax,1,-1
         WRITE(10,*) dsigu(1)*0.5_rsh +(k-1)*dsigu(1),ect(k,imax/2,2)
       END DO
       WRITE(10,'(/)') ! ecriture d une ligne blanche
       WRITE(10,'(/)') ! ecriture d une ligne blanche
     END IF

   END IF

! ----------------------------------------------------------
! duplicate nz along the boundaries or use cyclic boundaries
! ----------------------------------------------------------

   IF (l_obc_cycl_y) THEN


     CALL ex_peri_y_rsh(2,kmax+1,limin  ,limaxp1,ljmin  ,ljmaxp1,nz(:,limin :limaxp1,ljmin  :ljmaxp1))
!     CALL ex_peri_y_rsh(2,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,kz(:,liminm1:limaxp1,ljminm1:ljmaxp1))


   ELSE IF (.NOT. l_obc_cycl_y .AND. l_obc_south) THEN

     IF (ljmin.EQ.jmin) THEN
       DO i=liminp2,limaxm1
           nz(:,i,jmin+1+jshift_obc_south)=nz(:,i,jmin+2+jshift_obc_south)
!           kz(:,i,jmin+1)=kz(:,i,jmin+2)
       END DO
     END IF

   ELSE IF (.NOT. l_obc_cycl_y .AND. l_obc_north) THEN

     IF (ljmax.EQ.jmax) THEN
       DO i=liminp2,limaxm1
         nz(:,i,jmax-jshift_obc_north)=nz(:,i,jmax-1-jshift_obc_north)
!           kz(:,i,jmax)=kz(:,i,jmax-1)
       END DO
     END IF

   END IF

   IF (l_obc_cycl_x) THEN


     CALL ex_peri_x_rsh(2,kmax+1,limin  ,limaxp1,ljmin  ,ljmaxp1,nz(:,limin :limaxp1,ljmin  :ljmaxp1))
!     CALL ex_peri_x_rsh(2,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,kz(:,liminm1:limaxp1,ljminm1:ljmaxp1))


   ELSE IF (.NOT. l_obc_cycl_x .AND. l_obc_west) THEN

     IF (limin.EQ.imin) THEN
       DO j=ljminp2,ljmax
         nz(:,imin+1+ishift_obc_west,j)=nz(:,imin+2+ishift_obc_west,j)
       ENDDO
     END IF

   ELSE IF (.NOT. l_obc_cycl_x .AND. l_obc_east) THEN

     IF (limax.EQ.imax) THEN
       DO j=ljminp2,ljmax
         nz(:,imax-ishift_obc_east,j)=nz(:,imax-1-ishift_obc_east,j)
       ENDDO
     END IF

   END IF
!$OMP END SINGLE

!$OMP END PARALLEL

   CALL ex_i_rsh(0,1,kmax+1,limin,limaxp1,ljmin,ljmaxp1,nz(:,limin:limaxp1,ljmin:ljmaxp1))
   CALL ex_j_rsh(0,1,kmax+1,limin,limaxp1,ljmin,ljmaxp1,nz(:,limin:limaxp1,ljmin:ljmaxp1))


  END SUBROUTINE turbkl
