  SUBROUTINE dyn2dx_tls_coef(u,v,xev,ustar,vstar)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2dx_tls_coef  ***
!&E
!&E ** Purpose : Estimate of the coefficients of the tridiagonal
!&E              linear system and the RHS relative to the
!&E              momentum and continuity equations
!&E
!&E ** Description : estimate along x-direction
!&E
!&E ** Called by : dyn2dxy.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xev,u,v,ustar,vstar,igdu,hx,htx,hty,dx,dx_u,dx_v,
!&E                     dy_f,dy_u,fu,hey,hex,surf,phieau,
!&E                     tfiusr,vis2d,tosx,pa,potgene,buoy2d
!&E
!&E ** Modified variables : y2dxe,y2dv,axe,bxe,cxe,cv,bv,av
!&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        !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E        !  2009-01  (A. Mano) MPI
!&E        !  2009-01  (V. Garnier) OMP
!&E        !  2009-01  (A. Mano T. Odaka) MPI2D
!&E        !  2009-07  (N. Bruneau F. Dumas) coupling with swan
!&E        !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E        !  2014-07  (F. Ganthy) Obstruction-induced force on momentum
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d, ONLY  : botstressmax,tfiusr,ig,id,jb,jh,          &
                          htx,hty,hx,hey,hex,dx,dx_u,dy_f,surf,     &
                          dx_v,dx_f,dy_v,dy_u,fu,vis2d,halfdt,      &
                          hminfrot,dt,icon,h0,hminfrot,hmaxfrot,z0b,&
                          fwet,fwetp,frot_lin
   USE comvarp2d,  ONLY : v_utmp
   USE comfluxdo,  ONLY : nbfl,buoy2d,phieau
   USE comaree,    ONLY : l_genpot,potgene
   USE cometdyn,   ONLY : tosx,pa
   USE comvars3d,  ONLY : rhoref
   USE comcolig2d, ONLY : y2dxe,y2dv,axe,bxe,cxe,cv,bv,av,igdu,mn,ms,whx




   IMPLICIT NONE

!! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xev
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: u,v,ustar,vstar

!! * Local declarations
   INTEGER                                       :: i,j,k
   REAL(KIND=rlg)                                :: phieautot              !dyn2dx_cty.h
   REAL(KIND=rsh)  :: ft,gdt,vit,factft,vdusdy,vdudy,vdudy2,   &           !dyn2dx_mtm.h
                      udusdx,ududx,ududx2,dusdy,dusdy2,v_utmp2

   REAL(KIND=rsh)  :: uc,up,um,vp,vm                                       !dyn2dx_mtm.h


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

   !PRINT*, 'enter DYN2DX_TLS_COEF routine'


!$OMP DO SCHEDULE(RUNTIME)

   DO j=ljminp1,ljmax
     DO i =MAX0(ig(j)-1,limin),MIN0(id(j)-1,limax)
       v_utmp(i,j)=0.25_rsh*(v(i,j)+v(i+1,j)+v(i,j-1)+v(i+1,j-1))
     END DO
   END DO
!$OMP END DO




!$OMP DO SCHEDULE(RUNTIME)

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

