 MODULE subobc3dapply

#ifdef key_substance

   !!======================================================================
   !!                   ***  MODULE  subobc3dapply  ***
   !! 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_obc3d_apply               ! routine called by stepu.F90 stepv.F90

   !! * Interface
   INTERFACE sub_obc3d_apply
#ifdef key_lagoon
     MODULE PROCEDURE sub_obc3d_apply_lagoon
#else
     MODULE PROCEDURE sub_obc3d_apply_ref
#endif
   END INTERFACE

   !! * Shared module variables

   !! * Private variables

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

#if ! defined key_lagoon
  SUBROUTINE sub_obc3d_apply_ref(u,v,cv_wat)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_obc3d_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,u,obc_cv_wat_w,hx,dx_u,obc_cv_wat_e,dy_v,
   !&E                     hy,obc_cv_wat_s,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       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
   !&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,nv_adv,obc_cv_wat_e,           &
                            obc_cv_wat_s,obc_cv_wat_n,obc_cv_name
   USE obcogcm,      ONLY : obc_coefrel
#ifdef key_biolo
#if defined key_N_tracer
   USE bioloinit,     ONLY : nb_source_tracerN,p_marque_tracerN,iv_tracer_N,iv_signed_N,nb_source_river_tracerN
#endif
#if defined key_P_tracer
   USE bioloinit,     ONLY : nb_source_tracerP,p_marque_tracerP,iv_tracer_P,iv_signed_P,nb_source_river_tracerP
#endif
#endif

   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,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,ijgnul,iv
   REAL(KIND=rsh)            :: uu,vv,cou,dtrelax
   REAL(KIND=rlg)            :: dtm
#if defined key_N_tracer || defined key_P_tracer
   INTEGER                             :: is,ivtra,iv_tra,iv_sign
#endif

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

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

