  SUBROUTINE turb_smago(uz,vz,igd)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE turb_smago  ***
!&E
!&E ** Purpose : compute the horizontal viscosity and diffusion coefficients,
!&E              depending on the local shear conditions
!&E
!&E ** Description : use the method by Smagorinsky
!&E
!&E ** Called by : dyn3dx, dyn3dy
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : igd,vz,uz,dy,dx,tfiusr,vish_xe,dy_f,dx_f,tfivsr,vish_phi,h0_g
!&E
!&E ** Modified variables : vish_xe,vish_phi
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2004-08-20  Original code
!&E       !  2007-09    (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-08    (T. Odaka) MPI
!&E       !  2009-04    (A. Mano) MPI2D
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,ONLY : ig,id,jb,jh,h0_g
   USE comvars2d,ONLY : vismul,vismin,vismax,fvisc,cosmag, &
                       l_sponge,l_spongn,l_spongs,l_spongw,sponge_nbcell
   USE comvars2d,ONLY : dx,dy,dx_f,dy_f,tfivsr,tfiusr
   USE comvars2d,ONLY : vish_phi,vish_xe


   USE toolmpi,  ONLY : ex_i_rsh, ex_j_rsh

   IMPLICIT NONE

!! * Arguments
   INTEGER       ,DIMENSION(     liminm1:limaxp1,ljminm1:ljmaxp1),INTENT(in) :: igd
   REAL(kind=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: uz,vz

!! * Local declarations
   INTEGER        :: i,j,k
   REAL(kind=rsh) :: nu0,nub
   REAL(kind=rsh) :: vb,vc,vd,ve,vf,vg,vh,vi,vzuxe,uzvxe,vzuphi,uzvphi
   REAL(kind=rsh) :: ua,ub,uc,ud,ue,uf,ug,uh

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

   !PRINT*, 'enter TURB_SMAGO routine'


!$OMP do SCHEDULE(runtime ) &

!$OMP& private(i,j,k,vb,vc,vd,ve,vf,vg,vh,vi, &
!$OMP& vzuxe,uzvxe,vzuphi,uzvphi, &
!$OMP& ua,ub,uc,ud,ue,uf,ug,uh)

   DO j=ljminp1,ljmaxm1
     DO i=liminp1,limaxm1
       IF (igd(i,j)/=4 .OR. igd(i  ,j+1)/=4 .OR. igd(i  ,j-1)/=4 &
                       .OR. igd(i+1,j  )/=4 .OR. igd(i-1,j  )/=4) THEN
         DO k=1,kmax
           vg=vz(k,i-1,j-1)
           vh=vz(k,i,j-1)
           vi=vz(k,i+1,j-1)

           vd=vz(k,i-1,j)
           ve=vz(k,i,j)
           vf=vz(k,i+1,j)

           vb=vz(k,i,j+1)
           vc=vz(k,i+1,j+1)

           ug=uz(k,i-1,j-1)
           uh=uz(k,i,j-1)

           ud=uz(k,i-1,j)
           ue=uz(k,i,j)
           uf=uz(k,i+1,j)

           ua=uz(k,i-1,j+1)
           ub=uz(k,i,j+1)
           uc=uz(k,i+1,j+1)

           vzuphi=0.25_rsh*(vb+vc-vh-vi)
           uzvxe =0.25_rsh*(ua+ub-ug-uh)
           vzuxe =0.25_rsh*(vf+vi-vd-vg)
           uzvphi=0.25_rsh*(uc+uf-ua-ud)

           vish_xe(k,i,j)=cosmag*dy(i,j)*dx(i,j)*SQRT(((ue-ud)/dx(i,j)          &
                          -tfiusr(i,j)*0.5_rsh*(ve+vh))**2+((ve-vh)/dy(i,j))**2 &
                          +0.5_rsh*(uzvxe/dy(i,j)+tfiusr(i,j)*0.5_rsh*(ue+ud)   &
                          +vzuxe/dx(i,j))**2)
           vish_xe(k,i,j)=MIN(MAX(vish_xe(k,i,j),vismin),vismax)

           vish_phi(k,i,j)=cosmag*dy_f(i,j)*dx_f(i,j)*SQRT((uzvphi/dx_f(i,j)        &
                           -tfivsr(i,j)*0.5_rsh*(vf+ve))**2+(vzuphi/dy_f(i,j))**2  &
                           +0.5_rsh*((ub-ue)/dy_f(i,j)+tfivsr(i,j)*0.5_rsh*(ub+ue) &
                           +(vf-ve)/dx_f(i,j))**2)
           vish_phi(k,i,j)=MIN(MAX(vish_phi(k,i,j),vismin),vismax)
         END DO

       END IF
     END DO
   END DO

!$OMP END DO

!   a l ouest
!------------
!$OMP DO SCHEDULE(runtime) private(nu0)
   DO j=ljmin,ljmax
     IF (l_spongw) THEN
       IF (h0_g(ig(j),j) > -valmanq) THEN
         DO i=ig(j),MAX(MIN(ig(j)+sponge_nbcell,imax),imin)

           IF ((i.GE.limin).AND.(i.LE.limax)) THEN

             nu0=fvisc*0.01_rsh*dy(i,j)**1.15_rsh
             vish_phi(:,i,j)=nu0*vismul
             vish_xe (:,i,j)=nu0*vismul

           END IF

         END DO
       END IF
     END IF

!   a l est
!----------
     IF (l_sponge) THEN
       IF (h0_g(id(j),j) > -valmanq) THEN
         DO i=MIN(MAX(id(j)-sponge_nbcell,imin),imax),id(j)

           IF ((i.GE.limin).AND.(i.LE.limax)) THEN

             nu0=fvisc*0.01_rsh*dy(i,j)**1.15_rsh
             vish_phi(:,i,j)=nu0*vismul
             vish_xe (:,i,j)=nu0*vismul

           END IF

         END DO
       END IF
     END IF
   END DO
!$OMP END  DO

!   au nord
!----------
!$OMP DO SCHEDULE(runtime) PRIVATE(nu0)
   DO i=liminp1,limax
     IF (l_spongn) THEN
       IF (h0_g(i,jh(i)) > -valmanq) THEN
         DO j=MIN(MAX(jh(i)-sponge_nbcell,jmin),jmax),jh(i)

           IF ((j.GE.ljmin).AND.(j.LE.ljmax)) THEN

             nu0=fvisc*0.01_rsh*dy(i,j)**1.15_rsh
             vish_phi(:,i,j)=nu0*vismul
             vish_xe (:,i,j)=nu0*vismul

           END IF

         END DO
       END IF
     END IF

!   au sud
!----------
     IF (l_spongs) THEN
       IF (h0_g(i,jb(i)) > -valmanq) THEN
         DO j=jb(i),MAX(MIN(jb(i)+sponge_nbcell,jmax),jmin)

           IF ((j.GE.ljmin).AND.(j.LE.ljmax)) THEN

             nu0=fvisc*0.01_rsh*dy(i,j)**1.15_rsh
             vish_phi(:,i,j)=nu0*vismul
             vish_xe (:,i,j)=nu0*vismul

           END IF

         END DO
       END IF
     END IF
   END DO
!$OMP END  DO

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-1,0,kmax,liminm1,limax  ,ljminm1,ljmax  ,vish_phi(:,liminm1:limax  ,ljminm1:ljmax  ))
   CALL ex_j_rsh(-1,0,kmax,liminm1,limax  ,ljminm1,ljmax  ,vish_phi(:,liminm1:limax  ,ljminm1:ljmax  ))
   CALL ex_i_rsh( 0,1,kmax,limin  ,limaxp1,ljmin  ,ljmaxp1,vish_xe (:,limin  :limaxp1,ljmin  :ljmaxp1))
   CALL ex_j_rsh( 0,1,kmax,limin  ,limaxp1,ljmin  ,ljmaxp1,vish_xe (:,limin  :limaxp1,ljmin  :ljmaxp1))
!$OMP end master
!$OMP barrier
!$OMP flush(vish_phi,vish_xe)

   !PRINT*, 'exit TURB_SMAGO routine'

  END SUBROUTINE turb_smago
