!!======================================================================

  SUBROUTINE dyn2dyctymtm2(v,xev,xeuv,u,uv,vstar,ustar)

!&E---------------------------------------------------------------------
!&E                 ***  dyn2dyctymtm2  ***
!&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.
!&E
!&E ** Description :
!&E
!&E ** Called by : dyn2dxy.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xev,xeuv,u,v,uv,ustar,vstar,igdv,hy,htx,hty,
!&E                     dx_f,dx_v,dy,dy_u,dy_v,fv,hey,hex,surf,phieau,
!&E                     tfivsr,vis2d,tosy,pa,potgene,buoy2d
!&E                     me,mo,why
!&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) Commentaires et uniformisation dx,dy
!&E        !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E        !  2009-01  (A. Mano T. Odaka) MPI2D
!&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  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E        !  2011-11  (F. Vandermeirsch) correct spatial extension on which we can estimate quick and quickest scheme
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,  ONLY : jb,jh,surf,dx_u,dx_v,dx_f,dy,dy_u,dy_v,icon,halfdt, &
                          hy,hex,hey,htx,hty,fv,tfivsr,vis2d,dt,h0,frot_lin,  &
                          hminfrot,botstressmax,hmaxfrot,z0b,fwet,fwetp
   USE comvarp2d,  ONLY : u_vtmp
   USE comvars3d,  ONLY : rhoref
   USE comcolig2d, ONLY : y2dv,cv,bv,av,igdv,axe,bxe,cxe,y2dxe,me,mo,why
   USE comaree,    ONLY : l_genpot,potgene
   USE cometdyn,   ONLY : pa,tosy
   USE comfluxdo,  ONLY : phieau,buoy2d,nbfl



   IMPLICIT NONE

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

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

   REAL(KIND=rsh)                               :: vc,vm,vp,um,up


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

   gdt = grav*dt


!$OMP DO schedule(runtime)

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




!$OMP DO schedule(runtime)


   DO j=ljminp1,ljmax
     DO i=liminp2,limaxm1
       j1=jb(i)
       j2=jh(i)
       IF (j >= j1+1 .AND. j <= j2-1) THEN

! Evaluation of axe, bxe, cxe and y2dxe
!--------------------------------------

! 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)=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) )  &
                             +alphan*( hex(i  ,j)*uv(i  ,j)*dy_u(i  ,j)   &
                                      -hex(i-1,j)*uv(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) )  &
                            -phieautot )



       END IF

! Evaluation of av, bv, cv and y2dv
!----------------------------------

       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_vtmp(i,j)*halfdt                 &
! sea surface pressure gradient
                  -0.5_rlg*gdt/dy_v(i,j)*alphap*( xev(i,j+1)-xev(i,j) )  &

! NL advection : u*u*tan(phi)/R (spherical coordinates)
                  -halfdt*tfivsr(i,j)*u_vtmp(i,j)**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

     END DO    ! ends loop on i
   END DO     ! ends loop on j
!$OMP END DO

  END SUBROUTINE dyn2dyctymtm2
