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

  SUBROUTINE dyn2dx_tls_obc(i1,i2,j,tfg,tfd,u,xe,xeuv)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2dx_tls_obc  ***
!&E
!&E ** Purpose : Modify coefficients of matrix
!&E              to introduce implicit part of the sea surface pressure gradient
!&E
!&E ** Description : - Characteristics method (ssh and u)
!&E                  - or Pure prescription of sea surface height
!&E
!&E ** Called by : dyn2dxy, dyn2dyx, dyn3dxy, dyn3dyx
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xe(i1,j), xe(i1+1,j), xe(i2,j), xe(i2-1,j)
!&E                     u(i1,j), u(i1+1,j), u(i2-1,j), u(i2-2,j)
!&E
!&E ** Modified variables : y2dv
!&E
!&E ** Reference : Blayo & Debreu, 2005. Ocean modelling, 9, 231-252
!&E                Vandermeirsch, convention 03/2.210 458F, sept 2004
!&E
!&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
!&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
!&E
!&E
!&E ** History :
!&E       !  2007-05  (F. Vandermeirsch, J. Marin) Original code
!&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2009-01  (A. Mano T. Odaka) MPI2D
!&E       !  2009-01  (F. Vandermeirsch) introduction of obc_ssh_* and characteristics method
!&E       !  2009-07  (F. Vandermeirsch)  MC 2nd order scheme
!&E       !  2009-07  (V. Garnier) MPI2d for characteristics method
!&E       !  2010-04  (F. Vandermeirsch) OBC : introduction of Flather operator
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2011-11  (F. Vandermeirsch) correction : spatial average
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d
   USE comcolig2d, ONLY : y2dv,av,bv,cv
   USE comfluxdo
   USE obccombine, ONLY : obc_ue, obc_uw, obc_sshe, obc_sshw, &
                          l_obc_west,l_obc_east,              &
                          l_obc_char
   USE toolmpi,     only : ex_i_rsh

   IMPLICIT NONE

!! * Arguments
   INTEGER,                                            INTENT( in )     :: i1,i2,j
   REAL(KIND=rsh),                                     INTENT( in )     :: tfg,tfd
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: xe
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: xeuv
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: u

!! * local declarations
   REAL(KIND=rsh)               :: h,w1,w3

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

   IF (.NOT. l_obc_char) THEN

! introduce sea suface height at open boundaries

     IF (l_obc_west) THEN
!IF ( i1==imin+31+ishift_obc_west .AND. limin<=i1 .AND. i1<=limax) THEN
!IF ( i1==imin+31+ishift_obc_west .AND. limin==imin) THEN
     IF ( i1==imin+31 .AND. limin==imin) THEN
      !IF (Agrif_Root()) THEN
       y2dv(i1  ,j)=y2dv(i1  ,j)+alphan*grav*halfdt/dx_u(i1,j)*obc_sshw(j)*tfg
      !ELSE
       ! y2dv(i1  ,j)=y2dv(i1  ,j)+alphan*grav*halfdt/dx_u(i1,j)*xeuv(i1,j)*tfg
      !ENDIF
     ENDIF
     END IF

     IF (l_obc_east) THEN
     IF ( i2==imax-ishift_obc_east .AND. limin<=i2 .AND. i2<=limax) THEN
      !IF (Agrif_Root()) THEN
       y2dv(i2-1,j)=y2dv(i2-1,j)-alphan*grav*halfdt/dx_u(i2-1,j)*obc_sshe(j)*tfd
      !ELSE
       ! y2dv(i2-1,j)=y2dv(i2-1,j)-alphan*grav*halfdt/dx_u(i2-1,j)*xeuv(i2,j)*tfd
      !ENDIF
     ENDIF
     END IF

   ELSE    ! Characteristics method

! western boundary impact

     IF (l_obc_west .AND. i1==imin+1+ishift_obc_west) THEN
     IF (limin<=i1 .AND. i1<=limax) THEN
       IF (tfg.eq.1) THEN
         h = h0(i1,j)
         w3=0.5_rsh*(obc_uw(i1-1,j)+obc_uw(i1,j))+SQRT(grav/h)*obc_sshw(j)


! Methode des caracteristiques
         bv(i1,j)=1.0_rlg/(1.0_rlg/bv(i1,j)+ 0.5_rlg*alphan*al1/ac0*SQRT(grav*h)*halfdt/dx_u(i1,j))
         cv(i1,j)=cv(i1,j)-0.25_rlg*alphan*al1/ac0*grav*halfdt/dx_u(i1,j)
         y2dv(i1,j)=y2dv(i1,j)+0.5_rlg*alphan/ac0*SQRT(grav*h)*halfdt/dx_u(i1,j)*(w3    &
                                    -ac1*u(i1,j)-0.5_rlg*al3*u(i1+1,j)                  &
                                    +2.0_rlg*SQRT(grav/h)*(0.25_rlg*al2*xe(i1,j)+ac2*xe(i1+1,j)))

       END IF
     ENDIF
     END IF

! eastern boundary impact
     IF (l_obc_east .AND. i2==imax-ishift_obc_east) THEN
     IF (limin<=i2 .AND. i2<=limax) THEN
       IF (tfd.eq.1) THEN
         h = h0(i2,j)
! Methode des caracteristiques
! w1 caracteristique entrante, calculee avec variable OGCM a position xe
! w1_N+0.5 = u_ext - sqrt(g/h0)*h_ext
! on est en i2 -1/2 au point xe
         w1=0.5_rsh*(obc_ue(i2-1,j)+obc_ue(i2,j))-SQRT(grav/h)*obc_sshe(j)
! w3 caracteristique sortante, calculee avec variables diagnostiques
! schema d ordre 2 : (al1,al2,al3) = (1,1,-1)
! (w3_N+0.5)n+.5 = al1*(w3_N)n+.5 + al2*(w3_N)n + al3*(w3_N-.5)n
! issu du developpement xe_N+0.5 = 0.5*sqrt(g/h0)*(w3_N-w1_N+0.5)
         av(i2-1,j)=av(i2-1,j)+0.25_rlg*alphan*al1/ac0*grav*halfdt/dx_u(i2-1,j)
         bv(i2-1,j)=1.0_rlg/(1.0_rlg/bv(i2-1,j)+ 0.5_rlg*alphan*al1/ac0*SQRT(grav*h)*halfdt/dx_u(i2-1,j))
         y2dv(i2-1,j) = y2dv(i2-1,j)-0.5_rlg*alphan/ac0*SQRT(grav*h)*halfdt/dx_u(i2-1,j)*(-w1   &
                                    +ac1*u(i2-1,j)+0.5_rlg*al3*u(i2-2,j)                         &
                                    +2.0_rlg*SQRT(grav/h)*(0.25_rlg*al2*xe(i2,j)+ac2*xe(i2-1,j)))
       END IF
     ENDIF
     END IF

   END IF

   !PRINT*, 'exit DYN2DX_TLS_OBC routine'

  END SUBROUTINE dyn2dx_tls_obc
