  SUBROUTINE dyn2dy_tls_solve(xeuv)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2dy_tls_solve  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : y2dxe,y2dv,axe,bxe,cxe,cv,bv,av
!&E
!&E ** Modified variables : xeuv,uv (1D)
!&E
!&E ** Reference :
!&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
!&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
!&E
!&E ** History :
!&E       !  2010-07  (L. Debreu, F. Dumas) rewritting of ADI
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters
   USE comvars2d,   ONLY : jb,jh

   USE toolmpi,     ONLY : ciminp2,cimaxm1
   USE comcolig_mpi,ONLY : y2dxe=>cy2dxe,axe=>caxe,bxe=>cbxe,cxe=>ccxe


   IMPLICIT NONE

!! * Arguments

   REAL(kind=rsh),DIMENSION(ciminp2:cimaxm1,jmin:jmax),INTENT(inout)   :: xeuv


!! * Local declarations
   INTEGER                                       :: i,j,k,n
   REAL(kind=rlg)                                :: xk
   REAL(kind=rlg),DIMENSION(jmax-jmin-1)         :: p0,q0,u0

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


!$OMP DO SCHEDULE(RUNTIME)
   DO i=ciminp2,cimaxm1


! nothing to do if no wet point all along the column
     IF (jb(i) < jh(i)) THEN

       n=jh(i)-jb(i)-1

! factorisation du systeme lineaire tridiagonal en LU
!----------------------------------------------------

       q0(1)=-cxe(i,jb(i)+1)/bxe(i,jb(i)+1)
       DO k=jb(i)+2,jh(i)-1
         p0(k-jb(i)  )=1.0_rlg/(axe(i,k)*q0(k-jb(i)-1)+bxe(i,k))
         q0(k-jb(i)  )=-cxe(i,k)*p0(k-jb(i)  )
       END DO


! resolution du systeme triangulaire inferieur
!---------------------------------------------

       u0(1)=y2dxe(i,jb(i)+1)/bxe(i,jb(i)+1)
       DO k=jb(i)+2,jh(i)-1
         u0(k-jb(i)  )=(y2dxe(i,k)-axe(i,k)*u0(k-jb(i)-1))*p0(k-jb(i)  )
       END DO

! resolution du systeme triangulaire superieur
!---------------------------------------------

       xk=0.0_rlg
       k=n+1

       DO j=jh(i)-1,jb(i)+1,-1
         k=k-1
         xk=q0(k)*xk+u0(k)
         xeuv(i,j)=xk
       END DO

     END IF

   END DO
!$OMP END DO

   RETURN

  END SUBROUTINE dyn2dy_tls_solve
