  SUBROUTINE dyn2dy_cyc_solve(xeuv)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2dy_cyc_solve  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : xeuv
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2010-07  (L. Debreu) rewritting of ADI - substitution of Lapack routine
!&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,ktemp
   REAL(kind=rlg)                                            :: xk,p0
   REAL(kind=rlg),DIMENSION(2*(jmax-jmin)-3)                 :: q0,s0,u0

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


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


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

! introduction des conditions aux limites
!----------------------------------------
! factorisation du systeme lineaire tridiagonal en LU
!----------------------------------------------------

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


       k = jb(i)+1
       
       q0(1)=-cxe(i,k)/bxe(i,k)
       s0(1)=-axe(i,k)/bxe(i,k)
       u0(1)=y2dxe(i,k)/bxe(i,k)

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



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

       k=n
       q0(k)=0.0_rlg
       s0(k)=1.0_rlg
       
       DO j=jh(i)-2,jb(i)+1,-1
         k=k-1
         s0(k) = s0(k)  + q0(k)*s0(k+1)
         q0(k) = u0(k)  + q0(k)*q0(k+1)
       END DO

       k=jh(i)-1
       
       xk=(y2dxe(i,k)-cxe(i,k)*q0(1)-axe(i,k)*q0(n-1))/(cxe(i,k)*s0(1)+axe(i,k)*s0(n-1)+bxe(i,k))
       
       xeuv(i,k)=xk
       
       k = 0
       
       DO j=jb(i)+1,jh(i)-2
         k=k+1
         xeuv(i,j)=xk*s0(k)+q0(k)
       END DO

   END DO
!$OMP END DO

   RETURN

  END SUBROUTINE dyn2dy_cyc_solve
