  SUBROUTINE dyn2dy_tls_coef(v,u,xeu,vstar,ustar)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2dy_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 y-direction
!&E
!&E ** Called by : dyn2dyx.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xeu,u,v,ustar,vstar,igdv,hy,htx,hty,dy,dx_f,dx_v,
!&E                     dy_v,dy_u,fv,hey,hex,surf,phieau,
!&E                     tfivsr,vis2d,tosy,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        !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E        !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E        !  2008-07  (T. Odaka) MPI
!&E        !  2009-01  (A. Mano) MPI2D
!&E        !  2009-01  (V. Garnier) OMP
!&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,tfivsr,ig,id,jb,jh,          &
                          htx,hty,hy,hey,hex,dx_u,dx_v,dx_f,surf,   &
                          dy,dy_v,dy_u,fv,vis2d,halfdt,frot_lin,    &
                          dt,icon,h0,hminfrot,hmaxfrot,z0b,fwet,fwetp
   USE comvarp2d,  ONLY : u_vtmp
   USE comfluxdo,  ONLY : nbfl,buoy2d,phieau
   USE comaree,    ONLY : l_genpot,potgene
   USE cometdyn,   ONLY : tosy,pa
   USE comvars3d,  ONLY : rhoref
   USE comcolig2d, ONLY : y2dxe,y2dv,axe,bxe,cxe,cv,bv,av,igdv,me,mo,why




   IMPLICIT NONE

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

!! * Local declarations
   INTEGER                                       :: i,j,k,j1,j2
   REAL(KIND=rlg)                                :: phieautot              !dyn2dy_cty.h
   REAL(KIND=rsh)  :: ft,gdt,vit,factft,udvsdx,udvdx,udvdx2,   &           !dyn2dy_mtm.h
                      vdvsdy,vdvdy,vdvdy2,dvsdx,dvsdx2,u_vtmp2

   REAL(KIND=rsh)  :: vc,vp,vm,up,um                                       !dyn2dy_mtm.h


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

   !PRINT*, 'enter DYN2DY_TLS_COEF routine'


!$OMP DO SCHEDULE(RUNTIME)

   DO j=ljmin,ljmax
     DO i =liminp1,limax
       IF (j >= jb(i)-1 .and. j <= jh(i)-1) THEN
         u_vtmp(i,j)=0.25_rsh*(u(i,j)+u(i,j+1)+u(i-1,j)+u(i-1,j+1))
       END IF
     END DO
   END DO
!$OMP END DO




!$OMP DO SCHEDULE(RUNTIME) &

!$OMP& PRIVATE(i,j,ft,gdt,vit,factft,dvsdx,dvsdx2,u_vtmp2,phieautot)
   DO j=ljminp1,ljmax
     DO i =liminp2,limaxm1