!&E   ***************   include dyn2dx_mtm   **********************************
!&E
!&E ** Purpose : calcul des coefficients de la matrice tridiagonale
!&E              du systeme lineaire et du second membre relatifs
!&E              a l equation de conservation de la quantite de mouvement
!&E              pour le calcul en ligne. (etape de prediction dans
!&E              l algorithme de predicteur-correcteur)
!&E
!&E ** Description :
!&E
!&E ** Called by : dyn2dyx.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xev,u,v,ustar,vstar,igdu,hx,htx,hty,
!&E                     dx,dx_u,dy_f,dy_u,fu,ms,mn,whx,
!&E                     tfiusr,vis2d,tosx,pa,potgene,buoy2d
!&E
!&E ** Modified variables : y2dv,cv,bv,av
!&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       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E       !  2009-03  (V. Garnier) modification of vdusdy term (centered scheme)
!&E       !  2009-03  (T. Duhaut V. Garnier) QUICK and QUICKEST scheme (with masks)
!&E       !  2009-07  (N. Bruneau F. Dumas) coupling with swan
!&E       !  2010-07  (L. Debreu, F. Dumas) rewritting of ADI
!&E       !  2011-03  (F. Dumas)  bottom stress parametrisation + remove case 2 and 3
!&E       !  2011-10  (L. Pineau, F. Dumas, V. garnier) correction bottom stress parametrisation
!&E       !  2011-11  (F. Vandermeirsch) correct spatial extension on which we can estimate quick and quickest scheme
!&E       !  2014-07  (F. Ganthy) Obstruction-induced force on momentum
!&E       !  2015-03  (F. Ganthy) Bug correction for implicit part of obstruction induced friction
!&E       !  2017-12  (F. Ganthy) Bug corrected for obstruction-induced force
!&E
!&E
!&E  RQ V Garnier vis2d est homogene spatialement donc OK
!&E     neanmoins, il faudrait plutot utiliser vish_xe(1,:,:) et vis
!&E     car ainsi le coef pourrait varier spatialement
!&E     Aujourd hui, faux si cela varie spatialement, il faut definir
!&E     coef de viscosite en u et v
!&E
!&E-------------------------------------------------------------------------

     gdt = grav*dt

! cas d un point u sans eau de part ni d autre
! drying u-point : no water on the right and left

     IF (igdu(i,j) == 4) THEN
       av(i,j)=0.0_rlg
       bv(i,j)=1.0_rlg
       cv(i,j)=0.0_rlg
       y2dv(i,j)=0.0_rlg

! cas d un point u avec de l eau de part et d autre
! u-point surrounded by water

     ELSE IF (igdu(i,j) == 1) THEN

       v_utmp2=0.25_rsh*(vstar(i,j)+vstar(i+1,j)+vstar(i,j-1)+vstar(i+1,j-1))

! sea surface pressure gradient : implicit part
       av(i,j)=-alphan*0.5_rlg*gdt/dx_u(i,j)                ! i-1
       cv(i,j)= alphan*0.5_rlg*gdt/dx_u(i,j)                ! i



! bottom stress



       vit=sqrt(u(i,j)**2+v_utmp(i,j)**2)
       factft=1.0_rsh+vit-botstressmax+abs(vit-botstressmax)
       ft=factft*0.5_rsh*dt*vit/htx(i,j)*MIN((0.4_rsh/LOG(MAX(MIN(htx(i,j),hmaxfrot),hminfrot)/(exp1*z0b(i,j))))**2,0.1_rsh)





