!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#if defined SOLVE3D || defined WET_DRY

      subroutine set_depth (tile)
      implicit none
# include "param.h"
      integer tile, trd
C$    integer omp_get_thread_num
# include "compute_tile_bounds.h"
      call set_depth_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine set_depth_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
! Create S-coordinate system: based on model topography h(i,j),
! fast-time-averaged free-surface field and vertical coordinate
! transformation metrics compute evolving depths of of the three-
! dimensional model grid (z_r,z_w) and vertical heights of model
! grid boxes.
!----------------------------------------------------------------------
!
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k
      real    cff_r,cff1_r,cff_w,cff1_w, z_r0,z_w0, zetatmp,htmp
# include "param.h"
# include "grid.h"
# ifdef AGRIF
#  include "zoom.h"
# endif
# include "ocean2d.h"
# include "ocean3d.h"
# include "coupling.h"
# include "scalars.h"
# ifdef MRL_WCI
#  include "forces.h"
# endif
# ifdef M3FAST
#  include "nbq.h"
# endif
      real eps
      parameter (eps=1.e-14)
!
# include "compute_auxiliary_bounds.h"
!
# ifdef EW_PERIODIC
#  define I_RANGE Istr,Iend
# else
#  define I_RANGE IstrR,IendR
# endif
# ifdef NS_PERIODIC
#  define J_RANGE Jstr,Jend
# else
#  define J_RANGE JstrR,JendR
# endif

!
!----------------------------------------------------------------------
!  Initialization of free surface zeta arrays and wet/dry mask
!
!  During initialization and/or restart: copy initial free surface
!  field into array for holding fast-time averaged free surface.
!----------------------------------------------------------------------
!
      if (iic.eq.0) then
!$acc kernels if(compute_on_device) default(present)
        do j=J_RANGE
          do i=I_RANGE
# ifdef WET_DRY
            wetdry(i,j)=1.
#  ifdef AGRIF
            rmask_childs(i,j)=1.
#  endif
            IF (h(i,j).eq.0.) THEN
              h(i,j)=eps
            END IF
            IF (zeta(i,j,1) .lt. Dcrit(i,j)-h(i,j)) THEN
              zeta(i,j,1)=Dcrit(i,j)-h(i,j)
              zeta(i,j,2)=zeta(i,j,1)
              zeta(i,j,3)=zeta(i,j,1)
              zeta(i,j,4)=zeta(i,j,1)
              wetdry(i,j)=0.
            ENDIF
            rmask_wet(i,j)=wetdry(i,j)
# endif
# ifdef SOLVE3D
#  ifdef NEW_S_COORD
            hinv(i,j)=1./(abs(h(i,j))+hc)
#  else
            hinv(i,j)=1./h(i,j)
#  endif
            Zt_avg1(i,j)=zeta(i,j,knew)  ! knew=1 unless modified
                                         ! in ana_initial (JET)
# endif /* SOLVE3D */
          enddo
        enddo

# ifdef NBQ_MASS
        do k=1,N
          do j=J_RANGE
            do i=I_RANGE
              rho_nbq(i,j,k)=0.
            enddo
          enddo
        enddo
# endif
!$acc end kernels
# if defined BD_EXCHANGE
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend, hinv)
#  ifdef WET_DRY
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend, h)
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend, wetdry)
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend, rmask_wet)
#   ifdef AGRIF
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend, rmask_childs)
#   endif
#  endif
#  if defined M3FAST || defined WET_DRY
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          zeta(START_2D_ARRAY,1))
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          zeta(START_2D_ARRAY,2))
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          zeta(START_2D_ARRAY,3))
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          zeta(START_2D_ARRAY,4))
#  endif
# endif
      endif   !  iic==0  <-- end initialization
