 MODULE obcapplyuz

   !!=======================================================================
   !!                        ***  MODULE obcapplyuz  ***
   !! 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_uz,obc_apply_uz_mc    ! called by stepu, dyn3dx 

   !! * Shared module variables

   !! * Private variables

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

  SUBROUTINE obc_apply_uz(uz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_apply_uz  ***
   !&E
   !&E ** Purpose : Apply open boundary conditions for uz
   !&E
   !&E ** Description : Available schemes :
   !&E                     - periodic along x or y
   !&E                     - Dirichlet
   !&E                     - MARS standard (CMS)
   !&E
   !&E ** Called by : stepu, 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
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE obccombine, ONLY : obc_uze,obc_uzw,obc_uzs,obc_uzn,     &
                          l_obc_mars,l_obc_diri,l_obc_cycl
#ifdef key_wave_qeuler_dynamics
   USE wave,       ONLY : wave_usz
#endif

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

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

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

   PRINT_DBG*, 'enter OBC_APPLY_UZ routine'

   !! cyclic open boundary conditions
OMPNOAGRIF SINGLE
   IF (l_obc_cycl_x) THEN
     !IF (jmax.EQ.ljmax) uz(:,limin:limax,jmax  )=uz(:,limin:limax,jmax-1)
     !IF (jmin.EQ.ljmin) uz(:,limin:limax,jmin+1)=uz(:,limin:limax,jmin+2)
     !IF (jmin.EQ.ljmin) uz(:,limin:limax,jmin  )=uz(:,limin:limax,jmin+2)
     ! mars obc 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) uz(:,i,j)=uz(:,i,j+1)
         IF (hx(i,j-1) /= -valmanq-h0fond) uz(:,i,j-1)=uz(:,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) uz(:,i,j)=uz(:,i,j-1)
       ENDIF_MPI
     END DO

#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_x_rsh(3,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
       uz(:,ig(j)-1,j)=uz(:,id(j)-2,j)
       uz(:,ig(j)  ,j)=uz(:,id(j)-1,j)
       uz(:,id(j)  ,j)=uz(:,ig(j)+1,j)
#ifdef key_shiftobc
       uz(:,ig(j)-2,j)=uz(:,id(j)-3,j)
       uz(:,id(j)+1,j)=uz(:,ig(j)+2,j)
       uz(:,id(j)+2,j)=uz(:,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,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO i=limin,limax
       uz(:,i,jh(i))=uz(:,i,jb(i)+1)
       uz(:,i,jb(i))=uz(:,i,jh(i)-1)
       uz(:,i,jb(i)-1)=uz(:,i,jh(i)-2)
#ifdef key_shiftobc
       uz(:,i,jb(i)-2)=uz(:,i,jh(i)-3)
       uz(:,i,jh(i)+1)=uz(:,i,jb(i)+2)
       uz(:,i,jh(i)+2)=uz(:,i,jb(i)+3)
#endif
     END DO
#endif
     !IF (imin.EQ.limin) uz(:,imin,ljmin:ljmax)=uz(:,imin+1,ljmin:ljmax)
     !IF (imax.EQ.limax) uz(:,imax,ljmin:ljmax)=uz(:,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) uz(:,i,j)=uz(:,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) uz(:,i,j)=uz(:,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) uz(:,ii,j)=obc_uzw(:,ii,j)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uzw(:,i,j)
#ifdef key_wave_qeuler_dynamics
           uz(:,i,j)=uz(:,i,j)-wave_usz(:,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) uz(:,ii,j)=obc_uze(:,ii,j)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uze(:,i,j)
#ifdef key_wave_qeuler_dynamics
           uz(:,i,j)=uz(:,i,j)-wave_usz(:,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) uz(:,i,jj)=obc_uzs(:,i,jj)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uzs(:,i,j)
#ifdef key_wave_qeuler_dynamics
           uz(:,i,j)=uz(:,i,j)-wave_usz(:,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) uz(:,i,jj)=obc_uzn(:,i,jj)
         END DO
#else
         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uzn(:,i,j)
#ifdef key_wave_qeuler_dynamics
           uz(:,i,j)=uz(:,i,j)-wave_usz(:,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)
       
   ! 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) uz(:,i,j0)=uz(:,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) uz(:,i0,j0)=uz(:,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) uz(:,i,j)=uz(:,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) uz(:,i,j)=uz(:,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) uz(:,i,j)=uz(:,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) uz(:,i,j)=uz(:,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_UZ routine'

  END SUBROUTINE obc_apply_uz

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

  SUBROUTINE obc_apply_uz_mc(uz,vz,boz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_apply_uz_mc  ***
   !&E
   !&E ** Purpose : Apply open boundary conditions for u
   !&E
   !&E ** Description : - Compatible with characteristics method
   !&E
   !&E
   !&E ** Called by : dyn3dxy, dyn3dyx
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :  uz,vz
   !&E
   !&E ** Modified variables : uz
   !&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       !  2009-07  (V. Garnier) MPI2D
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvarp2d
   USE comvars3d
   USE obccombine, ONLY : obc_uzw,obc_uze,obc_uzn,obc_uzs

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

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

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

   PRINT_DBG*, 'enter OBC_APPLY_UZ_MC routine'

   ! eastern and western boundaries
   ! ------------------------------
OMPNOAGRIF DO SCHEDULE(RUNTIME)
   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-1,j) /= -valmanq-h0fond) uz(:,i1-1,j)=obc_uzw(:,i1-1,j)
     ENDIF_MPI

   ! Eastern boundary

     i2 = id(j)
     IF_MPI (i2==imax-ishift_obc_east .AND. limin<=i2 .AND. i2<=limax) THEN
       IF (hx(i2,j) /= -valmanq-h0fond) uz(:,i2,j)=obc_uze(:,i2,j)
     ENDIF_MPI

   END DO     
OMPNOAGRIF END DO

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


   ! Southern open boundary

     j1=jb(i)

     IF_MPI (j1==jmin+1+jshift_obc_south .AND. ljmin<=j1 .AND. j1<=ljmax) THEN
       IF (hy(i,j1) /= -valmanq-h0fond) THEN
!FVDM      IF (hx(i,j1) /= -valmanq-h0fond .and. hy(i,j1) /= -valmanq-h0fond) THEN
! car MC barocline n est pas completement developpee
         DO k=1,kmax
           IF ((vz(k,i,j1-1)+vz(k,i,j1)).GE.0.0_rsh) THEN
             uz(k,i,j1)=obc_uzs(k,i,j1)
           ELSE
             uz(k,i,j1)=uz(k,i,j1+1)
           END IF
         END DO
       END IF
     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-1) /= -valmanq-h0fond) THEN
!FVDM      IF (hx(i,j2) /= -valmanq-h0fond .and. hy(i,j2-1) /= -valmanq-h0fond) THEN
! car MC barocline n est pas completement developpee
         DO k=1,kmax
           IF ((vz(k,i,j2-1)+vz(k,i,j2)).LE.0.0_rsh) THEN
             uz(k,i,j2)=obc_uzn(k,i,j2)
           ELSE
             uz(k,i,j2)=uz(k,i,j2-1)
           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_uz_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
       !jm IF (j>jmax/3 .and. j<2*jmax/3) THEN
       CALL obcch3d_uz_exp(uz,vz,boz,ig(j)-1,j,4)
       !jm END IF
   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_uz_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_uz_exp(uz,vz,boz,j,jb(j)-1,2)
   END DO
#endif
#endif

END SUBROUTINE obc_apply_uz_mc


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

END MODULE obcapplyuz
