  SUBROUTINE turbadv_uzvzwz(xe2,xe3,uz,vz)

#ifdef key_turbadv
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turbadv_uzvzwz  ***
   !&E
   !&E ** Purpose : - re-sample uz and vz to vertical half-grid cells (uzdw,vzdw)
   !&E              - estimate the vertical velocity at vertical half-grid cells (wzdw)
   !&E                from the vertical integration of the local continuity equation
   !&E
   !&E ** Description : progonstic estimate for advection of turbulence
   !&E
   !&E ** Called by : step3dxy.F90, step3dyx.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E        !   2015-08 (F. Ganthy) Original code
   !&E
   !&E---------------------------------------------------------------------

# include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : imin,imax,jmin,jmax,kmax,ljmin,ljmax,rsh,            &
                          ljminm1,ljminp2,ljmaxm1,limin,limax,limaxp1,ljmaxp1, &
                          liminm2,limaxp2,ljminm2,ljmaxp2,liminm1,ljminm1,     &
                          ishift_obc_west,ishift_obc_east,jshift_obc_south,    &
                          jshift_obc_north
   USE comvars2d, ONLY  : ig,id,h0,dy_u,dx_v,l_obc_cycl_y,l_obc_cycl_x,jb,jh,  &
                          surf,halfdt,hm,fwet,fwetp,nbouc,hx,hy,dhx,dhy,       &
                          goulet,govlet,hex,hey
   USE comvars3d, ONLY  : sig,sigw,dsigw,dsigu,wz,uzsw,vzsw,wzsw,fxsw,fysw,phisw
   USE comfluxdo, ONLY  : fxeau,fyeau,phieau
   USE toolgeom, ONLY   : f_hzw_x,f_hzw_y,f_hzw,f_hzu

#ifdef key_MPI_2D
   USE toolmpi, ONLY : ex_j_rsh, ex_i_rsh,ex_peri_x_rsh,ex_peri_y_rsh
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT (in)        :: xe2,xe3
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT (in)   :: uz,vz


   !! * Local declarations
   INTEGER        :: k,i,j
   REAL(KIND=rsh) :: hp,xet,div,divt,hz3_w

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

   PRINT_DBG*, 'enter TURBADV_uzvzwz routine'

! Flo2val : verifier la parallelisation OMP et MPI

!$OMP PARALLEL DEFAULT(none) SHARED(dsigw,dsigu,ljminp2,ljmaxm1,ljmin,       &
!$OMP& ljmax,limin,limax,ig,id,jh,jb,uz,vz,hm,h0,xe2,xe3,uzsw,vzsw,wzsw,     &
!$OMP& fxeau,fyeau,halfdt,fwet,fwetp,phieau,fxsw,fysw,phisw,ishift_obc_west, &
!$OMP& ishift_obc_east,jshift_obc_south,jshift_obc_north,l_obc_cycl_x,l_obc_cycl_y)
!$OMP& PRIVATE(i,j,k,hp,div,divt,ext)

!$OMP DO SCHEDULE(runtime)
   DO j=ljminp2,ljmaxp1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       hp = h0(i,j)+xe3(i,j)
       IF (hp.LT.hm) THEN
         uzsw(:,i,j)  = 0.0_rsh
         vzsw(:,i,j)  = 0.0_rsh
         wzsw(:,i,j)  = 0.0_rsh
         fxsw(:,i,j)  = 0.0_rsh
         fysw(:,i,j)  = 0.0_rsh
         phisw(:,i,j) = 0.0_rsh
       ELSE
         !----------------------------------------------------------------
         ! Re-sampling of horizontal velocity components and computation of 
         ! fxeau,fyeau,phisw
         !----------------------------------------------------------------
         !****************
         phisw(:,i,j) = 0.0_rsh
         DO k=1,kmax-1 ! water column
           uzsw(k,i,j)  =  uz(k,i,j)   * (sig(k+1) - sigw(k)) / (sig(k+1) - sig(k)) + &
                           uz(k+1,i,j) * (sigw(k)  - sig(k))  / (sig(k+1) - sig(k))
           vzsw(k,i,j)  =  vz(k,i,j)   * (sig(k+1) - sigw(k)) / (sig(k+1) - sig(k)) + &
                           vz(k+1,i,j) * (sigw(k)  - sig(k))  / (sig(k+1) - sig(k))
#ifdef key_siggen || key_gencoord
           fxsw(k,i,j)  = uzsw(k,i,j) * halfdt * &
                          f_hzw_x(hx(i,j),0.5_rsh*(xe2(i,j)+xe2(i+1,j))+dhx(i,j),k,i,j)*goulet(i,j)
           fysw(k,i,j)  = vzsw(k,i,j) * halfdt * &
                          f_hzw_y(hy(i,j),0.5_rsh*(xe2(i,j)+xe2(i,j+1))+dhy(i,j),k,i,j)*govlet(i,j)
#else
           fxsw(k,i,j)  = uzsw(k,i,j) * halfdt * hex(i,j)
           fysw(k,i,j)  = vzsw(k,i,j) * halfdt * hey(i,j)
#endif
           phisw(k,i,j) =  phieau(k,i,j)   * (sig(k+1) - sigw(k)) / (sig(k+1) - sig(k)) + &
                           phieau(k+1,i,j) * (sigw(k)  - sig(k))  / (sig(k+1) - sig(k))
         ENDDO ! End loop on k=1,kmax-1
       ENDIF ! End test on hp<hm

       !----------------------------------------------------------------------
       ! Estimate of the vertical velocity from the local continuity equation
       !----------------------------------------------------------------------
       IF (hp.LT.0.0_rsh) THEN
         wzsw(:,i,j) = 0.0_rsh
       ELSE