!
!----------------------------------------------------------------------
!  Set vertical grid
!----------------------------------------------------------------------
!
# ifdef SOLVE3D
#  if defined BD_EXCHANGE
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, Zt_avg1)
#  endif
!$acc kernels if(compute_on_device) default(present)

      do j=J_RANGE              !!! WARNING: Setting must be consistent
        do i=I_RANGE            !!!          with omega.F
          z_w(i,j,0)=-h(i,j)
        enddo

        do k=1,N,+1
#  ifdef NEW_S_COORD
          cff_w =hc*sc_w(k)
          cff_r =hc*sc_r(k)
          cff1_w=Cs_w(k)
          cff1_r=Cs_r(k)
#  else
          cff_w =hc*(sc_w(k)-Cs_w(k))
          cff_r =hc*(sc_r(k)-Cs_r(k))
          cff1_w=Cs_w(k)
          cff1_r=Cs_r(k)
#  endif
          do i=I_RANGE
#  ifdef M3FAST
            zetatmp=zeta(i,j,knew)
#  else
            zetatmp=Zt_avg1(i,j)
#  endif
#  if defined MASKING && !defined WET_DRY
            zetatmp=zetatmp*rmask(i,j)
#  endif
#  ifdef NEW_S_COORD
            htmp=abs(h(i,j))
#  else
            htmp=h(i,j)
#  endif
            z_w0=cff_w+cff1_w*htmp
            z_r0=cff_r+cff1_r*htmp
#  ifdef NEW_S_COORD
            z_w(i,j,k)=z_w0*h(i,j)*hinv(i,j)+zetatmp
     &                                       *(1.+z_w0*hinv(i,j))
            z_r(i,j,k)=z_r0*h(i,j)*hinv(i,j)+zetatmp
     &                                       *(1.+z_r0*hinv(i,j))
#  else
            z_w(i,j,k)=z_w0+zetatmp*(1.+z_w0*hinv(i,j))
            z_r(i,j,k)=z_r0+zetatmp*(1.+z_r0*hinv(i,j))
#  endif
          enddo
        enddo
      enddo

      do k=1,N
        do j=J_RANGE
          do i=I_RANGE
#  ifdef NBQ
#   ifdef NBQ_GRID_SLOW
            if (LAST_FAST_STEP)  Hz_bak(i,j,k)=Hz(i,j,k)
#   else
            if (FIRST_FAST_STEP) Hz_bak(i,j,k)=Hz(i,j,k)
#   endif
#   ifdef NBQ_HZ_PROGNOSTIC
            Hz_bak2(i,j,k)=Hz(i,j,k)
#   endif
#  else
            Hz_bak(i,j,k)=Hz(i,j,k)
#  endif
            Hz(i,j,k)=z_w(i,j,k)-z_w(i,j,k-1)

#  ifdef NBQ_MASS
            Hzr(i,j,k)=Hz(i,j,k)
            Hz(i,j,k) =Hz(i,j,k)*(1.+rho_grd(i,j,k))
     &                              +rho_nbq(i,j,k)
#  endif
          enddo
        enddo
      enddo

#  if defined M3FAST && ! defined EXACT_RESTART
      if (iic.eq.0 .or. (FIRST_TIME_STEP .and. iif.eq.1)) then
        do k=1,N,+1
          do j=J_RANGE
            do i=I_RANGE
              Hz_bak(i,j,k)=Hz(i,j,k)
            enddo
          enddo
        enddo
      endif
#  endif
!$acc end kernels

#  undef I_RANGE
#  undef J_RANGE

#  if defined BD_EXCHANGE
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend, 
     &                        z_w(START_2D_ARRAY,0))
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        z_r(START_2D_ARRAY,1))
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hz(START_2D_ARRAY,1))
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hz_bak(START_2D_ARRAY,1))
#   ifdef NBQ_MASS
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hzr(START_2D_ARRAY,1))
#   endif
#  endif
      return
      end

#  ifdef MORPHODYN
!======================================================================
!
!                    subroutine set_depth_movbat
!
!======================================================================
!

      subroutine set_depth_morphodyn (tile)
      implicit none