! non-linear terms (udu/dy and vdu/dy) : explicit

       uc=ustar(i,j)
       up=0.5_rsh*(uc+abs(uc))
       um=0.5_rsh*(uc-abs(uc))
       vp=0.5_rsh*(v_utmp2+abs(v_utmp2))
       vm=0.5_rsh*(v_utmp2-abs(v_utmp2))
       IF (i >= imin+2.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+3) THEN
         ududx  =(up*u(i-2,j)*whx(i-1,j)+(-7.0_rsh*up-3.0_rsh*um)*u(i-1,j)+  &
                  (3.0_rsh*up-3.0_rsh*um)*u(i,j)+(3.0_rsh*up+7.0_rsh*um)*u(i+1,j)-um*u(i+2,j)*whx(i+1,j))/(8.0_rsh*dx_u(i,j))
         ududx2 =(up*ustar(i-2,j)*whx(i-1,j)+(-7.0_rsh*up-3.0_rsh*um)*ustar(i-1,j)+  &
                  (3.0_rsh*up-3.0_rsh*um)*ustar(i,j)+(3.0_rsh*up+7.0_rsh*um)*ustar(i+1,j)-um*ustar(i+2,j)*whx(i+1,j)) &
                 /(8.0_rsh*dx_u(i,j))
         udusdx = (alphan*ududx2 + alphap*ududx)
         vdudy  =(vp*u(i,j-2)*ms(i,j-1)*ms(i,j)+(-7.0_rsh*vp-3.0_rsh*vm)*u(i,j-1)+  &
                  (3.0_rsh*vp-3.0_rsh*vm)*u(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*u(i,j+1)-vm*u(i,j+2)*mn(i,j+1)*mn(i,j)) &
                 /(8.0_rsh*dy_u(i,j)) 
         vdudy2 =(vp*ustar(i,j-2)*ms(i,j-1)*ms(i,j)+(-7.0_rsh*vp-3.0_rsh*vm)*ustar(i,j-1)+  &
                  (3.0_rsh*vp-3.0_rsh*vm)*ustar(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*ustar(i,j+1)-vm*ustar(i,j+2)*mn(i,j+1)*mn(i,j)) &
                 /(8.0_rsh*dy_u(i,j)) 
         vdusdy = (alphan*vdudy2 + alphap*vdudy)
       ELSE
! upwind scheme
         ududx  = up*(u(i,j)-u(i-1,j))/dx(i,j)+um*(u(i+1,j)-u(i,j))/dx(i+1,j)
         ududx2 = up*(ustar(i,j)-ustar(i-1,j))/dx(i,j)+um*(ustar(i+1,j)-ustar(i,j))/dx(i+1,j)
         udusdx = (alphan*ududx2 + alphap*ududx)
         vdudy  = vp*(u(i,j)-u(i,j-1))/dy_f(i,j-1)+vm*(u(i,j+1)-u(i,j))/dy_f(i,j)
         vdudy2 = vp*(ustar(i,j)-ustar(i,j-1))/dy_f(i,j-1)+vm*(ustar(i,j+1)-ustar(i,j))/dy_f(i,j)
         vdusdy = (alphan*vdudy2 + alphap*vdudy)
       END IF



       bv(i,j)=1.0_rlg/(                                &
! time evolution (u at n+1/2)
               1.0_rlg                                  &

! bottom stress : implicit part
               +ft*cfroi)

! rhs : explicit terms
       y2dv(i,j)=                                       &
! time evolution (u at n-1/2)
                 u(i,j)                                 &
! bottom stress (u at n-1/2)
                 -ft*cfroe*u(i,j)                       &
! coriolis
                 +fu(i,j)*v_utmp2*halfdt                &
! sea surface pressure gradient
                 -0.5_rlg*gdt/dx_u(i,j)*alphap*(xev(i+1,j)-xev(i,j))      &

! NL advection : u*v*tan(phi)/R (spherical coordinates)
                 +halfdt*u(i,j)*v_utmp2*tfiusr(i,j)     &
! non-linear (udu/dx) : explicit
                 -halfdt*udusdx                         &
! non-linear (vdu/dy) : explicit
                 -halfdt*vdusdy                         &

! viscosity (along x)
                 +halfdt*( (vis2d(i+1,j)*(u(i+1,j)-u(i  ,j)))/dx(i+1,j)   &
                          -(vis2d(i  ,j)*(u(i  ,j)-u(i-1,j)))/dx(i  ,j) ) &
                           /dx_u(i,j)                                     &
! viscosity (along y)
                 +halfdt*( vis2d(i,j  )*mn(i,j)*(u(i,j+1)-u(i,j  ))/dy_f(i,j)     &
                          -vis2d(i,j-1)*ms(i,j)*(u(i,j  )-u(i,j-1))/dy_f(i,j-1) ) &
                           /dy_u(i,j)



! wind stress and atmospherical pressure gradient
       y2dv(i,j)=y2dv(i,j)                              &

! wind stress in x-direction : explicit
                 +halfdt*tosx(i,j)/max(htx(i,j),hminfrot)/rhoref &

! atmospherical pressure gradient : explicit
                 -halfdt*(pa(i+1,j)-pa(i,j))/(dx_u(i,j)*rhoref)

! generator potential : explicit
       IF (l_genpot) &
         y2dv(i,j)=y2dv(i,j)+halfdt/dx_u(i,j)*(potgene(i+1,j)-potgene(i,j))

! internal pressure gradient : explicit
       IF (nbfl /= 0.and.icon==2) &
         y2dv(i,j)=y2dv(i,j)+halfdt/dx_u(i,j)*htx(i,j)*(buoy2d(i+1,j)-buoy2d(i,j))



     END IF   ! ends test on mask for wetting and drying

!&E   ***************   end include dyn2dx_mtm   ******************************
!&E   ***************   include dyn2dx_cty   **********************************
!&E
!&E ** Purpose : remplissage des termes du systeme lineaire
!&E              provenant de la discretisation de l equation de
!&E              continuite dans le systeme d equation de Saint-
!&E              Venant sur l ensemble de la grille de calcul.
!&E
!&E ** Description : Equation couplee a celle de la conservation de la
!&E                  quantite de mouvement dans la direction Ouest-Est.
!&E
!&E
!&E ** Called by : dyn2dx
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xev,u,v,vstar,dx_v,dy_u,hey,hex,surf,phieau
!&E
!&E ** Modified variables : y2dxe,axe,bxe,cxe
!&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        !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E        !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E        !  2009-07  (N. Bruneau, F. Dumas) coupling with swan
!&E
!&E
!&E-------------------------------------------------------------------------

! 1- terme de diagonale de la matrice tridiagonale
!    du modele barotrope
!--------------------------------------------------

! 2- terme de sous-diagonale de la matrice tridiagonale
!    issu de l implicitation de la divergence de flux zonale
!    (flux a l ouest du point xe(i,j) considere). Remplissage
!    du tableau hex qui ressert dans les flux de traceurs
!------------------------------------------------------------

   axe(i,j)=-alphan*halfdt*dy_u(i-1,j)/surf(i,j)*hex(i-1,j)
   bxe(i,j)=fwet(i,j)

! 3- terme de sous-diagonale de la matrice tridiagonale
!    issu de l implicitation de la divergence de flux zonale
!    (flux a l est du point xe(i,j) considere). Remplissage
!    du tableau hex qui ressert dans les flux de traceurs
!-----------------------------------------------------------

   cxe(i,j)=alphan*halfdt*dy_u(i,j)/surf(i,j)*hex(i,j)

! 4- terme du vecteur second membre du systeme lineaire
!    issu de la divergence de flux meridienne
!    (flux au nord et sud  du point xe(i,j) considere)
!    Prise en compte de tous les apports d eau dans
!    la colonne d eau. Remplissage
!    du tableau hey qui ressert dans les flux de traceurs
!----------------------------------------------------------

   phieautot=0.0_rsh
   DO k=1,kmax
     phieautot=phieautot+phieau(k,i,j)
   END DO
   y2dxe(i,j)=xev(i,j)*fwetp(i,j)-halfdt/surf(i,j)*(                   &
                       alphap*( hex(i  ,j)*u(i  ,j)*dy_u(i  ,j)       &
                              - hex(i-1,j)*u(i-1,j)*dy_u(i-1,j) )     &
                      +alphap*( hey(i,j  )*v(i,j  )*dx_v(i,j  )       &
                              - hey(i,j-1)*v(i,j-1)*dx_v(i,j-1) )     &
                      +alphan*( hey(i,j  )*vstar(i,j  )*dx_v(i,j  )   &
                              - hey(i,j-1)*vstar(i,j-1)*dx_v(i,j-1) ) &
                      - phieautot )



!&E   ***************   end include dyn2dx_cty   ******************************

     END DO
   END DO
!$OMP END DO

   !PRINT*, 'exit DYN2DX_TLS_COEF routine'

  END SUBROUTINE dyn2dx_tls_coef
