 MODULE subobc2dapply

#ifdef key_substance

   !!======================================================================
   !!                   ***  MODULE  subobc2dapply  ***
   !! Ocean dynamics:  apply open boundary conditions of variables
   !!                  (tracers other than temperature or salinity)
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC sub_obc2d_apply               ! routine called by stepu.F90 stepv.F90

   !! * Shared module variables

   !! * Private variables

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

  SUBROUTINE sub_obc2d_apply(u,v,cv_wat)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_obc2d_apply  ***
   !&E
   !&E ** Purpose : open boundary conditions for tracers other than temperature and salinity
   !&E              For incoming conditions, used values are stored in the 
   !&E              obc_cv_wat_s, obc_cv_wat_n, obc_cv_wat_w and obc_cv_wat_e arrays,
   !&E              filled in the "subreaddat" routine.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : cv_wat,hx,u,obc_cv_wat_w,dx_u,obc_cv_wat_e,hy,v,
   !&E                     obc_cv_wat_s,dy_v,obc_cv_wat_n
   !&E 
   !&E ** Modified variables : cv_wat
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2004-06
   !&E       !  2006-11 (V. Garnier) module - mise en forme
   !&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
   !&E       !  2009-04 (A. Mano) MPI2D
   !&E       !  2010-04 (T. Odaka) MPI2D cyclic
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : dt,ig,hx,h0fond,dx_u,id,jb,hy,dy_v,jh,          &
                            l_obc_cycl_x,l_obc_cycl_y
   USE comsubstance, ONLY : obc_cv_wat_w,nvp,nv_adv,obc_cv_wat_e,           &
                            obc_cv_wat_s,obc_cv_wat_n
   USE obcogcm,      ONLY : obc_coefrel
   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh

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

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

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

   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
!!!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
     CALL_MPI ex_peri_x_rsh(3,(nv_adv-nvp),liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(nvp+1:nv_adv,1,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(1,(nv_adv-nvp),liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(nvp+1:nv_adv,1,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
   ! western boundary
       cv_wat(nvp+1:nv_adv,1,ig(j),j)=cv_wat(nvp+1:nv_adv,1,id(j)-1,j)
       cv_wat(nvp+1:nv_adv,1,ig(j)-1,j)=cv_wat(nvp+1:nv_adv,1,id(j)-2,j)
!!!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(nvp+1:nv_adv,1,ig(j)-2,j)=cv_wat(nvp+1:nv_adv,1,id(j)-3,j)
#endif
   ! eastern boundary
       cv_wat(nvp+1:nv_adv,1,id(j),j)=cv_wat(nvp+1:nv_adv,1,ig(j)+1,j)
!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(nvp+1:nv_adv,1,id(j)+1,j)=cv_wat(nvp+1:nv_adv,1,ig(j)+2,j)
       cv_wat(nvp+1:nv_adv,1,id(j)+2,j)=cv_wat(nvp+1:nv_adv,1,ig(j)+3,j)
#endif
     END DO
#endif
   END IF

   ! periodic conditions along y-axis
   IF (l_obc_cycl_y) THEN
#ifdef key_MPI_2D
!!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
     CALL_MPI ex_peri_y_rsh(3,(nv_adv-nvp),liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(nvp+1:nv_adv,1,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(1,(nv_adv-nvp),liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(nvp+1:nv_adv,1,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
     DO i=limin,limax
   ! southern boundary
       cv_wat(nvp+1:nv_adv,1,i,jb(i))=cv_wat(nvp+1:nv_adv,1,i,jh(i)-1)
       cv_wat(nvp+1:nv_adv,1,i,jb(i)-1)=cv_wat(nvp+1:nv_adv,1,i,jh(i)-2)
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(nvp+1:nv_adv,1,i,jb(i)-2)=cv_wat(nvp+1:nv_adv,1,i,jh(i)-3)
#endif
   ! northern boundary
       cv_wat(nvp+1:nv_adv,1,i,jh(i))=cv_wat(nvp+1:nv_adv,1,i,jb(i)+1)
!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(nvp+1:nv_adv,1,i,jh(i)+1)=cv_wat(nvp+1:nv_adv,1,i,jb(i)+2)
       cv_wat(nvp+1:nv_adv,1,i,jh(i)+2)=cv_wat(nvp+1:nv_adv,1,i,jb(i)+3)
#endif
     END DO
#endif
   END IF
!$OMP END SINGLE

   !----------------------------------------------
   ! on traite l advection quand le courant sort
   ! sur les limites gauche et droite
   !----------------------------------------------

   IF (.NOT.l_obc_cycl_x) 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
         uu=u(i,j)
         IF (uu > 0.0_rsh) THEN
           cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_w(nvp+1:nv_adv,1,1,j)*dtrelax
         ELSE
           cou=uu*dtm/dx_u(i,j)
           cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh+cou)-cou*cv_wat(nvp+1:nv_adv,1,i+1,j)
         END IF
       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
         uu=u(i-1,j)
         IF (uu < 0.0_rsh) THEN
           cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_e(nvp+1:nv_adv,1,imax,j)*dtrelax
         ELSE
           cou=uu*dtm/dx_u(i-1,j)
           cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh-cou)+cou*cv_wat(nvp+1:nv_adv,1,i-1,j)
         END IF
       END IF
       ENDIF_MPI

     END DO
   END IF

   !----------------------------------------------
   ! on traite l advection quand le courant sort
   ! sur les limites basse et haute
   !----------------------------------------------

   IF (.NOT.l_obc_cycl_y) 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
             vv=v(i,j)
             IF (vv > 0.0_rsh) THEN
               cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_s(nvp+1:nv_adv,1,i,1)*dtrelax
             ELSE
               cou=vv*dtm/dy_v(i,j)
               cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh+cou)-cou*cv_wat(nvp+1:nv_adv,1,i,j+1)
             END IF
         END IF
       END IF

   ! northern boundary
       IF (j == jh(i) .AND. j==jmax-jshift_obc_north) THEN
         IF (hy(i,j-1) /= -valmanq-h0fond) THEN
             vv=v(i,j-1)
             IF (vv < 0.0_rsh) THEN
               cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_n(nvp+1:nv_adv,1,i,jmax)*dtrelax
             ELSE
               cou=vv*dtm/dy_v(i,j-1)
               cv_wat(nvp+1:nv_adv,1,i,j)=cv_wat(nvp+1:nv_adv,1,i,j)*(1.0_rsh-cou)+cou*cv_wat(nvp+1:nv_adv,1,i,j-1)
             END IF
         END IF
       END IF

     END DO
     END DO
   END IF

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

  END SUBROUTINE sub_obc2d_apply


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

#endif

 END MODULE subobc2dapply