#   include "param.h"
      integer tile, trd
C$    integer omp_get_thread_num
#   include "compute_tile_bounds.h"
      call set_depth_morphodyn_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine set_depth_morphodyn_tile (Istr,Iend,Jstr,Jend)

      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k
      real    cff_r,cff1_r,cff_w,cff1_w, z_r0,z_w0, zetatmp,htmp
#   include "param.h"
#   include "grid.h"
#   ifdef AGRIF
#    include "zoom.h"
#   endif
#   include "ocean2d.h"
#   include "ocean3d.h"
#   include "coupling.h"
#   include "scalars.h"
#   ifdef MRL_WCI
#    include "forces.h"
#   endif
#   ifdef M3FAST
#    include "nbq.h"
#   endif
      real eps
      parameter (eps=1.e-14)
!
#   include "compute_auxiliary_bounds.h"
!
#   ifdef EW_PERIODIC
#    define I_RANGE Istr,Iend
#   else
#    define I_RANGE IstrR,IendR
#   endif
#   ifdef NS_PERIODIC
#    define J_RANGE Jstr,Jend
#   else
#    define J_RANGE JstrR,JendR
#   endif

      do j=J_RANGE              !!! WARNING: Setting must be consistent
        do i=I_RANGE            !!!          with omega.F
          h(i,j)=h(i,j)+dh(i,j)
          z_w(i,j,0)=-h(i,j)
        enddo

        do k=1,N,+1
#   ifdef NEW_S_COORD
          cff_w =hc*sc_w(k)
          cff_r =hc*sc_r(k)
          cff1_w=Cs_w(k)
          cff1_r=Cs_r(k)
#   else
          cff_w =hc*(sc_w(k)-Cs_w(k))
          cff_r =hc*(sc_r(k)-Cs_r(k))
          cff1_w=Cs_w(k)
          cff1_r=Cs_r(k)
#   endif
          do i=I_RANGE
#   ifdef M3FAST
            zetatmp=zeta(i,j,knew)
#   else
            zetatmp=Zt_avg1(i,j)
#   endif
#  ifdef NEW_S_COORD
            htmp=abs(h(i,j))
            hinv(i,j)=1./(htmp+hc)
#  else
            htmp=h(i,j)
            hinv(i,j)=1./htmp
#  endif
#   if defined MASKING && !defined WET_DRY
            zetatmp=zetatmp*rmask(i,j)
#   endif
            z_w0=cff_w+cff1_w*htmp
            z_r0=cff_r+cff1_r*htmp
#   ifdef NEW_S_COORD
            z_w(i,j,k)=z_w0*h(i,j)*hinv(i,j)+zetatmp
     &                                       *(1.+z_w0*hinv(i,j))
            z_r(i,j,k)=z_r0*h(i,j)*hinv(i,j)+zetatmp
     &                                       *(1.+z_r0*hinv(i,j))
#   else
            z_w(i,j,k)=z_w0+zetatmp*(1.+z_w0*hinv(i,j))
            z_r(i,j,k)=z_r0+zetatmp*(1.+z_r0*hinv(i,j))
#   endif
            Hz(i,j,k)=z_w(i,j,k)-z_w(i,j,k-1)
#   ifdef NBQ_MASS
            Hzr(i,j,k)=Hz(i,j,k)
            Hz(i,j,k) =Hz(i,j,k)*(1.+rho_grd(i,j,k))
     &                              +rho_nbq(i,j,k)
#   endif
          enddo
        enddo
      enddo

#   undef I_RANGE
#   undef J_RANGE

#   if defined BD_EXCHANGE
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend, 
     &                        z_w(START_2D_ARRAY,0))
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        z_r(START_2D_ARRAY,1))
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hz(START_2D_ARRAY,1))
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, h)
#    ifdef NBQ_MASS
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hzr(START_2D_ARRAY,1))
#    endif
#   endif
      return
      end
