 MODULE obcapplyuz

!!=======================================================================
!!                        ***  MODULE obcapplyuz  ***
!! 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_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


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

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


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

   !PRINT*, 'enter OBC_APPLY_UZ routine'

!! cyclic open boundary conditions
!$OMP 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 (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
! Northern open boundary
       j=jh(i)
       IF (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
     END DO



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


   ELSE IF (l_obc_cycl_y) THEN


     CALL ex_peri_y_rsh(2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))


!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 (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
! Eastern boundary
       i=id(j)
       IF (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
     END DO
   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)-1
       IF (i==imin+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN

         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uzw(:,i,j)

         END IF

       ENDIF

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

         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uze(:,i,j)

         END IF

       ENDIF
     END DO
!$OMP END DO

!$OMP SINGLE
     DO i=liminp1,limaxm1
       j=jb(i)
       IF (j==jmin+1+jshift_obc_south .AND. ljmin<=j .AND. j<=ljmax) THEN

         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uzs(:,i,j)

         END IF

       ENDIF

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

         IF (hx(i,j) /= -valmanq-h0fond) THEN
           uz(:,i,j)=obc_uzn(:,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=ljminp2,ljmaxm1

! Western boundary

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

! Eastern boundary

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

     END DO
!$OMP END DO


!$OMP SINGLE
     DO i=liminp1,limaxm1

! Southern open boundary

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

! Northern open boundary

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

     END DO
!$OMP END SINGLE

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


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

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

! Eastern boundary

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

   END DO     
!$OMP END DO

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


! Southern open boundary

     j1=jb(i)

     IF (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

! Northern open boundary

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

   END DO
!$OMP END DO

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


END SUBROUTINE obc_apply_uz_mc


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

END MODULE obcapplyuz
