 MODULE obcapplyvz

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

#include "toolcpp.h"

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

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC obc_apply_vz,obc_apply_vz_mc  ! called by stepv, dyn3dy

   !! * Shared module variables

   !! * Private variables

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

  SUBROUTINE obc_apply_vz(vz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_apply_vz  ***
   !&E
   !&E ** Purpose : Apply open boundary conditions for vz
   !&E
   !&E ** Description : Available schemes :
   !&E                     - periodic along x or y
   !&E                     - Dirichlet
   !&E                     - MARS standard (CMS)
   !&E                     - Caracteristics method
   !&E
   !&E ** Called by : stepv, dyn3dy
   !&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       !  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
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE obccombine, ONLY : obc_vze,obc_vzw,obc_vzs,obc_vzn,     &
                          l_obc_mars,l_obc_diri,l_obc_cycl
#ifdef key_wave_qeuler_dynamics
   USE wave,       ONLY : wave_vsz
#endif

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

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

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

   PRINT_DBG*, 'enter OBC_APPLY_VZ routine'

   !! cyclic open boundary conditions
OMPNOAGRIF SINGLE
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_x_rsh(3,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
       vz(:,ig(j)-1,j)=vz(:,id(j)-2,j)
       vz(:,ig(j)  ,j)=vz(:,id(j)-1,j)
       vz(:,id(j)  ,j)=vz(:,ig(j)+1,j)
#ifdef key_shiftobc
       vz(:,ig(j)-2,j)=vz(:,id(j)-3,j)
       vz(:,id(j)+1,j)=vz(:,ig(j)+2,j)
       vz(:,id(j)+2,j)=vz(:,ig(j)+3,j)
#endif
     END DO
#endif
     !IF (jmax.EQ.ljmax) vz(:,limin:limax,jmax  )=vz(:,limin:limax,jmax-1)
     !IF (jmin.EQ.ljmin) vz(:,limin:limax,jmin  )=vz(:,limin:limax,jmin+1)
     ! Mars obc are applied
     DO i=liminp2,limaxm1
   ! Southern boundary
       j=jb(i)-1
       IF_MPI (j==jmin+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) vz(:,i,j)=vz(:,i,j+1  )
       ENDIF_MPI
   ! Northern boundary
       j=jh(i)
       IF_MPI (j==jmax-jshift_obc_north .AND. ljmin<=j .AND. j<=ljmax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) vz(:,i,j)=vz(:,i,j-1)
       ENDIF_MPI
     END DO

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

#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_y_rsh(3,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO i=limin,limax
       vz(:,i,jh(i)  )=vz(:,i,jb(i)+1)
       vz(:,i,jb(i)  )=vz(:,i,jh(i)-1)
       vz(:,i,jb(i)-1)=vz(:,i,jh(i)-2)
#ifdef key_shiftobc
       vz(:,i,jb(i)-2)=vz(:,i,jh(i)-3)
       vz(:,i,jh(i)+1)=vz(:,i,jb(i)+2)
       vz(:,i,jh(i)+2)=vz(:,i,jb(i)+3)
#endif
     END DO
#endif
   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)
       IF_MPI (i==imin+1+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)
           IF (hy(ii,j) /= -valmanq-h0fond) vz(:,ii,j)=obc_vzw(:,ii,j)
         END DO
#else
         IF (hy(i,j) /= -valmanq-h0fond) THEN
           vz(:,i,j)=obc_vzw(:,i,j)
#ifdef key_wave_qeuler_dynamics
           vz(:,i,j)=vz(:,i,j)-wave_vsz(:,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 (hy(ii,j) /= -valmanq-h0fond) vz(:,ii,j)=obc_vze(:,ii,j)
         END DO
#else
         IF (hy(i,j) /= -valmanq-h0fond) THEN
           vz(:,i,j)=obc_vze(:,i,j)
#ifdef key_wave_qeuler_dynamics
           vz(:,i,j)=vz(:,i,j)-wave_vsz(:,i-1,j)
#endif
         END IF
#endif
       ENDIF_MPI
     END DO
OMPNOAGRIF END DO

OMPNOAGRIF SINGLE
     DO i=liminp2,limaxm1

       j=jb(i)-1
       IF_MPI (j==jmin+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)-1
           IF (hy(i,jj) /= -valmanq-h0fond) vz(:,i,jj)=obc_vzs(:,i,jj)
         END DO
#else
         IF (hy(i,j) /= -valmanq-h0fond) THEN
           vz(:,i,j)=obc_vzs(:,i,j)
#ifdef key_wave_qeuler_dynamics
           vz(:,i,j)=vz(:,i,j)-wave_vsz(:,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 (hy(i,jj) /= -valmanq-h0fond) vz(:,i,jj)=obc_vzn(:,i,jj)
         END DO
#else
         IF (hy(i,j) /= -valmanq-h0fond) THEN
           vz(:,i,j)=obc_vzn(:,i,j)
#ifdef key_wave_qeuler_dynamics
           vz(:,i,j)=vz(:,i,j)-wave_vsz(:,i,j-1)
#endif
         END IF
#endif
       ENDIF_MPI

     END DO
OMPNOAGRIF END SINGLE

   !! Ends Dirichlet
   END IF

#ifdef key_lagoon
   !!  open boundary conditions (for LAGOON only)
   !! l_obc_mars=.true always
     DO ng=1,nbgrau
       i0=igrau(ng)
       j0=jgrau(ng)
       
   ! Southern boundary
       IF (type_grau(ng)==3) THEN
       !! grau sur la frontiere Sud
         j=j0-1 
         IF_MPI ( (i0 >= liminm2) .AND. (i0 <= limaxp2) .AND.(j >= ljminm2) .AND. (j <= ljmaxp2) ) THEN
           IF (hy(i0,j) /= -valmanq-h0fond) vz(:,i0,j)=vz(:,i0,j0)
         ENDIF_MPI
       END IF

   ! Northern boundary
       IF (type_grau(ng)==1) THEN
       !! grau sur la frontiere Nord
         IF_MPI ( (i0 >= liminm2) .AND. (i0 <= limaxp2) .AND.(j0 >= ljminm2) .AND. (j0 <= ljmaxp2) ) THEN
           IF (hy(i0,j0) /= -valmanq-h0fond) vz(:,i0,j0)=vz(:,i0,j0-1)
         ENDIF_MPI 
       END IF

     END DO

#else
   !! Default open boundary conditions (MARS standard)

   IF (l_obc_mars) THEN

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

   ! Western boundary

       i=ig(j)
       IF_MPI (i==imin+1+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
         IF (hy(i,j) /= -valmanq-h0fond) vz(:,i,j  )=vz(:,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 (hy(i,j) /= -valmanq-h0fond) vz(:,i,j)=vz(:,i-1,j)
       ENDIF_MPI

     END DO
OMPNOAGRIF END DO


OMPNOAGRIF SINGLE
     DO i=liminp2,limaxm1

   ! Southern boundary

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

   ! Northern boundary

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

     END DO
OMPNOAGRIF END SINGLE

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

   PRINT_DBG*, 'exit OBC_APPLY_VZ routine'

  END SUBROUTINE obc_apply_vz

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

  SUBROUTINE obc_apply_vz_mc(uz,vz,boz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_apply_vz_mc  **
   !&E
   !&E ** Purpose : Apply open boundary conditions for vz
   !&E
   !&E ** Description : - Compatible with characteristics method
   !&E
   !&E ** Called by : dyn3dxy, dyn3dyx
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :  uz,vz
   !&E
   !&E ** Modified variables : vz
   !&E
   !&E ** Reference : Blayo & Debreu, 2005. Ocean modelling, 9, 231-252
   !&E
   !&E ** History :
   !&E       !  2007-05  (F. Vandermeirsch, J. Marin, V. Garnier) Original code
   !&E       !  2009-05  (F. Vandermeirsch) characteristics method with new time scheme
   !&E
   !&E---------------------------------------------------------------------
   USE comvarp2d
   USE comvars3d
   USE obccombine, ONLY : obc_vze,obc_vzw,obc_vzn,obc_vzs

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( inout )  :: vz
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: uz
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )  :: boz

   !! * local declarations
   INTEGER                     :: i,j,k,i1,i2,j1,j2

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

OMPNOAGRIF DO SCHEDULE(RUNTIME)
   DO i=liminp2,limaxm1

   ! Southern boundary

     j1 = jb(i)
     IF_MPI (j1==jmin+1+jshift_obc_south .AND. ljmin<=j1 .AND. j1<=ljmax) THEN
       IF (hy(i,j1-1) /= -valmanq-h0fond) vz(:,i,j1-1)=obc_vzs(:,i,j1-1)
     ENDIF_MPI

   ! Northern open boundary

     j2 = jh(i)
     IF_MPI (j2==jmax-jshift_obc_north .AND. ljmin<=j2 .AND. j2<=ljmax) THEN
       IF (hy(i,j2) /= -valmanq-h0fond)   vz(:,i,j2)=obc_vzn(:,i,j2) 
     ENDIF_MPI

   END DO 
OMPNOAGRIF END DO


OMPNOAGRIF DO SCHEDULE(RUNTIME)
!   DO j=ljminp1,ljmaxm1
   DO j=ljminp2,ljmaxm1

   ! western  boundary

     i1=ig(j)
     IF_MPI (i1==imin+1+ishift_obc_west .AND. limin<=i1 .AND. i1<=limax) THEN
       IF (hx(i1,j) /= -valmanq-h0fond) THEN
!FVDM     IF (hx(i1,j) /= -valmanq-h0fond .and. hy(i1,j) /= -valmanq-h0fond) THEN
! car MC barocline n est pas completement developpee
         DO k=1,kmax
           IF ((uz(k,i1-1,j)+uz(k,i1,j)).GE.0.0_rsh) THEN
             vz(k,i1,j)=obc_vzw(k,i1,j)
           ELSE       
             vz(k,i1,j)=vz(k,i1+1,j)
           END IF
         END DO
       END IF
     ENDIF_MPI

   ! eastern  boundary

     i2=id(j)
     IF_MPI (i2==imax-ishift_obc_east .AND. limin<=i2 .AND. i2<=limax) THEN
       IF (hx(i2-1,j) /= -valmanq-h0fond) THEN
!FVDM     IF (hx(i2-1,j) /= -valmanq-h0fond .and. hy(i2,j) /= -valmanq-h0fond) THEN
! car MC barocline n est pas completement developpee
         DO k=1,kmax
           IF ((uz(k,i2-1,j)+uz(k,i2,j)).LE.0.0_rsh) THEN 
             vz(k,i2,j)=obc_vze(k,i2,j)
           ELSE
             vz(k,i2,j)=vz(k,i2-1,j)
           END IF
         END DO
       END IF
     ENDIF_MPI

   END DO  
OMPNOAGRIF END DO

!---------------------------------------------------       
!Traitement des CFO par methode des caracteristiques
!----------------------------------------------------
!
#ifdef key_obc_char_barocline

   PRINT*, 'key_obc_char_barocline not ready yet'
   STOP
#if defined key_MPI_2D
   PRINT*, 'obcapplyuz mc part not MPIed yet'
   STOP
#endif

#ifdef key_obc_char_east
   DO j=jmin+2,jmax-1
       IF (h0(id(j),j-1)/=-valmanq .AND. h0(id(j),j+1)/=-valmanq) &
         CALL obcch3d_vz_exp(uz,vz,boz,id(j),j,3)
   END DO
#endif

#ifdef key_obc_char_west
   DO j=jmin+2,jmax-1
       IF (h0(ig(j),j-1)/=-valmanq .AND. h0(ig(j),j+1)/=-valmanq) &
         !jm on applique la methode barocline que sur 1/3 de la frontiere
         !jmif (j>jmax/3 .and. j<2*jmax/3) THEN
           CALL obcch3d_vz_exp(uz,vz,boz,ig(j)-1,j,4)
         !jmendif
   END DO
#endif

#ifdef key_obc_char_north
   DO j=imin+2,imax-1
       IF (h0(j-1,jh(j))/=-valmanq .AND. h0(j+1,jh(j))/=-valmanq) &
         CALL obcch3d_vz_exp(uz,vz,boz,j,jh(j),1)
   END DO
#endif

#ifdef key_obc_char_south
   DO j=imin+2,imax-1
       IF (h0(j-1,jb(j))/=-valmanq .AND. h0(j+1,jb(j))/=-valmanq) &
         CALL obcch3d_vz_exp(uz,vz,boz,j,jb(j)-1,2)
   END DO


   WRITE (*,*) 'routine limavz ... frontiere Sud non implementee'
   STOP
#endif
#endif  
   
END SUBROUTINE obc_apply_vz_mc

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

END MODULE obcapplyvz
