 MODULE turbadvobc3dapply

#ifdef key_turbadv
   !!======================================================================
   !!                   ***  MODULE  turbadvobc3dapply  ***
   !! Ocean dynamics : apply open boundary conditions of variables
   !!                  turbulence cinetic energy and turbulence dissipation
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC turbadv_obc3d_apply           ! routine called by stepu.F90 stepv.F90
!#ifdef key_agrif
!   PUBLIC ts_obc3d_apply_agrif          ! routine called by stepu.F90 stepv.F90
!#endif

   !! * Shared module variables

   !! * Private variables

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

  SUBROUTINE turbadv_obc3d_apply(u,v,ect,psi)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turbadv_obc3d_apply  ***
   !&E
   !&E ** Purpose : apply open boundary conditions for turbulent kinetic energy
   !&E              and turbulent dissiaption field
   !&E
   !&E ** Description : apply external value when water mass enters the domain
   !&E
   !&E ** Called by : stepu, stepv
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : ect,psi,hx,u,obc_sw,obc_tw,dx_u,obc_se,obc_te,
   !&E                     v,obc_ss,obc_ts,dy_v,obc_sn,obc_tn
   !&E 
   !&E ** Modified variables : ect,psi
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !                       Original code
   !&E       !  2014-08 (F. Ganthy) Addapted following ts_advobc3dapply.F90
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters, ONLY : imin,imax,jmin,jmax,kmax,                      &
                          ljminp2,ljmaxm1,ljmin,ljmax,ljminm1,ljmaxp1,   &
                          limin,limax,liminm1,limaxp2,liminp2,limaxm1,   &
                          ljmaxp2,liminm2,ljminm2,rsh,rlg,num_testcase,  &
                          ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north

   USE comvars2d,  ONLY : dt,ig,hx,h0fond,dx_u,id,jb,hy,dy_v,jh,l_obc_cycl_x,l_obc_cycl_y
   USE comvars3d,  ONLY : l_testcase,valmanq
   !USE obcogcm,    ONLY : obc_tw,obc_se,obc_te,obc_ss,obc_ts,obc_sn,obc_tn,obc_sw,obc_coefrel
   USE obccombine, ONLY : l_obc_char,l_obc_south,l_obc_north,l_obc_east,l_obc_west
   USE_MPI toolmpi,ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: ect,psi
   REAL(KIND=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: u,v

   !! * Local declarations
   INTEGER                   :: i,j,k
   REAL(KIND=rsh)            :: uu,vv,cou,dtrelax
   REAL(KIND=rlg)            :: dtm

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

   IF (l_testcase .AND. (num_testcase==1 &
       .OR. num_testcase==33 .OR. num_testcase == 37 &
! cas kelvin a virer
       .OR. num_testcase==3 .OR. num_testcase==38) &
      ) THEN

     ! None

   ELSE ! general

     dtm=0.5_rlg*dt
     !dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
     !IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh

!$OMP SINGLE
   ! periodic conditions along x-axis
     IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc
       CALL_MPI ex_peri_x_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,ect(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_peri_x_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,psi(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
       CALL_MPI ex_peri_x_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,ect(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_peri_x_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,psi(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
       DO j=ljmin,ljmax
   ! western boundary
         ect(:,ig(j),j)   = ect(:,id(j)-1,j)
         ect(:,ig(j)-1,j) = ect(:,id(j)-2,j)
         psi(:,ig(j),j)   = psi(:,id(j)-1,j)
         psi(:,ig(j)-1,j) = psi(:,id(j)-2,j)
#ifdef key_shiftobc
         ect(:,ig(j)-2,j) = ect(:,id(j)-3,j)
         psi(:,ig(j)-2,j) = psi(:,id(j)-3,j)
#endif
   ! eastern boundary
         ect(:,id(j),j) = ect(:,ig(j)+1,j)
         psi(:,id(j),j) = psi(:,ig(j)+1,j)
#ifdef key_shiftobc
         ect(:,id(j)+1,j) = ect(:,ig(j)+2,j)
         ect(:,id(j)+2,j) = ect(:,ig(j)+3,j)
         psi(:,id(j)+1,j) = psi(:,ig(j)+2,j)
         psi(:,id(j)+2,j) = psi(:,ig(j)+3,j)
#endif
       END DO
#endif
     END IF  ! ends periodic conditions along x-axis

   ! periodic conditions along y-axis
     IF (l_obc_cycl_y) THEN

#ifdef key_MPI_2D
#ifdef key_shiftobc
       CALL_MPI ex_peri_y_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,ect(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_peri_y_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,psi(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
       CALL_MPI ex_peri_y_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,ect(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_peri_y_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,psi(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
       DO i=limin,limax
   ! southern boundary
         ect(:,i,jb(i))   = ect(:,i,jh(i)-1)
         ect(:,i,jb(i)-1) = ect(:,i,jh(i)-2)
         psi(:,i,jb(i))   = psi(:,i,jh(i)-1)
         psi(:,i,jb(i)-1) = psi(:,i,jh(i)-2)
#ifdef key_shiftobc
         ect(:,i,jb(i)-2) = ect(:,i,jh(i)-3)
         psi(:,i,jb(i)-2) = psi(:,i,jh(i)-3)
#endif
   ! northern boundary
         ect(:,i,jh(i)) = ect(:,i,jb(i)+1)
         psi(:,i,jh(i)) = psi(:,i,jb(i)+1)
#ifdef key_shiftobc
         ect(:,i,jh(i)+1) = ect(:,i,jb(i)+2)
         ect(:,i,jh(i)+2) = ect(:,i,jb(i)+3)
         psi(:,i,jh(i)+1) = psi(:,i,jb(i)+2)
         psi(:,i,jh(i)+2) = psi(:,i,jb(i)+3)
#endif
       END DO
#endif

     END IF  ! ends periodic conditions along y-axis
!$OMP END SINGLE


     IF (.NOT.l_obc_cycl_x .AND. (l_obc_west .OR. l_obc_east)) THEN

       DO j=ljminp2,ljmaxm1

   ! Western boundary
         i=ig(j)
         IF_MPI (i==imin+1+ishift_obc_west .AND. limin<=i .AND. i<=limax) THEN
           IF (hx(i,j) /= -valmanq-h0fond) THEN
             DO k=1,kmax
               IF (l_obc_char) THEN
                 uu=u(k,i-1,j)+u(k,i,j)
               ELSE
                 uu=u(k,i,j)
               END IF
               IF (uu > 0.0_rsh) THEN
                 ! No external source
                 ect(k,i,j)=ect(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
                 psi(k,i,j)=psi(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
               ELSE
                 IF (l_obc_char) THEN
               ! Extrapolation d ordre 2
                   ect(k,i,j)=2.0_rsh*ect(k,i+1,j)-ect(k,i+2,j)
                   psi(k,i,j)=2.0_rsh*psi(k,i+1,j)-psi(k,i+2,j)
                 ELSE
                   cou=uu*dtm/dx_u(i,j)
                   ect(k,i,j)=ect(k,i,j)*(1.0_rsh+cou)-cou*ect(k,i+1,j)
                   psi(k,i,j)=psi(k,i,j)*(1.0_rsh+cou)-cou*psi(k,i+1,j)
                 END IF
               END IF
             END DO
           END IF
         ENDIF_MPI

   ! Eastern boundary
         i=id(j)
         IF_MPI (i==imax-ishift_obc_east .AND. limin<=i .AND. i<=limax) THEN
           IF (hx(i-1,j) /= -valmanq-h0fond) THEN
             DO k=1,kmax
               IF (l_obc_char) THEN
                 uu=u(k,i-1,j)+u(k,i,j)
               ELSE
                 uu=u(k,i-1,j)
               END IF
               IF (uu < 0.0_rsh) THEN
                 ! No external source
                 ect(k,i,j)=ect(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
                 psi(k,i,j)=psi(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
               ELSE
                 IF (l_obc_char) THEN
                   ect(k,i,j)=2.0_rsh*ect(k,i-1,j)-ect(k,i-2,j)
                   psi(k,i,j)=2.0_rsh*psi(k,i-1,j)-psi(k,i-2,j)
                 ELSE
                   cou=uu*dtm/dx_u(i-1,j)
                   ect(k,i,j)=ect(k,i,j)*(1.0_rsh-cou)+cou*ect(k,i-1,j)
                   psi(k,i,j)=psi(k,i,j)*(1.0_rsh-cou)+cou*psi(k,i-1,j)
                 END IF
               END IF
             END DO
           END IF
         ENDIF_MPI

       END DO

     END IF  ! ends test on .not.l_obc_cycl_x

     IF (.NOT.l_obc_cycl_y .AND. (l_obc_north .OR. l_obc_south)) THEN

       DO j=ljmin,ljmax
       DO i=liminp2,limaxm1

   ! Southern boundary
         IF (j == jb(i) .AND. j==jmin+1+jshift_obc_south) THEN
           IF (hy(i,j) /= -valmanq-h0fond) THEN
             DO k=1,kmax
               IF (l_obc_char) THEN
                 vv=v(k,i,j-1)+v(k,i,j)
               ELSE
                 vv=v(k,i,j)
               END IF
               IF (vv > 0.0_rsh) THEN
                 ! No external source
                 ect(k,i,j)=ect(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
                 psi(k,i,j)=psi(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
               ELSE
                 IF (l_obc_char) THEN
                   ect(k,i,j)=2.0_rsh*ect(k,i,j+1)-ect(k,i,j+2)
                   psi(k,i,j)=2.0_rsh*psi(k,i,j+1)-psi(k,i,j+2)
                 ELSE
                   cou=vv*dtm/dy_v(i,j)
                   ect(k,i,j)=ect(k,i,j)*(1.0_rsh+cou)-cou*ect(k,i,j+1)
                   psi(k,i,j)=psi(k,i,j)*(1.0_rsh+cou)-cou*psi(k,i,j+1)
                 END IF
               END IF
             END DO
           END IF

   ! Northern boundary
         ELSE IF (j == jh(i) .AND. j==jmax-jshift_obc_north) THEN
           IF (hy(i,j-1) /= -valmanq-h0fond) THEN
             DO k=1,kmax-1
               IF (l_obc_char) THEN
                 vv=v(k,i,j-1)+v(k,i,j)
               ELSE
                 vv=v(k,i,j-1)
               END IF
               IF (vv < 0.0_rsh) THEN
                 ! No external source
                 ect(k,i,j)=ect(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
                 psi(k,i,j)=psi(k,i,j)*(1.0_rsh-dtrelax)+0.0_rsh*dtrelax
               ELSE
                 IF (l_obc_char) THEN
                   ect(k,i,j)=2.0_rsh*ect(k,i,j-1)-ect(k,i,j-2)
                   psi(k,i,j)=2.0_rsh*psi(k,i,j-1)-psi(k,i,j-2)
                 ELSE
                   cou=vv*dtm/dy_v(i,j-1)
                   ect(k,i,j)=ect(k,i,j)*(1.0_rsh-cou)+cou*ect(k,i,j-1)
                   psi(k,i,j)=psi(k,i,j)*(1.0_rsh-cou)+cou*psi(k,i,j-1)
                 END IF
               END IF
             END DO
           END IF
         END IF

       END DO
       END DO

     END IF  ! ends test on .not.l_obc_cycl_y

   END IF   ! ends test on l_testcase

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,ect(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,ect(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,psi(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,psi(:,liminm1:limaxp2,ljminm1:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(ect,psi)

  END SUBROUTINE turbadv_obc3d_apply

   !!======================================================================
#endif

  END MODULE turbadvobc3dapply