#  endif /* MORPHODYN */
!
!======================================================================
!
!                       subroutine set_HUV
!
!======================================================================
!
      subroutine set_HUV (tile)
      implicit none
#  include "param.h"
      integer tile, trd
C$    integer omp_get_thread_num
#  include "compute_tile_bounds.h"
      call set_HUV_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine set_HUV_tile (Istr,Iend,Jstr,Jend)
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k
#  include "param.h"
#  include "grid.h"
#  include "ocean3d.h"
#  include "scalars.h"
!
#  ifdef MRL_WCI
#   include "forces.h"
#  endif
#  include "compute_auxiliary_bounds.h"
!
#  ifdef EW_PERIODIC
#   define IU_RANGE Istr,Iend
#   define IV_RANGE Istr,Iend
#  else
#   define IU_RANGE Istr,IendR
#   define IV_RANGE IstrR,IendR
#  endif

#  ifdef NS_PERIODIC
#   define JU_RANGE Jstr,Jend
#   define JV_RANGE Jstr,Jend
#  else
#   define JU_RANGE JstrR,JendR
#   define JV_RANGE Jstr,JendR
#  endif
!$acc kernels if(compute_on_device) default(present)
      do k=1,N
        do j=JU_RANGE
          do i=IU_RANGE
            Huon(i,j,k)=0.5*(Hz(i,j,k)+Hz(i-1,j,k))*on_u(i,j)
     &                                       *( u(i,j,k,nrhs)
#  ifdef MRL_WCI
     &                                           + ust(i,j,k)
#  endif
     &                                                      )
#  if defined MRL_WCI && defined MASKING
            Huon(i,j,k)=Huon(i,j,k)*umask(i,j)
#  endif

#  if defined UV_VIS4 && defined UV_MIX_GEO
            z_u(i,j,k)=0.5*(z_r(i,j,k)+z_r(i-1,j,k))
            dz_u(i,j,k)=z_r(i,j,k)-z_r(i-1,j,k)
#  endif
          enddo
        enddo
        do j=JV_RANGE
          do i=IV_RANGE
            Hvom(i,j,k)=0.5*(Hz(i,j,k)+Hz(i,j-1,k))*om_v(i,j)
     &                                       *( v(i,j,k,nrhs)
#  ifdef MRL_WCI
     &                                           + vst(i,j,k)
#  endif
     &                                                      )
#  if defined MRL_WCI && defined MASKING
            Hvom(i,j,k)=Hvom(i,j,k)*vmask(i,j)
#  endif
#  if defined UV_VIS4 && defined UV_MIX_GEO
            z_v(i,j,k)=0.5*(z_r(i,j,k)+z_r(i,j-1,k))
            dz_v(i,j,k)=z_r(i,j,k)-z_r(i,j-1,k)
#  endif
          enddo
        enddo
      enddo
!$acc end kernels
#  undef IU_RANGE
#  undef JU_RANGE
#  undef IV_RANGE
#  undef JV_RANGE
!
! Exchange periodic boundaries, if so prescribed.
!
#  if defined BD_EXCHANGE
#   ifdef THREE_GHOST_POINTS_UV
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             Huon(START_2D_ARRAY,1))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             Hvom(START_2D_ARRAY,1))
#   else
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        Huon(START_2D_ARRAY,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hvom(START_2D_ARRAY,1))
#   endif
#   if defined UV_VIS4 && defined UV_MIX_GEO
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        z_u(START_2D_ARRAY,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        z_v(START_2D_ARRAY,1))
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        dz_u(START_2D_ARRAY,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        dz_v(START_2D_ARRAY,1))
#   endif
#  endif
      return
      end
!
!======================================================================
!
!                       subroutine set_HUV1
!
!======================================================================
!
      subroutine set_HUV1 (tile)
      implicit none
#  include "param.h"
      integer tile, trd
