!======================================================================
! 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"

#ifdef M3FAST

      subroutine step3d_fast (tile)
      implicit none
      integer tile, trd
# include "param.h"
# include "private_scratch.h"
C$    integer omp_get_thread_num
# include "compute_tile_bounds.h"
      trd=0
C$    trd=omp_get_thread_num()
!      
!======================================================================
!                       step3d_fast
!======================================================================
!
      call step3d_fast_tile ( Istr,Iend,Jstr,Jend
     &                   , A2d(1, 1,trd), A2d(1, 2,trd), A2d(1, 3,trd)
     &                   , A2d(1, 4,trd), A2d(1, 5,trd), A2d(1, 6,trd)
     &                   , A2d(1, 7,trd), A2d(1, 8,trd)
     &                   , A2d(1, 9,trd), A2d(1,10,trd)
     &                   , A2d(1,11,trd), A2d(1,12,trd)
     &                   , A2d(1,13,trd), A2d(1,14,trd)
     &                   , A2d(1,15,trd), A2d(1,16,trd), A2d(1,17,trd)
# ifdef NBQ
! Next 3 arrays have a last dim of 2
     &                   , A2d(1,18,trd), A2d(1,20,trd), A2d(1,22,trd)
     &                   , A2d(1,24,trd), A2d(1,25,trd), A2d(1,26,trd)
     &                   , A2d(1,27,trd), A2d(1,28,trd), A2d(1,29,trd)
     &                   , A3d(1, 1,trd), A3d(1, 2,trd)
     &                   , A3d(1, 3,trd), A3d(1, 4,trd)
     &                   , A3d(1, 5,trd), A3d(1, 6,trd)
     &                   , A3d(1, 7,trd), A3d(1, 8,trd) 
# endif
     &                   )

      return
      end
! 
!======================================================================
!                      step3d_fast_tile
!======================================================================
!
      subroutine step3d_fast_tile (Istr,Iend,Jstr,Jend
     &                          ,Dnew,rubar,rvbar
     &                          ,Drhs, UFx,UFe
     &                          ,VFx,VFe
     &                          ,urhs,vrhs 
     &                          ,DUon,DVom
     &                          ,ru_ext_nbq_sum, rv_ext_nbq_sum
     &                          ,ru_ext_nbq_old, rv_ext_nbq_old,work
# ifdef NBQ
     &                          ,dthetadiv_nbqdz_w
     &                          ,dZdxq_u,dZdyq_v
     &                          ,FX,FY,FC,DC,CF,BC
     &                          ,Hzw_half_nbq_inv,  Hzr_half_nbq_inv 
     &                          ,Hzw_half_nbq_inv_u,Hzw_half_nbq_inv_v
     &                          ,Hzu_half_qdmu,     Hzv_half_qdmv
     &                          ,FX_3D, FY_3D
# endif
     &                          )
!
!
!***********************************************************************
!
!                    SOLVE FAST MODE 3D EQUATIONS
!
!***********************************************************************
!
!  This routines:
!  1- Computes non-NBQ RHS forcing terms of momentum equations. First
!     computes the barotropic (external) RHS forcing term (rubar,rvbar)
!     then adds it to the internal RHS forcing (computed in pre_step3d). 
!  2- Solves the 3D momentum conservation equations for fast-mode 
!     components (qdmu_nbq, qdmv_nbq, qdmw_nbq) by time integration of 
!     all forces:
!       Compressible pressure force + second viscosity + gravity
!       + NT Coriolis force + restoring force + non-NBQ RHS forces
!  3- Solves mass conservation equation, i.e., computes compressible 
!     density rho_nbq by time integration of momentum divergence
!
!  In this version, a first guest of zeta is derived from the surface 
!  vertical velocity (surface characteristic relation) instead of the 
!  depth-averaged conservation of mass. This satisfies dynamical coupling 
!  with the surface layer. After solving the 3D momentum equations, a 
!  final zeta field is diagnozed from mass conservation (then Hz is also
!  corrected for the internal time step).
!
!  W-momentum equation is solved with explicit or implicit methods:
!  - Explicit scheme: w-momentum is updated right after (and the same 
!                     way as) u- and v-momentum.
!  - Implicit scheme: horizontal component of divergence is first 
!                     precomputed (as required by fast-mode mass 
!                     conservation) before tridiagonal Gauss Elimination 
!                     is carried out for qdmw_nbq(m).
!
!  For all components, a Forward-Backward scheme is implemented:
!  - Explicit scheme: Forward:  zeta, qdmu_nbq, qdmw_nbq. 
!                     Backward:       rho_nbq.
!  - Implicit scheme: Forward:  zeta, qdmu_nbq. 
!                     Backward:       qdmw_nbq, rho_nbq.
!
!  In the NBQ_PERF option, the vertical grid is not evolving at fast
!  time step to gain computational time.
!
!***********************************************************************
!
# ifdef CVTK_DEBUG
      USE debug
# endif
      implicit none
# include "param.h"
      integer Istr,Iend,Jstr,Jend, i,j,k, kbak,kold,
     &        imin,imax,jmin,jmax,
     &        k1, k2, kp1
# ifdef MPI
     &       ,err
# endif
# ifdef PSOURCE
     &       ,is
# endif
      integer,dimension(2) :: ijmax
      real    mybeta,myalpha,myepsilon,mygamma,
     &        VMAX,VMAXL, cff,cff0,cff1,cff2,cff3, 
     &        DUnew,DVnew, dum_s

# ifdef NBQ_THETAIMP
      real,parameter :: thetaimp_nbq = 0.5 
# else
      real,parameter :: thetaimp_nbq = 1.
# endif
!
! Parameters for NBQ pressure gradient scheme
!
      real,parameter :: gammau=1.0  ! 1.0  --> second-order
                                    ! 0.45 --> increased stability range
      real,parameter :: gammau_2=(1./3.)*(1.-gammau)

      real
     &                Dnew(PRIVATE_2D_SCRATCH_ARRAY),
     &               rubar(PRIVATE_2D_SCRATCH_ARRAY),
     &               rvbar(PRIVATE_2D_SCRATCH_ARRAY),
     &                Drhs(PRIVATE_2D_SCRATCH_ARRAY),
     &                 UFx(PRIVATE_2D_SCRATCH_ARRAY),
     &                 UFe(PRIVATE_2D_SCRATCH_ARRAY),
     &                 VFx(PRIVATE_2D_SCRATCH_ARRAY),
     &                 VFe(PRIVATE_2D_SCRATCH_ARRAY),
     &                urhs(PRIVATE_2D_SCRATCH_ARRAY),  
     &                vrhs(PRIVATE_2D_SCRATCH_ARRAY),
     &                DUon(PRIVATE_2D_SCRATCH_ARRAY),
     &                DVom(PRIVATE_2D_SCRATCH_ARRAY)
      real
     &       ru_ext_nbq_sum(PRIVATE_2D_SCRATCH_ARRAY),
     &       rv_ext_nbq_sum(PRIVATE_2D_SCRATCH_ARRAY),
     &       ru_ext_nbq_old(PRIVATE_2D_SCRATCH_ARRAY),
     &       rv_ext_nbq_old(PRIVATE_2D_SCRATCH_ARRAY),
     &                 work(PRIVATE_2D_SCRATCH_ARRAY)
