  SUBROUTINE turbadv_wcen(xe2,xe3,ect,psi)

#ifdef key_turbadv
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turbadv_wcen  ***
   !&E
   !&E ** Purpose : transport computation for turbulent kinetic energy and 
   !&E              turbulent dissipation
   !&E
   !&E ** Description : advection of turbulence through a centered scheme
   !&E
   !&E ** Called by : step3dxy.F90, step3dyx.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : xe2,xe3,ect,psi
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E        !   201-08 (F. Ganthy) Original code
   !&E
   !&E---------------------------------------------------------------------

# include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(limin:limax,ljmin:ljmax),INTENT(in)         :: xe2
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe3
   
   REAL(kind=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax),INTENT(inout) :: ect,psi

   !! * Local declarations
   INTEGER         :: i,j,k,l
   REAL(kind=rsh)  :: dtm,d2,d3,unsd3
   REAL(kind=rlg)  :: xk,yk
   REAL(kind=rlg),DIMENSION(kmax)   :: a,b,c,u0,v0,p0,yect,ypsi
   REAL(kind=rlg),DIMENSION(0:kmax) :: q0

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

   PRINT_DBG*, 'enter TURBADV_wcen routine'

! Flo2val : verifier la parallelisation OMP et MPI

!$OMP PARALLEL DEFAULT(none) PRIVATE(i,j,k,l,dtm,d2,d3,unds3,xk,yk,a,b,c,u0,   &
!$OMP& v0,p0,yect,ypsi,q0)
!$OMP SHARED(halfdt,ect,psi,xe2,xe3,hm,ig,id,jh,jb,h0,limin,limax,ljmin,ljmax, &
!$OMP& liminm1,limaxp2,ljminp2,ljmaxm1,ljminm1,ljmaxp2,fwet,fwetp,wz,dsigw,    &
!$OMP& dx_v,dy_u,surf,fxect,fyect,fxpsi,fypsi,kz)

   dtm = halfdt

   !-------------------------------------------------------------
   ! Computation of horizontal fluxes (advection and diffusion)
   !-------------------------------------------------------------
#if ! defined key_1DV
   fxect(:,:,:) = 0.0_rlg
   fyect(:,:,:) = 0.0_rlg
   fxpsi(:,:,:) = 0.0_rlg
   fypsi(:,:,:) = 0.0_rlg
   CALL turbadv_difx(ect,psi)
   CALL turbadv_dify(ect,psi)
#endif

!$OMP DO SCHEDULE(RUNTIME)

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

       d3 = h0(i,j)+xe3(i,j)
       d2 = h0(i,j)+xe2(i,j)

       IF(d3.GT.hm) THEN
         unsd3 = 1.0_rsh/d3

         !***************
         k=1
         !***************
         c(k) = dtm*d3*0.5_rlg*wzsw(k,i,j)/dsigw(k)
         b(k) = d3*fwet(i,j) + dtm*d3*0.5_rlg*wzsw(k,i,j)/dsigw(k)

         yect(k) = d2*fwetp(i,j)*ect(k,i,j) - &
                   (dy_u(i,j)*fxect(k,i,j) - dy_u(i-1,j)*fxect(k,i-1,j) + &
                    dx_v(i,j)*fyect(k,i,j) - dx_v(i,j-1)*fyect(k,i,j-1))/surf(i,j)

         ypsi(k) = d2*fwetp(i,j)*psi(k,i,j) - &
                   (dy_u(i,j)*fxpsi(k,i,j) - dy_u(i-1,j)*fxpsi(k,i-1,j) + &
                    dx_v(i,j)*fypsi(k,i,j) - dx_v(i,j-1)*fypsi(k,i,j-1))/surf(i,j)

         q0(k) = -c(k)/b(k)    ! LU factorization
         u0(k) = yect(k)/b(k)  ! TKE
         v0(k) = ypsi(k)/b(k)  ! PSI

         !***************
         DO k = 2,kmax-1
         !***************
         c(k) =  dtm*d3*0.5_rlg*wzsw(k,i,j)/dsigw(k)
         a(k) = -dtm*d3*0.5_rlg*wzsw(k-1,i,j)/dsigw(k)
         b(k) = d3*fwet(i,j) +                          &
                dtm * (d3*0.5_rlg*wzsw(k,i,j) - d3*0.5_rlg*wzsw(k-1,i,j))/dsigw(k)

         yect(k) = d2*fwetp(i,j)*ect(k,i,j) - &
                   (dy_u(i,j)*fxect(k,i,j) - dy_u(i-1,j)*fxect(k,i-1,j) + &
                    dx_v(i,j)*fyect(k,i,j) - dx_v(i,j-1)*fyect(k,i,j-1))/surf(i,j)

         ypsi(k) = d2*fwetp(i,j)*psi(k,i,j) - &
                   (dy_u(i,j)*fxpsi(k,i,j) - dy_u(i-1,j)*fxpsi(k,i-1,j) + &
                    dx_v(i,j)*fypsi(k,i,j) - dx_v(i,j-1)*fypsi(k,i,j-1))/surf(i,j)

         p0(k) =  a(k)*q0(k-1)+b(k) ! LU factorization
         q0(k) = -c(k)/p0(k)
         u0(k) = (yect(k)-a(k)*u0(k-1))/p0(k) ! TKE
         v0(k) = (ypsi(k)-a(k)*v0(k-1))/p0(k) ! PSI 
         !***************
         END DO
         !***************


         !-------------------------
         ! Estimate TKE and PSI
         !-------------------------
         xk = u0(kmax-1)
         ect(kmax-1,i,j) = xk
         yk = v0(kmax-1)
         psi(kmax-1,i,j) = yk

         DO k = kmax-2,1,-1
           xk = q0(k)*xk+u0(k)
           ect(k,i,j) = xk
           yk = q0(k)*yk+v0(k)
           psi(k,i,j) = yk
         ENDDO

       ENDIF ! End test on d3>hm
     ENDDO ! End loop on i
   ENDDO ! End loop on j
!$OMP END DO

!$OMP END PARALLEL

   PRINT_DBG*, 'exit TURBADV_wcen routine'
#endif

   END SUBROUTINE turbadv_wcen