C$    integer omp_get_thread_num
#  include "compute_tile_bounds.h"
      call set_HUV1_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine set_HUV1_tile (Istr,Iend,Jstr,Jend)
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k
#  include "param.h"
#  include "grid.h"
#  include "ocean3d.h"
#  include "scalars.h"
#  ifdef MRL_WCI
#   include "forces.h"
#  endif
!
#  include "compute_auxiliary_bounds.h"
!
#  ifdef EW_PERIODIC
#   define IU_RANGE Istr,Iend
#   define IV_RANGE Istr,Iend
#  else
#   define IU_RANGE Istr,IendR
#   define IV_RANGE IstrR,IendR
#  endif

#  ifdef NS_PERIODIC
#   define JU_RANGE Jstr,Jend
#   define JV_RANGE Jstr,Jend
#  else
#   define JU_RANGE JstrR,JendR
#   define JV_RANGE Jstr,JendR
#  endif

      do k=1,N
        do j=JU_RANGE
          do i=IU_RANGE
            Huon(i,j,k)=0.25*(3.*(Hz(i,j,k)+Hz(i-1,j,k))
     &                   -Hz_bak(i,j,k)-Hz_bak(i-1,j,k))
     &                        *on_u(i,j)*( u(i,j,k,nrhs)
#  ifdef MRL_WCI
     &                                      + ust(i,j,k)
#  endif
     &                                                 )
#  if defined MRL_WCI && defined MASKING
            Huon(i,j,k)=Huon(i,j,k)*umask(i,j)
#  endif

          enddo
        enddo
        do j=JV_RANGE
          do i=IV_RANGE
            Hvom(i,j,k)=0.25*( 3.*(Hz(i,j,k)+Hz(i,j-1,k))
     &                    -Hz_bak(i,j,k)-Hz_bak(i,j-1,k))
     &                         *om_v(i,j)*( v(i,j,k,nrhs)
#  ifdef MRL_WCI
     &                                       + vst(i,j,k)
#  endif
     &                                                  )
#  if defined MRL_WCI && defined MASKING
            Hvom(i,j,k)=Hvom(i,j,k)*vmask(i,j)
#  endif
          enddo
        enddo
      enddo
#  undef IU_RANGE
#  undef JU_RANGE
#  undef IV_RANGE
#  undef JV_RANGE
!
! Exchange periodic boundaries, if so prescribed.
!
#  if defined BD_EXCHANGE
#   ifdef THREE_GHOST_POINTS_UV
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             Huon(START_2D_ARRAY,1))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             Hvom(START_2D_ARRAY,1))
#   else
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        Huon(START_2D_ARRAY,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hvom(START_2D_ARRAY,1))
#   endif
#  endif
      return
      end
!
!======================================================================
!
!                       subroutine set_HUV2
!
!======================================================================
!
      subroutine set_HUV2 (tile)
      implicit none
      integer tile, trd, omp_get_thread_num
#  include "param.h"
#  include "private_scratch.h"
#  include "compute_tile_bounds.h"
      trd=omp_get_thread_num()
      call set_HUV2_tile (Istr,Iend,Jstr,Jend, A2d(1,1,trd),
     &                                         A2d(1,2,trd))
      return
      end

      subroutine set_HUV2_tile (Istr,Iend,Jstr,Jend,DC,FC)
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k
#  include "param.h"
      real DC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     FC(PRIVATE_1D_SCRATCH_ARRAY,0:N)
#  include "grid.h"
#  include "ocean3d.h"
#  include "scalars.h"
#  ifdef MRL_WCI
#   include "forces.h"
#  endif
#  include "coupling.h"
#ifdef OPENACC
#  include "openacc_work.h"
#endif
!
#  include "compute_auxiliary_bounds.h"
!
#  ifdef EW_PERIODIC
#   define IU_RANGE Istr,Iend
#   define IV_RANGE Istr,Iend
#  else
#   define IU_RANGE Istr,IendR
#   define IV_RANGE IstrR,IendR
#  endif