#ifdef key_siggen || key_gencoord
         divt = 0.0_rsh
         xet=(xe3(i,j)-xe2(i,j)*fwetp(i,j)/fwet(i,j))/halfdt
         !******************
         k = 1
         hz3_w = f_hzw(h0(i,j),xe3(i,j),k,i,j)
         IF(h0(i,j).LE.hm) hz3_w = hm/REAL(kmax,rsh)
         div = (fxeau(k,i,j)*dy_u(i,j)-fxeau(k,i-1,j)*dy_u(i-1,j))/surf(i,j)/halfdt &
               +(fyeau(k,i,j)*dx_v(i,j)-fyeau(k,i,j-1)*dx_v(i,j-1))/surf(i,j)/halfdt &
               -phieau(k,i,j)/(surf(i,j)*fwet(i,j))/dsigu(k)
         divt = divt + 0.5_rsh * (div*dsigu(k))
         wzsw(k,i,j) = -(xet*(sig(k)+1.0_rsh)+divt)/hz3_w
         !******************
         DO k=2,kmax
           hz3_w = f_hzu(h0(i,j),xe3(i,j),k,i,j)
           IF(h0(i,j).LE.hm) hz3_w = hm/REAL(kmax,rsh)
           div = (fxsw(k-1,i,j)*dy_u(i,j)-fxsw(k-1,i-1,j)*dy_u(i-1,j))/surf(i,j)/halfdt &
                +(fysw(k-1,i,j)*dx_v(i,j)-fysw(k-1,i,j-1)*dx_v(i,j-1))/surf(i,j)/halfdt &
                -phisw(k-1,i,j)/(surf(i,j)*fwet(i,j))/dsigw(k-1)
           divt = divt + div*dsigw(k-1)
           wzsw(k,i,j) = -(xet*(sig(k)+1.0_rsh)+divt)/hz3_w
         ENDDO ! End loop on k=2,kmax
#else
         divt = 0.0_rsh
         xet=(xe3(i,j)-xe2(i,j)*fwetp(i,j)/fwet(i,j))/halfdt
         !******************
         k = 1
         div = (fxeau(k,i,j)*dy_u(i,j)-fxeau(k,i-1,j)*dy_u(i-1,j))/surf(i,j)/halfdt &
               +(fyeau(k,i,j)*dx_v(i,j)-fyeau(k,i,j-1)*dx_v(i,j-1))/surf(i,j)/halfdt &
               -phieau(k,i,j)/(surf(i,j)*fwet(i,j))/dsigu(k)
         divt = divt + 0.5_rsh * (div*dsigu(k))
         wzsw(k,i,j) = -(xet*(sig(k)+1.0_rsh)+divt)/hp
         !******************
         DO k=2,kmax
           div = (fxsw(k-1,i,j)*dy_u(i,j)-fxsw(k-1,i-1,j)*dy_u(i-1,j))/surf(i,j)/halfdt &
                +(fysw(k-1,i,j)*dx_v(i,j)-fysw(k-1,i,j-1)*dx_v(i,j-1))/surf(i,j)/halfdt &
                -phisw(k-1,i,j)/(surf(i,j)*fwet(i,j))/dsigw(k-1)
           divt = divt + div*dsigw(k-1)
           wzsw(k,i,j) = -(xet*(sig(k)+1.0_rsh)+divt)/hp
         ENDDO ! End loop on k=2,kmax
#endif
       ENDIF ! End test on hp<0
     ENDDO ! End loop on i
   ENDDO ! End loop on j
!$OMP END DO

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
     CALL ex_peri_x_rsh(1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
#else
     DO j=ljmin,ljmax
       wzsw(:,ig(j),j) = wzsw(:,id(j)-1,j)
       wzsw(:,id(j),j) = wzsw(:,ig(j)+1,j)
     ENDDO
#endif
   ELSE
     IF (limax==imax) wzsw(:,imax-ishift_obc_east,:)=wzsw(:,imax-1-ishift_obc_east,:)
     IF (limin==imin) wzsw(:,imin+1+ishift_obc_west,:)=wzsw(:,imin+2+ishift_obc_west,:)
   END IF
   IF (l_obc_cycl_y) THEN
#ifdef key_MPI_2D
     CALL ex_peri_y_rsh(1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
#else
     DO i=limin,limax
       wzsw(:,i,jb(i))=wzsw(:,i,jh(i)-1)
       wzsw(:,i,jh(i))=wzsw(:,i,jb(i)+1)
     END DO
#endif
   ELSE
     IF (ljmin==jmin) wzsw(:,limin:limax,jmin+1+jshift_obc_south)=wzsw(:,limin:limax,jmin+2+jshift_obc_south)
     IF (ljmax==jmax) wzsw(:,limin:limax,jmax-jshift_obc_north)=wzsw(:,limin:limax,jmax-1-jshift_obc_north)
   END IF
!$OMP END SINGLE

!$OMP END PARALLEL

   CALL_MPI ex_i_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wzsw(:,liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL_MPI ex_j_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wzsw(:,liminm1:limaxp1,ljminm1:ljmaxp1))


   PRINT_DBG*, 'exit TURBADV_uzvzwz routine'
#endif

   END SUBROUTINE turbadv_uzvzwz
