  SUBROUTINE ts_advz_compact(xe,sal,temp)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ts_advz_compact  ***
!&E
!&E ** Purpose : compute vertical fluxes of temperature and salinity
!&E              using an Upwind Compact and Conservative scheme
!&E              adapted to an euler time scheme
!&E
!&E ** Description : For very low water heights ( < hm), no vertical fluxes are computed
!&E              because the scheme is degenerated vertically and vertical fluxes are useless.
!&E              This routine is called in "adv" only if the key
!&E              1 is activated (see tssub_eq_wcompact.h).
!&E
!&E ** Called by : adv
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xe,h0,fzeau,wz,sal,temp
!&E
!&E ** Modified variables : ftempz,fsalz
!&E
!&E ** Reference : Debreu, 2008 (note : UCC : Upwind Compact and Conservative)
!&E
!&E ** History :
!&E       !  2010-05 (L. Debreu, F. Dumas) Original code
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2012-08  (V. Garnier) Use of w(k,i,j) to estimate cfl (Courant–Friedrichs–Lewy) variable
!&E                                rewritting of the coefficients
!&E                                review of the OBC for the estimate of interpolated concentrations
!&E
!&E---------------------------------------------------------------------
!! * Modules used

   USE parameters
   USE comvars2d,  ONLY : id,ig,h0,hm,halfdt
   USE comvars3d,  ONLY : dsigu,wz
   USE comfluxdo,  ONLY : fzeau
   USE comfluvert, ONLY : fsalz,ftempz


   USE comvars3d,  ONLY : f_alpha_T, theta_wz


   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)      :: xe
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: sal,temp

!! * Local declarations
   INTEGER                                  :: i,j,k,n
   REAL(kind=rsh)                           :: h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2
   REAL(kind=rsh),DIMENSION(kmax)           :: ats,bts,cts,y2dt,y2ds,salhat,temphat

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

!$OMP DO SCHEDULE(RUNTIME) &
!$OMP& PRIVATE(i,j,k,n,h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2 &
!$OMP&, ats,bts,cts,y2dt,y2ds,salhat,temphat)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)

! coefficients of Upwind Compact and Conservative scheme
!-------------------------------------------------------

       h3=h0(i,j)+xe(i,j)
       IF (h3 > hm) THEN
       DO k=1,kmax-1
!cfl=(wz(k,i,j)+wz(k-1,i,j))*0.5_rsh*REAL(halfdt,rsh)/dsigu(k)
         cfl=wz(k,i,j)*REAL(halfdt,rsh)/dsigu(k)

         IF(wz(k,i,j) >= 0) THEN
           cfl = cfl * f_alpha_T(k  ,i,j)
         ELSE
           cfl = cfl * f_alpha_T(k+1,i,j)
         ENDIF

         IF (cfl > 0.0_rsh) THEN
           alpha1=3.0_rsh*(cfl-3.0_rsh)*(cfl-2.0_rsh)
           alpha2=3.0_rsh*(cfl+1.0_rsh)*(cfl+2.0_rsh)
           alpha0=-6.0_rsh*(cfl-3.0_rsh)*(cfl+2.0_rsh)
           gamma1=38.0_rsh-2.0_rsh*cfl**2
           gamma2=20.0_rsh+(cfl-3.0_rsh)*cfl
           mu1=(cfl+1.0_rsh)*(cfl+2.0_rsh)
           mu2=0.0_rsh
         ELSE
           alpha1=3.0_rsh*(cfl-1.0_rsh)*(cfl-2.0_rsh)
           alpha2=3.0_rsh*(cfl+3.0_rsh)*(cfl+2.0_rsh)
           alpha0=-6.0_rsh*(cfl-2.0_rsh)*(cfl+3.0_rsh)
           gamma1=20.0_rsh+(cfl+3.0_rsh)*cfl
           gamma2=38.0_rsh-2.0_rsh*cfl**2
           mu1=0.0_rsh
           mu2=(cfl-1.0_rsh)*(cfl-2.0_rsh)
         END IF

         ats(k)=alpha1
         bts(k)=alpha0
         cts(k)=alpha2

         IF (k >= 2 .AND. k <= kmax-2) THEN
           y2ds(k)=mu1*sal(k-1,i,j) +gamma1*sal(k,i,j) +gamma2*sal(k+1,i,j) +mu2*sal(k+2,i,j)
           y2dt(k)=mu1*temp(k-1,i,j)+gamma1*temp(k,i,j)+gamma2*temp(k+1,i,j)+mu2*temp(k+2,i,j)
         END IF
         IF (k==1) THEN
           y2ds(k)=(mu1+gamma1-alpha1)*sal(k,i,j) +gamma2*sal(k+1,i,j) +mu2*sal(k+2,i,j)
           y2dt(k)=(mu1+gamma1-alpha1)*temp(k,i,j)+gamma2*temp(k+1,i,j)+mu2*temp(k+2,i,j)
         END IF
         IF (k==kmax-1) THEN
           y2ds(k)= mu1*sal(k-1,i,j)+gamma1*sal(k,i,j) +(mu2+gamma2-alpha2)*sal(k+1,i,j)
           y2dt(k)=mu1*temp(k-1,i,j)+gamma1*temp(k,i,j)+(mu2+gamma2-alpha2)*temp(k+1,i,j)
         END IF


       END DO

       n=kmax-1
       CALL  inversion_tridiag(ats(1:kmax-1),bts(1:kmax-1),cts(1:kmax-1), salhat(1:kmax-1),y2ds(1:kmax-1),n)
       CALL  inversion_tridiag(ats(1:kmax-1),bts(1:kmax-1),cts(1:kmax-1),temphat(1:kmax-1),y2dt(1:kmax-1),n)