# ifdef NBQ
      real    Hzu_half_qdmu(PRIVATE_2D_SCRATCH_ARRAY,0:N), 
     &        Hzv_half_qdmv(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &     Hzw_half_nbq_inv(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &     Hzr_half_nbq_inv(PRIVATE_2D_SCRATCH_ARRAY,  N),
     &   Hzw_half_nbq_inv_u(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &   Hzw_half_nbq_inv_v(PRIVATE_2D_SCRATCH_ARRAY,0:N)
      real 
     &    dthetadiv_nbqdz_w(PRIVATE_2D_SCRATCH_ARRAY,2),
     &              dZdxq_u(PRIVATE_2D_SCRATCH_ARRAY,2),
     &              dZdyq_v(PRIVATE_2D_SCRATCH_ARRAY,2),
     &                   FX(PRIVATE_2D_SCRATCH_ARRAY),
     &                   FY(PRIVATE_2D_SCRATCH_ARRAY),
     &                   FC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &                   DC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &                   CF(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &                   BC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &         FX_3D(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &         FY_3D(PRIVATE_2D_SCRATCH_ARRAY,0:N)
#  define qdmw_nbq_old FX_3D
#  define dthetadiv_nbqdz_uv FY_3D
# endif
# ifdef UV_COR_NT
      real
     &               ntcoru(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &               ntcorv(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &               ntcorw(PRIVATE_2D_SCRATCH_ARRAY,0:N)
# endif

# include "grid.h"
# include "ocean2d.h"
# include "ocean3d.h"
# include "coupling.h"
# include "forces.h"
# ifdef NBQ
      real sum_nbq
#  if !defined VARYING_CSOUND
      real cff1_cst
#  endif
# endif
# include "openacc_work.h"
# ifdef MRL_WCI
      real vstu,ustv,dudx,dvdx,dude,dvde
# endif
# include "mixing.h"
# include "climat.h"
# include "scalars.h"
# include "sources.h"
# ifdef AGRIF
#  include "zoom.h"
        integer irhox, irhoy, irhot
# endif
# if defined INTERNAL || defined BODYTIDE
      real,parameter :: U0 = 0.02
      real,parameter :: omega = 2.*pi/(12.4*3600.)
# endif
# ifdef ACOUSTIC
      real dist_d  
# endif
# ifdef WET_DRY
      real cff1_WD,cff2_WD
# endif
# include "nbq.h"
!
# ifdef MASKING
#  define SWITCH *
# else
#  define SWITCH !
# endif
!
# ifdef MPI
#  include "mpi_cpl.h"
      include 'mpif.h'
#  define LOCALLM Lmmpi
#  define LOCALMM Mmmpi
# else
#  define LOCALLM Lm
#  define LOCALMM Mm
# 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
!
!*********************************************************************
! Set several approximations to speed-up computation
!*********************************************************************
!
# ifdef NBQ
!----------------------------------------------------------------------
! d./ds terms update frequency
!----------------------------------------------------------------------
!
#  ifdef NBQ_GRID_SLOW
#   define NSTEP_DS mod(iic,1).eq.1 .and. iif.eq.1
#  else
#   define NSTEP_DS .true. 
#  endif
!
!----------------------------------------------------------------------
! Grid update frequency (Hz ...)
!----------------------------------------------------------------------
!
#  ifdef NBQ_GRID_SLOW
#   define NSTEP_GRID .or. iif.eq.nfast
#  else
#   define NSTEP_GRID .or. mod(iif,1).eq.0
#  endif

# endif /* NBQ */
!
!*********************************************************************
!
!                 EXTERNAL (2D) FAST MODE PROCESSING
!
!  Compute external (2D) rhs terms rubar,rvbar of 2D fast-mode equations
!  using Generalized Forward-Backward AB3-AM4 algorithm, and update 
!  internal and external forcing terms for 3D fast mode equations: 
!  ru_int, ru_ext~(rufrc+rubar).
!
!  1- AB3 forward step for D,ubar,vbar
!  2- Advance zeta(m+1)
!  3- AM4 backward step for rubar,rvbar 
!  4- update rufrc,rvfrc
!  5- update ru_int,rv_int
!  6- make some backups
!
! Reference for Generalized FB scheme:
! ------------------------------------
! Shchepetkin, A.F., and J.C. McWilliams, 2009: Computational kernel 
! algorithms for fine-scale, multiprocess, longtime oceanic simulations. 
! Pp. 119–182 in Handbook of Numerical Analysis: Computational Methods 
! for the Atmosphere and Oceans. R.M. Teman and J.J. Tribbia, eds, 
! Elsevier Science.
!
!*********************************************************************

!
!=====================================================================
! Store rho.h at first slow and fast time-step
!=====================================================================
!
!$acc kernels if(compute_on_device) default(present)

# ifdef NBQ
      if (FIRST_TIME_STEP.and.FIRST_FAST_STEP) then
         do k=1,N
           do j=JstrV-2,Jend+1
             do i=IstrU-2,Iend+1
               rho_bak(i,j,k)=rho(i,j,k) 
             enddo
            enddo
          enddo
      endif

      ! Extrapolation in time:
      do j=JstrR,JendR
          do k=1,N
            do i=IstrR,IendR
              rho_grd(i,j,k)=(rho(i,j,k)*1.5-rho_bak(i,j,k)*0.5)/rho0
            ! rho_grd(i,j,k)=rho(i,j,k)/rho0
            enddo
          enddo
      enddo
# endif
 
# ifdef M3FAST_REINIT
!
!  Reinitialise fast mode at each first fast step
!
      if (FIRST_FAST_STEP) then
        do k=1,N
          do j=Jstr-1,Jend+1
            do i=IstrU-1,Iend+1
              qdmu_nbq(i,j,k)=0.5*u(i,j,k,nstp)*(Hz(i,j,k)+Hz(i-1,j,k))
            enddo
          enddo 
        enddo
        do k=1,N
          do j=JstrV-1,Jend+1
            do i=Istr-1,Iend+1
              qdmv_nbq(i,j,k)=0.5*v(i,j,k,nstp)*(Hz(i,j,k)+Hz(i,j-1,k))
            enddo
          enddo 
        enddo
#  ifdef NBQ
        do j=JstrV-1,Jend
          do k=1,N-1
            do i=IstrU-1,Iend
              qdmw_nbq(i,j,k)=0.5*wz(i,j,k,nstp)*(Hz(i,j,k)+Hz(i,j,k+1))
            enddo
          enddo
          do i=IstrU-1,Iend
            qdmw_nbq(i,j,0)=0.5*wz(i,j,0,nstp)*Hz(i,j,1)
          enddo
          do i=IstrU-1,Iend
            qdmw_nbq(i,j,N)=0.5*wz(i,j,N,nstp)*Hz(i,j,N)
          enddo
        enddo
#  endif
      endif    ! FIRST_FAST_STEP
# endif
!$acc end kernels

!
!=====================================================================
! AB3 Forward Step: compute total depth of water column and vertically
! --- ------- ----  integrated mass fluxes which are needed to compute 
! rhs terms of the barotropic momentum equations (rubar,rvbar).
!=====================================================================
!
!----------------------------------------------------------------------
!  Set indices to extrapolate (D,ubar,vbar) at m+1/2 (AB3)
!----------------------------------------------------------------------
!
      mybeta=0.281105 ! parameter for AB3 extrapolation

       if (FIRST_FAST_STEP) then     
                                      ! Meaning of temporal indices
        kbak=kstp                     ! ------- -- -------- -------
        kold=kstp                     ! m-2     m-1      m      m+1
        cff1= 1.0                     ! kold    kbak     kstp   knew
        cff2= 0.0
        cff3= 0.0
       elseif (FIRST_FAST_STEP+1) then  
        kbak=kstp-1                   ! AB2 forward scheme
        if (kbak.lt.1) kbak=4
        kold=kbak
        cff1= 1.5
        cff2=-0.5
        cff3= 0.0
      else                             ! AB3 forward scheme
        kbak=kstp-1 
        if (kbak.lt.1) kbak=4
        kold=kbak-1
        if (kold.lt.1) kold=4
        cff1= 1.5+mybeta
        cff2=-2.0*mybeta-0.5
        cff3= mybeta
      endif
!
# ifdef CVTK_DEBUG
C$OMP BARRIER
C$OMP MASTER
       call check_tab2d(zeta(:,:,kstp),'zeta step3d_fast #0','r',
     &    ondevice=.TRUE.)
C$OMP END MASTER
# endif  

!$acc kernels if(compute_on_device) default(present)
!
!----------------------------------------------------------------------
!  Extrapolate (D,ubar,vbar) at m+1/2
!----------------------------------------------------------------------
!
! Total depth/mass at m+1/2
!
      do j=JstrV-2,Jend+1
        do i=IstrU-2,Iend+1
# ifdef NBQ_MASS
          Drhs(i,j)=cff1*(zeta(i,j,kstp)+h(i,j))*rhobar_nbq(i,j,kstp)
     &             +cff2*(zeta(i,j,kbak)+h(i,j))*rhobar_nbq(i,j,kbak)
     &             +cff3*(zeta(i,j,kold)+h(i,j))*rhobar_nbq(i,j,kold)
# else
          Drhs(i,j)=cff1*(zeta(i,j,kstp)+h(i,j))
     &             +cff2*(zeta(i,j,kbak)+h(i,j))
     &             +cff3*(zeta(i,j,kold)+h(i,j))
# endif /* NBQ_MASS */
        enddo
      enddo
      
!
! Depth-average ubar velocity at m+1/2
!
      do j=Jstr-1,Jend+1
        do i=IstrU-1,Iend+1
          urhs(i,j)=cff1*ubar(i,j,kstp) 
     &             +cff2*ubar(i,j,kbak)
     &             +cff3*ubar(i,j,kold)
# if defined MRL_WCI && defined MASKING
          urhs(i,j)=urhs(i,j)*umask(i,j)+ust2d(i,j)*(umask(i,j)-1.0)
# endif
          DUon(i,j)=0.5*(Drhs(i,j)+Drhs(i-1,j))*on_u(i,j)*( urhs(i,j)
# ifdef MRL_WCI
     &                                                   + ust2d(i,j)
# endif
     &                                                              )
        enddo
      enddo
!
! Depth-average ubar velocity at m+1/2
!
      do j=JstrV-1,Jend+1
        do i=Istr-1,Iend+1
          vrhs(i,j)=cff1*vbar(i,j,kstp)
     &             +cff2*vbar(i,j,kbak)
     &             +cff3*vbar(i,j,kold)
# if defined MRL_WCI && defined MASKING
          vrhs(i,j)=vrhs(i,j)*vmask(i,j)+vst2d(i,j)*(vmask(i,j)-1.0)
# endif
          DVom(i,j)=0.5*(Drhs(i,j)+Drhs(i,j-1))*om_v(i,j)*( vrhs(i,j)
# ifdef MRL_WCI
     &                                                   + vst2d(i,j)
# endif
     &                                                              )
        enddo
      enddo
!$acc end kernels

# ifdef OBC_VOLCONS
      call set_DUV_bc_tile (Istr,Iend,Jstr,Jend, Drhs, DUon,DVom)
# endif

# ifdef CVTK_DEBUG_ADVANCED
C$OMP BARRIER
C$OMP MASTER
      call check_tab2d(zeta(:,:,kstp),'zeta step3d_fast #1','r',
     &   ondevice=.TRUE.)
C$OMP END MASTER
# endif 
!
!-----------------------------------------------------------------------
! Load RHS values into additional AGRIF shared array for nesting
!-----------------------------------------------------------------------
! 
# ifdef AGRIF
      irhox = Agrif_Irhox()
      irhoy = Agrif_Irhoy()
      irhot = Agrif_Irhot()

      if (FIRST_FAST_STEP) then
        do j=Jstr-1,Jend+1
          do i=Istr-1,Iend+1
            Zt_avg3(i,j,0)=zeta(i,j,kstp)       
          enddo
        enddo 
        do j=JstrR,JendR
          do i=Istr,IendR
            DU_avg3(i,j,0)=DUon(i,j)
          enddo
        enddo 
        do j=Jstr,JendR
          do i=IstrR,IendR
            DV_avg3(i,j,0)=DVom(i,j)
          enddo
        enddo 
      endif

#  if defined BD_EXCHANGE
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                        Zt_avg3(START_2D_ARRAY,0))
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        DU_avg3(START_2D_ARRAY,0))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        DV_avg3(START_2D_ARRAY,0))
#  endif

#  ifdef CVTK_DEBUG_ADVANCED
      if (.not.agrif_Root()) then
C$OMP BARRIER
C$OMP MASTER
        call check_tab2d(Zt_avg3(:,:,0),'Zt_avg3 (0) step3d_fast','r')
        call check_tab2d(DU_avg3(:,:,0),'DU_avg3 (0) step3d_fast','u')
        call check_tab2d(DV_avg3(:,:,0),'DV_avg3 (0) step3d_fast','v')
C$OMP END MASTER  
      endif
#  endif  
# endif /* AGRIF */
!
!=====================================================================
!
!  Compute zeta(m+1) -- first guest -- and update grid
!  
!=====================================================================
!
!----------------------------------------------------------------------
!  Get derived z grid variables at first fast step
!----------------------------------------------------------------------
!

# ifdef NBQ
      if (FIRST_FAST_STEP) then
        call grid_nbq_tile(Istr,Iend,Jstr,Jend,
     &                     Hzw_half_nbq_inv,   Hzr_half_nbq_inv,
     &                     Hzw_half_nbq_inv_u, Hzw_half_nbq_inv_v,
     &                     Hzu_half_qdmu,      Hzv_half_qdmv)
      endif
!
!-------------------------------------------------------------------
!  zeta(m+1) is derived here from the surface vertical velocity 
!  (surface characteristic relation) instead of depth-averaged 
!  conservation of mass. This satisfies dynamical coupling with 
!  the surface layer.
!
!  The surface kinematic relation is:
!
!  zeta(m+1) = zeta(m) + dtfast * [ wsurf_nbq(m)
!                                  - usurf_nbq(m)*dzeta/dx(m+0.5)
!                                  - vsurf_nbq(m)*dzeta/dy(m+0.5) ]
!
!  with zeta(m+0.5)=(1.5+beta)*zeta(m)-(0.5+2*beta)*zeta(m-1)
!                                             +beta*zeta(m-2) --> AB3
!-------------------------------------------------------------------
!  Computes surface velocities
!-------------------------------------------------------------------
!
!$acc kernels if(compute_on_device) default(present)    
      if (IstrU.le.Iend) then
        do j=Jstr,Jend
          do i=IstrU-1,Iend+1     
            usurf_nbq(i,j)=(qdmu_nbq(i,j,N)  
     &                     *2./(Hz(i,j,N)+Hz(i-1,j,N))
#  ifdef MRL_WCI
     &                     +ust(i,j,N)
#   ifdef WET_DRY
     &                      *umask_wet(i,j)
#   endif
#  endif
     &                     )
#  ifdef MASKING
     &                     *umask(i,j) 
#  endif
          enddo 
        enddo 
      endif

      if (JstrV.le.Jend) then
        do j=JstrV-1,Jend+1
          do i=Istr,Iend     
            vsurf_nbq(i,j)=(qdmv_nbq(i,j,N)              
     &                     *2./(Hz(i,j,N)+Hz(i,j-1,N))
#  ifdef MRL_WCI
     &                     +vst(i,j,N)
#   ifdef WET_DRY
     &                      *vmask_wet(i,j)
#   endif
#  endif
     &                      )
#  ifdef MASKING
     &                     *vmask(i,j) 
#  endif
          enddo
        enddo 
      endif

      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
#  ifdef NBQ
          wsurf_nbq(i,j)=(qdmw_nbq(i,j,N)         
#   ifdef NBQ_MASS     
     &                       /(1.+rho_grd(i,j,N))
#   endif
     &                       *Hzw_half_nbq_inv(i,j,N)
#  else
          wsurf_nbq(i,j)=(We(i,j,N)*pm(i,j)*pn(i,j)
#  endif
#  ifdef MRL_WCI
     &                     +wst(i,j,N)
#   ifdef WET_DRY
     &                       *rmask_wet(i,j)
#   endif
#  endif
     &                    )
#  ifdef MASKING
     &                       *rmask(i,j) 
#  endif
        enddo
      enddo 

!
!-------------------------------------------------------------------
! Apply point sources for river runoff simulations
!-------------------------------------------------------------------
!
#  ifdef PSOURCE
      do is=1,Nsrc 
#   ifdef MPI
        i=Isrc_mpi(is,mynode)
        j=Jsrc_mpi(is,mynode)
#   else
        i=Isrc(is)
        j=Jsrc(is)
#   endif
        if (IstrR.le.i .and. i.le.IendR .and.
     &      JstrR.le.j .and. j.le.JendR) then
          if (Dsrc(is).eq.0) then
            urhs(i,j)=2.*Qbar(is)/( on_u(i,j)
     &                             *(Drhs(i-1,j)+Drhs(i,j)) )
            DUon(i,j)=Qbar(is)
            usurf_nbq(i,j)=2.*Qsrc(is,N)
     &                     /( on_u(i,j)*(Hz(i,j,N)+Hz(i-1,j,N)) )
          else
            vrhs(i,j)=2.*Qbar(is)/( om_v(i,j)
     &                             *(Drhs(i,j-1)+Drhs(i,j)) )
            DVom(i,j)=Qbar(is)
            vsurf_nbq(i,j)=2.*Qsrc(is,N)
     &                     /( om_v(i,j)*(Hz(i,j,N)+Hz(i,j-1,N)) )
          endif
        endif
      enddo
#  endif
!
!-------------------------------------------------------------------
! Advance zeta at m+1
!-------------------------------------------------------------------
!
#  define zab3 UFx
      if (FIRST_FAST_STEP) then
        cff1 = 1.0
        cff2 = 0.0
        cff3 = 0.0
      elseif (FIRST_FAST_STEP+1) then  ! AB2
        cff1 = 1.5
        cff2 =-0.5
        cff3 = 0.0
      else                             ! AB3
        cff1 = 1.5+mybeta
        cff2 =-2.0*mybeta-0.5
        cff3 = mybeta
      endif

      do j=JstrV-2,Jend+1
        do i=IstrU-2,Iend+1
          zab3(i,j) =   cff1 * zeta(i,j,kstp)
     &                + cff2 * zeta(i,j,kbak)
     &                + cff3 * zeta(i,j,kold)
        enddo
      enddo

      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zeta(i,j,knew)=(zeta(i,j,kstp) + dtfast*( wsurf_nbq(i,j) 
     &                     -0.5*(usurf_nbq(i  ,j)
     &                                     *(zab3(i  ,j)
     &                                      -zab3(i-1,j))*pm_u(i,j)
     &                          +usurf_nbq(i+1,j)
     &                                     *(zab3(i+1,j)
     &                                      -zab3(i  ,j))*pm_u(i+1,j) 
     &                          )
     &                     -0.5*(vsurf_nbq(i  ,j)
     &                                     *(zab3(i,j  ) 
     &                                      -zab3(i,j-1))*pn_v(i,j)
     &                          +vsurf_nbq(i,j+1)
     &                                     *(zab3(i,j+1)
     &                                      -zab3(i,j  ))*pn_v(i,j+1) 
     &                          )
     &                                                ) )
#  ifdef MASKING
     &                                                  *rmask(i,j)
#  endif
        enddo
      enddo
#  undef zab3



# else /* NBQ */
!
!********************************************************************
!
!  Advance zeta at m+1 from continuity for the hydrostatic case *****
!
!********************************************************************
!
!  First, apply point sources
!
#  ifdef PSOURCE
      do is=1,Nsrc 
#   ifdef MPI
        i=Isrc_mpi(is,mynode)
        j=Jsrc_mpi(is,mynode)
#   else
        i=Isrc(is)
        j=Jsrc(is)
#   endif
        if (IstrR.le.i .and. i.le.IendR .and.
     &      JstrR.le.j .and. j.le.JendR) then
          if (Dsrc(is).eq.0) then
            urhs(i,j)=2.*Qbar(is)/( on_u(i,j)
     &                             *(Drhs(i-1,j)+Drhs(i,j)) )
            DUon(i,j)=Qbar(is)
          else
            vrhs(i,j)=2.*Qbar(is)/( om_v(i,j)
     &                             *(Drhs(i,j-1)+Drhs(i,j)) )
            DVom(i,j)=Qbar(is)
          endif
        endif
      enddo
#  endif
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zeta(i,j,knew)=(zeta(i,j,kstp) + dtfast*pm(i,j)*pn(i,j)
     &                                   *(DUon(i,j)-DUon(i+1,j  )
     &                                    +DVom(i,j)-DVom(i  ,j+1)))
#  ifdef MASKING
     &                                                   *rmask(i,j)
#  endif
        enddo
      enddo
      

# endif /* NBQ */
!
!-----------------------------------------------------------------------
!  Add nudging terms
!-----------------------------------------------------------------------
!
# ifdef ZNUDGING
#  ifdef ZONAL_NUDGING
      if (FIRST_TIME_STEP .or. mod(iic,10).eq.0) then
        if (FIRST_FAST_STEP) then
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   zeta(START_2D_ARRAY,knew),zetazon)
        endif
      endif
      if (FIRST_TIME_STEP) then
        if (FIRST_FAST_STEP) then
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   ssh(START_2D_ARRAY),sshzon)
        endif
      endif
#  endif  /* ZONAL_NUDGING */
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zeta(i,j,knew)=zeta(i,j,knew) + dtfast*Znudgcof(i,j)
#  ifdef ZONAL_NUDGING
     &                                 *(sshzon(j)-zetazon(j))
#  else
     &                              *(ssh(i,j)-zeta(i,j,knew))
#  endif /* ZONAL_NUDGING */
#  ifdef MASKING
     &                                             *rmask(i,j)
#  endif
        enddo
      enddo
# endif /* ZNUDGING */
!
!-----------------------------------------------------------------------
!  Compute wet/dry masks
!-----------------------------------------------------------------------
!
! First: modify new free-surface to ensure that depth 
!        is > Dcrit in masked cells.
!
# if defined WET_DRY && defined MASKING
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          cff=0.5+SIGN(0.5,Dcrit(i,j)-h(i,j))
          zeta(i,j,knew)=zeta(i,j,knew)+ 
     &                   cff*(Dcrit(i,j)-h(i,j))*(1.-rmask(i,j))
        enddo
      enddo 
# endif

!$acc end kernels

!
! Then compute wet/dry masks
!
# ifdef WET_DRY
      call wetdry_tile (Istr,Iend,Jstr,Jend)
# endif
!
!-----------------------------------------------------------------------
!  Set boundary conditions for the free-surface
!-----------------------------------------------------------------------
!
      call zetabc_tile (Istr,Iend,Jstr,Jend)
!
!-----------------------------------------------------------------------
!  Perform exchanges
!-----------------------------------------------------------------------
!
# if defined BD_EXCHANGE
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                        zeta(START_2D_ARRAY,knew))
# endif
!
!-----------------------------------------------------------------------
!  Update vertical grid
!
!  As soon as zeta(m+1) is known, the grid can be updated at fast
!  step m+1. In PERF option, the grid is updated at lower frequency 
!  and only during re-evaluation of zeta(m+1).
!
!  Caution: Hz_bak must be assigned only once in set_depth. The
!  following code must thus be consistant with set_depth routine.
!-----------------------------------------------------------------------
!
# if defined NBQ && !defined NBQ_GRID_SLOW
      if ( (FIRST_TIME_STEP .and. FIRST_FAST_STEP) NSTEP_GRID ) then
!
! Update main grid parameters
!
        call set_depth_tile(Istr,Iend,Jstr,Jend)
!
! Update derived grid variables
!
        call grid_nbq_tile(Istr,Iend,Jstr,Jend,
     &                     Hzw_half_nbq_inv,   Hzr_half_nbq_inv,
     &                     Hzw_half_nbq_inv_u, Hzw_half_nbq_inv_v,
     &                     Hzu_half_qdmu,      Hzv_half_qdmv )
      endif
# endif /* NBQ_GRID_SLOW */

!$acc kernels if(compute_on_device) default(present)
# ifndef NOT_NBQ_AM4
! 
!=====================================================================
! AM4 backward step: compute depth-averaged RHS of 2D fast-mode 
! --- -------- ----  momentum equations
!=====================================================================
!
!----------------------------------------------------------------------
! Surface pressure gradient
!----------------------------------------------------------------------
!
! Interpolate zeta fields half-step backward (AM4) for the subsequent 
! computation of barotropic pressure-gradient
!
      myalpha   = 0.25
      myepsilon = 0.00976186 - 0.13451357*myalpha
      mygamma   = 0.08344500 - 0.51358400*myalpha
!
      if (FIRST_FAST_STEP) then
        cff0=0.                  !---> Compute pressure-gradient
        cff1=1.                  !     terms using just zeta(:,:,kstp)
        cff2=0.
        cff3=0.
      elseif (FIRST_FAST_STEP+1) then
        cff0= 1.0833333333333    ! AM3 backward scheme
        cff1=-0.1666666666666    ! with coefficients chosen for
        cff2= 0.0833333333333    ! maximum stability, while maintaining
        cff3= 0.                 ! third-accuracy; alpha_max=1.73
      else
        cff0=0.5+2.*myepsilon+mygamma+2.*myalpha  ! AM4 backward scheme
        cff1=1.-cff0-mygamma-myepsilon            ! with implicit diffusion
        cff2=mygamma                              ! given by myalpha
        cff3=myepsilon
      endif
# endif
!
# define zwrk    UFx
# define rzeta   UFe
# define rzeta2  VFe
# define rzetaSA VFx

      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
# ifdef NOT_NBQ_AM4
          zwrk(i,j)=zeta(i,j,kstp)
# else
          zwrk(i,j)=(cff0*zeta(i,j,knew)+cff1*zeta(i,j,kstp)
     &              +cff2*zeta(i,j,kbak)+cff3*zeta(i,j,kold))
# endif
# ifdef NBQ_MASS
     &                                *rhobar_nbq(i,j,kstp) 
# endif
# ifdef VAR_RHO_2D
          rzeta(i,j)=(1.+rhoS(i,j))*zwrk(i,j)
          rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
          rzetaSA(i,j)=zwrk(i,j)*(rhoS(i,j)-rhoA(i,j))
# else
          rzeta(i,j)=zwrk(i,j)
          rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
# endif
        enddo
      enddo
!
! Compute surface pressure gradient
!
      cff=0.5*g
      do j=Jstr,Jend
        do i=Istr,Iend
          rubar(i,j)=cff*on_u(i,j)*( (h(i-1,j)+h(i,j))*(rzeta(i-1,j)
     &                        -rzeta(i,j)) +rzeta2(i-1,j)-rzeta2(i,j)
# ifdef VAR_RHO_2D
     &              +(h(i-1,j)-h(i,j))*( rzetaSA(i-1,j)+rzetaSA(i,j)
     &                        +0.333333333333*(rhoA(i-1,j)-rhoA(i,j))
     &                                      *(zwrk(i-1,j)-zwrk(i,j)))
# endif
# ifdef MRL_WCI
     &                  + ( h(i-1,j)+h(i,j)+rzeta(i-1,j)+rzeta(i,j) )
     &                                       *( sup(i,j)-sup(i-1,j) )
# endif
     &                                                              )
 
          rvbar(i,j)=cff*om_v(i,j)*( (h(i,j-1)+h(i,j))*(rzeta(i,j-1)
     &                        -rzeta(i,j)) +rzeta2(i,j-1)-rzeta2(i,j)
# ifdef VAR_RHO_2D
     &              +(h(i,j-1)-h(i,j))*( rzetaSA(i,j-1)+rzetaSA(i,j)
     &                        +0.333333333333*(rhoA(i,j-1)-rhoA(i,j))
     &                                      *(zwrk(i,j-1)-zwrk(i,j)))
# endif
# ifdef MRL_WCI
     &                  + ( h(i,j-1)+h(i,j)+rzeta(i,j-1)+rzeta(i,j) )
     &                                       *( sup(i,j)-sup(i,j-1) )
# endif
     &                                                              )
        enddo
      enddo            !--> discard  zwrk, rzeta, rzeta2, rzetaSA

# undef rzetaSA
# undef rzeta2
# undef rzeta
# undef zwrk
!
# ifdef UV_ADV
!----------------------------------------------------------------------
! Compute horizontal advection terms for momentum equations (2D only)
!-------- ---------- --------- ----- --- -------- --------- --- -----
!
! Centered second order advection scheme
!
! Numerical diffusion of momentum is implicitely added through 3D
! forcing of advection in rufrc and rvfrc (i.e., diffusion is
! at slow time scale)
!
! NOTE: mathematically necessary (minimal) index ranges for momentum-
! flux components are 
!
!      UFx(IstrU-1:Iend,Jstr:Jend)   VFx(Istr:Iend+1,JstrV:Jend)
!      UFe(IstrU:Iend,Jstr:Jend+1)   VFe(Istr,Iend,JstrV-1,Jend)
!
! however, for computational efficiency, these ranges are
! unified by suppressing U,V-suffices in order to allow fusion of the
! consecutive loops. This leads to slight increase of the redundant
! computations near western and southern boundaries in non-periodic
! directions.
!----------------------------------------------------------------------
!
      do j=Jstr,Jend
        do i=Istr-1,Iend
          UFx(i,j)=0.25*(DUon(i,j)+DUon(i+1,j))
     &                 *(urhs(i,j)+urhs(i+1,j))

          VFx(i+1,j)=0.25*(DUon(i+1,j)+DUon(i+1,j-1))
     &                   *(vrhs(i+1,j)+vrhs(i,j))
#  ifdef MASKING
     &                                 *pmask(i+1,j)
#  endif
        enddo
      enddo
      do j=Jstr-1,Jend
        do i=Istr,Iend
          VFe(i,j)=0.25*(DVom(i,j)+DVom(i,j+1))
     &                 *(vrhs(i,j)+vrhs(i,j+1))

          UFe(i,j+1)=0.25*(DVom(i,j+1)+DVom(i-1,j+1))
     &                   *(urhs(i,j+1)+urhs(i,j))
#  ifdef MASKING
     &                                 *pmask(i,j+1)
#  endif
        enddo
      enddo
      do j=Jstr,Jend
        do i=Istr,Iend
          rubar(i,j)=rubar(i,j)-UFx(i,j)+UFx(i-1,j)
     &                         -UFe(i,j+1)+UFe(i,j)

          rvbar(i,j)=rvbar(i,j)-VFx(i+1,j)+VFx(i,j)
     &                         -VFe(i,j)+VFe(i,j-1)
        enddo
      enddo    !--> discard UFx,VFe,UFe,VFx, DUon,DVom
# endif /* UV_ADV */
!
!-----------------------------------------------------------------------
! Compute Coriolis (2D and 3D) term and advective curvilinear metric
! terms (2D only).
!-----------------------------------------------------------------------
!
# if defined UV_COR || (defined CURVGRID && defined UV_ADV)
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          cff=Drhs(i,j)*(
#  ifdef UV_COR
     &                   fomn(i,j)
#  endif
#  if (defined CURVGRID && defined UV_ADV)
     &          +0.5*( dndx(i,j)*(vrhs(i,j)+vrhs(i,j+1))
     &                -dmde(i,j)*(urhs(i,j)+urhs(i+1,j)))
#  endif
     &                   )
#  ifdef MRL_WCI
#   if defined CURVGRID && defined UV_ADV
          cff1 = Drhs(i,j)*(
     &    0.5*( dndx(i,j)*(vst2d(i,j)+vst2d(i,j+1))
     &         -dmde(i,j)*(ust2d(i,j)+ust2d(i+1,j)) ))
#   else
          cff1 = 0.0
#   endif
          UFx(i,j)=(cff+cff1)*(vrhs(i,j)+vrhs(i,j+1))
     &                 +cff*(vst2d(i,j)+vst2d(i,j+1))
          VFe(i,j)=(cff+cff1)*(urhs(i,j)+urhs(i+1,j))
     &                 +cff*(ust2d(i,j)+ust2d(i+1,j))
#  else
          UFx(i,j)=cff*(vrhs(i,j)+vrhs(i,j+1))
          VFe(i,j)=cff*(urhs(i,j)+urhs(i+1,j))
#  endif 
        enddo
      enddo
      do j=Jstr,Jend
        do i=IstrU,Iend
          rubar(i,j)=rubar(i,j)+0.25*(UFx(i,j)+UFx(i-1,j))
        enddo
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
          rvbar(i,j)=rvbar(i,j)-0.25*(VFe(i,j)+VFe(i,j-1))
        enddo 
      enddo
# endif /* UV_COR */
!
!-----------------------------------------------------------------------
! Linear and/or quadratic bottom stress.
!-----------------------------------------------------------------------
!
# ifndef BSTRESS_FAST
#  ifdef BBL
        do j=Jstr,Jend
          do i=IstrU,Iend
            rubar(i,j)=rubar(i,j) - bustr(i,j)
     &                             *om_u(i,j)*on_u(i,j)
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            rvbar(i,j)=rvbar(i,j) - bvstr(i,j)
     &                             *om_v(i,j)*on_v(i,j)
          enddo
        enddo
#  else
      if (rdrg2.gt.0.) then
        do j=Jstr,Jend
          do i=IstrU,Iend
            cff=0.25*( vbar(i  ,j,kstp)+vbar(i  ,j+1,kstp)
     &                +vbar(i-1,j,kstp)+vbar(i-1,j+1,kstp))
 
            rubar(i,j)=rubar(i,j) - ubar(i,j,kstp)*( rdrg+rdrg2
     &              *sqrt(ubar(i,j,kstp)*ubar(i,j,kstp)+cff*cff)
     &                               )*om_u(i,j)*on_u(i,j)
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            cff=0.25*( ubar(i,j  ,kstp)+ubar(i+1,j  ,kstp)
     &                +ubar(i,j-1,kstp)+ubar(i+1,j-1,kstp))
 
            rvbar(i,j)=rvbar(i,j) - vbar(i,j,kstp)*( rdrg+rdrg2
     &              *sqrt(cff*cff+vbar(i,j,kstp)*vbar(i,j,kstp))
     &                               )*om_v(i,j)*on_v(i,j)
          enddo
        enddo
      else if (rdrg.gt.0.0) then
        do j=Jstr,Jend
          do i=IstrU,Iend
            rubar(i,j)=rubar(i,j) - rdrg*ubar(i,j,kstp)
     &                             *om_u(i,j)*on_u(i,j)
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            rvbar(i,j)=rvbar(i,j) - rdrg*vbar(i,j,kstp)
     &                             *om_v(i,j)*on_v(i,j)
          enddo
        enddo
      endif
#  endif
# endif /* BSTRESS_FAST */
!
!-----------------------------------------------------------------------
! Add 2D vortex-force terms combined with advection terms
!-----------------------------------------------------------------------
!
# ifdef MRL_WCI
      do j=Jstr,Jend
        do i=IstrU,Iend
          vstu = 0.25*( vst2d(i  ,j)+vst2d(i  ,j+1)
     &                 +vst2d(i-1,j)+vst2d(i-1,j+1) )
          dudx = 0.50*( urhs(i+1,j)-urhs(i-1,j  ) )
          dvdx = 0.50*( vrhs(i  ,j)-vrhs(i-1,j  )
     &                 +vrhs(i,j+1)-vrhs(i-1,j+1) )
          rubar(i,j) = rubar(i,j) + 0.5*on_u(i,j)*
     &                      ( Drhs(i-1,j)+Drhs(i,j) )
     &               *( ust2d(i,j)*dudx + vstu*dvdx )
        enddo
      enddo

      do j=JstrV,Jend
        do i=Istr,Iend
          ustv = 0.25*( ust2d(i,j  )+ust2d(i+1,j  )
     &                 +ust2d(i,j-1)+ust2d(i+1,j-1) )
          dude = 0.50*( urhs(i,j  )-urhs(i  ,j-1)
     &                 +urhs(i+1,j)-urhs(i+1,j-1) )
          dvde = 0.50*( vrhs(i,j+1)-vrhs(i  ,j-1) )
          rvbar(i,j) = rvbar(i,j) + 0.5*om_v(i,j)*
     &                      ( Drhs(i,j-1)+Drhs(i,j) )
     &               *( ustv*dude + vst2d(i,j)*dvde )
        enddo
      enddo
# endif
!
!=====================================================================
! Compute rufrc & rvfrc: internal mode forcing for barotropic fast mode
!
! During the first fast time step convert rufrc & fvfrc into forcing
! terms by subtracting the fast-time "rubar" and "rvbar" from them;
! These forcing terms are then extrapolated forward in time using
! optimized Adams-Bashforth weights, so that the resultant rufrc
! and rvfrc are centered effectively at time n+1/2. From now on,
! these newly computed forcing terms will remain constant during
! the fast time stepping and will be added to "rubar" and "rvbar"
! during all subsequent fast time steps.
!=====================================================================
!
      if (FIRST_FAST_STEP) then
!
        if (FIRST_TIME_STEP) then
          cff3=0.                        ! This version is designed
          cff2=0.                        ! for coupling during 3D
          cff1=1.                        ! predictor sub-step: here
        elseif (FIRST_TIME_STEP+1) then  ! forcing term "rufrc" is
          cff3= 0.                       ! computed as instantaneous
          cff2=-0.5                      ! value at 3D time step
          cff1= 1.5                      ! "nstp" first, and then
        else                             ! extrapolated half-step
          cff3= 0.281105                 ! forward using  AM3-like
          cff2=-0.5-2.*cff3              ! weights optimized for
          cff1= 1.5+cff3                 ! maximum stability (with
        endif                            ! special care for startup)
        do j=Jstr,Jend
         do i=IstrU,Iend
           cff=rufrc(i,j)-rubar(i,j)
           rufrc(i,j)=cff1*cff + cff2*rufrc_bak(i,j,3-nstp)
     &                         + cff3*rufrc_bak(i,j,nstp)
           rufrc_bak(i,j,nstp)=cff
         enddo
        enddo
        do j=JstrV,Jend
         do i=Istr,Iend
           cff=rvfrc(i,j)-rvbar(i,j)
           rvfrc(i,j)=cff1*cff + cff2*rvfrc_bak(i,j,3-nstp)
     &                         + cff3*rvfrc_bak(i,j,nstp)
           rvfrc_bak(i,j,nstp)=cff
         enddo
        enddo
      endif   !<-- FIRST_FAST_STEP
!
!=====================================================================
!   Update internal and external forcing terms for NBQ mode
!
!   Compute external forcing terms ru_ext_nbq and updated internal 
!   forcing terms ru_int_nbq for NBQ equations
!
!   ru_int_nbq     : RHS (3D) ( *mask & 2D correction)
!   ru_ext_nbq     : RHS (2D)
!   ru_ext_nbq_old : RHS (2D) at previous time-step
!   ru_ext_nbq_sum : time-integrated RHS (2D)
!=====================================================================
!
!-----------------------------------------------------------------------
!  First fast time step only
!-----------------------------------------------------------------------
!
      if (FIRST_FAST_STEP) then
!
        do j=Jstr,Jend
          do i=IstrU,Iend
            ru_ext_nbq_sum(i,j)=0.
            ru_ext_nbq_old(i,j)=0.
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            rv_ext_nbq_sum(i,j)=0.
            rv_ext_nbq_old(i,j)=0.
          enddo
        enddo 
# ifdef MASKING
        do k=1,N
          do j=Jstr,Jend
            do i=IstrU,Iend
              ru_int_nbq(i,j,k)=ru_int_nbq(i,j,k)*umask(i,j)
            enddo
          enddo
        enddo
        do k=1,N
          do j=JstrV,Jend
            do i=Istr,Iend
              rv_int_nbq(i,j,k)=rv_int_nbq(i,j,k)*vmask(i,j)
            enddo
          enddo
        enddo
# endif
      endif ! FIRST_FAST_STEP
!
!-----------------------------------------------------------------------
!  All fast time steps
!-----------------------------------------------------------------------
!
# ifndef NBQ
      if (mod(iif-1,inc_faststep).eq.0) then
        nb_faststep = 0
        do j=Jstr,Jend
          do i=IstrU,Iend
            rubar_sum(i,j) = 0.
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            rvbar_sum(i,j) = 0.
          enddo
        enddo
      endif
      nb_faststep = nb_faststep + 1
# endif

# define ru_ext_nbq UFx
      do j=Jstr,Jend
        do i=IstrU,Iend
# ifdef NBQ
          ru_ext_nbq(i,j)=(rufrc(i,j)+rubar(i,j))
     &                                   *pm_u(i,j)*pn_u(i,j)
     &                                  /(Drhs(i,j)+Drhs(i-1,j))
#  ifdef MASKING
     &                                                *umask(i,j)
#  endif
          ru_ext_nbq_old(i,j)=ru_ext_nbq(i,j)-ru_ext_nbq_old(i,j)
# else
          rubar_nbq(i,j)=(rufrc(i,j)+rubar(i,j))*pm_u(i,j)*pn_u(i,j)
#  ifdef MASKING
     &                                                *umask(i,j)
#  endif
          ru_ext_nbq(i,j)=rubar_nbq(i,j)/(Drhs(i,j)+Drhs(i-1,j))
          rubar_sum(i,j)=rubar_sum(i,j)+ru_ext_nbq(i,j)
# endif
          ru_ext_nbq_sum(i,j)=ru_ext_nbq_sum(i,j)+ru_ext_nbq(i,j)
        enddo
      enddo
# ifdef NBQ     
      do k=1,N
        do j=Jstr,Jend
          do i=IstrU,Iend
            ru_int_nbq(i,j,k)=ru_int_nbq(i,j,k)+
     &                        ru_ext_nbq_old(i,j)
     &                        *(Hz(i-1,j,k)+Hz(i,j,k))
          enddo
        enddo
      enddo
      do j=Jstr,Jend
        do i=IstrU,Iend
          ru_ext_nbq_old(i,j)=ru_ext_nbq(i,j)
        enddo
      enddo
# endif
# undef ru_ext_nbq

# define rv_ext_nbq UFx
      do j=JstrV,Jend
        do i=Istr,Iend
# ifdef NBQ
          rv_ext_nbq(i,j)=(rvfrc(i,j)+rvbar(i,j))
     &                                   *pm_v(i,j)*pn_v(i,j)
     &                                   /(Drhs(i,j)+Drhs(i,j-1))
#  ifdef MASKING
     &                                                *vmask(i,j)
#  endif
          rv_ext_nbq_old(i,j)=rv_ext_nbq(i,j)-rv_ext_nbq_old(i,j)
# else
          rvbar_nbq(i,j)=(rvfrc(i,j)+rvbar(i,j))*pm_v(i,j)*pn_v(i,j)
#  ifdef MASKING
     &                                                *vmask(i,j)
#  endif
          rv_ext_nbq(i,j)=rvbar_nbq(i,j)/(Drhs(i,j)+Drhs(i,j-1))
          rvbar_sum(i,j)=rvbar_sum(i,j)+rv_ext_nbq(i,j)
# endif
          rv_ext_nbq_sum(i,j)=rv_ext_nbq_sum(i,j)+rv_ext_nbq(i,j)
        enddo
      enddo
# ifdef NBQ        
      do k=1,N
        do j=JstrV,Jend
          do i=Istr,Iend
            rv_int_nbq(i,j,k)=rv_int_nbq(i,j,k)+
     &                        rv_ext_nbq_old(i,j)
     &                        *(Hz(i,j-1,k)+Hz(i,j,k))
          enddo
        enddo
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
          rv_ext_nbq_old(i,j)=rv_ext_nbq(i,j)
        enddo
      enddo
# endif
# undef rv_ext_nbq 

# if defined CVTK_DEBUG_ADVANCED && defined NBQ
!$acc end kernels
C$OMP BARRIER
C$OMP MASTER
       call check_tab3d(ru_int_nbq,'ru_int_nbq (A)','u',
     &   ondevice=.TRUE.)
       call check_tab3d(rv_int_nbq,'rv_int_nbq (A)','v',
     &   ondevice=.TRUE.)
C$OMP END MASTER     
!$acc kernels if(compute_on_device) default(present)
# endif  
!
!=====================================================================
!  Initializations and backups for 3D NBQ equations
!=====================================================================
!
!-----------------------------------------------------------------------
!  Initializations of ru_nbq_avg2 (as qdmu_nbq) at first fast step
!  It is used at last fast step for computation of new ru_nbq_avg2
!-----------------------------------------------------------------------
!
      if (FIRST_FAST_STEP) then
        do k=1,N
          do j=Jstr,Jend
            do i=IstrU,Iend
              ru_nbq_avg2(i,j,k)=qdmu_nbq(i,j,k)
            enddo
          enddo 
        enddo
        do k=1,N
          do j=JstrV,Jend
            do i=Istr,Iend
              rv_nbq_avg2(i,j,k)=qdmv_nbq(i,j,k)
            enddo
          enddo 
        enddo
# ifdef NBQ
        do k=0,N
          do j=Jstr,Jend
            do i=Istr,Iend
              rw_nbq_avg2(i,j,k)=qdmw_nbq(i,j,k)
            enddo
          enddo 
        enddo
# endif
      endif    ! FIRST_FAST_STEP
!
!-----------------------------------------------------------------------
!  zw_nbq backups
!-----------------------------------------------------------------------
!
# ifdef NBQ
      if (FIRST_TIME_STEP .and. FIRST_FAST_STEP) then
        do j=JstrR,JendR
          do k=0,N
            do i=IstrR,IendR
              zw_nbq (i,j,k,:)=z_w(i,j,k)
            enddo
          enddo
        enddo
      endif
# endif
!
!-------------------------------------------------------------------
!  Acoustic wave emission
!-------------------------------------------------------------------
!
# ifdef ACOUSTIC 
      if (FIRST_TIME_STEP .and. FIRST_FAST_STEP) then
        period_exp = 0.025/2.
        for_a_exp  = 2.5
        amp_exp    = 1.e-2
        hmax_exp   = 128.
        dg_exp     = 2.
        time_nbq   = 0. 
      endif
# endif
!
!------------------------------------------------------------------
!  Fast bottom friction
!
!  Set bottom stress using logarithmic or linear 
!  and/or quadratic formulation. 
!-------------------------------------------------------------------
!
# if defined BSTRESS_FAST && !defined BBL
#  define UBOT UFx
#  define VBOT VFe

      if (mod(iif-1,inc_faststep).eq.0) then

        if (maxval(Zob).ne.0.) then
          do j=JstrV-1,Jend+1
            do i=IstrU-1,Iend+1
              UBOT(i,j)=2.*qdmu_nbq(i,j,1)/(Hz(i,j,1)+Hz(i-1,j,1))
              VBOT(i,j)=2.*qdmv_nbq(i,j,1)/(Hz(i,j,1)+Hz(i,j-1,1))
            enddo
          enddo
          do j=JstrV-1,Jend
            do i=IstrU-1,Iend
              !cff1=MAX(z_r(i,j,1)-z_w(i,j,0),Zob(i,j)+1.e-4)
              cff1=z_r(i,j,1)-z_w(i,j,0)
              cff=vonKar/LOG(cff1/Zob(i,j))
              work(i,j)=MIN(Cdb_max,MAX(Cdb_min,cff*cff))
            enddo
          enddo
          do j=Jstr,Jend
            do i=IstrU,Iend
              cff=0.25*(VBOT(i  ,j)+VBOT(i  ,j+1)+
     &                  VBOT(i-1,j)+VBOT(i-1,j+1))
              bustr(i,j)=0.5*(work(i-1,j)+work(i,j))*UBOT(i,j)*
     &                   SQRT(UBOT(i,j)*UBOT(i,j)+cff*cff)
            enddo
          enddo
          do j=JstrV,Jend
            do i=Istr,Iend
              cff=0.25*(UBOT(i,j  )+UBOT(i+1,j  )+
     &                  UBOT(i,j-1)+UBOT(i+1,j-1))
              bvstr(i,j)=0.5*(work(i,j-1)+work(i,j))*VBOT(i,j)*
     &                   SQRT(VBOT(i,j)*VBOT(i,j)+cff*cff)
            enddo
          enddo
          
        elseif (rdrg2.gt.0.) then
          do j=JstrV-1,Jend+1
            do i=IstrU-1,Iend+1
              UBOT(i,j)=2.*qdmu_nbq(i,j,1)/(Hz(i,j,1)+Hz(i-1,j,1))
              VBOT(i,j)=2.*qdmv_nbq(i,j,1)/(Hz(i,j,1)+Hz(i,j-1,1))
            enddo
          enddo
          do j=JstrV,Jend
            do i=Istr,Iend
              cff=0.25*(VBOT(i  ,j)+VBOT(i  ,j+1)+
     &                  VBOT(i-1,j)+VBOT(i-1,j+1))
              bustr(i,j)=rdrg2*UBOT(i,j)*
     &                   SQRT(UBOT(i,j)*UBOT(i,j)+cff*cff)
            enddo
          enddo
          do j=Jstr,Jend
            do i=IstrU,Iend
              cff=0.25*(UBOT(i,j  )+UBOT(i+1,j  )+
     &                  UBOT(i,j-1)+UBOT(i+1,j-1))
              bvstr(i,j)=rdrg2*VBOT(i,j)*
     &                    SQRT(VBOT(i,j)*VBOT(i,j)+cff*cff)
            enddo
          enddo
        else
          do j=Jstr,Jend
            do i=IstrU,Iend
              bustr(i,j)=rdrg*2.*qdmu_nbq(i,j,1)/(Hz(i,j,1)+Hz(i-1,j,1))
            enddo
          enddo
          do j=JstrV,Jend
            do i=Istr,Iend
              bvstr(i,j)=rdrg*2.*qdmv_nbq(i,j,1)/(Hz(i,j,1)+Hz(i,j-1,1))
            enddo
          enddo
        endif

!
! Set boundary conditions
!
#  ifndef EW_PERIODIC
      IF (EASTERN_EDGE) THEN
        DO j=Jstr,Jend
          bustr(Iend+1,j)=bustr(Iend,j)
        END DO
        DO j=JstrV,Jend
          bvstr(Iend+1,j)=bvstr(Iend,j)
        END DO
      END IF
      IF (WESTERN_EDGE) THEN
        DO j=Jstr,Jend
          bustr(IstrU-1,j)=bustr(IstrU,j)
        END DO
        DO j=JstrV,Jend
          bvstr(Istr-1,j)=bvstr(Istr,j)
        END DO
      END IF
#  endif
#  ifndef NS_PERIODIC
      IF (NORTHERN_EDGE) THEN
        DO i=IstrU,Iend
          bustr(i,Jend+1)=bustr(i,Jend)
        END DO
        DO i=Istr,Iend
          bvstr(i,Jend+1)=bvstr(i,Jend)
        END DO
      END IF
      IF (SOUTHERN_EDGE) THEN
        DO i=IstrU,Iend
          bustr(i,Jstr-1)=bustr(i,Jstr)
        END DO
        DO i=Istr,Iend
          bvstr(i,JstrV-1)=bvstr(i,JstrV)
        END DO
      END IF
#  endif
#  if !defined EW_PERIODIC && !defined NS_PERIODIC
      IF (SOUTHERN_EDGE.and.WESTERN_EDGE) THEN
        bustr(Istr,Jstr-1)=0.5*(bustr(Istr+1,Jstr-1)+bustr(Istr,Jstr))
        bvstr(Istr-1,Jstr)=0.5*(bvstr(Istr,Jstr)+bvstr(Istr-1,Jstr+1))
      END IF
      IF (SOUTHERN_EDGE.and.EASTERN_EDGE) THEN
        bustr(Iend+1,Jstr-1)=0.5*(bustr(Iend+1,Jstr)+bustr(Iend,Jstr-1))
        bvstr(Iend+1,Jstr)=0.5*(bvstr(Iend+1,Jstr+1)+bvstr(Iend,Jstr))
      END IF
      IF (NORTHERN_EDGE.and.WESTERN_EDGE) THEN
        bustr(Istr,Jend+1)=0.5*(bustr(Istr,Jend)+bustr(Istr+1,Jend+1))
        bvstr(Istr-1,Jend+1)=0.5*(bvstr(Istr-1,Jend)+bvstr(Istr,Jend+1))
      END IF
      IF (NORTHERN_EDGE.and.EASTERN_EDGE) THEN
        bustr(Iend+1,Jend+1)=0.5*(bustr(Iend+1,Jend)+bustr(Iend,Jend+1))
        bvstr(Iend+1,Jend+1)=0.5*(bvstr(Iend+1,Jend)+bvstr(Iend,Jend+1))
      END IF
#  endif
#  if defined BD_EXCHANGE
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,bustr(START_2D_ARRAY))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,bvstr(START_2D_ARRAY))
#  endif


      endif
