 MODULE obcapplyvz

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





























































!! * Modules used
   USE parameters
   USE comvars2d
   USE 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


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

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


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

   !PRINT*, 'enter OBC_APPLY_VZ routine'

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


     CALL ex_peri_x_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))


!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 (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
! Northern boundary
       j=jh(i)
       IF (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
     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 (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
! Eastern boundary
       i=id(j)
       IF (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
     END DO



     CALL ex_peri_y_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,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
           vz(:,i,j)=obc_vzw(:,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
           vz(:,i,j)=obc_vze(:,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
           vz(:,i,j)=obc_vzs(:,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
           vz(:,i,j)=obc_vzn(:,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) vz(:,i,j  )=vz(:,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) vz(:,i,j)=vz(:,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) vz(:,i,j)=vz(:,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) vz(:,i,j)=vz(:,i,j-1)
       ENDIF

     END DO
!$OMP END SINGLE

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


   !PRINT*, '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

!$OMP DO SCHEDULE(RUNTIME)
   DO i=liminp2,limaxm1

! Southern boundary

     j1 = jb(i)
     IF (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

! Northern open boundary

     j2 = jh(i)
     IF (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

   END DO 
!$OMP END DO


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

! western  boundary

     i1=ig(j)
     IF (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

! eastern  boundary

     i2=id(j)
     IF (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

   END DO  
!$OMP END DO

!---------------------------------------------------
!Traitement des CFO par methode des caracteristiques
!----------------------------------------------------
!

   
END SUBROUTINE obc_apply_vz_mc

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

END MODULE obcapplyvz
