  SUBROUTINE turb(xeuv,uz,vz,pn2)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE turb  ***
!&E
!&E ** Purpose : compute the vertical viscosity and diffusion coefficients
!&E
!&E ** Description :
!&E
!&E ** Called by : stepu, stepv
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : z0b, uz, vz, xeuv, h0, hm, kzinit, nzinit, turb_0eq_option,
!&E                     bgdiff, sig, pn2, ig, id, jb, jh
!&E
!&E ** Modified variables : kz, nz
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2004-09-15  Original code
!&E       !  2008-12    (V. Garnier) Utilisation de sig (au lieu de sigw) pour calcul de u*
!&E       !  2008-12    (V. Garnier) shift of sigma values : -1 < sig < 0
!&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
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters,ONLY : rsh,rlg,imin,imax,jmin,jmax,kmax,epsilon,  &
                         ljmin,ljmax,ljminm1,ljminp2,ljmaxp1,       &
                         limin,limax,liminm1,limaxp1,liminp2,       &
                         limaxm1,limaxp2,ljmaxp2,liminm2,ljminm2
   USE parameters,ONLY : ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north
   USE comvars2d, ONLY : ig,id,h0,z0b,icon,hm,jb,jh,l_testcase
   USE comvars2d, ONLY : l_obc_cycl_x,l_obc_cycl_y,frot_lin
   USE comvars3d, ONLY : turb_0eq_option,nz,nzinit,kz,kzinit,sig,sigw,bgdiff
   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,limin:limax,ljmin:ljmax),INTENT(in) :: pn2

!! * Local declarations
   INTEGER        :: k,i,j
   REAL(kind=rsh) :: l,h,ufond,ufond2,nzprov,kzprov,filt
   REAL(kind=rsh) :: ri,dvdz,usup,vsup,z,uinf,vinf,c1,c2,n0



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

!-----------------------------------------
! viscosite constante
!-----------------------------------------

   IF (turb_0eq_option == 1 .AND. .NOT.l_testcase) THEN
     nz(:,limin  :limaxp1,ljmin  :ljmaxp1)=nzinit
     kz(:,liminm1:limaxp1,ljminm1:ljmaxp1)=kzinit
     RETURN
   END IF
   
!-----------------------------------------
! modele de Prandt  (not really) :
!                  Nz = n0 * sigma*sqrt(1-sigma)
!                  with n0 = u* H/15
!-----------------------------------------
! viscosite=f(u*) longueur de melange constante sur la hauteur d eau
! n0 coefficient de viscosite
!            - dependant de u*
!            - fonction de la profondeur
!            - * profondeur selon  Csanady l=H/20. ou  Elder=> l=H/15
! u*=sqrt(tau)=sqrt(Cd*u2)=0.4/log((sig(1)+1.0_rsh)*h/z0b(i,j))*ufond
! u*=0.4/log(dz_fond/z0b)*h*|ufond|
! nz=u* z_fond/h sqrt(z_surf/h)
   
   IF (turb_0eq_option == 2) THEN
     k=1
     DO j=ljminp2,ljmax
       DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

         h=h0(i,j)+xeuv(i,j)
         IF (h >= hm) THEN
           k=1
           ufond2=(uz(k,i,j)+uz(k,i-1,j))**2+(vz(k,i,j)+vz(k,i,j-1))**2
           ufond=0.5_rsh*SQRT(ufond2)

! n0 (ustar) is estimated in the middle of the cell

           n0=0.4_rsh/LOG((sig(1)+1.0_rsh)*h/z0b(i,j))*ufond*h/15.0_rsh

! nz is estimated at each interface
           nz(1:kmax,i,j)=n0*(sigw(1:kmax)+1.0_rsh)*SQRT(ABS(sigw(1:kmax)))+bgdiff


         END IF
       ENDDO
     ENDDO

   ELSE IF (turb_0eq_option == 3) THEN

!-----------------------------------------
! modele de Quetin :
!                    Nz =       l*l dU/dz
!                    Kz = f(Ri) l*l dU/dz  ; f(Ri)=exp(-0.6*Ri)
!                    with mixing length l = Kz(1-0.7z/D)
!-----------------------------------------

     filt=0.9_rsh
     c1=0.4_rsh
     c2=0.7_rsh

     DO j=ljminp2,ljmax
       DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

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

           k=0
           uinf=0.5_rsh*(uz(k+1,i,j)+uz(k+1,i-1,j))
           vinf=0.5_rsh*(vz(k+1,i,j)+vz(k+1,i,j-1))
           DO k=1,kmax-1
             usup=0.5_rsh*(uz(k+1,i,j)+uz(k+1,i-1,j))
             vsup=0.5_rsh*(vz(k+1,i,j)+vz(k+1,i,j-1))

! water depth from the bottom
             z=(1.0_rsh+sigw(k))*h