#  undef UBOT
#  undef VBOT
# endif
!
!------------------------------------------------------------------
!  Implicit part: system setup
!-------------------------------------------------------------------
!    
      do j=Jstr,Jend
        do i=Istr,Iend
# ifdef NBQ
          work(i,j)=pm(i,j)*pn(i,j)
# else
          if (LAST_FAST_STEP) work(i,j)=pm(i,j)*pn(i,j)
# endif
          DU_nbq(i,j)=0.
          DV_nbq(i,j)=0.
        enddo
      enddo
!	
!------------------------------------------------------------------
!  Store qdmw_nbq into working array
!------------------------------------------------------------------
!
# ifdef NBQ
      if (LAST_FAST_STEP) then
        do k=0,N 
          do j=Jstr,Jend             
            do i=Istr,Iend
               rw_nbq(i,j,k)=qdmw_nbq(i,j,k) 
            enddo
          enddo
        enddo
      endif
# endif
!	
!------------------------------------------------------------------
!  Recover Hz at first fast step (if final Hz correction needed)
!------------------------------------------------------------------
!
# ifdef NBQ_HZCORRECT
      if (iic.gt.ntstart .and. FIRST_FAST_STEP) then
        do k=1,N
          do j=JstrV-2,Jend+1
            do i=IstrU-2,Iend+1
              Hz(i,j,k) = Hz_correct(i,j,k)
            enddo
          enddo
        enddo
      endif
