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





























































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

   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC obc_apply_v,obc_apply_v_mc      ! called by step2dyx, stepv,
! dyn2dy, dyn3dy

!! * Shared module variables

!! * Private variables

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

  SUBROUTINE obc_apply_v(v)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE obc_apply_v  ***
!&E
!&E ** Purpose : Apply open boundary conditions for v
!&E
!&E ** Description : Available schemes :
!&E                     - periodic along x or y
!&E                     - Dirichlet
!&E                     - MARS standard (CMS)
!&E
!&E ** Called by : stepv, step2dyx, dyn2dy, dyn3dy
!&E
!&E ** External calls :
!&E
!&E ** Reference : Blayo Debreu
!&E
!&E ** History :
!&E                    Original code
!&E       !  2007-03  (V. Garnier) introduction of obc periodic along y
!&E       !  2007-05  (F. Vandermeirsch, J. Marin, V. Garnier) introduction of caracteristics method
!&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_ve,obc_vw,obc_vs,obc_vn,     &
                          l_obc_mars,l_obc_diri,l_obc_cycl


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

!! * local declarations
   INTEGER                      :: i,j,ii,jj


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

   !PRINT*, 'enter OBC_APPLY_V routine'

!! cyclic open boundary conditions
!$OMP SINGLE
   IF (l_obc_cycl_x) THEN


     CALL ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))


!IF (jmax.EQ.ljmax) v(limin:limax,jmax  )=v(limin:limax,jmax-1)
!IF (jmin.EQ.ljmin) v(limin:limax,jmin  )=v(limin:limax,jmin+1)
! mars obc are applied
     DO i=liminp2,limaxm1
! Southern boundary
       j=jb(i)-1
       IF (j==jmin+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i,j+1)
       ENDIF
! Northern boundary
       j=jh(i)
       IF (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i,j-1)
       ENDIF
     END DO

   ELSE IF (l_obc_cycl_y) THEN
!IF (imin.EQ.limin)  v(imin  ,:)=v(imin+2,:)
!IF (imin.EQ.limin)  v(imin+1,:)=v(imin+2,:)
!IF (imax.EQ.limax)  v(imax  ,:)=v(imax-1,:)
!l_obc_mars are applied
     DO j=ljminp1,ljmaxm1
! Western boundary
       i=ig(j)
       IF (i==imin+1+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
!cvalprov         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i,j-1)
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i+2,j)
         IF (hy(i-1,j) /= -valmanq-h0fond) v(i-1,j)=v(i+2,j)
       ENDIF
! Eastern boundary
       i=id(j)
       IF (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i-1,j)
       ENDIF
     END DO


     CALL ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))


   END IF   ! ends periodic part
!$OMP END SINGLE


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