!$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)*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(1,(nv_adv)*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
   ! western boundary
       cv_wat(1:nv_adv,:,ig(j),j)=cv_wat(1:nv_adv,:,id(j)-1,j)
       cv_wat(1:nv_adv,:,ig(j)-1,j)=cv_wat(1:nv_adv,:,id(j)-2,j)
!!!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(1:nv_adv,:,ig(j)-2,j)=cv_wat(1:nv_adv,:,id(j)-3,j)
#endif
   ! eastern boundary
       cv_wat(1:nv_adv,:,id(j),j)=cv_wat(1:nv_adv,:,ig(j)+1,j)
!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(1:nv_adv,:,id(j)+1,j)=cv_wat(1:nv_adv,:,ig(j)+2,j)
       cv_wat(1:nv_adv,:,id(j)+2,j)=cv_wat(1:nv_adv,:,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*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(1,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
     DO i=limin,limax
   ! southern boundary
       cv_wat(1:nv_adv,:,i,jb(i))=cv_wat(1:nv_adv,:,i,jh(i)-1)
       cv_wat(1:nv_adv,:,i,jb(i)-1)=cv_wat(1:nv_adv,:,i,jh(i)-2)
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(1:nv_adv,:,i,jb(i)-2)=cv_wat(1:nv_adv,:,i,jh(i)-3)
#endif
   ! northern boundary
       cv_wat(1:nv_adv,:,i,jh(i))=cv_wat(1:nv_adv,:,i,jb(i)+1)
!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
       cv_wat(1:nv_adv,:,i,jh(i)+1)=cv_wat(1:nv_adv,:,i,jb(i)+2)
       cv_wat(1:nv_adv,:,i,jh(i)+2)=cv_wat(1:nv_adv,:,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
         DO k=1,kmax
           uu=u(k,i,j)
           IF (uu > 0.0_rsh) THEN
#if defined key_zerograd
             DO iv=1,nv_adv
               IF (obc_cv_name(iv)=='zerograd') THEN
                   ! zero gradient
                   ijgnul=1
                   dtrelax=0.0_rsh
               ELSE
                   dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                   IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                   ijgnul=0
               ENDIF
               cv_wat(iv,k,i,j)=cv_wat(iv,k,i+ijgnul,j)*(1.0_rsh-dtrelax) +obc_cv_wat_w(iv,k,1,j)*dtrelax
            ENDDO
#if defined key_N_tracer
            DO is=nb_source_river_tracerN+1,nb_source_tracerN
              DO ivtra=1,nb_var_tracerN
                 iv_tra=iv_tracer_N(ivtra,is)
                 iv_sign=iv_signed_N(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
              ENDDO
            ENDDO
#endif
#if defined key_P_tracer
            DO is=nb_source_river_tracerP+1,nb_source_tracerP
              DO ivtra=1,nb_var_tracerP
                 iv_tra=iv_tracer_P(ivtra,is)
                 iv_sign=iv_signed_P(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
              ENDDO
            ENDDO
#endif

#else

             !cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i+ijgnul,j)*(1.0_rsh-dtrelax)+obc_cv_wat_w(1:nv_adv,k,1,j)*dtrelax
             cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_w(1:nv_adv,k,1,j)*dtrelax
#endif
           ELSE
             cou=uu*dtm/dx_u(i,j)
             cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh+cou)-cou*cv_wat(1:nv_adv,k,i+1,j)
           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
           uu=u(k,i-1,j)
           IF (uu < 0.0_rsh) THEN
#if defined key_zerograd
             DO iv=1,nv_adv
               IF (obc_cv_name(iv)=='zerograd') THEN
                   ! zero gradient
                   ijgnul=1
                   dtrelax=0.0_rsh
               ELSE
                   dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                   IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                   ijgnul=0
               ENDIF
               cv_wat(iv,k,i,j)=cv_wat(iv,k,i-ijgnul,j)*(1.0_rsh-dtrelax) +obc_cv_wat_e(iv,k,imax,j)*dtrelax
            ENDDO
#if defined key_N_tracer
            DO is=nb_source_river_tracerN+1,nb_source_tracerN
              DO ivtra=1,nb_var_tracerN
                 iv_tra=iv_tracer_N(ivtra,is)
                 iv_sign=iv_signed_N(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
              ENDDO
            ENDDO
#endif
#if defined key_P_tracer
            DO is=nb_source_river_tracerP+1,nb_source_tracerP
              DO ivtra=1,nb_var_tracerP
                 iv_tra=iv_tracer_P(ivtra,is)
                 iv_sign=iv_signed_P(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
              ENDDO
            ENDDO
#endif

#else
             !cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i-ijgnul,j)*(1.0_rsh-dtrelax)+obc_cv_wat_e(1:nv_adv,k,imax,j)*dtrelax
             cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_e(1:nv_adv,k,imax,j)*dtrelax
#endif
           ELSE
             cou=uu*dtm/dx_u(i-1,j)
             cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-cou)+cou*cv_wat(1:nv_adv,k,i-1,j)
           END IF
         END DO
       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
           DO k=1,kmax
             vv=v(k,i,j)
             IF (vv > 0.0_rsh) THEN
#if defined key_zerograd
             DO iv=1,nv_adv
               IF (obc_cv_name(iv)=='zerograd') THEN
                   ! zero gradient
                   ijgnul=1
                   dtrelax=0.0_rsh
               ELSE
                   dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                   IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                   ijgnul=0
               ENDIF
               cv_wat(iv,k,i,j)=cv_wat(iv,k,i,j+ijgnul)*(1.0_rsh-dtrelax) +obc_cv_wat_s(iv,k,i,1)*dtrelax
            ENDDO
#if defined key_N_tracer
            DO is=nb_source_river_tracerN+1,nb_source_tracerN
              DO ivtra=1,nb_var_tracerN
                 iv_tra=iv_tracer_N(ivtra,is)
                 iv_sign=iv_signed_N(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
              ENDDO
            ENDDO
#endif
#if defined key_P_tracer
            DO is=nb_source_river_tracerP+1,nb_source_tracerP
              DO ivtra=1,nb_var_tracerP
                 iv_tra=iv_tracer_P(ivtra,is)
                 iv_sign=iv_signed_P(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
              ENDDO
            ENDDO
#endif

#else
               !cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j+ijgnul)*(1.0_rsh-dtrelax)+obc_cv_wat_s(1:nv_adv,k,i,1)*dtrelax
               cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_s(1:nv_adv,k,i,1)*dtrelax
#endif
             ELSE
               cou=vv*dtm/dy_v(i,j)
               cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh+cou)-cou*cv_wat(1:nv_adv,k,i,j+1)
             END IF
           END DO
         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
           DO k=1,kmax
             vv=v(k,i,j-1)
             IF (vv < 0.0_rsh) THEN
#if defined key_zerograd
             DO iv=1,nv_adv
               IF (obc_cv_name(iv)=='zerograd') THEN
                   ! zero gradient
                   ijgnul=1
                   dtrelax=0.0_rsh
               ELSE
                   dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                   IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                   ijgnul=0
               ENDIF
               cv_wat(iv,k,i,j)=cv_wat(iv,k,i,j-ijgnul)*(1.0_rsh-dtrelax)+obc_cv_wat_n(iv,k,i,jmax)*dtrelax
            ENDDO
#if defined key_N_tracer
            DO is=nb_source_river_tracerN+1,nb_source_tracerN
               DO ivtra=1,nb_var_tracerN
                 iv_tra=iv_tracer_N(ivtra,is)
                 iv_sign=iv_signed_N(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
               ENDDO
            ENDDO
#endif
#if defined key_P_tracer
            DO is=nb_source_river_tracerP+1,nb_source_tracerP
               DO ivtra=1,nb_var_tracerP
                 iv_tra=iv_tracer_P(ivtra,is)
                 iv_sign=iv_signed_P(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
               ENDDO
            ENDDO
#endif

#else
               !cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j-ijgnul)*(1.0_rsh-dtrelax)+obc_cv_wat_n(1:nv_adv,k,i,jmax)*dtrelax
               cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_n(1:nv_adv,k,i,jmax)*dtrelax
#endif
             ELSE
               cou=vv*dtm/dy_v(i,j-1)
               cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-cou)+cou*cv_wat(1:nv_adv,k,i,j-1)
             END IF
           END DO
         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_obc3d_apply_ref
#endif

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

#if defined key_lagoon
  SUBROUTINE sub_obc3d_apply_lagoon(u,v,cv_wat)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_obc3d_apply_lagoon  ***
   !&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,u,obc_cv_wat_w,hx,dx_u,obc_cv_wat_e,dy_v,
   !&E                     hy,obc_cv_wat_s,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       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
   !&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 comvars2d,    ONLY : igrau,jgrau,nbgrau,type_grau
   USE comsubstance, ONLY : obc_cv_wat_w,nv_adv,obc_cv_wat_e,           &
                            obc_cv_wat_s,obc_cv_wat_n,obc_cv_name
   USE comsubstance, ONLY : obc_cv_wat_grau
   USE obcogcm,      ONLY : obc_coefrel
#ifdef key_biolo
#if defined key_N_tracer
   USE bioloinit,     ONLY : nb_source_tracerN,p_marque_tracerN,iv_tracer_N,iv_signed_N,nb_source_river_tracerN
#endif
#if defined key_P_tracer
   USE bioloinit,     ONLY : nb_source_tracerP,p_marque_tracerP,iv_tracer_P,iv_signed_P,nb_source_river_tracerP
#endif
#endif

   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,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,ijgnul,iv,ng
   REAL(KIND=rsh)            :: uu,vv,cou,dtrelax
   REAL(KIND=rlg)            :: dtm
#if defined key_N_tracer || defined key_P_tracer
   INTEGER                             :: is,ivtra,iv_tra,iv_sign
#endif

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

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

! for lagoons : l_obc_cyclx ou y =.False. always

   DO ng=1,nbgrau
      i=igrau(ng)
      j=jgrau(ng)

     ! western boundary
      IF (type_grau(ng)==4) THEN
      !! grau sur la frontiere Ouest
        IF_MPI ( (i >= limin) .AND. (i <= limax) .AND.(j >= ljmin) .AND. (j <= ljmax) ) THEN
          DO k=1,kmax
            uu=u(k,i,j)
            IF (uu > 0.0_rsh) THEN
#if defined key_zerograd
              DO iv=1,nv_adv
                IF (obc_cv_name(iv)=='zerograd') THEN
                  ! zero gradient
                  ijgnul=1
                  dtrelax=0.0_rsh
                ELSE
                  dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                  IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                  ijgnul=0
                ENDIF
                cv_wat(iv,k,i,j)=cv_wat(iv,k,i+ijgnul,j)*(1.0_rsh-dtrelax) +obc_cv_wat_grau(iv,ng)*dtrelax
              ENDDO
#if defined key_N_tracer
              DO is=nb_source_river_tracerN+1,nb_source_tracerN
                DO ivtra=1,nb_var_tracerN
                   iv_tra=iv_tracer_N(ivtra,is)
                   iv_sign=iv_signed_N(ivtra,is)
                   cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
                ENDDO
              ENDDO
#endif
#if defined key_P_tracer
              DO is=nb_source_river_tracerP+1,nb_source_tracerP
                DO ivtra=1,nb_var_tracerP
                  iv_tra=iv_tracer_P(ivtra,is)
                  iv_sign=iv_signed_P(ivtra,is)
                  cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
                ENDDO
              ENDDO
#endif

#else
!  no key_zerograd
              cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_grau(1:nv_adv,ng)*dtrelax
#endif
            ELSE
              cou=uu*dtm/dx_u(i,j)
              cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh+cou)-cou*cv_wat(1:nv_adv,k,i+1,j)
            END IF
          END DO
        ENDIF_MPI
      END IF

    ! eastern boundary
      IF (type_grau(ng)==2) THEN
        !! grau sur la frontiere Est
        IF_MPI ( (i >= limin) .AND. (i <= limax) .AND.(j >= ljmin) .AND. (j <= ljmax) ) THEN

          DO k=1,kmax
            uu=u(k,i-1,j)
            IF (uu < 0.0_rsh) THEN
#if defined key_zerograd
              DO iv=1,nv_adv
                IF (obc_cv_name(iv)=='zerograd') THEN
                  ! zero gradient
                  ijgnul=1
                  dtrelax=0.0_rsh
                ELSE
                  dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                  IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                  ijgnul=0
                ENDIF
                cv_wat(iv,k,i,j)=cv_wat(iv,k,i-ijgnul,j)*(1.0_rsh-dtrelax) +obc_cv_wat_grau(iv,ng)*dtrelax
              ENDDO
#if defined key_N_tracer
              DO is=nb_source_river_tracerN+1,nb_source_tracerN
               DO ivtra=1,nb_var_tracerN
                iv_tra=iv_tracer_N(ivtra,is)
                iv_sign=iv_signed_N(ivtra,is)
                cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
               ENDDO
              ENDDO
#endif
#if defined key_P_tracer
              DO is=nb_source_river_tracerP+1,nb_source_tracerP
               DO ivtra=1,nb_var_tracerP
                 iv_tra=iv_tracer_P(ivtra,is)
                 iv_sign=iv_signed_P(ivtra,is)
                 cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
               ENDDO
              ENDDO
#endif

#else
              cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_grau(1:nv_adv,ng)*dtrelax
#endif
            ELSE
              cou=uu*dtm/dx_u(i-1,j)
              cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-cou)+cou*cv_wat(1:nv_adv,k,i-1,j)
            END IF
          END DO
        ENDIF_MPI
      END IF

     ! southern boundary
     IF (type_grau(ng)==3) THEN
      !! grau sur la frontiere Sud
        IF_MPI ( (i >= limin) .AND. (i <= limax) .AND.(j >= ljmin) .AND. (j <= ljmax) ) THEN
           DO k=1,kmax
             vv=v(k,i,j)
             IF (vv > 0.0_rsh) THEN
#if defined key_zerograd
               DO iv=1,nv_adv
                 IF (obc_cv_name(iv)=='zerograd') THEN
                   ! zero gradient
                   ijgnul=1
                   dtrelax=0.0_rsh
                 ELSE
                   dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                   IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                   ijgnul=0
                 ENDIF
                 cv_wat(iv,k,i,j)=cv_wat(iv,k,i,j+ijgnul)*(1.0_rsh-dtrelax) +obc_cv_wat_grau(iv,ng)*dtrelax
               ENDDO
#if defined key_N_tracer
               DO is=nb_source_river_tracerN+1,nb_source_tracerN
                DO ivtra=1,nb_var_tracerN
                  iv_tra=iv_tracer_N(ivtra,is)
                  iv_sign=iv_signed_N(ivtra,is)
                  cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
                ENDDO
               ENDDO
#endif
#if defined key_P_tracer
               DO is=nb_source_river_tracerP+1,nb_source_tracerP
                DO ivtra=1,nb_var_tracerP
                  iv_tra=iv_tracer_P(ivtra,is)
                  iv_sign=iv_signed_P(ivtra,is)
                  cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
                ENDDO
               ENDDO
#endif

#else
               cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_grau(1:nv_adv,ng)*dtrelax
#endif
             ELSE
               cou=vv*dtm/dy_v(i,j)
               cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh+cou)-cou*cv_wat(1:nv_adv,k,i,j+1)
             END IF
           END DO
        ENDIF_MPI
     END IF

    ! northern boundary
     IF (type_grau(ng)==1) THEN
       !! grau sur la frontiere Nord
       IF_MPI ( (i >= limin) .AND. (i <= limax) .AND.(j >= ljmin) .AND. (j <= ljmax) ) THEN
          DO k=1,kmax
            vv=v(k,i,j-1)
            IF (vv < 0.0_rsh) THEN
#if defined key_zerograd
              DO iv=1,nv_adv
                IF (obc_cv_name(iv)=='zerograd') THEN
                  ! zero gradient
                  ijgnul=1
                  dtrelax=0.0_rsh
                ELSE
                  dtrelax=dtm*obc_coefrel  ! dt / time scale of relaxation
                  IF (obc_coefrel==0.0_rsh) dtrelax=1.0_rsh
                  ijgnul=0
                ENDIF
                cv_wat(iv,k,i,j)=cv_wat(iv,k,i,j-ijgnul)*(1.0_rsh-dtrelax)+obc_cv_wat_grau(iv,ng)*dtrelax
              ENDDO
#if defined key_N_tracer
              DO is=nb_source_river_tracerN+1,nb_source_tracerN
                DO ivtra=1,nb_var_tracerN
                  iv_tra=iv_tracer_N(ivtra,is)
                  iv_sign=iv_signed_N(ivtra,is)
                  cv_wat(iv_tra,k,i,j)=p_marque_tracerN(ivtra,is)*cv_wat(iv_sign,k,i,j)
                ENDDO
              ENDDO
#endif
#if defined key_P_tracer
              DO is=nb_source_river_tracerP+1,nb_source_tracerP
                DO ivtra=1,nb_var_tracerP
                  iv_tra=iv_tracer_P(ivtra,is)
                  iv_sign=iv_signed_P(ivtra,is)
                  cv_wat(iv_tra,k,i,j)=p_marque_tracerP(ivtra,is)*cv_wat(iv_sign,k,i,j)
                ENDDO
              ENDDO
#endif

#else
              cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-dtrelax)+obc_cv_wat_grau(1:nv_adv,ng)*dtrelax
#endif
            ELSE
              cou=vv*dtm/dy_v(i,j-1)
              cv_wat(1:nv_adv,k,i,j)=cv_wat(1:nv_adv,k,i,j)*(1.0_rsh-cou)+cou*cv_wat(1:nv_adv,k,i,j-1)
            END IF
          END DO
        ENDIF_MPI
      END IF

   END DO   ! ends nb_grau 

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_obc3d_apply_lagoon
#endif
   !!======================================================================

#endif
 END MODULE subobc3dapply