# endif
!
!----------------------------------------------------------------------
!  Store boundary values of nbq variables at previous 
!  time-step for use in radiation boundary conditions
!----------------------------------------------------------------------
!
# if defined OBC_NBQ && defined OBC_NBQORLANSKI
      call nbq_bry_store_tile (Istr,Iend,Jstr,Jend)
# endif
!
!
!***********************************************************************
!
!                  SOLVE FAST MODE 3D NBQ EQUATIONS
!
!***********************************************************************
!
!  Fast mode 3D momentum and mass-conservation equations can now be solved.
!
!  W-momentum equation is solved with explicit or implicit methods:
!  - Explicit scheme: w-momentum is updated right after (and the same 
!                     way as) u- and v-momentum.
!  - Implicit scheme: horizontal component of divergence is first 
!                     precomputed (as required by fast-mode mass 
!                     conservation) before tridiagonal Gauss Elimination 
!                     is carried out for qdmw_nbq(m).
!
!  Compressible pressure-force and second viscosity are calculated 
!  in thetadiv_nbq. Caution: this variable contains theta in the 
!  first part of the algorithm and momentum divergence in the remaining.
!
!  A Forward-backward scheme is implemented:
!  - Explicit scheme: Forward:  zeta, qdmu_nbq, qdmw_nbq. 
!                     Backward:       rho_nbq.
!  - Implicit scheme: Forward:  zeta, qdmu_nbq. 
!                     Backward:       qdmw_nbq, rho_nbq.
!
!***********************************************************************
!
# ifdef UV_COR_NT
!=====================================================================
!  Compute Non-traditional Coriolis pseudo-force
!
!    ntcoru = -e W cos(a) * Hz      [m2/s2]
!    ntcorv = +e W sin(a) * Hz
!    ntcorw = +e (U cos(a) - V sin(a) ) * Hz
!  
!  with e = 2 Omega cos(Phi)
!       a = angle between North and meridional grid axis
!
!=====================================================================
!
#  define Huw ntcoru
#  define Hvw ntcorv
      do j=Jstr,Jend
        do i=Istr,Iend+1
          Huw(i,j,0)=0.5*qdmu_nbq(i,j,1)
          Huw(i,j,N)=0.5*qdmu_nbq(i,j,N)
        enddo
      enddo
      do j=Jstr,Jend+1
        do i=Istr,Iend
          Hvw(i,j,0)=0.5*qdmv_nbq(i,j,1)
          Hvw(i,j,N)=0.5*qdmv_nbq(i,j,N)
        enddo
      enddo
      do k=1,N-1
        do j=Jstr,Jend
          do i=Istr,Iend+1
            Huw(i,j,k)=0.5*(qdmu_nbq(i,j,k)+qdmu_nbq(i,j,k+1))
          enddo
        enddo
        do j=Jstr,Jend+1
          do i=Istr,Iend
            Hvw(i,j,k)=0.5*(qdmv_nbq(i,j,k)+qdmv_nbq(i,j,k+1))
          enddo
        enddo
      enddo
      do k=0,N
        do j=Jstr,Jend
          do i=Istr,Iend
            ntcorw(i,j,k) = 0.5*e(i,j)*(
     &          cosa(i,j)*(Huw(i,j,k)+Huw(i+1,j,k))
     &        - sina(i,j)*(Hvw(i,j,k)+Hvw(i,j+1,k))
     &                                 )
          enddo
        enddo
      enddo
#  undef Huw
#  undef Hvw
!
      do k=1,N
        do j=JstrV-1,Jend
          do i=IstrU-1,Iend
            UFx(i,j)=-0.5*e(i,j)*cosa(i,j)
     &                            *( qdmw_nbq(i,j,k-1)
     &                              +qdmw_nbq(i,j,k  )) 
            VFe(i,j)=+0.5*e(i,j)*sina(i,j)
     &                            *( qdmw_nbq(i,j,k-1)
     &                              +qdmw_nbq(i,j,k)  )
          enddo
        enddo
        do j=Jstr,Jend
          do i=IstrU,Iend
            ntcoru(i,j,k)=0.5*(UFx(i,j)+UFx(i-1,j))
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            ntcorv(i,j,k)=0.5*(VFe(i,j)+VFe(i,j-1))
          enddo
        enddo
      enddo
# endif
!
!=====================================================================
!  Compute "Pressure - Viscosity" component (theta)
!=====================================================================
!
# ifdef NBQ
      do k=1,N
        do j=JstrV-2,Jend+1
          do i=IstrU-2,Iend+1
            thetadiv_nbq(i,j,k)=( -visc2_nbq*( thetadiv_nbq(i,j,k)
     &                                       +thetadiv3_nbq(i,j,k))
     &                         +soundspeed2_nbq(i,j)*rho_nbq(i,j,k) )
     &                                     *Hzr_half_nbq_inv(i,j,k)  
          enddo
        enddo
      enddo
!
!=====================================================================
!  Integrate fast-mode momentum equations
!
!  dqdm/dt=dtfast*(Compressible pressure force + second viscosity
!                  + gravity + non-NBQ forces + nudging)
!
!  If explicit scheme: all (x,y,z) mom components are done here
!  If implicit scheme: z-mom component is done after
!=====================================================================
!
!-----------------------------------------------------------------------
!  Pressure-Viscosity forces in XI- and ETA-Directions
!-----------------------------------------------------------------------
!
! U - component

#  ifdef NBQ_GRID_SLOW
        if (NSTEP_DS) then
#  endif
! Top Boundary conditions
      do j=Jstr-1,Jend
        do i=Istr-1,Iend
#  ifdef NBQ_GRID_SLOW
          dthetadiv_nbqdz(i,j,N,1)  = - thetadiv_nbq(i,j,N)
#  else
          dthetadiv_nbqdz(i,j,N)    = - thetadiv_nbq(i,j,N)
#  endif
        enddo
      enddo

      do k=1,N-1
        do j=Jstr-1,Jend
          do i=Istr-1,Iend
#  ifdef NBQ_GRID_SLOW
                  dthetadiv_nbqdz(i,j,k,1) = thetadiv_nbq(i,j,k+1)
     &                                   - thetadiv_nbq(i,j,k)
#  else
                  dthetadiv_nbqdz(i,j,k)   = thetadiv_nbq(i,j,k+1) 
     &                                   - thetadiv_nbq(i  ,j,k)
#  endif
             enddo
          enddo
      enddo

! Bottom boundary condition
      do j=Jstr,Jend
        do i=Istr-1,Iend
          dthetadiv_nbqdz_uv(i,j,0)=0. 
        enddo
      enddo

      do k=1,N
        do j=Jstr,Jend
          do i=Istr,Iend
#  ifdef NBQ_GRID_SLOW
                dthetadiv_nbqdz_uv(i,j,k)=Hzw_half_nbq_inv_u(i,j,k)*(
     &                                    dthetadiv_nbqdz(i,j,k,1)
     &                                   +dthetadiv_nbqdz(i-1,j,k,1))
#  else
                dthetadiv_nbqdz_uv(i,j,k)=Hzw_half_nbq_inv_u(i,j,k)*(
     &                                    dthetadiv_nbqdz(i,j,k)
     &                                   +dthetadiv_nbqdz(i-1,j,k))
#  endif          
          enddo
        enddo
      enddo

#  ifdef NBQ_GRID_SLOW
      do k=1,N
        do j=Jstr,Jend
          do i=IstrU,Iend
            if (k.gt.1.and.k.lt.N) then 
                dthetadiv_nbqdz(i,j,k,1)=(z_r(i  ,j,k)
     &                                     -z_r(i-1,j,k))  
     &              *(dthetadiv_nbqdz_uv(i,j,k)+
     &                dthetadiv_nbqdz_uv(i,j,k-1)) ! dZdx * (d(delta p)dz)_u
          elseif (k.gt.1) then
              dthetadiv_nbqdz(i,j,k,1)=(z_r(i ,j,k)
     &                                     -z_r(i-1,j,k))       
     &                 *dthetadiv_nbqdz_uv(i,j,k-1) ! dZdx * (d(delta p)dz)_u
     &                 +(z_w(i,j,N)-z_w(i-1,j,N))
     &                 *dthetadiv_nbqdz_uv(i,j,k)
          else
              dthetadiv_nbqdz(i,j,k,1)=(z_r(i ,j,k)
     &                                     -z_r(i-1,j,k))       
     &             *2.*dthetadiv_nbqdz_uv(i,j,k) ! dZdx * (d(delta p)dz)_u
          endif
          enddo 
        enddo
      enddo
      
      endif ! NSTEP_DS
#  endif

!-----------------------------------------------------------------------
!  Fast-mode U-momentum: qdmu_nbq
!-----------------------------------------------------------------------
!
      do k=1,N
        do j=Jstr,Jend
          do i=IstrU,Iend
            if (k.gt.1.and.k.lt.N) then 
#  ifdef NBQ_GRID_SLOW
              dum_s=dthetadiv_nbqdz(i,j,k,1)
#  else
              dum_s=(z_r(i,j,k)-z_r(i-1,j,k))                   
     &              *(dthetadiv_nbqdz_uv(i,j,k)+
     &                dthetadiv_nbqdz_uv(i,j,k-1)) ! dZdx * (d(delta p)dz)_u
#  endif
          elseif (k.gt.1) then
#  ifdef NBQ_GRID_SLOW
            dum_s=dthetadiv_nbqdz(i,j,k,1)
#  else
            dum_s=(z_r(i,j,k)-z_r(i-1,j,k))                      
     &                *dthetadiv_nbqdz_uv(i,j,k-1) ! dZdx * (d(delta p)dz)_u
     &                +(z_w(i,j,N)-z_w(i-1,j,N))   
     &                *dthetadiv_nbqdz_uv(i,j,k)
#  endif
          else
#  ifdef NBQ_GRID_SLOW
            dum_s=dthetadiv_nbqdz(i,j,k,1)
#  else
            dum_s=(z_r(i,j,k)-z_r(i-1,j,k))                      
     &              *2.*dthetadiv_nbqdz_uv(i,j,k) ! dZdx * (d(delta p)dz)_u
#  endif
          endif
#  ifdef MASKING
          if (umask(i-1,j)*umask(i+1,j) .eq. 0.) then
            dum_s=dum_s-(thetadiv_nbq(i  ,j,k)-
     &                       thetadiv_nbq(i-1,j,k))  ! - d(delta p)dx
          else
#  endif
            dum_s=dum_s
     &                -(gammau  *thetadiv_nbq(i  ,j,k)+
     &                  gammau_2*thetadiv_nbq(i+1,j,k)-
     &                  gammau  *thetadiv_nbq(i-1,j,k)-
     &                  gammau_2*thetadiv_nbq(i-2,j,k)) ! - d(delta p)dx
#  ifdef MASKING
          endif
#  endif
          dum_s=dum_s*Hzu_half_qdmu(i,j,k)
#  ifdef UV_COR_NT
          dum_s=dum_s+ntcoru(i,j,k)
#  endif
#  if defined INTERNAL || defined BODYTIDE
          dum_s=dum_s+0.5*omega*U0*cos(omega*time)
     &                       *(Hz(i-1,j,k)+Hz(i,j,k))
#  endif
#  ifdef BSTRESS_FAST
          if (k.eq.1) dum_s=dum_s-bustr(i,j)
#  endif
          qdmu_nbq(i,j,k)=qdmu_nbq(i,j,k) + dtnbq * (
     &                        dum_s + ru_int_nbq(i,j,k) )  
#  ifdef MASKING
          qdmu_nbq(i,j,k)=qdmu_nbq(i,j,k)*umask(i,j)
#  endif
          ru_nbq(i,j,k)=dum_s/work(i,j) 

#  if defined NBQ_NUDGING && defined NBQCLIMATOLOGY
          qdmu_nbq(i,j,k)=qdmu_nbq(i,j,k)*(1.-NBQnudgcof(i,j))
     &                        +u(i,j,k,nrhs)*Hzu_half_qdmu(i,j,k)
     &                                           *NBQnudgcof(i,j)
#  endif
          enddo 
        enddo
      enddo

      do j=Jstr,Jend
        do i=IstrU,Iend
          sum_nbq=0.
!$acc loop reduction(+:sum_nbq)
          do k=1,N
             sum_nbq=sum_nbq+qdmu_nbq(i,j,k)
          enddo
          DU_nbq(i,j)=DU_nbq(i,j)+sum_nbq
        enddo
      enddo

#  ifdef NBQ_GRID_SLOW
        if (NSTEP_DS) then
#  endif

! Bottom boundary condition
      do j=Jstr-1,Jend
        do i=Istr,Iend
          dthetadiv_nbqdz_uv(i,j,0)=0.
        enddo
      enddo 

      do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
#  ifdef NBQ_GRID_SLOW
                 dthetadiv_nbqdz_uv(i,j,k)=Hzw_half_nbq_inv_v(i,j,k)*(
     &                                     dthetadiv_nbqdz(i,j,k,1)
     &                                    +dthetadiv_nbqdz(i,j-1,k,1))
#  else
                 dthetadiv_nbqdz_uv(i,j,k)=Hzw_half_nbq_inv_v(i,j,k)*(
     &                                     dthetadiv_nbqdz(i,j,k)
     &                                    +dthetadiv_nbqdz(i,j-1,k))
#  endif
            enddo
          enddo  
      enddo

#  ifdef NBQ_GRID_SLOW

        do k=1,N
          do j=Jstr,Jend
            if (j.ge.JstrV) then
            do i=Istr,Iend
              if (k.gt.1.and.k.lt.N) then 
                  dthetadiv_nbqdz(i,j,k,2)=(z_r(i,j  ,k)
     &                                     -z_r(i,j-1,k)) 
     &                *(dthetadiv_nbqdz_uv(i,j,k)+
     &                  dthetadiv_nbqdz_uv(i,j,k-1)) ! dZdy * (d(delta p)dz)_v
              elseif (k.gt.1) then
                  dthetadiv_nbqdz(i,j,k,2)=(z_r(i,j  ,k)
     &                                     -z_r(i,j-1,k)) 
     &                *dthetadiv_nbqdz_uv(i,j,k-1) ! dZdy * (d(delta p)dz)_v
     &                +(z_w(i,j,N)-z_w(i,j-1,N))        
     &                *dthetadiv_nbqdz_uv(i,j,k)
              else
                  dthetadiv_nbqdz(i,j,k,2)=(z_r(i,j  ,k)
     &                                     -z_r(i,j-1,k)) 
     &                *2.*dthetadiv_nbqdz_uv(i,j,k) ! dZdy * (d(delta p)dz)_v
              endif
            enddo
           endif !<-- j.ge.JstrV
          enddo
        enddo

        endif ! NSTEP_DS
#  endif
!
!-----------------------------------------------------------------------
!  Fast-mode V-momentum: qdmv_nbq
!-----------------------------------------------------------------------
!
        do k=1,N
          do j=Jstr,Jend
            if (j.ge.JstrV) then
            do i=Istr,Iend
              if (k.gt.1.and.k.lt.N) then 
#  ifdef NBQ_GRID_SLOW
                dum_s=dthetadiv_nbqdz(i,j,k,2)
#  else
                dum_s=(z_r(i,j,k)-z_r(i,j-1,k)) 
     &                *(dthetadiv_nbqdz_uv(i,j,k)+
     &                  dthetadiv_nbqdz_uv(i,j,k-1)) ! dZdy * (d(delta p)dz)_v
#  endif
              elseif (k.gt.1) then
#  ifdef NBQ_GRID_SLOW
                dum_s=dthetadiv_nbqdz(i,j,k,2)
#  else
                dum_s=(z_r(i,j,k)-z_r(i,j-1,k))            
     &                *dthetadiv_nbqdz_uv(i,j,k-1) ! dZdy * (d(delta p)dz)_v
     &                +(z_w(i,j,N)-z_w(i,j-1,N))      
     &                *dthetadiv_nbqdz_uv(i,j,k)
#  endif
              else
#  ifdef NBQ_GRID_SLOW
                dum_s=dthetadiv_nbqdz(i,j,k,2) 
#  else
                dum_s=(z_r(i,j,k)-z_r(i,j-1,k)) 
     &                *2.*dthetadiv_nbqdz_uv(i,j,k) ! dZdy * (d(delta p)dz)_v
#  endif
              endif
                
#  ifdef MASKING
              if (vmask(i,j-1)*vmask(i,j+1) .eq. 0.) then
              dum_s=dum_s
     &             -(thetadiv_nbq(i,j  ,k)-
     &               thetadiv_nbq(i,j-1,k)) ! - d(delta p)dy
              else
#  endif
              dum_s=dum_s
     &             -(gammau  *thetadiv_nbq(i,j  ,k)+
     &               gammau_2*thetadiv_nbq(i,j+1,k)-
     &               gammau  *thetadiv_nbq(i,j-1,k)-
     &               gammau_2*thetadiv_nbq(i,j-2,k)) ! - d(delta p)dy
