  SUBROUTINE turblepskznz(xeuv,pn2,psi,ect)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE turblepskznz  ***
!&E
!&E ** Purpose : compute the vertical viscosity and diffusion coefficients,
!&E              using the k GLS 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,ect,psi,eps,pn2,kz,nz
!&E
!&E ** Modified variables : eps,nz,kz
!&E
!&E ** History :
!&E       !  2004-09
!&E       !  2009-03  (T. Duhaut) corrections
!&E       !  2009-04  (A. Mano) MPI2D
!&E       !  2009-12  (T. Duhaut) stability function
!&E       !  2011-10  (V. Garnier) OMP
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2012-02  (F. Dumas) Initialization of nz, kz at kmax
!&E       !  2013-05  (S. Petton) Reformulation of stability functions using a unifying notation
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d
   USE cometdyn
   USE comturb
   USE comvarp3d,   ONLY : uz,vz
   USE toolgeom,    ONLY : f_zw,f_hzw
   USE obccombine,  ONLY : l_obc_south, l_obc_north,l_obc_east,l_obc_west
   USE toolmpi, ONLY : ex_j_rsh,ex_i_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,limin  :limax  ,ljmin  :ljmax  ),INTENT(inout):: psi,ect
   REAL(kind=rsh),DIMENSION(kmax,limin  :limax  ,ljmin  :ljmax  ),INTENT(in)   :: pn2

!! * Local declarations
   INTEGER        :: i,j,k
   REAL(kind=rsh) :: hp,hz_w
   REAL(kind=rsh) :: M2, N2
   REAL(kind=rsh) :: alpha_m,alpha_n,c_mu,c_mu_prim,Denom,alpha_m_max
! testcase
   INTEGER        :: kk
   REAL(kind=rsh) :: zlim

!!----------------------------------------------------------------------
!! * Executable part
   
!------------------------------------------------------------------------------
! calcul du taux de dissipation (eps) de l ect en utilisant la formulation de la GLS
! des fonctions de stabilite (c_mu, c_mu_prim) et des kz et nz
!------------------------------------------------------------------------------

!$OMP PARALLEL DEFAULT(none) SHARED(h0,hm,xeuv,ect,eps,psi,pn2,uz,vz,kz,nz,   &
!$OMP& ljminp2,ljmaxm1,limin,limax,ljmin,ljmax,id,ig,dsigu,dsigw,nbouc, &
!$OMP& cmu,pturb,mturb,nturb,bgdiff,epsmin,l_obc_cycl_x,l_obc_cycl_y,         &
!$OMP& l_obc_west,l_obc_east,l_obc_north, l_obc_south,                        &
!$OMP& sf_d0,sf_d1,sf_d2,sf_d3,sf_d4,sf_d5,                                   &
!$OMP& sf_n0,sf_n1,sf_n2,sf_nb0,sf_nb1,sf_nb2,liminp2,limaxm1,                &

!$OMP& limaxp1,ljmaxp1,                                                       &


!$OMP& lim_am0,lim_am1,lim_am2,lim_am3,lim_am4,lim_am5,lim_am6,alpha_n_min)   &

!$OMP& PRIVATE(i,j,k,hp,hz_w,M2,N2,alpha_m,alpha_n,c_mu,c_mu_prim,Denom,      &
!$OMP& alpha_m_max,zlim,kk)

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


! Dissipation rate
        eps(:,i,j)=cmu**(3+pturb/nturb)*ect(:,i,j)**(1.5_rsh+mturb/nturb)*psi(:,i,j)**(-1./nturb)
        eps(:,i,j) = MAX(eps(:,i,j),epsmin)

        DO k=1,kmax-1

! Shear and Buoyancy

          M2= 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

          N2=pn2(k,i,j)
! Dimensionless coefficients
          alpha_m = M2*ect(k,i,j)**2/(eps(k,i,j)**2)
          alpha_n = N2*ect(k,i,j)**2/(eps(k,i,j)**2)

! Stability and Realisability
          alpha_n = MIN(MAX(alpha_n_min , alpha_n),1.0e10)
          alpha_m_max = (lim_am0 + lim_am1 * alpha_n + lim_am2 * alpha_n**2 + lim_am3 * alpha_n**3) / &
                        (lim_am4 + lim_am5 * alpha_n + lim_am6 * alpha_n**2)
          alpha_m = MIN(alpha_m , alpha_m_max) 

! Functions
          Denom = sf_d0  +  sf_d1*alpha_n +  sf_d2*alpha_m + sf_d3*alpha_n*alpha_m + sf_d4*alpha_n**2.0_rsh + sf_d5*alpha_m**2.0_rsh
          c_mu      = (sf_n0  +  sf_n1*alpha_n +  sf_n2*alpha_m)/Denom
          c_mu_prim = (sf_nb0 + sf_nb1*alpha_n + sf_nb2*alpha_m)/Denom

! Vertical viscosity Nz and diffusivity Kz
          nz(k,i,j)=ect(k,i,j)**2.0_rsh/eps(k,i,j)*c_mu
          kz(k,i,j)=ect(k,i,j)**2.0_rsh/eps(k,i,j)*c_mu_prim

          nz(k,i,j)=MIN(MAX(nz(k,i,j),bgdiff),0.9_rsh)
          kz(k,i,j)=MIN(MAX(kz(k,i,j),bgdiff),0.9_rsh)

        END DO ! end loop on k

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

       ELSE
! RAZ de eps pour banc decouvrant
         eps(:,i,j)=epsmin
       END IF ! end test on h0+xe



     END DO 
   END DO
!$OMP END DO

!$OMP SINGLE
   IF ( l_testcase .AND. num_testcase==22 ) THEN
     zlim=50.0_rsh
     DO k=kmax-1,1,-1
       IF ( ect(k,2,2) > 0.00001_rsh ) THEN
         zlim=50*(dsigu(1)*0.5_rsh +(k-1)*dsigu(1))
       ELSE
         IF ( ect(k,2,2)-ect(k+1,2,2) < 0.0_rsh ) THEN
           WRITE(153,*) nbouc*15_rsh/3600.0_rsh,(50.0_rsh-zlim)
           IF (MOD(nbouc,54) ==0) THEN
             DO kk=kmax,1,-1
               WRITE(13,*) dsigu(1)*0.5_rsh +(kk-1)*dsigu(1),ect(kk,2,2)
             END DO
             WRITE(13,'(/)') ! ecriture d une ligne blanche
           END IF
         ELSE
           zlim=dsigu(1)*0.5_rsh +(k-1)*dsigu(1)
         END IF
       END IF
     END DO

   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 turblepskznz