#  ifdef NS_PERIODIC
#   define JU_RANGE Jstr,Jend
#   define JV_RANGE Jstr,Jend
#  else
#   define JU_RANGE JstrR,JendR
#   define JV_RANGE Jstr,JendR
#  endif

!$acc kernels if(compute_on_device) default(present)

 !     do j=JU_RANGE
        DOLOOP2D_R(IU_RANGE,JU_RANGE)
        do i=IU_RANGE
          DC(i,0)=0.
          FC(i,0)=0.
        enddo
        do k=1,N
          do i=IU_RANGE
            DC(i,k)=0.5*(Hz(i,j,k)+Hz(i-1,j,k))*on_u(i,j)
            DC(i,0)=DC(i,0)+DC(i,k)
            FC(i,0)=FC(i,0)+DC(i,k)*u(i,j,k,nrhs)
          enddo
        enddo
        do i=IU_RANGE
          FC(i,0)=(FC(i,0)-DU_avg2(i,j))/DC(i,0)
#  ifdef MRL_WCI
     &                               +ust2d(i,j)
#  endif
        enddo
        do k=1,N
          do i=IU_RANGE
            u(i,j,k,nrhs)=(u(i,j,k,nrhs)-FC(i,0))
#  ifdef MASKING
     &                                *umask(i,j)
#   ifdef MRL_WCI
     &               +ust(i,j,k)*(umask(i,j)-1.0)
#   endif
#  endif
            Huon(i,j,k)=DC(i,k)*( u(i,j,k,nrhs)
#  ifdef MRL_WCI
     &                             + ust(i,j,k)
#  endif
     &                                        )
          enddo
        enddo
!      enddo
        ENDDOLOOP2D


!      do j=JV_RANGE
        DOLOOP2D_R(IV_RANGE,JV_RANGE)
        do i=IV_RANGE
          DC(i,0)=0.
          FC(i,0)=0.
        enddo
        do k=1,N
          do i=IV_RANGE
            DC(i,k)=0.5*(Hz(i,j,k)+Hz(i,j-1,k))*om_v(i,j)
            DC(i,0)=DC(i,0)+DC(i,k)
            FC(i,0)=FC(i,0)+DC(i,k)*v(i,j,k,nrhs)
          enddo
        enddo
        do i=IV_RANGE
          FC(i,0)=(FC(i,0)-DV_avg2(i,j))/DC(i,0)
#  ifdef MRL_WCI
     &                               +vst2d(i,j)
#  endif
        enddo
        do k=1,N
          do i=IV_RANGE
            v(i,j,k,nrhs)=(v(i,j,k,nrhs)-FC(i,0))
#  ifdef MASKING
     &                                *vmask(i,j)
#   ifdef MRL_WCI
     &               +vst(i,j,k)*(vmask(i,j)-1.0)
#   endif
#  endif
            Hvom(i,j,k)=DC(i,k)*( v(i,j,k,nrhs)
#  ifdef MRL_WCI
     &                             + vst(i,j,k)
#  endif
     &                                        )
          enddo
        enddo
!     enddo
        ENDDOLOOP2D

!$acc end kernels

#  undef IU_RANGE
#  undef JU_RANGE
#  undef IV_RANGE
#  undef JV_RANGE
!
! Exchange periodic boundaries, if so prescribed.
!
#  if defined BD_EXCHANGE
#   ifdef THREE_GHOST_POINTS_UV  
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             Huon(START_2D_ARRAY,1))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             Hvom(START_2D_ARRAY,1))
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             u(START_2D_ARRAY,1,nrhs))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             v(START_2D_ARRAY,1,nrhs))
#   else
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        Huon(START_2D_ARRAY,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        Hvom(START_2D_ARRAY,1))
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        u(START_2D_ARRAY,1,nrhs))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        v(START_2D_ARRAY,1,nrhs))
#   endif
#  endif
# endif /* SOLVE3D */
      return
      end
#else
      subroutine set_depth_empty
      return
      end
#endif /* SOLVE3D || WET_DRY */