#  ifdef MASKING
              endif
#  endif
              dum_s=dum_s*Hzv_half_qdmv(i,j,k)
#  ifdef UV_COR_NT
              dum_s=dum_s+ntcorv(i,j,k)
#  endif
#  if defined INTERNAL || defined BODYTIDE
              dum_s=dum_s+0.25*(f(i,j)+f(i,j-1))*
     &                       U0*sin(omega*time)*
     &                       (Hz(i,j,k)+Hz(i,j-1,k))
#  endif
#  ifdef BSTRESS_FAST
              if (k.eq.1) dum_s=dum_s-bvstr(i,j)
#  endif
              qdmv_nbq(i,j,k)=qdmv_nbq(i,j,k) + dtnbq * (
     &                        dum_s + rv_int_nbq(i,j,k) )
#  ifdef MASKING
              qdmv_nbq(i,j,k)=qdmv_nbq(i,j,k)*vmask(i,j)
#  endif
              rv_nbq(i,j,k)=dum_s/work(i,j)  

#  if defined NBQ_NUDGING && defined NBQCLIMATOLOGY
              qdmv_nbq(i,j,k)=qdmv_nbq(i,j,k)*(1.-NBQnudgcof(i,j))
     &                        +v(i,j,k,nrhs)*Hzv_half_qdmv(i,j,k)
     &                                           *NBQnudgcof(i,j)
#  endif
            enddo
           endif !<-- j.ge.JstrV
          enddo
        enddo

        do j=JstrV,Jend
          do i=Istr,Iend
            sum_nbq=0.
!$acc loop reduction(+:sum_nbq)
            do k=1,N
               sum_nbq=sum_nbq+qdmv_nbq(i,j,k)
            enddo
            DV_nbq(i,j)=DV_nbq(i,j)+sum_nbq
          enddo
        enddo

# else /* NBQ */

!
!********************************************************************
!
!  Advance fast 3D u,v for the hydrostatic case *****
!
!********************************************************************
!
#  ifdef BSTRESS_FAST
      do j=Jstr,Jend
        do i=IstrU,Iend
          rubar_nbq(i,j)=rubar_nbq(i,j)-bustr(i,j)
        enddo
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend      
          rvbar_nbq(i,j)=rvbar_nbq(i,j)-bvstr(i,j)
        enddo
      enddo
#  endif
      if ((mod(iif-1,inc_faststep) .eq. inc_faststep-1) .OR.
     &    (LAST_FAST_STEP)) then
        do k=1,N
          do j=Jstr,Jend
            do i=IstrU,Iend
              dum_s=0.
#  ifdef BSTRESS_FAST
              if (k.eq.1) dum_s=dum_s-bustr(i,j)
#  endif
              qdmu_nbq(i,j,k)=qdmu_nbq(i,j,k)+dtfast * (
     &                      nb_faststep*(dum_s + ru_int_nbq(i,j,k))
     &                     +rubar_sum(i,j)*(Hz(i-1,j,k)+Hz(i,j,k)))    
#  ifdef MASKING
     &                                        *umask(i,j)
#  endif
              DU_nbq(i,j)=DU_nbq(i,j)+qdmu_nbq(i,j,k)
              if (LAST_FAST_STEP) ru_nbq(i,j,k)=dum_s/work(i,j) 
            enddo 
          enddo
          do j=JstrV,Jend
            do i=Istr,Iend
              dum_s=0.
#  ifdef BSTRESS_FAST
              if (k.eq.1) dum_s=dum_s-bvstr(i,j)
#  endif
              qdmv_nbq(i,j,k)=qdmv_nbq(i,j,k)+dtfast * (
     &                       nb_faststep*(dum_s + rv_int_nbq(i,j,k))
     &                      +rvbar_sum(i,j)*(Hz(i,j-1,k)+Hz(i,j,k)))
#  ifdef MASKING
     &                                        *vmask(i,j)
#  endif
              DV_nbq(i,j)=DV_nbq(i,j)+qdmv_nbq(i,j,k)
              if (LAST_FAST_STEP) rv_nbq(i,j,k)=dum_s/work(i,j)
            enddo
          enddo
        enddo
      endif

# endif /* NBQ */

!
!-----------------------------------------------------------------------
! Apply point sources for river runoff simulations
!-----------------------------------------------------------------------
!
# ifdef PSOURCE
      do is=1,Nsrc 
#  ifdef MPI
        i=Isrc_mpi(is,mynode)
        j=Jsrc_mpi(is,mynode)
#  else
        i=Isrc(is)
        j=Jsrc(is)
#  endif
        if (IstrR.le.i .and. i.le.IendR .and.
     &      JstrR.le.j .and. j.le.JendR) then
          if (Dsrc(is).eq.0) then
            DU_nbq(i,j)=0.
#  ifdef WET_DRY
            if (Qbar(is).gt.0) then
              umask_wet(i,j)= + 1
            else
              umask_wet(i,j)= - 1
            endif
#  endif
            do k=1,N
              qdmu_nbq(i,j,k)=Qsrc(is,k)*pn_u(i,j)
              DU_nbq(i,j)=DU_nbq(i,j)+qdmu_nbq(i,j,k)
            enddo
          else
            DV_nbq(i,j)=0.
#  ifdef WET_DRY
            if (Qbar(is).gt.0) then
              vmask_wet(i,j)= + 1
            else
              vmask_wet(i,j)= - 1
            endif
#  endif
            do k=1,N
              qdmv_nbq(i,j,k)=Qsrc(is,k)*pm_v(i,j)
              DV_nbq(i,j)=DV_nbq(i,j)+qdmv_nbq(i,j,k)
            enddo
          endif
        endif
      enddo
# endif
!
!-----------------------------------------------------------------------
!  U & V momentum wet mask
!-----------------------------------------------------------------------
!
# ifdef WET_DRY
      do j=Jstr,Jend
        do i=IstrU,Iend
          cff1_WD=ABS(ABS(umask_wet(i,j))-1.)
          cff2_WD=0.5+SIGN(0.5,DU_nbq(i,j))*umask_wet(i,j)
          umask_wet(i,j)=0.5*umask_wet(i,j)*cff1_WD
     &                         +cff2_WD*(1.-cff1_WD)
          DU_nbq(i,j)=DU_nbq(i,j)*umask_wet(i,j)
#  ifdef MRL_WCI
          ust2d(i,j)=ust2d(i,j)*umask_wet(i,j)
#  endif
          do k=1,N
            qdmu_nbq(i,j,k)=qdmu_nbq(i,j,k)*umask_wet(i,j)
          enddo
        enddo 
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
          cff1_WD=ABS(ABS(vmask_wet(i,j))-1.)
          cff2_WD=0.5+SIGN(0.5,DV_nbq(i,j))*vmask_wet(i,j)
          vmask_wet(i,j)=0.5*vmask_wet(i,j)*cff1_WD
     &                         +cff2_WD*(1.-cff1_WD)
          DV_nbq(i,j)=DV_nbq(i,j)*vmask_wet(i,j)
#  ifdef MRL_WCI
          vst2d(i,j)=vst2d(i,j)*vmask_wet(i,j)
#  endif
          do k=1,N
            qdmv_nbq(i,j,k)=qdmv_nbq(i,j,k)*vmask_wet(i,j)
          enddo
        enddo 
      enddo
# endif

# ifdef CVTK_DEBUG_ADVANCED
!$acc end kernels
C$OMP BARRIER
C$OMP MASTER
      call check_tab3d(qdmu_nbq,'qdmu_nbqint','uint',
     &   ondevice=.TRUE.)
      call check_tab3d(qdmv_nbq,'qdmv_nbqint','vint',
     &   ondevice=.TRUE.)
C$OMP END MASTER
!$acc kernels if(compute_on_device) default(present)
# endif
!
!-----------------------------------------------------------------------
!  U & V momentum open boundary conditions
!-----------------------------------------------------------------------
!
      do j=Jstr,Jend
        do i=IstrU,Iend
          ubar(i,j,knew)=urhs(i,j)
        enddo 
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
          vbar(i,j,knew)=vrhs(i,j)
        enddo 
      enddo

!$acc end kernels

# ifdef NBQ
      M2bc_nbq_flag=.true. ! apply boundary wet/dry conditions
                           ! and compute DU_nbq
      call u2dbc_tile   (Istr,Iend,Jstr,Jend, work)
      call v2dbc_tile   (Istr,Iend,Jstr,Jend, work)
# endif
      call unbq_bc_tile (Istr,Iend,Jstr,Jend, work)
      call vnbq_bc_tile (Istr,Iend,Jstr,Jend, work)
!
!-----------------------------------------------------------------------
! Exchange periodic boundaries and computational margins
!-----------------------------------------------------------------------
!
# if defined BD_EXCHANGE
#  ifndef NBQ
      if ((mod(iif-1,inc_faststep) .eq. inc_faststep-1) .OR.
     &     (LAST_FAST_STEP)) then
#  endif
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        qdmu_nbq(START_2D_ARRAY,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        qdmv_nbq(START_2D_ARRAY,1))
#  ifndef NBQ
      endif
#  endif
!      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
!     &                        DU_nbq(START_2D_ARRAY))
!      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
!     &                        DV_nbq(START_2D_ARRAY))
# endif
!$$$# ifdef CVTK_DEBUG_ADVANCED
!$$$C$OMP BARRIER
!$$$C$OMP MASTER
!$$$      call check_tab3d(qdmu_nbq,'qdmu_nbq','u')
!$$$      call check_tab3d(qdmv_nbq,'qdmv_nbq','v')
!$$$C$OMP END MASTER
!$$$# endif  

!$acc kernels if(compute_on_device) default(present)
# ifdef NBQ
!
!-----------------------------------------------------------------------
!      Fast mode W-Momentum equation: qdmw_nbq
!                 ==> EXPLICIT scheme
!-----------------------------------------------------------------------
!
#  ifndef NBQ_IMP
      do j=Jstr,Jend
        do k=1,N-1
          do i=Istr,Iend   
            dum_s = thetadiv_nbq(i,j,k) - thetadiv_nbq(i,j,k+1)   
#   ifdef UV_COR_NT
            dum_s = dum_s + ntcorw(i,j,k)
#   endif
            qdmw_nbq(i,j,k)=qdmw_nbq(i,j,k)   
     &                      + dtnbq * ( dum_s + rw_int_nbq(i,j,k) )
#   ifdef NBQ_GRAV
     &               -0.25*(rho_nbq(i,j,k  )*Hzr_half_nbq_inv(i,j,k  )
     &                     +rho_nbq(i,j,k+1)*Hzr_half_nbq_inv(i,j,k+1))
     &                    *(Hzr(i,j,k)+Hzr(i,j,k+1))
     &                    *g*dtnbq
#   endif
#   ifdef MASKING
          qdmw_nbq(i,j,k)=qdmw_nbq(i,j,k) * rmask(i,j)
#   endif
          enddo             
        enddo

        k=N
        do i=Istr,Iend
          dum_s =   thetadiv_nbq(i,j,N)
#   ifdef UV_COR_NT
          dum_s = dum_s + ntcorw(i,j,k)
#   endif
          qdmw_nbq(i,j,N)=qdmw_nbq(i,j,N)   
     &                   + dtnbq * ( dum_s + rw_int_nbq(i,j,N) )
#   ifdef NBQ_GRAV
     &                   -rho_nbq(i,j,N)*0.5
     &                   *g*dtnbq
#   endif
#   ifdef MASKING
          qdmw_nbq(i,j,k)=qdmw_nbq(i,j,k) * rmask(i,j)
#   endif
        enddo
      enddo !<-- j loop
!
!-----------------------------------------------------------------------
! Vertical momentum open boundary conditions
!-----------------------------------------------------------------------
!
      call wnbq_bc_tile (Istr,Iend,Jstr,Jend, work)

#   ifdef UV_COR_NT
#    if defined BD_EXCHANGE
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
     &                        qdmw_nbq(START_2D_ARRAY,0))
#    endif
#   endif

#  endif  /* ! NBQ_IMP */
!
!
!=====================================================================
!      Fast-mode conservation of mass
!  
!  ... and integration of W-momentum with IMPLICIT scheme
!
!  From now on, thetadiv_nbq array is used for divergence (not theta):
!
!=====================================================================
!
      if (IstrU.gt.Iend) then
        do j=Jstr,Jend
          do i=Istr,Iend+1
            FX(i,j)=0.
          enddo
        enddo
      endif
      if (JstrV.gt.Jend) then
        do j=Jstr,Jend+1
          do i=Istr,Iend 
            FY(i,j)=0.
          enddo
        enddo
      endif
!
!-----------------------------------------------------------------------
!  X-component dZdx*qdmu for horizontal divergence
!-----------------------------------------------------------------------
!

#  if !defined NBQ_GRID_SLOW
#   define dZdxq_u dthetadiv_nbqdz
!#   define dZdxq_w FX_3D
#  endif

#  ifdef NBQ_GRID_SLOW
        if (NSTEP_DS) then
#  endif
      do k=1,N
        do j=Jstr,Jend
          do i=Istr,Iend+1
            dZdxq_u(i,j,k)=(z_r(i,j,k)-z_r(i-1,j,k)) 
     &                           *qdmu_nbq(i,j,k)  ! (dZdx * (rho u))_u
          enddo
        enddo        
      enddo

! Bottom boundary conditions
#  ifdef NBQ_FREESLIP
            do j=Jstr,Jend
              do i=Istr,Iend+1 
                dZdxq_w(i,j,0 ) = (z_w(i,j,0)-z_w(i-1,j,0))
     &                             *qdmu_nbq(i,j,1)  
     &                             /(Hzr(i,j,1)+Hzr(i-1,j,1))
              enddo
            enddo

            do j=Jstr,Jend
              do i=Istr,Iend    
                qdmw_nbq(i,j,0)=0.5*(dZdxq_w(i  ,j,0)*pm_u(i  ,j)
     &                              +dZdxq_w(i+1,j,0)*pm_u(i+1,j) )
     &                                                 * Hzr(i,j,1)
#   ifdef MASKING
                qdmw_nbq(i,j,0)=qdmw_nbq(i,j,0)*rmask(i,j)
#   endif 
              enddo
            enddo 

#  else /* NBQ_FREESLIP */

#   ifndef NBQ_GRID_SLOW
            do j=Jstr,Jend
              do i=Istr,Iend+1
                dZdxq_w(i,j,0)=0.  
                qdmw_nbq(i,j,0)=0.     
              enddo
            enddo 
#   endif
 
#  endif /* NBQ_FREESLIP */

! Top boundary conditions
            do j=Jstr,Jend
              do i=Istr,Iend+1
                dZdxq_w(i,j,N )= (z_w(i,j,N)-z_w(i-1,j,N))
     &                            *qdmu_nbq(i,j,N)
     &                            /(Hzr(i,j,N)+Hzr(i-1,j,N))
              enddo
            enddo  

! Inner domain   
      do k=1,N-1
            do j=Jstr,Jend
              do i=Istr,Iend+1
                dZdxq_w(i,j,k )=Hzw_half_nbq_inv_u(i,j,k)
     &                                 *(dZdxq_u(i,j,k)+
     &                                   dZdxq_u(i,j,k+1))
              enddo 
            enddo
      enddo

#  ifdef NBQ_GRID_SLOW
        else   ! NSTEP_DS: Update d./ds terms
#  endif  /* NBQ_GRID_SLOW */

! Bottom boundary conditions
#  ifdef NBQ_FREESLIP
            do j=Jstr,Jend
              do i=Istr,Iend+1 
                dZdxq_w(i,j,0)=(z_w(i,j,0)-z_w(i-1,j,0))
     &                          *qdmu_nbq(i,j,1)  
     &                          /(Hzr(i,j,1)+Hzr(i-1,j,1))
              enddo
            enddo

            do j=Jstr,Jend
              do i=Istr,Iend    
                qdmw_nbq(i,j,0)=0.5*(dZdxq_w(i  ,j,0)*pm_u(i  ,j) 
     &                              +dZdxq_w(i+1,j,0)*pm_u(i+1,j) ) 
     &                                                  *Hzr(i,j,1)    
#   ifdef MASKING
                qdmw_nbq(i,j,0)=qdmw_nbq(i,j,0)*rmask(i,j)
#   endif 
              enddo
            enddo 

#  endif /* NBQ_FREESLIP */

#  ifdef NBQ_GRID_SLOW
        endif ! NSTEP_DS
#  endif  /* NBQ_GRID_SLOW */

#  if !defined NBQ_GRID_SLOW
#   undef dZdxq_u
#   undef dZdxq_w
#  endif

#  define FX dthetadiv_nbqdz_uv
      do k=1,N
          if (IstrU.le.Iend) then
            do j=Jstr,Jend
              do i=Istr,Iend+1
                FX(i,j,k)=-pm_u(i,j)*
     &       (dZdxq_w(i,j,k )-dZdxq_w(i,j,k-1))
#  ifdef MASKING
         FX(i,j,k)=FX(i,j,k)*umask(i,j)
#  endif                
              enddo
            enddo

            do j=Jstr,Jend
              do i=Istr,Iend
                thetadiv_nbq(i,j,k)=
     &  FX(i,j,k)+FX(i+1,j,k)
              enddo
            enddo

          else ! IstrU.gt.Iend

            do j=Jstr,Jend
              do i=Istr,Iend
                thetadiv_nbq(i,j,k)=0.
              enddo
            enddo

          endif 
      enddo
#  undef FX

!
!-----------------------------------------------------------------------
!  Y-component dZdy*qdmv for horizontal divergence
!-----------------------------------------------------------------------
!
#  if !defined NBQ_GRID_SLOW
#   define dZdyq_v dthetadiv_nbqdz
!#   define dZdyq_w dthetadiv_nbqdz_uv
#  endif

#  ifdef NBQ_GRID_SLOW
        if (NSTEP_DS) then
#  endif
      do k=1,N
        do j=Jstr,Jend+1
          do i=Istr,Iend
            dZdyq_v(i,j,k)=(z_r(i,j,k)-z_r(i,j-1,k)) 
     &                       *qdmv_nbq(i,j,k)  ! (dZdy * (rho v))_v
          enddo
        enddo
      enddo

! Bottom boundary conditions
#  ifdef NBQ_FREESLIP
            do j=Jstr,Jend+1
              do i=Istr,Iend
                 dZdyq_w(i,j,0 )= (z_w(i,j,0)-z_w(i,j-1,0))
     &                             *qdmv_nbq(i,j,1) 
     &                             /(Hzr(i,j,1)+Hzr(i,j-1,1))
              enddo
            enddo
            do j=Jstr,Jend
              do i=Istr,Iend   
                qdmw_nbq(i,j,0)=qdmw_nbq(i,j,0)
     &                          +0.5*(dZdyq_w(i,j  ,0)*pn_v(i,j  )  
     &                               +dZdyq_w(i,j+1,0)*pn_v(i,j+1) )
     &                                                  * Hzr(i,j,1)
#   ifdef MASKING
                qdmw_nbq(i,j,0)=qdmw_nbq(i,j,0)*rmask(i,j)
#   endif 
              enddo
            enddo

#  else  /* NBQ_FREESLIP */

#   ifndef NBQ_GRID_SLOW   
            do j=Jstr,Jend +1
              do i=Istr,Iend  
                dZdyq_w(i,j,0)=0.  
                qdmw_nbq(i,j,0)=0.
              enddo
            enddo 
#   endif
#  endif /* NBQ_FREESLIP */

