 MODULE subobc3dapply



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





























































!! * Modules used

   IMPLICIT NONE
   PRIVATE

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

!! * Interface
   INTERFACE sub_obc3d_apply

     MODULE PROCEDURE sub_obc3d_apply_ref

   END INTERFACE

!! * Shared module variables

!! * Private variables

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


  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


   USE 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


!!----------------------------------------------------------------------
!! * 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

!!!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder

     CALL ex_peri_x_rsh(1,(nv_adv)*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))


   END IF

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

!!!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder

     CALL ex_peri_y_rsh(1,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))


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


!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

           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

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

!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

           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

     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

!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

             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

!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

             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

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(cv_wat)

  END SUBROUTINE sub_obc3d_apply_ref


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


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


 END MODULE subobc3dapply



