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

  SUBROUTINE dyn2dy_tls_obc(j1,j2,i,tfb,tfh,v,xe,xeuv)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2dy_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(i,j1), xe(i,j1+1), xe(i,j2), xe(i,j2-1)
!&E                     v(i,j1), v(i,j1+1), v(i,j2-1), v(i,j2-2)
!&E
!&E ** Reference : Blayo & Debreu, 2005. Ocean modelling, 9, 231-252
!&E                Vandermeirsch, convention 03/2.210 458F, sept 2004
!&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-05  (F. Vandermeirsch, J. Marin) Original code
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E       !  2008-08  (T. Odaka) MPI
!&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_vn, obc_vs, obc_sshn, obc_sshs, &
                          l_obc_south,l_obc_north,            &
                          l_obc_char
   IMPLICIT NONE

!! * Arguments
   INTEGER,                                            INTENT( in )     :: j1,j2,i
   REAL(KIND=rsh),                                     INTENT( in )     :: tfh,tfb
   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 ) :: v
!! * 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_south) THEN
     IF ( j1==jmin+1+jshift_obc_south .AND. ljmin<=j1 .AND. j1<=ljmax) THEN
      !IF (Agrif_Root()) THEN
       y2dv(i,j1  )=y2dv(i,j1  )+alphan*grav*halfdt/dy_v(i,j1)*obc_sshs(i)*tfb
      !ELSE
       ! y2dv(i,j1  )=y2dv(i,j1  )+alphan*grav*halfdt/dy_v(i,j1)*xeuv(i,j1)*tfb
      !ENDIF
     ENDIF
     END IF

     IF (l_obc_north) THEN
     IF ( j2==jmax-jshift_obc_north .AND. ljmin<=j2 .AND. j2<=ljmax) THEN
      !IF (Agrif_Root()) THEN
       y2dv(i,j2-1)=y2dv(i,j2-1)-alphan*grav*halfdt/dy_v(i,j2-1)*obc_sshn(i)*tfh
      !ELSE
       ! y2dv(i,j2-1)=y2dv(i,j2-1)-alphan*grav*halfdt/dy_v(i,j2-1)*xeuv(i,j2)*tfh
      !ENDIF
     ENDIF
     END IF


   ELSE    ! Characteristics method

! southern boundary impact
     IF (l_obc_south .AND. j1==jmin+1+jshift_obc_south) THEN
     IF (ljmin<=j1 .AND. j1<=ljmax) THEN
       IF (tfb.eq.1) THEN
         h = h0(i,j1)
!  Methode des caracteristiques + extrapolation
! on est en j1 -1/2 au point xe
         w3=0.5_rsh*(obc_vs(i,j1-1)+obc_vs(i,j1))+SQRT(grav/h)*obc_sshs(i)   
         bv(i,j1)=1.0_rlg/(1.0_rlg/bv(i,j1)+ 0.5_rlg*alphan*al1/ac0*SQRT(grav*h)*halfdt/dy_v(i,j1))
         cv(i,j1)=cv(i,j1)-0.25_rlg*alphan*al1/ac0*grav*halfdt/dy_v(i,j1)
         y2dv(i,j1)= y2dv(i,j1)+0.5_rlg*alphan/ac0*SQRT(grav*h)*halfdt/dy_v(i,j1)*(w3   &
                                    -ac1*v(i,j1)-0.5_rlg*al3*v(i,j1+1)                  &
                                    +2.0_rlg*SQRT(grav/h)*(0.25_rlg*al2*xe(i,j1)+ac2*xe(i,j1+1)))       
       END IF 
     ENDIF
     END IF

! northern boundary impact
     IF (l_obc_north .AND. j2==jmax-jshift_obc_north) THEN
     IF (ljmin<=j2 .AND. j2<=ljmax) THEN
       IF (tfh.eq.1) THEN
         h = h0(i,j2)
         w1=0.5_rsh*(obc_vn(i,j2-1)+obc_vn(i,j2))-SQRT(grav/h)*obc_sshn(i)      

!  Methode des caractéristiques + extrapolation
         av(i,j2-1)=av(i,j2-1)+0.25_rlg*alphan*al1/ac0*grav*halfdt/dy_v(i,j2-1)
         bv(i,j2-1)=1.0_rlg/(1.0_rlg/bv(i,j2-1)+0.5_rlg*alphan*al1/ac0*SQRT(grav*h)*halfdt/dy_v(i,j2-1))
         y2dv(i,j2-1) = y2dv(i,j2-1)-0.5_rlg*alphan/ac0*SQRT(grav*h)*halfdt/dy_v(i,j2-1)*(-w1    &
                                    +ac1*v(i,j2-1)+0.5_rlg*al3*v(i,j2-2)                         &
                                    +2.0_rlg*SQRT(grav/h)*(0.25_rlg*al2*xe(i,j2)+ac2*xe(i,j2-1)))

       END IF
     ENDIF
     END IF

   END IF

  END SUBROUTINE dyn2dy_tls_obc 