! Top boundary conditions

        do j=Jstr,Jend+1
          do i=Istr,Iend
                dZdyq_w(i,j,N )= (z_w(i,j,N)-z_w(i,j-1,N))
     &                           *qdmv_nbq(i,j,N)
     &                           /(Hzr(i,j,N)+Hzr(i,j-1,N))
          enddo
        enddo

      do k=1,N-1
        do j=Jstr,Jend+1
          do i=Istr,Iend
                dZdyq_w(i,j,k )=Hzw_half_nbq_inv_v(i,j,k)
     &                           *(dZdyq_v(i,j,k)+
     &                             dZdyq_v(i,j,k+1)) ! (dZdy * (rho v))_uw/Hzw_v
          enddo 
        enddo
      enddo

#  ifdef NBQ_GRID_SLOW
        else  ! NSTEP_DS
#  endif /* NBQ_GRID_SLOW */      

! Bottom boundary conditions

#  ifdef NBQ_FREESLIP
            do j=Jstr,Jend+1
              do i=Istr,Iend
                dZdyq_w(i,j,0)= (z_w(i,j,0)-z_w(i,j-1,0))
     &                           *qdmv_nbq(i,j,1) 
     &                           /(Hzr(i,j,1)+Hzr(i,j-1,1))
              enddo
            enddo
            do j=Jstr,Jend
              do i=Istr,Iend   
                qdmw_nbq(i,j,0)=qdmw_nbq(i,j,0)
     &                         +0.5*(dZdyq_w(i,j  ,k)*pn_v(i,j  )  
     &                              +dZdyq_w(i,j+1,k)*pn_v(i,j+1) )
     &                                                 * Hzr(i,j,1) 
#   ifdef MASKING
                qdmw_nbq(i,j,0)=qdmw_nbq(i,j,0)*rmask(i,j)
#   endif 
              enddo
            enddo
#  endif
#  ifdef NBQ_GRID_SLOW
        endif ! NSTEP_DS
#  endif /* NBQ_GRID_SLOW */
#  if !defined NBQ_GRID_SLOW
#   undef dZdyq_v
#   undef dZdyq_w
#  endif

#  define FY dthetadiv_nbqdz_uv
          if (JstrV.le.Jend) then
          do k=1,N
            do j=Jstr,Jend+1
              do i=Istr,Iend 
                FY(i,j,k)=-pn_v(i,j)*
     &   (dZdyq_w(i,j,k)-dZdyq_w(i,j,k-1))
#  ifdef MASKING
        FY(i,j,k)=FY(i,j,k)*vmask(i,j)
#  endif                 
              enddo
            enddo
            do j=Jstr,Jend
              do i=Istr,Iend
                thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)
     &    +FY(i,j,k)+FY(i,j+1,k)
              enddo
            enddo
          enddo
          endif
#  undef FY


!
!-----------------------------------------------------------------------
!  Compute total horizontal Divergence divH(qdmH(m+1))
!-----------------------------------------------------------------------
!
#  ifdef OPENACC
      DOEXTEND(k,1,N,FX,FY)
#  else
      do k=1,N !<-- k loop
#  endif

        if (IstrU.le.Iend) then
          do j=Jstr,Jend
            do i=Istr,Iend+1
              FX(i,j)=on_u(i,j)*qdmu_nbq(i,j,k)
            enddo
          enddo
        endif
        if (JstrV.le.Jend) then
          do j=Jstr,Jend+1  
            do i=Istr,Iend
              FY(i,j)=om_v(i,j)*qdmv_nbq(i,j,k)
            enddo
          enddo
        endif
        if (IstrU.gt.Iend) then
          do j=Jstr,Jend
            do i=Istr,Iend  
              thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)
     &                           +pm(i,j)*pn(i,j)*(FY(i,j+1)-FY(i,j))  
#  ifdef MASKING
              thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)*rmask(i,j)
#  endif                              
            enddo
          enddo
        elseif (JstrV.gt.Jend) then
          do j=Jstr,Jend  
            do i=Istr,Iend   
              thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)
     &                            +pm(i,j)*pn(i,j)*(FX(i+1,j)-FX(i,j))  
#  ifdef MASKING
              thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)*rmask(i,j)
#  endif                              
            enddo
          enddo
        else
          do j=Jstr,Jend
            do i=Istr,Iend
              thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)
     &                            +pm(i,j)*pn(i,j)*(FX(i+1,j)-FX(i,j)+
     &                                              FY(i,j+1)-FY(i,j))  
#  ifdef MASKING
              thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)*rmask(i,j)
#  endif                              
            enddo
          enddo
        endif  
#  ifdef OPENACC
      ENDDOEXTEND
#  else
      enddo ! <-- k=1,N
#  endif
!
!-----------------------------------------------------------------------
!  thetadiv2_nbq: complet time-corrective term  (dh/dt included) 
!  thetadiv3_nbq: reduced time-corrective term  (no dh/dt)
!-----------------------------------------------------------------------
!
#  ifdef NBQ_GRID_SLOW
       if (FIRST_FAST_STEP) then
#  endif

        do k=0,N
          do j=JstrV-2,Jend+1
            do i=IstrU-2,Iend+1
              zw_nbq(i,j,k,knew)=z_w(i,j,k)
            enddo
          enddo
        enddo

#  define FC qdmw_nbq_old
#  define CF dthetadiv_nbqdz_uv

        do j=JstrV-2,Jend+1 !<-- j loop
          do i=IstrU-2,Iend+1
#  ifdef NBQ_HZ_PROGNOSTIC
            FC(i,j,0)=0.    ! Bottom boundary condition
#  endif
            CF(i,j,0)=0.
          enddo
        enddo

        do k=1,N-1
          do j=JstrV-2,Jend+1
            do i=IstrU-2,Iend+1
              FC(i,j,k) =
     &          -(zw_nbq(i,j,k,knew)-zw_nbq(i,j,k,kstp))/dtgrid_nbq
              CF(i,j,k)=   
     &          FC(i,j,k) 
     &          *0.5*(    rho_grd(i,j,k  ) 
     &                   +rho_nbq(i,j,k  )*Hzr_half_nbq_inv(i,j,k  ) 
     &                   +rho_grd(i,j,k+1)
     &                   +rho_nbq(i,j,k+1)*Hzr_half_nbq_inv(i,j,k+1) )

#  ifdef NBQ_HZ_PROGNOSTIC
             FC(i,j,k) = CF(i,j,k)-FC(i,j,k)
#  endif
            enddo
          enddo
        enddo

        do j=JstrV-2,Jend+1
          do i=IstrU-2,Iend+1
#  ifdef NBQ_HZ_PROGNOSTIC
             FC(i,j,N)=   
     &        -(zw_nbq(i,j,N,knew)-zw_nbq(i,j,N,kstp))/dtgrid_nbq
     &         *( 1.+rho_grd(i,j,N)
     &          )
#  endif
            CF(i,j,N)=  
     &        -(zw_nbq(i,j,N,knew)-zw_nbq(i,j,N,kstp))/dtgrid_nbq
     &         *( rho_grd(i,j,N)
     &          )
          enddo
        enddo !<-- j loop

        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
              thetadiv3_nbq(i,j,k)=
     &    CF(i,j,k)-CF(i,j,k-1) 
#  ifdef NBQ_HZ_PROGNOSTIC
              thetadiv2_nbq(i,j,k)=
     &    FC(i,j,k)-FC(i,j,k-1) 
#  endif
            enddo
          enddo
        enddo
#  undef FC
#  undef CF

#  ifdef NBQ_GRID_SLOW
       endif !<-- FIRST_FAST_STEP
#  endif
!
!-------------------------------------------------------------------
!  Solve implicit W-momentum equation
!-------------------------------------------------------------------
!
#  ifdef NBQ_IMP
      do k=0,N
        do j=Jstr,Jend
          do i=Istr,Iend
            qdmw_nbq_old(i,j,k)=qdmw_nbq(i,j,k)
          enddo
        enddo
      enddo

      !do j=Jstr,Jend ! <-- j loop

      DOLOOP2D(Istr,Iend,Jstr,Jend)

        do k=1,N
          do i=Istr,Iend
            FC(i,k)= soundspeed2_nbq(i,j)*rho_nbq(i,j,k)
     &              -(thetaimp_nbq*soundspeed2_nbq(i,j)*dtnbq+visc2_nbq)
     &                       *(thetadiv_nbq(i,j,k)+thetadiv3_nbq(i,j,k))
    
#   ifdef NBQ_THETAIMP 
            FC(i,k)=FC(i,k)
     &       -thetaimp_nbq*(1.-thetaimp_nbq)*soundspeed2_nbq(i,j)*dtnbq*
     &                   ( Hzw_half_nbq_inv(i,j,k  )*qdmw_nbq(i,j,k)
     &                    -Hzw_half_nbq_inv(i,j,k-1)*qdmw_nbq(i,j,k-1) )
#   endif

            FC(i,k)=FC(i,k)*Hzr_half_nbq_inv(i,j,k) 
          enddo
        enddo
     
!.........Inner layers

        do k=1,N-1
          do i=Istr,Iend   
            dum_s = FC(i,k) - FC(i,k+1)
#   ifdef UV_COR_NT
            dum_s = dum_s + ntcorw(i,j,k)
#   endif
            qdmw_nbq(i,j,k) = qdmw_nbq(i,j,k)   
     &                        + dtnbq * ( dum_s + rw_int_nbq(i,j,k) )
#   ifdef MASKING
            qdmw_nbq(i,j,k) = qdmw_nbq(i,j,k) * rmask(i,j)
#   endif               
          enddo             
        enddo

!.........Surface BC

        k=N
        do i=Istr,Iend
          dum_s =   FC(i,k)
#   ifdef UV_COR_NT
          dum_s = dum_s + ntcorw(i,j,k)
#   endif
          qdmw_nbq(i,j,k) = qdmw_nbq(i,j,k)   
     &                      + dtnbq * ( dum_s + rw_int_nbq(i,j,k) )
#   ifdef MASKING
          qdmw_nbq(i,j,k) = qdmw_nbq(i,j,k) * rmask(i,j)
#   endif              
        enddo   

    !  enddo !<-- j loop

        ENDDOLOOP2D

!---------------------
! Gaussian Elimination
!---------------------

!.......Compute coef.

#   if !defined VARYING_CSOUND
      cff1_cst = 1./(dtnbq*(thetaimp_nbq**2*csound_nbq**2
     &            *dtnbq+visc2_nbq))
#   endif
      !do j=Jstr,Jend  !<-- j loop
      DOLOOP2D(Istr,Iend,Jstr,Jend)

      do k=1,N
        do i=Istr,Iend
          FC(i,k) =
     &      Hzw_half_nbq_inv(i,j,k)*Hzr_half_nbq_inv(i,j,k)
          BC(i,k) =
     &      Hzw_half_nbq_inv(i,j,k-1)*Hzr_half_nbq_inv(i,j,k)
        enddo
      enddo

!..........Bottom BC
        k=1
        do i=Istr,Iend
#   if !defined VARYING_CSOUND
          cff1 = cff1_cst
#   else
          cff1=1./(dtnbq*(thetaimp_nbq**2*soundspeed2_nbq(i,j)
     &            *dtnbq+visc2_nbq)) 
#   endif
          cff=(cff1+FC(i,1)+BC(i,2))
          CF(i,1)=(-FC(i,2))
     &                                                          /cff
          DC(i,1)=qdmw_nbq(i,j,1)*cff1/cff   
     &           +qdmw_nbq(i,j,0)/cff*BC(i,1)
        enddo

!..........Inner layers

        do k=2,N-1
          do i=Istr,Iend
#   if !defined VARYING_CSOUND
          cff1 = cff1_cst
#   else
          cff1=1./(dtnbq*(thetaimp_nbq**2*soundspeed2_nbq(i,j)
     &            *dtnbq+visc2_nbq)) 
#   endif
            cff=(cff1+
     &               FC(i,k)+BC(i,k+1)
     &              +BC(i,k)*CF(i,k-1))
            CF(i,k)=(-FC(i,k+1))/cff
            DC(i,k)=(qdmw_nbq(i,j,k)*cff1+
     &               BC(i,k)*DC(i,k-1)) /cff           
          enddo            
        enddo

!..........Surface BC
        k=N
        do i=Istr,Iend
#   if !defined VARYING_CSOUND
          cff1 = cff1_cst
#   else
          cff1=1./(dtnbq*(thetaimp_nbq**2*soundspeed2_nbq(i,j)
     &            *dtnbq+visc2_nbq)) 
#   endif
          cff=(cff1+FC(i,N)
     &             +BC(i,N)*CF(i,N-1))  
          CF(i,N)=0. 
          DC(i,k)=(qdmw_nbq(i,j,N)*cff1+
     &             BC(i,N)*DC(i,N-1))/cff
        enddo 

!..........Solves tri-diag system
        do i=Istr,Iend
          qdmw_nbq(i,j,N)=DC(i,k)   
#   ifdef NBQ_GRAV
     &                    -rho_nbq(i,j,N)*0.5*g*dtnbq
#   endif     
        enddo

        do k=N-1,1,-1
          do i=Istr,Iend
            qdmw_nbq(i,j,k)=DC(i,k)-CF(i,k)*qdmw_nbq(i,j,k+1)
#   ifdef NBQ_GRAV
     &            -0.25*(rho_nbq(i,j,k)*Hzr_half_nbq_inv(i,j,k)
     &                 +rho_nbq(i,j,k+1)*Hzr_half_nbq_inv(i,j,k+1))
     &                *(Hzr(i,j,k)+Hzr(i,j,k+1))
     &                 *g*dtnbq
#   endif
          enddo            
        enddo                        

      !enddo   !<-- j loop 
      ENDDOLOOP2D
!
!-------------------------------------------------------------------
!  W-momentum open boundary conditions
!-------------------------------------------------------------------
!
#   ifdef OBC_NBQ
      call wnbq_bc_tile (Istr,Iend,Jstr,Jend, work)
#   endif

#   ifdef UV_COR_NT
#    if defined BD_EXCHANGE
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
     &                        qdmw_nbq(START_2D_ARRAY,0))
#    endif
#   endif

#  endif /* NBQ_IMP */
!
!-------------------------------------------------------------------
!  Mass conservation equation
!
!      div(m+1) = div(m+1) + divZ(qdmw_nbq(m+1))
!-------------------------------------------------------------------
!		
!.......Computes fluxes:  
!
    !  do j=Jstr,Jend  !<-- j loop

        DOLOOP2D(Istr,Iend,Jstr,Jend)
            
#  ifdef NBQ_FREESLIP
        do i=Istr,Iend
          FC(i,0)=Hzw_half_nbq_inv(i,j,0) *
     &            (thetaimp_nbq*qdmw_nbq(i,j,0)
     &           +(1.-thetaimp_nbq)*qdmw_nbq_old(i,j,0)) ! Bottom BC
        enddo
#  else
        do i=Istr,Iend
          FC(i,0)=0.   ! Bottom BC
        enddo
#  endif
        do k=1,N-1
          do i=Istr,Iend
            FC(i,k)=Hzw_half_nbq_inv(i,j,k) *
     &              (thetaimp_nbq*qdmw_nbq(i,j,k)
     &              +(1.-thetaimp_nbq)*qdmw_nbq_old(i,j,k))
            thetadiv_nbq(i,j,k)=thetadiv_nbq(i,j,k)
     &                         +FC(i,k)-FC(i,k-1)  
          enddo
        enddo
        do i=Istr,Iend
          FC(i,N)=Hzw_half_nbq_inv(i,j,N) *
     &            (thetaimp_nbq*qdmw_nbq(i,j,N)
     &           +(1.-thetaimp_nbq)*qdmw_nbq_old(i,j,N))
          thetadiv_nbq(i,j,N)=thetadiv_nbq(i,j,N)
     &                        +FC(i,N)-FC(i,N-1)    
        enddo

      !enddo  !<-- j loop
        ENDDOLOOP2D
!                 
#  if defined BD_EXCHANGE
!$acc end kernels
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        thetadiv_nbq(START_2D_ARRAY,1))
#   ifdef NBQ_HZ_PROGNOSTIC
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        thetadiv2_nbq(START_2D_ARRAY,1))
#   endif
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        thetadiv3_nbq(START_2D_ARRAY,1))
!$acc kernels if(compute_on_device) default(present)
#  endif
!
!=====================================================================
!
!      Solve rho_nbq from mass conservation 
!
!      using a Forward-Backward scheme:
!          rho_nbq(m+1) = rho_nbq(m) - dtfast*DIV(m+1)
!
!=====================================================================
!
      do k=1,N
        do j=JstrV-2,Jend+1
          do i=IstrU-2,Iend+1
            rho_nbq(i,j,k) = rho_nbq(i,j,k)  
     &                       - dtfast*(thetadiv_nbq(i,j,k)+
     &                                 thetadiv3_nbq(i,j,k))
          enddo
        enddo
      enddo

!     call rnbq_bc_tile(Istr,Iend,Jstr,Jend, work) !!! leave commented
!
!# if defined BD_EXCHANGE
!      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
!     &                       rho_nbq(START_2D_ARRAY,1))
!# endif
!
!-------------------------------------------------------------------
!  Acoustic wave emission
!-------------------------------------------------------------------
!
#  ifdef ACOUSTIC
        time_nbq = time_nbq + dtfast
        do k=1,N
          do j=JstrV-2,Jend+1
            do i=IstrU-2,Iend+1
              dist_d=sqrt((xr(i,j)-xl/2.)**2+(0.*(yr(i,j)-el/2.))**2
     &                             +(abs(z_w(i,j,k))-hmax_exp/2.)**2)

              rho_nbq(i,j,k)=rho_nbq(i,j,k)+amp_exp
     &                         *sin(2*acos(-1.)*time_nbq/period_exp)
     &                                 *exp(-dist_d**2/for_a_exp**2)
     &                                                   *Hzr(i,j,k)
            enddo
          enddo
        enddo
#  endif
!
!-------------------------------------------------------------------
!  rhobar_nbq: depth-mean density (/rho0)
!-------------------------------------------------------------------
!		
#  ifdef NBQ_MASS
!
! Compute rhobar_nbq(m+1) used in zeta diagnostic from 
! depth-integrated continuity equation
!
      do j=Jstr-1,Jend+1
        do i=Istr-1,Iend+1
          rhobar_nbq(i,j,knew)=0.
        enddo  
      enddo
      do k=1,N
        do j=Jstr-1,Jend+1
          do i=Istr-1,Iend+1
            rhobar_nbq(i,j,knew)= rhobar_nbq(i,j,knew)
     &                             +rho_nbq(i,j,k)
     &                             +rho_grd(i,j,k)*Hzr(i,j,k)
          enddo  
        enddo  
      enddo
      do j=Jstr-1,Jend+1
        do i=Istr-1,Iend+1
          rhobar_nbq(i,j,knew) = 1.+(rhobar_nbq(i,j,knew)) 
     &                            / (z_w(i,j,N)-z_w(i,j,0))
        enddo
      enddo
!      
! LAURENT: the next exchange should not be needed
!
#   if defined BD_EXCHANGE
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend, 
     &                           rhobar_nbq(START_2D_ARRAY,knew))
#   endif
#   ifdef CVTK_DEBUG_ADVANCED
!      call check_tab2d(rhobar_nbq(:,:,knew),'rhobar_nbq','r')
#   endif    
#  endif /* NBQ_MASS */

# endif /* NBQ */
!
!
!***********************************************************************
!
!
!                     FAST-MODE POST-PROCESSING
!
!
!***********************************************************************
!
!=====================================================================
!  Update total mass of water volume Dnew 
!=====================================================================
!
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          Dnew(i,j)=(zeta(i,j,knew)+h(i,j))
# ifdef NBQ_MASS
     &               *rhobar_nbq(i,j,knew)
