 MODULE obcapplyu
   !!=======================================================================
   !!                        ***  MODULE obcapplyu  ***
   !! Open boundary treatment : apply open boundary schemes
   !!                           for 2d velocities
   !!                           at eastern and western boundaries.
   !!=======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC obc_apply_u,obc_apply_u_mc    ! called by step2dxy, stepu,
                                                ! dyn2dx, dyn3dx

   !! * Shared module variables

   !! * Private variables

 CONTAINS
   !!======================================================================

  SUBROUTINE obc_apply_u(u)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_apply_u  ***
   !&E
   !&E ** Purpose : Apply open boundary conditions for u
   !&E
   !&E ** Description : Available schemes :
   !&E                     - periodic along x or y
   !&E                     - Dirichlet
   !&E                     - MARS standard (CMS)
   !&E
   !&E ** Called by : stepu, step2dxy, dyn2dx, dyn3dx
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !           Original code
   !&E       !  2007-03  (V. Garnier) introduction of obc periodic along y
   !&E       !  2008-07  (V. Garnier T. Odaka) MPI
   !&E       !  2009-03  (V. Garnier) MPI2D
   !&E       !  2010-04 (T. Odaka) MPI2D cyclic
   !&E       !  2010-04 (F. Vandermeirsch) review of MC and introduction of Flather
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE obccombine, ONLY : obc_ue,obc_uw,obc_us,obc_un,     &
                          l_obc_mars,l_obc_diri,l_obc_cycl
#ifdef key_wave_qeuler_dynamics
   USE wave,       ONLY : usz_zint
#endif

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(inout) :: u

   !! * local declarations
   INTEGER                      :: i,j,ii,jj
#ifdef key_lagoon
   INTEGER              :: ng,i0,j0
#endif

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

   PRINT_DBG*, 'enter OBC_APPLY_U routine'

   !! cyclic open boundary conditions
OMPNOAGRIF SINGLE
   IF (l_obc_cycl_x) THEN

     !IF (jmax.EQ.ljmax) u(limin:limax,jmax  )=u(limin:limax,jmax-1)
     !IF (jmin.EQ.ljmin) u(limin:limax,jmin+1)=u(limin:limax,jmin+2)
     !IF (jmin.EQ.ljmin) u(limin:limax,jmin  )=u(limin:limax,jmin+2)
     ! obc mars are applied
     DO i=liminp1,limaxm1
   ! Southern open boundary
       j=jb(i)
       IF_MPI (j==jmin+1+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i,j+1)
         IF (hx(i,j-1) /= -valmanq-h0fond) u(i,j-1)=u(i,j+1)
       ENDIF_MPI
   ! Northern open boundary
       j=jh(i)
       IF_MPI (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i,j-1)
       ENDIF_MPI
     END DO

#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_x_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
       u(ig(j)-1,j)=u(id(j)-2,j)
       u(ig(j)  ,j)=u(id(j)-1,j)
       u(id(j)  ,j)=u(ig(j)+1,j)
#ifdef key_shiftobc
       u(ig(j)-2,j)=u(id(j)-3,j)
       u(id(j)+1,j)=u(ig(j)+2,j)
       u(id(j)+2,j)=u(ig(j)+3,j)
#endif
     END DO
#endif
   ELSE IF (l_obc_cycl_y) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_y_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO i=limin,limax
       u(i,jh(i))=u(i,jb(i)+1)
       u(i,jb(i))=u(i,jh(i)-1)
       u(i,jb(i)-1)=u(i,jh(i)-2)
#ifdef key_shiftobc
       u(i,jb(i)-2)=u(i,jh(i)-3)
       u(i,jh(i)+1)=u(i,jb(i)+2)
       u(i,jh(i)+2)=u(i,jb(i)+3)
#endif
     END DO
#endif
     !IF (imin.EQ.limin) u(imin,ljmin:ljmax)=u(imin+1,ljmin:ljmax)
     !IF (imax.EQ.limax) u(imax,ljmin:ljmax)=u(imax-1,ljmin:ljmax)
     ! l_obc_mars are applied
     DO j=ljminp2,ljmaxm1
   ! Western boundary
       i=ig(j)-1
       IF_MPI (i==imin+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i+1,j)
       ENDIF_MPI
   ! Eastern boundary
       i=id(j)
       IF_MPI (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i-1,j)
       ENDIF_MPI
     END DO
   END IF   ! ends periodic part
OMPNOAGRIF END SINGLE

   !! Dirichlet
   IF (l_obc_diri) THEN