!&E   ***************   include dyn2dy_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 colonne. (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 : xeu,u,v,ustar,vstar,igdv,hy,htx,hty,
!&E                     dx_f,dx_v,dy,dy_v,fv,me,mo,why,
!&E                     tfivsr,vis2d,tosy,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 udvsdx 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-------------------------------------------------------------------------

     gdt=grav*dt

     IF (j >= jb(i) .and. j <= jh(i)-1) THEN

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

      IF (igdv(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 (igdv(i,j) == 1) THEN

        u_vtmp2=0.25_rsh*(ustar(i,j)+ustar(i,j+1)+ustar(i-1,j)+ustar(i-1,j+1))

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



! bottom stress



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





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

        vc=vstar(i,j)
        vp=0.5_rsh*(vc+abs(vc))
        vm=0.5_rsh*(vc-abs(vc))
        up=0.5_rsh*(u_vtmp2+abs(u_vtmp2))
        um=0.5_rsh*(u_vtmp2-abs(u_vtmp2))
        IF (i >= imin+3.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+2) THEN
          vdvdy=(vp*v(i,j-2)*why(i,j-1)+(-7.0_rsh*vp-3.0_rsh*vm)*v(i,j-1)+  &
                 (3.0_rsh*vp-3.0_rsh*vm)*v(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*v(i,j+1)-vm*v(i,j+2)*why(i,j+1))/(8.0_rsh*dy_v(i,j)) 
          vdvdy2=(vp*vstar(i,j-2)*why(i,j-1)+(-7.0_rsh*vp-3.0_rsh*vm)*vstar(i,j-1)+  &
                  (3.0_rsh*vp-3.0_rsh*vm)*vstar(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*vstar(i,j+1)-vm*vstar(i,j+2)*why(i,j+1)) &
                 /(8.0_rsh*dy_v(i,j)) 
          vdvsdy = (alphan*vdvdy2 + alphap*vdvdy)
          udvdx=(up*v(i-2,j)*mo(i-1,j)*mo(i,j)+(-7.0_rsh*up-3.0_rsh*um)*v(i-1,j)+  &
                 (3.0_rsh*up-3.0_rsh*um)*v(i,j)+(3.0_rsh*up+7.0_rsh*um)*v(i+1,j)-um*v(i+2,j)*me(i+1,j)*me(i,j))/(8.0_rsh*dx_v(i,j))
          udvdx2=(up*vstar(i-2,j)*mo(i-1,j)*mo(i,j)+(-7.0_rsh*up-3.0_rsh*um)*vstar(i-1,j)+  &
                  (3.0_rsh*up-3.0_rsh*um)*vstar(i,j)+(3.0_rsh*up+7.0_rsh*um)*vstar(i+1,j)-um*vstar(i+2,j)*me(i+1,j)*me(i,j)) &
                 /(8.0_rsh*dx_v(i,j))
          udvsdx = (alphan*udvdx2 + alphap*udvdx)
        ELSE
! upwind scheme
          vdvdy  = vp*(v(i,j)-v(i,j-1))/dy(i,j)+vm*(v(i,j+1)-v(i,j))/dy(i,j+1)
          vdvdy2 = vp*(vstar(i,j)-vstar(i,j-1))/dy(i,j)+vm*(vstar(i,j+1)-vstar(i,j))/dy(i,j+1)
          vdvsdy = (alphan*vdvdy2 + alphap*vdvdy)
          udvdx  = up*(v(i,j)-v(i-1,j))/dx_f(i-1,j)+um*(v(i+1,j)-v(i,j))/dx_f(i,j)
          udvdx2 = up*(vstar(i,j)-vstar(i-1,j))/dx_f(i-1,j)+um*(vstar(i+1,j)-vstar(i,j))/dx_f(i,j)
          udvsdx = (alphan*udvdx2 + alphap*udvdx)
        END IF



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

! bottom stress : implicit part
                +ft*cfroi)

! rhs : explicit terms
        y2dv(i,j)=                                       &
! time evolution (v at n)
                  v(i,j)                                 &
! bottom stress (v at n)
                  -ft*cfroe*v(i,j)                       &
! coriolis
                  -fv(i,j)*u_vtmp2*halfdt                &
! sea surface pressure gradient
                  -0.5_rlg*gdt/dy_v(i,j)*alphap*( xeu(i,j+1)-xeu(i,j) )  &

! NL advection : u*u*tan(phi)/R (spherical coordinates)
                  -halfdt*tfivsr(i,j)*u_vtmp2**2         &
! non-linear (vdv/dy) : explicit
                  -halfdt*vdvsdy                         &
! non-linear (udv/dx) : explicit
                  -halfdt*udvsdx                         &

! viscosity (along y)
                  +halfdt*( vis2d(i,j+1)*(v(i,j+1)-v(i,j  ))/dy(i,j+1)   &
                           -vis2d(i,j  )*(v(i,j  )-v(i,j-1))/dy(i,j  ) ) &
                           /dy_v(i,j)                                    &
! viscosity (along x)
                  +halfdt*( vis2d(i  ,j)*me(i,j)*(v(i+1,j)-v(i  ,j))/dx_f(i  ,j)   &
                           -vis2d(i-1,j)*mo(i,j)*(v(i  ,j)-v(i-1,j))/dx_f(i-1,j) ) &
                           /dx_v(i,j)



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

! wind stress in y-direction : explicit
                  +halfdt*tosy(i,j)/max(hty(i,j),hminfrot)/rhoref &

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

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

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



      END IF  ! ends test on mask for wetting and drying
     END IF    ! fin de test sur l emprise de la  colonne

!&E   ***************   end include dyn2dy_mtm   ******************************

!&E   ***************   include dyn2dy_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 Nord-Sud.
!&E
!&E
!&E ** Called by : dyn2dy_tls_coef
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xeu,u,v,ustar,dx_v,dy_u,hey,hex,surf,phieau
!&E
!&E ** Modified variables : y2dxe,axe,bxe,cxe
!&E
!&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-------------------------------------------------------------------------

       j1=jb(i)
       j2=jh(i)
       IF (j >= j1+1.and.j <= j2-1) THEN

! 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 meridienne
!    (flux au sud  du point xe(i,j) considere). Remplissage
!    du tableau hey qui ressert dans les flux de traceurs
!------------------------------------------------------------

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

! 3- terme de sous-diagonale de la matrice tridiagonale
!    issu de l implicitation de la divergence de flux meridienne
!    (flux au nord du point xe(i,j) considere). Remplissage
!    du tableau hey qui ressert dans les flux de traceurs
!------------------------------------------------------------
!
         cxe(i,j)=alphan*halfdt*dx_v(i,j)/surf(i,j)*hey(i,j)
!
! 4- terme du vecteur second membre du systeme lineaire
!    issu de la divergence de flux zonale
!    (flux a l est et a l ouest du point xe(i,j) considere)
!    Prise en compte de tous les apports d eau dans
!    la colonne d eau. Remplissage
!    du tableau hex 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)=xeu(i,j)*fwetp(i,j)-halfdt/surf(i,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) )     &
                             +alphap*( hex(i  ,j)*u(i  ,j)*dy_u(i  ,j)       &
                                     - hex(i-1,j)*u(i-1,j)*dy_u(i-1,j) )     &
                             +alphan*( hex(i  ,j)*ustar(i  ,j)*dy_u(i  ,j)   &
                                     - hex(i-1,j)*ustar(i-1,j)*dy_u(i-1,j) ) &
                                     - phieautot )



       END IF

!&E   ***************   end include dyn2dy_cty   ******************************

     END DO
   END DO
!$OMP END DO

   !PRINT*, 'exit DYN2DY_TLS_COEF routine'

  END SUBROUTINE dyn2dy_tls_coef