# endif
        enddo
      enddo
            
# if defined CVTK_DEBUG_ADVANCED && defined NBQ
!$acc end kernels
C$OMP BARRIER
C$OMP MASTER
      call check_tab3d(rho_nbq,'rho_nbq','r',
     &   ondevice=.TRUE.)
C$OMP END MASTER
!$acc kernels if(compute_on_device) default(present)
# endif    

!=====================================================================
!  Update W NBQ forcing for internal mode (rw_nbq ~ rw_nbq_avg1)
!  Note: here rw_nbq contains qdmw_nbq(m)
!=====================================================================
!
# ifdef NBQ
      if (LAST_FAST_STEP) then
        do k=0,N 
          do j=Jstr,Jend              
            do i=Istr,Iend
              rw_nbq(i,j,k)=((qdmw_nbq(i,j,k)-rw_nbq(i,j,k))
     &                              /dtnbq-ndtnbq*rw_int_nbq(i,j,k))
     &                                            /(pm(i,j)*pn(i,j))
            enddo
          enddo
        enddo
      endif

!# if defined CVTK_DEBUG && defined NBQ
!      call check_tab3d(rw_nbq,'rw_nbq','rint')
!#  endif  
# endif
!
!=====================================================================
! Get filtered RHS terms
! and multiply by dx*dy to get units of rho*Hz*dx*dy*ru
!=====================================================================
!
      if (LAST_FAST_STEP) then
!
!-----------------------------------------------------------------------
! Store average fields AVG1 of rho and rhobar
!-----------------------------------------------------------------------
!     
# ifdef NBQ_MASS
        do j=Jstr,Jend
          do i=Istr,Iend
            rhobar_nbq_avg1(i,j)=rhobar_nbq(i,j,knew)
          enddo
        enddo 
        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
              rho_nbq_avg1(i,j,k)=1.d0
     &                            + ( rho_nbq(i,j,k)/Hzr(i,j,k)
     &                               +rho_grd(i,j,k) )  
            enddo
          enddo 
        enddo
# endif /* NBQ_MASS */
!
!-----------------------------------------------------------------------
!  Compute average fields AVG2 of RHS NBQ forcing
!  Note: here ru_nbq_avg2, ru_nbq_2d_old ... are working arrays
!-----------------------------------------------------------------------
!
        do k=1,N
          do j=Jstr,Jend
            do i=IstrU,Iend
!             ru_nbq_avg1(i,j,k)= ru_nbq(i,j,k)
# ifdef NBQ
              ru_int_nbq(i,j,k) = ru_int_nbq(i,j,k)
     &                      -ru_ext_nbq_old(i,j)*(Hz(i-1,j,k)+Hz(i,j,k))
# endif
              ru_nbq_avg2(i,j,k)=
     &                          ((qdmu_nbq(i,j,k)-ru_nbq_avg2(i,j,k))/dt
     &                   -ru_int_nbq(i,j,k)-(ru_ext_nbq_sum(i,j)/nfast)*
     &                      (Hz(i,j,k)+Hz(i-1,j,k)))*on_u(i,j)*om_u(i,j)
            enddo
          enddo 
        enddo          

        do k=1,N
          do j=JstrV,Jend
            do i=Istr,Iend             
!             rv_nbq_avg1(i,j,k)= rv_nbq(i,j,k) 
# ifdef NBQ
              rv_int_nbq(i,j,k) = rv_int_nbq(i,j,k)
     &                      -rv_ext_nbq_old(i,j)*(Hz(i,j-1,k)+Hz(i,j,k))
# endif
              rv_nbq_avg2(i,j,k)=
     &                          ((qdmv_nbq(i,j,k)-rv_nbq_avg2(i,j,k))/dt
     &                   -rv_int_nbq(i,j,k)-(rv_ext_nbq_sum(i,j)/nfast)*
     &                      (Hz(i,j,k)+Hz(i,j-1,k)))*on_v(i,j)*om_v(i,j)
            enddo
          enddo 
        enddo

# ifdef NBQ
        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
!             rw_nbq_avg1(i,j,k)= rw_nbq(i,j,k)
              rw_nbq_avg2(i,j,k)=
     &                          ((qdmw_nbq(i,j,k)-rw_nbq_avg2(i,j,k))/dt
     &                           -rw_int_nbq(i,j,k))*on_r(i,j)*om_r(i,j)
            enddo
          enddo 
        enddo
# endif

      endif ! LAST_FAST_STEP
!
!-----------------------------------------------------------------------
! Dismiss coupling of NBQ, NBQ2EXT & NBQ2INT for testing
!-----------------------------------------------------------------------
!
# ifdef NBQ_NOCOUPLING 
      ru_nbq      =0.   ! 3D
      rv_nbq      =0.
      ru_nbq_avg2 =0.
      rv_nbq_avg2 =0.
#  ifdef NBQ
      rw_nbq      =0.
      rw_nbq_avg2 =0.
#  endif
#  ifdef NBQ_MASS
      rhobar_nbq  =1.
      rho_nbq     =0.
      rho_nbq_avg1=1.
#  endif
# endif
!
!-----------------------------------------------------------------------
! Exchange NBQ coupling
!-----------------------------------------------------------------------
!
# if defined BD_EXCHANGE
!$acc end kernels
#  ifdef NBQ_MASS
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                        rhobar_nbq_avg1(START_2D_ARRAY))
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        rho_nbq_avg1(START_2D_ARRAY,1))
#  endif
#  ifndef NBQ
      IF (LAST_FAST_STEP) then
#  endif
!     call exchange_u3d_tile (Istr,Iend,Jstr,Jend,  
!    &                        ru_nbq_avg2(START_2D_ARRAY,1))
!     call exchange_v3d_tile (Istr,Iend,Jstr,Jend,  
!    &                        rv_nbq_avg2(START_2D_ARRAY,1))
#  ifndef NBQ
      endif
#  endif
#  ifdef NBQ
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend,  
     &                        rw_nbq_avg2(START_2D_ARRAY,0))
#  endif 
!$acc kernels if(compute_on_device) default(present)
# endif
!
!=====================================================================
!  Depth-averaged velocity & forcing from fast mode
!=====================================================================
!
! Output: (ubar,vbar), (DU_avg1,DV_avg1)
!
# define Dstp DUon

      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          Dstp(i,j)=zeta(i,j,kstp)+h(i,j)
        enddo
      enddo

!$acc end kernels

      cff=2*dtfast
      cff1=0.5*weight(1,iif)
      cff2=0.5*weight(2,iif)

!$acc kernels if(compute_on_device) default(present)
      if (iif.eq.1) then
         DU_avg1=0.
         DV_avg1=0.
         DU_avg2=0.
         DV_avg2=0.
      endif

      do j=Jstr,Jend
        do i=IstrU,Iend
# ifdef NBQ
          DUnew=DU_nbq(i,j) *2.
# else
          DUnew=( (Dstp(i,j)+Dstp(i-1,j))*ubar(i,j,kstp)
     &                             + cff * rubar_nbq(i,j) )
#  ifdef MASKING
     &                                     *umask(i,j)
#  endif
# endif
# ifdef WET_DRY
     &                                     *umask_wet(i,j)
# endif
          ubar(i,j,knew)=DUnew/(Dnew(i,j)+Dnew(i-1,j))
          DU_avg1(i,j,nnew)=DU_avg1(i,j,nnew) 
     &                             +cff1*on_u(i,j)*( DUnew
# ifdef MRL_WCI
     &                +(Dnew(i,j)+Dnew(i-1,j))*ust2d(i,j)
#  ifdef WET_DRY
     &  *umask_wet(i,j)
#  endif
# endif
     &                                                   )
          DU_avg2(i,j)=DU_avg2(i,j)+cff2*on_u(i,j)*( DUnew 
# ifdef MRL_WCI
     &                +(Dnew(i,j)+Dnew(i-1,j))*ust2d(i,j)
#  ifdef WET_DRY
     & *umask_wet(i,j)
#  endif
# endif
     &                                                   )
        enddo
      enddo 
  
      do j=JstrV,Jend
        do i=Istr,Iend
# ifdef NBQ
          DVnew=DV_nbq(i,j) *2.
# else
          DVnew=( (Dstp(i,j)+Dstp(i,j-1))*vbar(i,j,kstp)
     &                             + cff * rvbar_nbq(i,j) )
#  ifdef MASKING
     &                                     *vmask(i,j)
#  endif
# endif
# ifdef WET_DRY
     &                                     *vmask_wet(i,j)
# endif
          vbar(i,j,knew)=DVnew/(Dnew(i,j)+Dnew(i,j-1))
          DV_avg1(i,j,nnew)=DV_avg1(i,j,nnew) 
     &                              +cff1*om_v(i,j)*(DVnew
# ifdef MRL_WCI
     &                +(Dnew(i,j)+Dnew(i,j-1))*vst2d(i,j)
#  ifdef WET_DRY
     & *vmask_wet(i,j)
#  endif
# endif
     &                                                   )
          DV_avg2(i,j)=DV_avg2(i,j)+cff2*om_v(i,j)*( DVnew
# ifdef MRL_WCI
     &                +(Dnew(i,j)+Dnew(i,j-1))*vst2d(i,j)
#  ifdef WET_DRY
     & * vmask_wet(i,j)
#  endif
# endif
     &                                                   )
        enddo
      enddo
!
!-----------------------------------------------------------------------
! Apply point sources for hydrostatic case
!-----------------------------------------------------------------------
!
# if defined PSOURCE && !defined NBQ
      do is=1,Nsrc 
#  ifdef MPI
        i=Isrc_mpi(is,mynode)
        j=Jsrc_mpi(is,mynode)
#  else
        i=Isrc(is)
        j=Jsrc(is)
#  endif
        if (IstrR.le.i .and. i.le.IendR .and.
     &      JstrR.le.j .and. j.le.JendR) then
          if (Dsrc(is).eq.0) then
            ubar(i,j,knew)=2.*Qbar(is)/( on_u(i,j)
     &                       *(Dnew(i-1,j)+Dnew(i,j)) )
            DU_avg1(i,j,nnew)=Qbar(is)
          else
            vbar(i,j,knew)=2.*Qbar(is)/( om_v(i,j)
     &                       *(Dnew(i,j-1)+Dnew(i,j)) )
            DV_avg1(i,j,nnew)=Qbar(is)
          endif
        endif
      enddo
# endif
!
!-----------------------------------------------------------------------
!  Set 2D Momemtum nudging
!-----------------------------------------------------------------------
!
# if defined M2NUDGING && defined M2CLIMATOLOGY

#  ifdef ZONAL_NUDGING
      if (FIRST_TIME_STEP .or. mod(iic,10).eq.0) then
        if (FIRST_FAST_STEP) then
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   ubar(START_2D_ARRAY,knew),ubzon)
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   vbar(START_2D_ARRAY,knew),vbzon)
        endif
      endif
      if (FIRST_TIME_STEP) then
        if (FIRST_FAST_STEP) then 
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   ubclm(START_2D_ARRAY),ubclmzon)
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   vbclm(START_2D_ARRAY),vbclmzon)
        endif
      endif
#  endif /* ZONAL_NUDGING */

      do j=Jstr,Jend
        do i=IstrU,Iend
#  ifdef ZONAL_NUDGING        
          DUnew = dtfast*M2nudgcof(i,j)*(ubclmzon(j)-ubzon(j))
#  else          
          DUnew = dtfast*M2nudgcof(i,j)*(ubclm(i,j)-ubar(i,j,knew))
#  endif
#  ifdef MASKING
     &                 * umask(i,j)
#  endif 
#  ifdef WET_DRY
     &                 * umask_wet(i,j)
#  endif
          ubar(i,j,knew)=ubar(i,j,knew) + DUnew
          DU_avg1(i,j,nnew)=DU_avg1(i,j,nnew) +cff1*DUnew*
     &                         (Dnew(i,j)+Dnew(i-1,j))*on_u(i,j)
          DU_avg2(i,j)     =DU_avg2(i,j)      +cff2*DUnew*
     &                         (Dnew(i,j)+Dnew(i-1,j))*on_u(i,j)
        enddo
      enddo
      
      do j=JstrV,Jend
        do i=Istr,Iend
#  if defined ZONAL_NUDGING  
          DVnew = dtfast*M2nudgcof(i,j)*(vbclmzon(j)-vbzon(j)) 
#  else
          DVnew = dtfast*M2nudgcof(i,j)*(vbclm(i,j)-vbar(i,j,knew))    
#  endif      
#  ifdef MASKING
     &                 * vmask(i,j)
#  endif 
#  ifdef WET_DRY
     &                 * vmask_wet(i,j)
#  endif
          vbar(i,j,knew)=vbar(i,j,knew) + DVnew
          DV_avg1(i,j,nnew)=DV_avg1(i,j,nnew) +cff1*DVnew*
     &                         (Dnew(i,j)+Dnew(i,j-1))*om_v(i,j)
          DV_avg2(i,j)     =DV_avg2(i,j)      +cff2*DVnew*
     &                         (Dnew(i,j)+Dnew(i,j-1))*om_v(i,j)
        enddo
      enddo
# endif /* M2NUDGING */
!
!----------------------------------------------------------------------- 
! Set boundary conditions and compute integral mass flux accross
! all open boundaries, if any.
!-----------------------------------------------------------------------
!
# ifdef NBQ
      M2bc_nbq_flag=.false.  ! skip wet/dry conditions
                             ! and DU_nbq computation
# endif
!$acc end kernels
      call u2dbc_tile (Istr,Iend,Jstr,Jend, work) 
      call v2dbc_tile (Istr,Iend,Jstr,Jend, work)

!$acc kernels if(compute_on_device) default(present)
# ifdef WET_DRY
#  ifndef EW_COM_PERIODIC
      if (WESTERN_EDGE) then
        DO j=Jstr,Jend
          ubar(Istr,j,knew)=ubar(Istr,j,knew)*umask_wet(Istr,j)
#   ifdef MRL_WCI
          ust2d(Istr,j)=ust2d(Istr,j)*umask_wet(Istr,j)
#   endif
        END DO
        DO j=JstrV,Jend
          vbar(Istr-1,j,knew)=vbar(Istr-1,j,knew)*vmask_wet(Istr-1,j)
#   ifdef MRL_WCI
          vst2d(Istr-1,j)=vst2d(Istr-1,j)*vmask_wet(Istr-1,j)
#   endif
        END DO
      END IF
      if (EASTERN_EDGE) then
        DO j=Jstr,Jend
          ubar(Iend+1,j,knew)=ubar(Iend+1,j,knew)*umask_wet(Iend+1,j)
#   ifdef MRL_WCI
          ust2d(Iend+1,j)=ust2d(Iend+1,j)*umask_wet(Iend+1,j)
#   endif
        END DO
        DO j=JstrV,Jend
          vbar(Iend+1,j,knew)=vbar(Iend+1,j,knew)*vmask_wet(Iend+1,j)
#   ifdef MRL_WCI
          vst2d(Iend+1,j)=vst2d(Iend+1,j)*vmask_wet(Iend+1,j)
#   endif
        END DO
      END IF
#  endif
#  ifndef NS_COM_PERIODIC
      if (SOUTHERN_EDGE) then
        DO i=IstrU,Iend
          ubar(i,Jstr-1,knew)=ubar(i,Jstr-1,knew)*umask_wet(i,Jstr-1)
#   ifdef MRL_WCI
          ust2d(i,Jstr-1)=ust2d(i,Jstr-1)*umask_wet(i,Jstr-1)
#   endif
        END DO
        DO i=IstrU,Iend
          vbar(i,Jstr,knew)=vbar(i,Jstr,knew)*vmask_wet(i,Jstr)
#   ifdef MRL_WCI
          vst2d(i,Jstr)=vst2d(i,Jstr)*vmask_wet(i,Jstr)
#   endif
        END DO
      END IF
      if (NORTHERN_EDGE) then
        DO i=Istr,Iend
          ubar(i,Jend+1,knew)=ubar(i,Jend+1,knew)*umask_wet(i,Jend+1)
#   ifdef MRL_WCI
          ust2d(i,Jend+1)=ust2d(i,Jend+1)*umask_wet(i,Jend+1)
#   endif
        END DO
        DO i=Istr,Iend
          vbar(i,Jend+1,knew)=vbar(i,Jend+1,knew)*vmask_wet(i,Jend+1)
#   ifdef MRL_WCI
          vst2d(i,Jend+1)=vst2d(i,Jend+1)*vmask_wet(i,Jend+1)
#   endif
        END DO
      END IF
#  endif
# endif
!
! zeta vill be recomputed via depth-integrated continuity equation
!
# if defined BD_EXCHANGE
!$acc end kernels
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        ubar(START_2D_ARRAY,knew))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        vbar(START_2D_ARRAY,knew))
!$acc kernels if(compute_on_device) default(present)
# endif
          
# ifdef OBC_VOLCONS
      call obc_flux_tile (Istr,Iend,Jstr,Jend)
# endif

!      call unbq_bc_tile (Istr,Iend,Jstr,Jend, work)
!      call vnbq_bc_tile (Istr,Iend,Jstr,Jend, work)
!#  if defined BD_EXCHANGE
!      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
!     &                        qdmu_nbq(START_2D_ARRAY,1))
!      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
!     &                        qdmv_nbq(START_2D_ARRAY,1))
!#  endif
!
!-----------------------------------------------------------------------
! Compute fast-time averaged barotropic mass fluxes along physical
! boundaries.
!-----------------------------------------------------------------------
!
# ifndef EW_PERIODIC
      if (WESTERN_EDGE) then
        do j=Jstr-1,JendR
          Dnew(Istr-1,j)=(h(Istr-1,j)+zeta(Istr-1,j,knew))
#  ifdef NBQ_MASS
     &                          *rhobar_nbq(Istr-1,j,knew)
#  endif
        enddo
      endif
      if (EASTERN_EDGE) then
        do j=Jstr-1,JendR
          Dnew(Iend+1,j)=(h(Iend+1,j)+zeta(Iend+1,j,knew))
#  ifdef NBQ_MASS
     &                          *rhobar_nbq(Iend+1,j,knew)
#  endif
        enddo
      endif
# endif /* !EW_PERIODIC */

# ifndef NS_PERIODIC
      if (SOUTHERN_EDGE) then
        do i=Istr-1,IendR
          Dnew(i,Jstr-1)=(h(i,Jstr-1)+zeta(i,Jstr-1,knew))
#  ifdef NBQ_MASS
     &                          *rhobar_nbq(i,Jstr-1,knew)
#  endif
        enddo
      endif
      if (NORTHERN_EDGE) then
        do i=Istr-1,IendR
          Dnew(i,Jend+1)=(h(i,Jend+1)+zeta(i,Jend+1,knew))
#  ifdef NBQ_MASS
     &                          *rhobar_nbq(i,Jend+1,knew)
#  endif
        enddo
      endif
# endif /* !NS_PERIODIC */
!$acc end kernels

      cff1=0.5*weight(1,iif)
      cff2=0.5*weight(2,iif)

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