! for straight open boundary only

OMPNOAGRIF DO SCHEDULE(static)
     DO j=ljmin,ljmax

       i=ig(j)-1
       IF_MPI (i==imin+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
#ifdef key_shiftobc_doTSobc
         DO ii=MAX0(imin,imin+ishift_obc_west-1),ig(j)-1
           IF (hx(ii,j) /= -valmanq-h0fond) u(ii,j)=obc_uw(ii,j)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           u(i,j)=obc_uw(i,j)
#ifdef key_wave_qeuler_dynamics
           u(i,j)=u(i,j)-usz_zint(i+1,j)
#endif
         END IF
#endif
       ENDIF_MPI

       i=id(j)
       IF_MPI (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN
#ifdef key_shiftobc_doTSobc
         DO ii=id(j),MIN0(imax,imax-ishift_obc_east+2)
           IF (hx(ii,j) /= -valmanq-h0fond) u(ii,j)=obc_ue(ii,j)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           u(i,j)=obc_ue(i,j)
#ifdef key_wave_qeuler_dynamics
           u(i,j)=u(i,j)-usz_zint(i-1,j)
#endif
         END IF
#endif
       ENDIF_MPI
     END DO
OMPNOAGRIF END DO

OMPNOAGRIF SINGLE
     DO i=liminp1,limaxm1
       j=jb(i)
       IF_MPI (j==jmin+1+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN
#ifdef key_shiftobc_doTSobc
         DO jj=MAX0(jmin,jmin+jshift_obc_south-1),jb(i)
           IF (hx(i,jj) /= -valmanq-h0fond) u(i,jj)=obc_us(i,jj)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           u(i,j)=obc_us(i,j)
#ifdef key_wave_qeuler_dynamics
           u(i,j)=u(i,j)-usz_zint(i,j+1)
#endif
         END IF
#endif
       ENDIF_MPI

       j=jh(i)
       IF_MPI (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN
#ifdef key_shiftobc_doTSobc
         DO jj=jh(i),MIN0(jmax,jmax-jshift_obc_north+2)
           IF (hx(i,jj) /= -valmanq-h0fond) u(i,jj)=obc_un(i,jj)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           u(i,j)=obc_un(i,j)
#ifdef key_wave_qeuler_dynamics
           u(i,j)=u(i,j)-usz_zint(i,j-1)
#endif
         END IF
#endif
       ENDIF_MPI
     END DO
OMPNOAGRIF END SINGLE

   !! Ends Dirichlet
   END IF

#ifdef key_lagoon
   !! open boundary conditions (LAGOON only) l_obc_mars=.true.

  ! OBC threw graus

     DO ng=1,nbgrau
       i0=igrau(ng)
       j0=jgrau(ng)
       
   ! Western boundary
       IF (type_grau(ng)==4) THEN
       !! grau sur la frontiere Ouest
         i=i0-1 
         IF_MPI ( (i >= liminm2) .AND. (i <= limaxp2) .AND.(j0 >= ljminm2) .AND. (j0 <= ljmaxp2) ) THEN
            IF (hx(i,j0) /= -valmanq-h0fond) u(i,j0)=u(i0,j0)
         ENDIF_MPI 
       END IF

   ! Eastern boundary
       IF (type_grau(ng)==2) THEN
       !! grau sur la frontiere Est
         IF_MPI ( (i0 >= liminm2) .AND. (i0 <= limaxp2) .AND.(j0 >= ljminm2) .AND. (j0 <= ljmaxp2) ) THEN
            IF (hx(i0,j0) /= -valmanq-h0fond) u(i0,j0)=u(i0-1,j0)
         ENDIF_MPI 
       END IF
       
     END DO

#else

   !! Default open boundary conditions (MARS standard)

   IF (l_obc_mars) THEN

OMPNOAGRIF DO SCHEDULE(static)
     DO j=ljminp2,ljmaxm1

   ! Western boundary

       i=ig(j)-1
       IF_MPI (i==imin+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i+1,j)
       ENDIF_MPI

   ! Eastern boundary

       i=id(j)
       IF_MPI (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i-1,j)
       ENDIF_MPI

     END DO
OMPNOAGRIF END DO


OMPNOAGRIF SINGLE
     DO i=liminp1,limaxm1

   ! Southern open boundary

       j=jb(i)
       IF_MPI (j==jmin+1+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i,j+1)
       ENDIF_MPI

   ! Northern open boundary

       j=jh(i)
       IF_MPI (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hx(i,j) /= -valmanq-h0fond) u(i,j)=u(i,j-1)
       ENDIF_MPI

     END DO
OMPNOAGRIF END SINGLE

   !! Ends MARS standard open boundary conditions  (straight boundaries)
   END IF
#endif

OMPMPINOAGRIF barrier
OMPMPINOAGRIF master
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
!part really usefull call ex_j_rsh(-2,2,1,liminm1,limax,ljminm2,ljmaxp2,u(liminm1:limax,ljminm2:ljmaxp2))
!larger exchanges to make transfer in step2d* and step3d* easier
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPINOAGRIF end master
OMPMPINOAGRIF barrier
OMPMPINOAGRIF flush(u)

   PRINT_DBG*, 'exit OBC_APPLY_U routine'

  END SUBROUTINE obc_apply_u

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

  SUBROUTINE obc_apply_u_mc(ua,u,va,v,xea,xe,l_star)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_apply_u_mc  ***
   !&E
   !&E ** Purpose : Apply open boundary conditions for u
   !&E
   !&E ** Description : - Characteristics method
   !&E
   !&E ** Called by : dyn2dxy, dyn2dyx, dyn3dxy, dyn3dyx
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :  ua,u,va,v,xea,xe
   !&E
   !&E ** Modified variables : u,xe
   !&E
   !&E ** Reference : Blayo & Debreu, 2005. Ocean modelling, 9, 231-252
   !&E                Vandermeirsch, convention 03/2.210 458F, sept 2004
   !&E
   !&E ** History :
   !&E       !  2007-05  (F. Vandermeirsch, J. Marin, V. Garnier) Original code
   !&E       !  2009-07  (F. Vandermeirsch)  MC 2nd order scheme
   !&E       !  2009-07  (V. Garnier) MPI2D
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2011-11  (F. Vandermeirsch) evolution
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE obccombine, ONLY : obc_ue,obc_uw,obc_us,obc_un,obc_vs,obc_vn
   USE obccombine, ONLY : obc_sshs,obc_sshn,obc_sshe,obc_sshw, &
                          l_obc_south,l_obc_north,l_obc_west,l_obc_east

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( inout ) :: u
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( inout ) :: xe
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: xea
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: ua,v,va
   LOGICAL,                                        INTENT( in )     :: l_star

   !! * local declarations
   INTEGER                      :: i,j,i1,i2,j1,j2
   REAL(KIND=rsh)               :: h,w3,w1u,w1ua,w1xea,w1,w3u,w3ua,w3xea, &
                                   w1v,w1va,w3v,w3va

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

   PRINT_DBG*, 'enter OBC_APPLY_U_MC routine'

   ! eastern and western boundaries
   ! ------------------------------

OMPNOAGRIF DO SCHEDULE(RUNTIME) 
!FVDM
! cohérence des "coins" en xe 
   !DO j=ljminp1,ljmaxm1 
   DO j=ljminp2,ljmaxm1 

     i1 = ig(j)
     i2 = id(j)

   ! western boundary

     IF (l_obc_west) THEN
       IF_MPI (i1==imin+1+ishift_obc_west .AND. limin<=i1 .AND. i1<=limax) THEN
!         IF (htx(i1,j)>0.0_rsh) THEN
         IF (h0(i1,j)>0.0_rsh) THEN
           h = h0(i1,j)
           w3=0.5_rsh*(obc_uw(i1-1,j)+obc_uw(i1,j))+SQRT(grav/h)*obc_sshw(j)
           w1ua=ua(i1,j)-0.5_rsh*SQRT(grav/h)*(xea(i1,j)+xea(i1+1,j))
           w1xea =0.5_rsh*(ua(i1,j)+ua(i1+1,j))-SQRT(grav/h)*xea(i1+1,j)

           IF (l_star) THEN 
#ifdef key_char_flather	     
   ! Flather [w1 = (w1)ext et u0 = u1]
             xe(i1,j)=SQRT(h/grav)*(w3-u(i1,j))
             u(i1-1,j)= u(i1,j)
#else
   ! Characteristics method
             xe(i1,j)=(0.5_rsh*SQRT(h/grav)*(w3-al1*u(i1,j)-al2*w1ua-al3*w1xea) &
                                             +0.25_rsh*al1*xe(i1+1,j))/ac0
             w1u=u(i1,j)-0.5_rsh*SQRT(grav/h)*(xe(i1,j)+xe(i1+1,j))
             u(i1-1,j)= w3+al1*w1u+al2*w1ua+al3*w1xea-u(i1,j)
#endif
           ELSE
             u(i1-1,j)= obc_uw(i1,j) 
             !xe(i1,j)=xea(i1,j) 
           ENDIF
         END IF
       ENDIF_MPI

   END IF 

   ! eastern boundary

     IF (l_obc_east) THEN
       IF_MPI (i2==imax-ishift_obc_east .AND. limin<=i2 .AND. i2<=limax) THEN
!         IF (htx(i2-1,j)>0.0_rsh) THEN
         IF (h0(i2,j)>0.0_rsh) THEN
           h = h0(i2,j)
         ! w1 caracteristique entrante 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)
           w3ua=ua(i2-1,j)+0.5_rlg*SQRT(grav/h)*(xea(i2,j)+xea(i2-1,j))
           w3xea=0.5_rlg*(ua(i2-1,j)+ua(i2-2,j))+SQRT(grav/h)*xea(i2-1,j)

         ! xe_N+0.5 = 0.5*sqrt(g/h0)*(w3_N+0.5-w1_N+0.5)
           IF (l_star) THEN 
             xe(i2,j)=(0.5_rsh*SQRT(h/grav)*(-w1+al1*u(i2-1,j)+al2*w3ua+al3*w3xea) &
                                             +0.25_rsh*al1*xe(i2-1,j))/ac0
           ELSE
             !xe(i2,j)=xea(i2,j) 
             u(i2,j)=obc_ue(i2,j)
           ENDIF
         ! 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
           w3u=u(i2-1,j)+0.5_rsh*SQRT(grav/h)*(xe(i2,j)+xe(i2-1,j))
         ! u_N+1 = w1_N+0.5 + w3_N+0.5 - u_N
           u(i2,j)= w1+al1*w3u+al2*w3ua+al3*w3xea-u(i2-1,j)
         END IF
       ENDIF_MPI
     END IF

   END DO
OMPNOAGRIF END DO


OMPNOAGRIF DO SCHEDULE(RUNTIME) 
!FVDM
! cohérence des "coins" en vitesse 
!DO i=liminp1,limaxm1
   DO i=limin,limax

     j1 = jb(i)
     j2 = jh(i)

   ! southern boundary

     IF (l_obc_south) THEN
       IF_MPI (j1==jmin+1+jshift_obc_south .AND. ljmin<=j1 .AND. j1<=ljmax) THEN
!         IF (hty(i,j1)>0.0_rsh) THEN
         IF (h0(i,j1)>0.0_rsh) THEN
           h = h0(i,j1)

           ! no update of xe at the southern boundary because not strong obc for calculation in row

           IF (hx(i,j1) /= -valmanq-h0fond) THEN
             IF ((v(i,j1-1)+v(i,j1)).GE.0.0_rsh) THEN
               u(i,j1)= obc_us(i,j1)
             ELSE
               u(i,j1)=u(i,j1+1)
             END IF
           END IF
         END IF
       ENDIF_MPI
     END IF

   ! northern boundary

     IF (l_obc_north) THEN
       IF_MPI (j2==jmax-jshift_obc_north .AND. ljmin<=j2 .AND. j2<=ljmax) THEN
!         IF (hty(i,j2-1)>0.0_rsh) THEN
         IF (h0(i,j2)>0.0_rsh) THEN
           h = h0(i,j2)

           ! no update of xe at the northern boundary because not strong obc for calculation in row
           IF (hx(i,j2) /= -valmanq-h0fond) THEN
             IF ((v(i,j2-1)+v(i,j2)).LE.0.0_rsh) THEN
               u(i,j2)= obc_un(i,j2)
             ELSE
               u(i,j2)=u(i,j2-1)
             END IF
           END IF
         END IF
       ENDIF_MPI
     END IF
   END DO
OMPNOAGRIF END DO

OMPMPINOAGRIF barrier
OMPMPINOAGRIF master
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
!part really usefull call
!ex_j_rsh(-2,2,1,liminm1,limax,ljminm2,ljmaxp2,u(liminm1:limax,ljminm2:ljmaxp2))
!larger exchanges to make transfer in step2d* and step3d* easier
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xe(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xe(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPINOAGRIF end master
OMPMPINOAGRIF barrier
OMPMPINOAGRIF flush(u,xe)

   PRINT_DBG*, 'exit OBC_APPLY_U_MC routine'

  END SUBROUTINE obc_apply_u_mc

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

 END MODULE obcapplyu