!       unsl=1.0_rsh/(c1*(z))+1.0_rsh/(c2*(h-z))
!       l=1.0_rsh/unsl
! l = 0.4 z_bot (1-0.7z_bot/D)
             l=0.4_rsh*z*(1.0_rsh-0.7_rsh*(1.0_rsh+sigw(k))) ! l = Kz(1-0.7z/D)
             dvdz=SQRT((vsup-vinf)**2+(usup-uinf)**2)/(h*(sig(k+1)-sig(k)))

             dvdz=MAX(REAL(dvdz,rlg),epsilon)
             IF (icon /=0 ) THEN
               ri=MAX(pn2(k,i,j)/(dvdz*dvdz),-0.09_rsh)
             ELSE
               ri=0.0_rsh
             END IF
             nzprov=l*l*dvdz/SQRT(1.0_rsh+10.0_rsh*ri)+bgdiff
             nz(k,i,j)=filt*nzprov+(1.0_rsh-filt)*nz(k,i,j)
             kzprov=l*l*dvdz*EXP(-0.6_rsh*ri)+bgdiff
             kz(k,i,j)=filt*kzprov+(1.0_rsh-filt)*kz(k,i,j)
             vinf=vsup
             uinf=usup
           ENDDO
         END IF
       ENDDO
     ENDDO

   ELSE IF (turb_0eq_option == 4) THEN

!-----------------------------------------
! Pacanovski et Philander (1981)
!-----------------------------------------
   
!$OMP parallel default(none)                &
!$OMP& private(i,j,filt,c1,c2,h,k,uinf,vinf &
!$OMP& ,z,usup,vsup,l,dvdz,nzprov,kzprov,   &

!$OMP& ri) shared(xeuv,h0,uz,vz,nz,kz,      &

!$OMP& sig,sigw,                            &

!$OMP& limin,limax,ljmin,ljmax,limaxm1,     &
!$OMP& liminp2,ljminp2,                     &
!$OMP& hm,bgdiff,icon,ig,id,pn2)
     filt=0.9_rsh
     c1=0.4_rsh
     c2=0.7_rsh


!$OMP do schedule(runtime)

     DO j=ljminp2,ljmax
       DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

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

           k=0
           uinf=0.5_rsh*(uz(k+1,i,j)+uz(k+1,i-1,j))
           vinf=0.5_rsh*(vz(k+1,i,j)+vz(k+1,i,j-1))
           IF (icon == 0) THEN
             DO k=1,kmax-1


! water depth from the bottom
               z=(1.0_rsh+sigw(k))*h
!       unsl=1.0_rsh/(c1*z)+1.0_rsh/(c2*(h-z))
!       l=1.0_rsh/unsl
               l=0.4_rsh*z*(1.0_rsh-0.7_rsh*(1.0_rsh+sigw(k)))

               usup=0.5_rsh*(uz(k+1,i,j)+uz(k+1,i-1,j))
               vsup=0.5_rsh*(vz(k+1,i,j)+vz(k+1,i,j-1))

               dvdz=SQRT((vsup-vinf)**2+(usup-uinf)**2)/(h*(sig(k+1)-sig(k)))

               nzprov=l*l*dvdz+bgdiff
               nz(k,i,j)=filt*nzprov+(1.0_rsh-filt)*nz(k,i,j)
               kzprov=l*l*dvdz+bgdiff
               kz(k,i,j)=filt*kzprov+(1.0_rsh-filt)*kz(k,i,j)
               vinf=vsup
               uinf=usup
             ENDDO
           ELSE
             DO k=1,kmax-1
               usup=0.5_rsh*(uz(k+1,i,j)+uz(k+1,i-1,j))
               vsup=0.5_rsh*(vz(k+1,i,j)+vz(k+1,i,j-1))

               dvdz=SQRT((vsup-vinf)**2+(usup-uinf)**2)/(h*(sig(k+1)-sig(k)))

               dvdz=MAX(dvdz,0.0001_rsh)
               ri=MAX(pn2(k,i,j)/(dvdz*dvdz),-0.19_rsh)
               nzprov=0.01_rsh/(1.0_rsh+5.0_rsh*ri)**2+bgdiff
               nz(k,i,j)=filt*nzprov+(1.0_rsh-filt)*nz(k,i,j)
               kzprov= nz(k,i,j)/(1.0_rsh+5.0_rsh*ri)+bgdiff
               kz(k,i,j)=filt*kzprov+(1.0_rsh-filt)*kz(k,i,j)
               vinf=vsup
               uinf=usup
             ENDDO
           END IF
         END IF
       ENDDO

     ENDDO
!$OMP end do

!$OMP end parallel

   END IF

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

   IF (turb_0eq_option >= 2) THEN

     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


     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 IF

  END SUBROUTINE turb