# ifndef EW_PERIODIC
      if (WESTERN_EDGE) then
        do j=JstrR,JendR
          cff=(Dnew(IstrU-1,j)+Dnew(IstrU-2,j))*(ubar(IstrU-1,j,knew)
#  ifdef MRL_WCI
     &                                              +ust2d(IstrU-1,j)
#  endif
     &                                              )*on_u(IstrU-1,j)
          DU_avg1(IstrU-1,j,nnew)=DU_avg1(IstrU-1,j,nnew)+cff1*cff
          DU_avg2(IstrU-1,j)=DU_avg2(IstrU-1,j)+cff2*cff
        enddo
        do j=JstrV,Jend
          cff=(Dnew(Istr-1,j)+Dnew(Istr-1,j-1) )*(vbar(Istr-1,j,knew)
#  ifdef MRL_WCI
     &                                               +vst2d(Istr-1,j)
#  endif
     &                                               )*om_v(Istr-1,j)
          DV_avg1(Istr-1,j,nnew)=DV_avg1(Istr-1,j,nnew)+cff1*cff
          DV_avg2(Istr-1,j)=DV_avg2(Istr-1,j)+cff2*cff
        enddo
      endif
      if (EASTERN_EDGE) then
        do j=JstrR,JendR
          cff=(Dnew(Iend+1,j)+Dnew(Iend,j))*(ubar(Iend+1,j,knew)
#  ifdef MRL_WCI
     &                                          +ust2d(Iend+1,j)
#  endif
     &                                          )*on_u(Iend+1,j)
          DU_avg1(Iend+1,j,nnew)=DU_avg1(Iend+1,j,nnew)+cff1*cff
          DU_avg2(Iend+1,j)=DU_avg2(Iend+1,j)+cff2*cff
        enddo
        do j=JstrV,Jend
          cff=(Dnew(Iend+1,j)+Dnew(Iend+1,j-1))*(vbar(Iend+1,j,knew)
#  ifdef MRL_WCI
     &                                              +vst2d(Iend+1,j)
#  endif
     &                                              )*om_v(Iend+1,j)
          DV_avg1(Iend+1,j,nnew)=DV_avg1(Iend+1,j,nnew)+cff1*cff
          DV_avg2(Iend+1,j)=DV_avg2(Iend+1,j)+cff2*cff
        enddo
      endif
# endif /* !EW_PERIODIC */
# ifndef NS_PERIODIC
      if (SOUTHERN_EDGE) then
        do i=IstrU,Iend
          cff=(Dnew(i,Jstr-1)+Dnew(i-1,Jstr-1))*(ubar(i,Jstr-1,knew)
#  ifdef MRL_WCI
     &                                              +ust2d(i,Jstr-1)
#  endif
     &                                              )*on_u(i,Jstr-1)
          DU_avg1(i,Jstr-1,nnew)=DU_avg1(i,Jstr-1,nnew)+cff1*cff
          DU_avg2(i,Jstr-1)=DU_avg2(i,Jstr-1)+cff2*cff
        enddo
        do i=IstrR,IendR
          cff=(Dnew(i,JstrV-1)+Dnew(i,JstrV-2))*(vbar(i,JstrV-1,knew)
#  ifdef MRL_WCI
     &                                              +vst2d(i,JstrV-1)
#  endif
     &                                              )*om_v(i,JstrV-1)
          DV_avg1(i,JstrV-1,nnew)=DV_avg1(i,JstrV-1,nnew)+cff1*cff
          DV_avg2(i,JstrV-1)=DV_avg2(i,JstrV-1)+cff2*cff
        enddo
      endif
      if (NORTHERN_EDGE) then
        do i=IstrU,Iend
          cff=(Dnew(i,Jend+1)+Dnew(i-1,Jend+1))*(ubar(i,Jend+1,knew)
#  ifdef MRL_WCI
     &                                              +ust2d(i,Jend+1)
#  endif
     &                                              )*on_u(i,Jend+1)
          DU_avg1(i,Jend+1,nnew)=DU_avg1(i,Jend+1,nnew)+cff1*cff
          DU_avg2(i,Jend+1)=DU_avg2(i,Jend+1)+cff2*cff
        enddo
        do i=IstrR,IendR
          cff=(Dnew(i,Jend+1)+Dnew(i,Jend))*(vbar(i,Jend+1,knew)
#  ifdef MRL_WCI
     &                                          +vst2d(i,Jend+1)
#  endif
     &                                          )*om_v(i,Jend+1)
          DV_avg1(i,Jend+1,nnew)=DV_avg1(i,Jend+1,nnew)+cff1*cff
          DV_avg2(i,Jend+1)=DV_avg2(i,Jend+1)+cff2*cff
        enddo
      endif
# endif /* !NS_PERIODIC */
!
!=====================================================================
!
!  Adjust ZETA using depth-integrated continuity equation 
!  and update grid
!
!  Once rhobar_nbq and depth-averaged momentum is updated, surface anomalies
!  can be adjusted to satisfy the low-frequency mass conservation 
!  equation. However, this adjustment do not satisfy conservation at machine 
!  precision and a final correction is needed.
!
!  This recomputation of zeta(m+1) using div(ubar) is only done for PRECISE
!  option. In PERF option, zeta(m+1) is computed only once based on the 
!  surface characteristic relation. But the final numerical correction 
!  is applied in all cases.
!
!=====================================================================
!
# ifdef NBQ
!-----------------------------------------------------------------------
!  Update zeta(m+1)
!-----------------------------------------------------------------------
!
        do j=Jstr-1,Jend+1
          do i=Istr-1,Iend+1
            Dnew(i,j)=(zeta(i,j,knew)+h(i,j))
#  ifdef NBQ_MASS
     &                                          *rhobar_nbq(i,j,knew)
#  endif
          enddo
        enddo

        do j=Jstr,Jend
          do i=Istr,Iend
            zeta(i,j,knew)=( (h(i,j)+zeta(i,j,kstp))
#  ifdef NBQ_MASS
     &                                          *rhobar_nbq(i,j,kstp)
#  endif
     &                                + (dtfast*pm(i,j)*pn(i,j)*0.5*(
     &      (Dnew(i  ,j)+Dnew(i-1,j))*(ubar(i  ,j,knew)
#  ifdef MRL_WCI
     &                               +ust2d(i  ,j)
#  endif  
     &                                                 )*on_u(i  ,j)
     &     -(Dnew(i+1,j)+Dnew(i  ,j))*(ubar(i+1,j,knew)
#  ifdef MRL_WCI
     &                               +ust2d(i+1,j)
#  endif 
     &                                                 )*on_u(i+1,j)
     &     +(Dnew(i,j  )+Dnew(i,j-1))*(vbar(i,j  ,knew)
#  ifdef MRL_WCI
     &                               +vst2d(i,j  )
#  endif
     &                                                 )*om_v(i,j  )
     &     -(Dnew(i,j+1)+Dnew(i,j  ))*(vbar(i,j+1,knew)
#  ifdef MRL_WCI
     &                               +vst2d(i,j+1)
#  endif
     &                                                 )*om_v(i,j+1))) )
#  ifdef NBQ_MASS
     &                                             /rhobar_nbq(i,j,knew)
#  endif
     &                                                          - h(i,j)
          enddo
        enddo
        

!
!-----------------------------------------------------------------------
!  Set masking for zeta, including wet/dry conditions
!-----------------------------------------------------------------------
!
#  ifdef MASKING
      do j=Jstr,Jend
        do i=Istr,Iend
          zeta(i,j,knew)=zeta(i,j,knew)*rmask(i,j)
#   ifdef WET_DRY
!    modify new free-surface to ensure that depth 
!    is > Dcrit in masked cells.
          cff=0.5+SIGN(0.5,Dcrit(i,j)-h(i,j))
          zeta(i,j,knew)=zeta(i,j,knew)+ 
     &                   cff*(Dcrit(i,j)-h(i,j))*(1.-rmask(i,j))
#   endif
        enddo
      enddo 
#  endif /* MASKING */
!
!-----------------------------------------------------------------------
!  Set boundary conditions for the free-surface
!  --> ensure closed boundaries
!-----------------------------------------------------------------------
!
!$acc end kernels

#  ifndef OBC_NBQ
      call zetabc_tile (Istr,Iend,Jstr,Jend)
#  endif

# endif /* NBQ */
!
!-----------------------------------------------------------------------
!  Update Zt_avg1 at last fast step
!-----------------------------------------------------------------------
!
      if (LAST_FAST_STEP) then
!$acc parallel loop if(compute_on_device) default(present)
        do j=JstrR,JendR
          do i=IstrR,IendR
            Zt_avg1(i,j)=zeta(i,j,knew)
          enddo
        enddo
!#  if defined BD_EXCHANGE
!        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
!     &                          Zt_avg1(START_2D_ARRAY))
!#  endif
      endif
!
!-----------------------------------------------------------------------
! Update grid parameters at m+1: Hz, z_r, z_w
! in prognostic or diagnostic way
!-----------------------------------------------------------------------
!
# ifdef NBQ
#  ifdef NBQ_GRID_SLOW
      if (LAST_FAST_STEP) then
#  endif

#  ifdef NBQ_HZ_PROGNOSTIC
!
!  Prognostic evaluation using momentum divergence
!
!$acc parallel loop if(compute_on_device) default(present)
        do j=JstrV-1,Jend
          do i=IstrU-1,Iend
            do k=1,N
              Hz(i,j,k)=Hz_bak2(i,j,k) - dtfast*(thetadiv_nbq(i,j,k)
     &                                         +thetadiv2_nbq(i,j,k))
              Hzr(i,j,k)=(Hz(i,j,k)-rho_nbq(i,j,k))/(1.+rho(i,j,k)/rho0)
              z_w(i,j,k)=z_w(i,j,k-1)+Hzr(i,j,k)
              z_r(i,j,k)=0.5*(z_w(i,j,k)+z_w(i,j,k-1))
            enddo
          enddo
        enddo
#   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,
     &                          Hzr(START_2D_ARRAY,1))
#   endif
#  else
!
!  Diagnostic evaluation from zeta(m+1)
!
        call set_depth_tile(Istr,Iend,Jstr,Jend)
#  endif
!
! Compute derived grid parameters if fast update
!
#  ifndef NBQ_GRID_SLOW
        call grid_nbq_tile(Istr,Iend,Jstr,Jend,
     &                     Hzw_half_nbq_inv,   Hzr_half_nbq_inv,
     &                     Hzw_half_nbq_inv_u, Hzw_half_nbq_inv_v,
     &                     Hzu_half_qdmu,      Hzv_half_qdmv)
#  endif

#  ifdef NBQ_GRID_SLOW
      endif !<-- LAST_FAST_STEP
#  endif
# endif /* NBQ */
!
!-----------------------------------------------------------------------
!  Exchange boundary information.!$acc kernels if(compute_on_device) default(present)
!-----------------------------------------------------------------------
!
# if defined BD_EXCHANGE
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                        zeta(START_2D_ARRAY,knew))
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        DU_avg1(START_2D_ARRAY,nnew))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        DV_avg1(START_2D_ARRAY,nnew))
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        DU_avg2(START_2D_ARRAY))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        DV_avg2(START_2D_ARRAY))
#  if defined MRL_WCI && defined WET_DRY
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        ust2d(START_2D_ARRAY))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        vst2d(START_2D_ARRAY))
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Apply conservation requirements for nesting
!-----------------------------------------------------------------------
!
# ifdef AGRIF
      if (.NOT.Agrif_Root()) THEN
        do j=Jstr-1,Jend+1
          do i=Istr-1,Iend+1
            Zt_avg3(i,j,iif)=zeta(i,j,knew)
          enddo
        enddo 
        do j=JstrR,JendR
          do i=Istr,IendR
            DU_avg3(i,j,iif) = 0.5*(h(i  ,j)+zeta(i  ,j,knew)+
     &                              h(i-1,j)+zeta(i-1,j,knew)) 
     &                               *on_u(i,j)*ubar(i,j,knew)
          enddo
        enddo 
        do j=Jstr,JendR
          do i=IstrR,IendR
            DV_avg3(i,j,iif) = 0.5*(h(i,j  )+zeta(i,j  ,knew)+
     &                              h(i,j-1)+zeta(i,j-1,knew)) 
     &                               *om_v(i,j)*vbar(i,j,knew)
          enddo
        enddo
 
#  if defined BD_EXCHANGE
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          Zt_avg3(START_2D_ARRAY,iif))
        call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                          DU_avg3(START_2D_ARRAY,iif))
        call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                          DV_avg3(START_2D_ARRAY,iif))
#  endif

#  ifdef CVTK_DEBUG_ADVANCED
C$OMP BARRIER
C$OMP MASTER
        call check_tab2d(Zt_avg3(:,:,iif),'Zt_avg3 step3d_fast','r')
        call check_tab2d(DU_avg3(:,:,iif),'DU_avg3 step3d_fast','u')
        call check_tab2d(DV_avg3(:,:,iif),'DV_avg3 step3d_fast','v')
C$OMP END MASTER
#  endif   
      endif

#  ifdef AGRIF_CONSERV_VOL
      if (iif.eq.nfast) then
        if (agrif_root()) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              DU_avg1(i,j,5) = dt * DU_avg2(i,j)
              DV_avg1(i,j,5) = dt * DV_avg2(i,j)
            enddo
          enddo
        else
          do j=JstrR,JendR
            do i=IstrR,IendR
              DU_avg1(i,j,5) = dt * DU_avg2(i,j)
              DV_avg1(i,j,5) = dt * DV_avg2(i,j)
              DU_avg1(i,j,4) = DU_avg1(i,j,4) + DU_avg1(i,j,5)
              DV_avg1(i,j,4) = DV_avg1(i,j,4) + DV_avg1(i,j,5)
            enddo
          enddo       
        endif
      endif
#  endif
# endif /* AGRIF */
!
!-----------------------------------------------------------------------
!  Copy density for extrapolation: 
!-----------------------------------------------------------------------
!

# ifdef NBQ_MASS
      if (LAST_FAST_STEP) then
!$acc parallel loop if(compute_on_device) default(present)
        do k=1,N
          do j=JstrV-2,Jend+1
            do i=IstrU-2,Iend+1
              rho_bak(i,j,k)=rho(i,j,k)
            enddo
           enddo
         enddo
      endif
# endif
!
!-----------------------------------------------------------------------
!  Correct Hz(m+1) for internal mode
!  by inverting internal continuity equation
!-----------------------------------------------------------------------
!
# ifdef NBQ_HZCORRECT
     
      if (LAST_FAST_STEP) then
!$acc parallel loop if(compute_on_device) default(present)
        do k=1,N
          do j=JstrV-2,Jend+1
            do i=IstrU-2,Iend+1
              Hz_correct(i,j,k)=Hz(i,j,k)
            enddo
          enddo
        enddo

#  if defined BD_EXCHANGE
        call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                          DU_avg2(START_2D_ARRAY))
        call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                          DV_avg2(START_2D_ARRAY))
#  endif

!$acc kernels if(compute_on_device) default(present)
        do j=Jstr,Jend
          do i=Istr,Iend
            dum_s=0.
            do k=1,N
              dum_s=dum_s+Hz(i,j,k)-Hz_bak(i,j,k)
            enddo
            do k=1,N
              Hz(i,j,k)=Hz(i,j,k)
     &                -(dum_s
     &             + (DU_avg2(i+1,j)-DU_avg2(i,j)
     &               +DV_avg2(i,j+1)-DV_avg2(i,j)
     &               )*pm(i,j)*pn(i,j)
     &               *dt)
     &               /(z_w(i,j,N)-z_w(i,j,0))
     &               *(z_w(i,j,k)-z_w(i,j,k-1))
#  ifdef MASKING
     &               *rmask(i,j)
#  endif

#  ifdef NBQ_HZCORR_DEBUG
              Hz_corr(i,j,k)=Hz(i,j,k)-Hz_correct(i,j,k)
#  endif

!# ifdef NBQ_MASS
!              Hzr(i,j,k)=(Hz(i,j,k)-rho_nbq(i,j,k)) ! adjust Zw,Zr?
!     &                        /(1.+rho(i,j,k)/rho0)
!# endif
!              z_w(i,j,k)=z_w(i,j,k-1)+Hzr(i,j,k)
!              z_r(i,j,k)=0.5*(z_w(i,j,k)+z_w(i,j,k-1))

            enddo
          enddo
        enddo
!$acc end kernels
#  if defined BD_EXCHANGE
        call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                          Hz(START_2D_ARRAY,1))
!# ifdef NBQ_MASS
!         call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
!     &                           Hzr(START_2D_ARRAY,1))
!# endif
!        call exchange_r3d_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))
#  endif 

      endif ! iif.eq.nfast

# endif  /* NBQ_HZCORRECT */

!
!-----------------------------------------------------------------------
!  TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP
!-----------------------------------------------------------------------
!

!$acc kernels if(compute_on_device) default(present)
      VMAXL=100.
      VMAX=0.
      do j=Jstr,Jend
        do i=Istr,Iend
          VMAX=max(VMAX,abs(ubar(i,j,knew)))
        enddo
      enddo
      do j=Jstr,Jend
        do i=Istr,Iend
          VMAX=max(VMAX,abs(vbar(i,j,knew)))
        enddo
      enddo
!$acc end kernels

# ifdef CVTK_DEBUG_ADVANCED
C$OMP BARRIER
C$OMP MASTER
      call check_tab2d(DU_avg1(:,:,nnew),'DU_avg1','u',
     &   ondevice=.TRUE.)
      call check_tab2d(DV_avg1(:,:,nnew),'DV_avg1','v',
     &   ondevice=.TRUE.)
#  ifdef NBQ_HZCORRECT
      if (LAST_FAST_STEP) then
      call check_tab3d(Hz,'Hz last fast step','r',
     &   ondevice=.TRUE.)
      endif
#  endif
C$OMP END MASTER
# endif 

      IF (VMAX.GT.VMAXL) THEN
!$acc update if(compute_on_device) host(ubar(Istr:Iend,Jstr:Jend,knew))
!$acc update if(compute_on_device) host(vbar(Istr:Iend,Jstr:Jend,knew))
        cff1=maxval(abs(ubar(Istr:Iend,Jstr:Jend,knew)))
        cff2=maxval(abs(vbar(Istr:Iend,Jstr:Jend,knew)))
        IF (cff1>cff2) THEN
          ijmax= maxloc(abs(ubar(Istr:Iend,Jstr:Jend,knew)))
        ELSE
          ijmax= maxloc(abs(vbar(Istr:Iend,Jstr:Jend,knew)))
        ENDIF
        imax=ijmax(1)
        jmax=ijmax(2)
# ifdef MPI
        imax=imax+iminmpi-1
        jmax=jmax+jminmpi-1
# endif
        write(stdout,'(9(A/))')
     &     '                                         ',
     &     '                                         ',
     &     ' ======================================= ',
     &     ' =                                     = ',
     &     ' =  STEP3D_FAST:   ABNORMAL JOB END    = ',
     &     ' =                 BLOW UP             = ',
     &     ' =                                     = ',
     &     ' ======================================= ',
     &     '                                         '

# ifdef MPI
        write(stdout,'(A,I4)') '      mynode =',mynode
# endif

        if (VMAX.eq.666.) then
          write(stdout,'(A,F10.2)')
     &                                            '  VMAX (M/S) =   NaN'
        else
          write(stdout,'(A,F10.2)') 
     &                                            '  VMAX (M/S) =',VMAX
        endif
        write(stdout,'(A,2I6)')  
     &                                       '  IMAX JMAX   =',imax,jmax
        write(stdout,'(A,2I6/)')
     &                                       '  IINT IFAST  =',iic,iif

!     GC quik fix Comment this because it blocks the correct mpi_abort
!     croco do not hang out  
!!      call wrt_his            ! get output during blow-up

        may_day_flag=1
# ifdef MPI
        call mpi_abort (MPI_COMM_WORLD, err)
# else
        stop                            !-->  EXIT
# endif
      ENDIF

# undef zwrk 
# undef rzeta  
# undef rzeta2  
# undef rzetaSA 

      return
      end
#else
      subroutine step3d_fast_empty
      end
#endif /* M3FAST */