! pour faire simple : decentre amont pres des limites (fond et surface)
!----------------------------------------------------------------------

       ftempz(0,i,j)=0.0_rsh
       ftempz(kmax,i,j)=0.0_rsh
       fsalz(0,i,j)=0.0_rsh
       fsalz(kmax,i,j)=0.0_rsh

       DO k=1,kmax-1

         IF(fzeau(k,i,j) >=0) THEN
           fsalz(k,i,j) =fzeau(k,i,j)*salhat(k)*f_alpha_T(k,i,j)+ &
                       (max(fzeau(k,i,j),0.0_rlg)*sal(k,i,j)+min(fzeau(k,i,j),0.0_rlg)*sal(k+1,i,j))* &
                       (1.0_rsh-f_alpha_T(k,i,j))*(1.0_rsh-theta_wz)
           ftempz(k,i,j)=fzeau(k,i,j)*temphat(k)*f_alpha_T(k,i,j)+ &
                       (max(fzeau(k,i,j),0.0_rlg)*temp(k,i,j)+min(fzeau(k,i,j),0.0_rlg)*temp(k+1,i,j))* &
                       (1.0_rsh-f_alpha_T(k,i,j))*(1.0_rsh-theta_wz)
         ELSE
           fsalz(k,i,j) =fzeau(k,i,j)*salhat(k)*f_alpha_T(k+1,i,j)+ &
                       (max(fzeau(k,i,j),0.0_rlg)*sal(k,i,j)+min(fzeau(k,i,j),0.0_rlg)*sal(k+1,i,j))* &
                       (1.0_rsh-f_alpha_T(k+1,i,j))*(1.0_rsh-theta_wz)
           ftempz(k,i,j)=fzeau(k,i,j)*temphat(k)*f_alpha_T(k+1,i,j)+ &
                       (max(fzeau(k,i,j),0.0_rlg)*temp(k,i,j)+min(fzeau(k,i,j),0.0_rlg)*temp(k+1,i,j))* &
                       (1.0_rsh-f_alpha_T(k+1,i,j))*(1.0_rsh-theta_wz)
         ENDIF

       END DO

       ENDIF     ! fin du test sur la hauteur d eau
     END DO     ! fin de la boucle en i
   END DO     ! fin de la boucle en j
!$OMP end do

  END SUBROUTINE ts_advz_compact

!****************************************
  SUBROUTINE inversion_tridiag(a,b,c,x,y,n)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE inversion_tridiag  ***
!&E
!&E ** Purpose : inverse tridiagonal system
!&E
!&E ** Description : a   : sous-diagonale
!&E                  b   :      diagonale
!&E                  c   : sur-diagonale
!&E                  y2d : second membre
!&E                  x   : inconnue
!&E
!&E ** Called by : ts_advz_compact
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xe,h0,fzeau,wz,sal,temp
!&E
!&E ** Modified variables : ftempz,fsalz
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2010-05 (L. Debreu, F. Dumas) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used

   USE parameters
   IMPLICIT NONE

!! * Arguments
   INTEGER,                    INTENT(IN) :: n
   REAL(KIND=rsh),DIMENSION(n),INTENT(IN) :: a,b,c,y
   REAL(KIND=rsh),DIMENSION(n),INTENT(OUT) :: x

!! * Local declarations
   INTEGER :: k
   REAL(KIND=rsh),DIMENSION(n) :: u0,p0,q0
   REAL(KIND=rsh)              :: xk

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

   q0(1)=-c(1)/b(1)
   u0(1)=y(1)/b(1)
   DO k=2,n-1
     p0(k)=a(k)*q0(k-1)+b(k)
     q0(k)=-c(k)/p0(k)
     u0(k)=(y(k)-a(k)*u0(k-1))/p0(k)
   END DO
   p0(n)=a(n)*q0(n-1)+b(n)
   u0(n)=(y(n)-a(n)*u0(n-1))/p0(n)

   xk=u0(n)
   x(n)=xk
   DO k=n-1,1,-1
     xk=q0(k)*xk+u0(k)
     x(k)=xk
   END DO

  END SUBROUTINE inversion_tridiag