!$OMP DO SCHEDULE(static)
     DO j=ljmin,ljmax

       i=ig(j)
       IF (i==imin+1+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN

         IF (hy(i,j) /= -valmanq-h0fond) THEN
           v(i,j)=obc_vw(i,j)

         END IF

       ENDIF

       i=id(j)
       IF (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN

         IF (hy(i,j) /= -valmanq-h0fond) THEN
           v(i,j)=obc_ve(i,j)

         END IF

       ENDIF
     END DO
!$OMP END DO

!$OMP SINGLE
     DO i=liminp2,limaxm1

       j=jb(i)-1
       IF (j==jmin+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN

         IF (hy(i,j) /= -valmanq-h0fond) THEN
           v(i,j)=obc_vs(i,j)

         END IF

       ENDIF

       j=jh(i)
       IF (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN

         IF (hy(i,j) /= -valmanq-h0fond) THEN
           v(i,j)=obc_vn(i,j)

         END IF

       ENDIF

     END DO
!$OMP END SINGLE

!! Ends Dirichlet
   END IF




!! Default open boundary conditions (MARS standard)

   IF (l_obc_mars) THEN

!$OMP DO SCHEDULE(static)
     DO j=ljminp1,ljmaxm1

! Western boundary

       i=ig(j)
       IF (i==imin+1+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i+1,j)
       ENDIF

! Eastern boundary

       i=id(j)
       IF (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i-1,j)
       ENDIF

     END DO
!$OMP END DO


!$OMP SINGLE
     DO i=liminp2,limaxm1

! Southern boundary

       j=jb(i)-1
       IF (j==jmin+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i,j+1)
       ENDIF

! Northern boundary

       j=jh(i)
       IF (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) v(i,j)=v(i,j-1)
       ENDIF

     END DO
!$OMP END SINGLE

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


!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
!part really usefull call ex_j_rsh(-2,2,1,limin,limaxp1,ljminm2,ljmaxp2,v(limin:limaxp1,ljminm2:ljmaxp2))
!larger exchanges to make transfer in step2d* and step3d* easier
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(v)

   !PRINT*, 'exit OBC_APPLY_V routine'

  END SUBROUTINE obc_apply_v

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

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

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE obc_apply_v_mc  ***
!&E
!&E ** Purpose : Apply open boundary conditions for v
!&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_vs,obc_vn,obc_ve,obc_vw
   USE obccombine, ONLY : obc_sshe,obc_sshw,obc_sshs,obc_sshn, &
                          l_obc_south,l_obc_north,l_obc_west,l_obc_east

!! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( inout ) :: v
   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 ) :: u,ua,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*, 'enter OBC_APPLY_V_MC routine'

! southern and northern boundaries
! --------------------------------
 
!$OMP DO SCHEDULE(RUNTIME)
!FVDM
! cohérence des "coins" en xe
!   DO i=liminp1,limaxm1
   DO i=liminp2,limaxm1

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

! southern boundary

     IF (l_obc_south) THEN
       IF (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)
! 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)
           w1va=va(i,j1)-0.5_rsh*SQRT(grav/h)*(xea(i,j1)+xea(i,j1+1))
           w1xea =0.5_rsh*(va(i,j1)+va(i,j1+1))-SQRT(grav/h)*xea(i,j1+1)
           IF (l_star) THEN
             xe(i,j1)=(0.5_rsh*SQRT(h/grav)*(w3-al1*v(i,j1)-al2*w1va-al3*w1xea) &
                                             +0.25_rsh*al1*xe(i,j1+1))/ac0
           ELSE
!xe(i,j1)=xea(i,j1)
             v(i,j1-1)=obc_vs(i,j1)
           ENDIF
           w1v=v(i,j1)-0.5_rsh*SQRT(grav/h)*(xe(i,j1)+xe(i,j1+1))
           v(i,j1-1)= w3+al1*w1v+al2*w1va+al3*w1xea-v(i,j1)
         END IF
       ENDIF
     END IF

! northern boundary
     IF (l_obc_north) THEN
       IF (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)
           w1=0.5_rsh*(obc_vn(i,j2-1)+obc_vn(i,j2))-SQRT(grav/h)*obc_sshn(i)
           w3va=va(i,j2-1)+0.5_rsh*SQRT(grav/h)*(xea(i,j2)+xea(i,j2-1))
           w3xea=0.5_rsh*(va(i,j2-1)+va(i,j2-2))+SQRT(grav/h)*xea(i,j2-1)
           IF (l_star) THEN

             xe(i,j2)=(0.5_rsh*SQRT(h/grav)*(-w1+al1*v(i,j2-1)+al2*w3va+al3*w3xea)    &
                                             +0.25_rsh*al1*xe(i,j2-1))/ac0

           ELSE
!xe(i,j2)=xea(i,j2)
             v(i,j2)=obc_vn(i,j2-1)
           ENDIF
           w3v=v(i,j2-1)+0.5_rsh*SQRT(grav/h)*(xe(i,j2)+xe(i,j2-1))
           v(i,j2)= w1+al1*w3v+al2*w3va+al3*w3xea-v(i,j2-1)
         END IF
       ENDIF
     END IF

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


!$OMP DO SCHEDULE(RUNTIME)
!FVDM
! cohérence des "coins" en vitesse
!   DO j=ljminp1,ljmaxm1
   DO j=ljmin,ljmax

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

! western boundary

     IF (l_obc_west) THEN
       IF (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)

! no update of xe at the westhern boundary because not strong obc
! for calculation along column

           IF (hy(i1,j) /= -valmanq-h0fond) THEN
             IF ((u(i1-1,j)+u(i1,j)).GE.0.0_rsh) THEN
               v(i1,j)= obc_vw(i1,j)
             ELSE
               v(i1,j)=v(i1+1,j)
             END IF
           END IF
         END IF
       ENDIF
     END IF

! eastern boundary

     IF (l_obc_east) THEN
       IF (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)

! no update of xe at the westhern boundary because not strong obc
! for calculation along column

           IF (hy(i2,j) /= -valmanq-h0fond) THEN
             IF ((u(i2-1,j)+u(i2,j)).LE.0.0_rsh) THEN
! caract. entrante v_N+0.5(t+1) = v_ext
               v(i2,j)= obc_ve(i2,j)
             ELSE
! caract. sortante v_N+0.5(t+1) = v_N-0.5(t)
! solution de dv/dt+u*dv/dx=0, u=dx/dt
               v(i2,j)=v(i2-1,j)
            END IF
          END IF
        END IF
      ENDIF
     END IF

   END DO
!$OMP END DO

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(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 ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xe(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xe(liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(v,xe)

   !PRINT*, 'exit OBC_APPLY_V_MC routine'

  END SUBROUTINE obc_apply_v_mc

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

 END MODULE obcapplyv
