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

#if !defined M3FAST || (defined AGRIF && defined NBQ_CHILD_ONLY)

# define PRED_COUPLED_MODE

      subroutine step2d (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()

      call step2D_FB_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)
# ifdef M2_HADV_UP3
     &                    , A2d(1,10,trd), A2d(1,11,trd)
# else
     &                    , A2d(1,10,trd), A2d(1,11,trd)
     &                    , A2d(1,12,trd), A2d(1,13,trd)
# endif
     &                    )

      return
      end

      subroutine step2D_FB_tile (Istr,Iend,Jstr,Jend, zeta_new,
     &                           Dnew,rubar,rvbar,
     &                           Drhs, UFx,UFe,
     &                           VFx,VFe
# ifdef M2_HADV_UP3
     &                          ,wrk1,wrk2
# else
     &                          ,urhs,vrhs
     &                          ,DUon,DVom
# endif
     &                          )

!
!======================================================================
! Perform one time step for barotropic mode (free-surface and baro-
! tropic 2D momentum equations) using Generalized Forward-Backward
! AB3-AM4 algorithm. Also calculate fast-time averages to interact
! with baroclinic mode.
!
! Forward-Backward schemes for coupled equations (here free-surface zeta
! and 2D momentum ubar and vbar) are stable for longer time-steps
! (alpha_max>1) than synchronous schemes while only needing evaluation
! of each r.h.s. term ones. Rather than the original Euler FB schemes,
! a generalized FB scheme with an AB3-AM4 step is used here for its added
! robustness (for Coriolis and advection terms). It is a combination of
! AB3-like step for zeta and AM4-like step for ubar and vbar but with
! coefficients chosen for maximum stability (e.g. beta=5/12 in original
! AB3 is replaced by 0.281105):
!
! AB3 Forward step:
! ----------------
! ubarrhs = (3/2+beta)*ubar(n) - (1/2+2*beta)*ubar(n-1) + beta*ubar(n-2)
! vbarrhs = (3/2+beta)*vbar(n) - (1/2+2*beta)*vbar(n-1) + beta*vbar(n-2)
! zeta(n+1) = zeta(n) + RZETA(ubarhs,vbarrhs)
!
! AM4 Backward step:
! -----------------
! zetarhs = (1/2+gamma+2*epsilon) * zeta(n+1) +
!           (1/2-2*gamma-3*epsilon) * zeta(n) +
!           gamma * zeta(n-1) + epsilon * zeta(n-2)
! ubar(n+1) = ubar(n) + RUBAR(zetarhs,ubarrhs,vbarrhs)
! vbar(n+1) = vbar(n) + RVBAR(zetarhs,ubarrhs,vbarrhs)
!
! beta=0.281105   gamma=0.0880    epsilon=0.013
!
! Note that the first 2D step is performed with the original Euler FB
! scheme and the second 2D step with a AB2-AM3 FB scheme with coefficients
! again chosen for maximum stability.
!
! Reference:
! ---------
! 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.
!
!======================================================================
!
      USE debug
      implicit none
# include "param.h"
      integer Istr,Iend,Jstr,Jend, i,j, kbak, kold,
# ifdef MPI
     &         err,
# endif
# if defined PSOURCE || defined PSOURCE_MASS
     &        is,
# endif
     &        imin,imax,jmin,jmax
      integer,dimension(2) :: ijmax
      real sum_c
      real    VMAX,VMAXL
      real zeta_new(PRIVATE_2D_SCRATCH_ARRAY),  cff,
     &         Dnew(PRIVATE_2D_SCRATCH_ARRAY),  cff0,
     &        rubar(PRIVATE_2D_SCRATCH_ARRAY),  cff1,
     &        rvbar(PRIVATE_2D_SCRATCH_ARRAY),  cff2,
     &         Drhs(PRIVATE_2D_SCRATCH_ARRAY),  cff3,
     &          UFx(PRIVATE_2D_SCRATCH_ARRAY),
     &          UFe(PRIVATE_2D_SCRATCH_ARRAY),  DUnew,
     &          VFx(PRIVATE_2D_SCRATCH_ARRAY),  DVnew,
     &          VFe(PRIVATE_2D_SCRATCH_ARRAY)
# ifdef M2_HADV_UP3
      real     wrk1(PRIVATE_2D_SCRATCH_ARRAY),
     &         wrk2(PRIVATE_2D_SCRATCH_ARRAY)
      real     curvX,curvE,cffE,cffX, gamma
      parameter (gamma = -0.25)
# else
      real     urhs(PRIVATE_2D_SCRATCH_ARRAY),
     &         vrhs(PRIVATE_2D_SCRATCH_ARRAY),
     &         DUon(PRIVATE_2D_SCRATCH_ARRAY),
     &         DVom(PRIVATE_2D_SCRATCH_ARRAY)
# endif
# ifdef M2FILTER_NONE
      real :: myepsilon,mybeta,myalpha,mygamma
# endif
# include "grid.h"
# include "ocean2d.h"
# include "ocean3d.h"
# include "coupling.h"
# include "forces.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 POT_TIDES
#  include "tides.h"
# endif
# ifdef AGRIF
#  include "zoom.h"
        integer irhox, irhoy, irhot
# endif
# if defined INTERNAL || defined BODYTIDE
      real U0, omega
# endif
# ifdef WET_DRY
      real cff1_WD,cff2_WD
# endif
!
# 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"
!
!======================================================================
! AB3 Forward Step: Computes zeta(n+1) and zetarhs
!======================================================================
!
! Preliminary step: compute total depth (meters) of the water
! ----------------- column and vertically integrated mass fluxes
! which are needed to compute divergence in rhs_zeta and as input
! data to compute nonlinear advection terms for the barotropic
! momentum equations.
!
# ifdef M2FILTER_NONE
#  ifdef TANK
         myalpha   = 0.01
#  else
         myalpha   = 0.25
#  endif
         mybeta    = 0.281105
         myepsilon = 0.00976186 - 0.13451357*myalpha
         mygamma   = 0.08344500 - 0.51358400*myalpha
# endif
# ifdef AGRIF
        irhox = Agrif_Irhox()
        irhoy = Agrif_Irhoy()
        irhot = Agrif_Irhot()
# endif
      if (FIRST_2D_STEP) then         ! Meaning of temporal indices
        kbak=kstp                     ! ------- -- -------- -------
        kold=kstp                     ! m-2     m-1      m      m+1
        cff1=1.                       ! kold    kbak    kstp    knew
        cff2=0.
        cff3=0.
      elseif (FIRST_2D_STEP+1) then
        kbak=kstp-1
        if (kbak.lt.1) kbak=4
        kold=kbak
        cff1=1.                  ! Logically AB2-AM3 forward-backward
        cff2=0.                  ! scheme with coefficients chosen for
        cff3=0.                  ! maximum stability ... (see below)
      else
        kbak=kstp-1
        if (kbak.lt.1) kbak=4
        kold=kbak-1
        if (kold.lt.1) kold=4
# ifdef M2FILTER_NONE
        cff1= 1.5+mybeta
        cff2=-2.0*mybeta-0.5
        cff3= mybeta
# else
        cff1= 1.781105
        cff2=-1.06221
        cff3= 0.281105
# endif
      endif
!
# if defined EW_PERIODIC || defined MPI || !defined M2_HADV_UP3
      imin=IstrU-2
      imax=Iend+1
#  define IV_EXT_RANGE Istr-1,Iend+1
# else
      imin=max(IstrU-3,0)
      imax=min(Iend+2,Lm+1)
#  define IV_EXT_RANGE max(Istr-2,0),min(Iend+2,Lm+1)
# endif
# if defined NS_PERIODIC || defined MPI || !defined M2_HADV_UP3
      jmin=JstrV-2
      jmax=Jend+1
#  define JU_EXT_RANGE Jstr-1,Jend+1
# else
      jmin=max(JstrV-3,0)
      jmax=min(Jend+2,Mm+1)
#  define JU_EXT_RANGE max(Jstr-2,0),min(Jend+2,Mm+1)
# endif
!
!$acc kernels if(compute_on_device) default(present) async(1)

      do j=jmin,jmax
        do i=imin,imax
          Drhs(i,j)=cff1*zeta(i,j,kstp)+cff2*zeta(i,j,kbak)
     &                                 +cff3*zeta(i,j,kold)
     &                                             + h(i,j)
        enddo
      enddo

      do j=JU_EXT_RANGE
        do i=imin+1,imax
          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
# undef JU_EXT_RANGE

      do j=jmin+1,jmax
        do i=IV_EXT_RANGE
          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
# undef IV_EXT_RANGE

!
! 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 ( (imin+1) .le.i .and. i.le.imax .and.
     &       (jmin+1) .le.j .and. j.le.jmax) 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)
          elseif (Dsrc(is).eq.1) then
            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

# ifdef M2_HADV_UP3
#  if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                   urhs(START_2D_ARRAY))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                   vrhs(START_2D_ARRAY))
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                   Duon (START_2D_ARRAY))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                   Dvom(START_2D_ARRAY))
#  endif
# endif

# 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
!$acc end kernels
!$acc wait
       call check_tab2d(zeta(:,:,kstp),'zeta step2d #1','r',
     &   ondevice=.TRUE.)
       call check_tab2d(ubar(:,:,kstp),'ubar step2d #1','u',
     &   ondevice=.TRUE.)
       call check_tab2d(vbar(:,:,kstp),'vbar step2d #1','v',
     &   ondevice=.TRUE.)
!$acc kernels if(compute_on_device) default(present) async(1)
C$OMP END MASTER
# endif
!
!-----------------------------------------------------------------------
! Advance free-surface:   Compute zeta_new, which is at new time
!-------- ---- --------   step, and interpolate half-step backward
! for the subsequent computation of barotropic pressure-gradient
! terms. It should be noted that because Forward Euler step is used
! to update zeta during the first barotropic step, the pressure term
! must be computed via Backward step to keep it numerically stable.
! However, this would interfere with the computation of forcing terms
! "rufrc,rvfrc" because computation of pressure gradient in 3D mode
! uses exactly the initial value of "zeta", rather than value changed
! by one barotropic time step.  To resolve this conflict, the
! pressure gradient term computation during the first barotropic
! step is computed in two stages: first use just zeta(:,:,kstp) to
! insure exact consistency with 3D mode; then, after "rufrc,rvfrc"
! are finalized, add a correction term based on the difference
! zeta_new(:,:)-zeta(:,:,kstp) to "rubar,rvbar" to make them
! consistent with Backward step for pressure gradient terms.
!-----------------------------------------------------------------------
!
      if (FIRST_2D_STEP) then
# ifdef SOLVE3D
        cff0=0.                !---> Compute pressure-gradient
        cff1=1.                !  terms using just zeta(:,:,kstp)
# else
        cff0=1.
        cff1=0.
# endif
        cff2=0.
        cff3=0.
      elseif (FIRST_2D_STEP+1) then
        cff0= 1.0833333333333    ! Logically AB2-AM3 forward-backward
        cff1=-0.1666666666666    ! scheme with coefficients chosen for
        cff2= 0.0833333333333    ! maximum stability, while maintaining
        cff3= 0.                 ! third-accuracy; alpha_max=1.73
      else
# ifdef M2FILTER_NONE
        cff0=0.5+2.*myepsilon+mygamma+2.*myalpha
        cff1=1.-cff0-mygamma-myepsilon
        cff2=mygamma
        cff3=myepsilon
# else
        cff0=0.614
        cff1=0.285
        cff2=0.088
        cff3=0.013
# endif
      endif

# define zwrk UFx
# define rzeta  UFe
# define rzeta2  VFe
# define rzetaSA VFx

!
!-----------------------------------------------------------------------
! Computes zeta(n+1)
!-----------------------------------------------------------------------
!
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zeta_new(i,j)=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 ZONAL_NUDGING
          zeta(i,j,knew)=zeta_new(i,j)
# endif
        enddo
      enddo
!
!-----------------------------------------------------------------------
! Add nudging terms
!-----------------------------------------------------------------------
!
# ifdef ZNUDGING
#  ifdef ZONAL_NUDGING
!$acc end kernels
      if (iic.eq.ntstart .or. mod(iic,10).eq.0) then
        if (FIRST_2D_STEP) then
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   zeta(START_2D_ARRAY,knew),zetazon)
        endif
      endif
      if (iic.eq.ntstart) then
        if (FIRST_2D_STEP) then
          call zonavg_2d(Istr,Iend,Jstr,Jend,
     &                   ssh(START_2D_ARRAY),sshzon)
        endif
      endif
!$acc kernels if(compute_on_device) default(present) async(1)
#  endif  /* ZONAL_NUDGING */
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zeta_new(i,j)=zeta_new(i,j) + dtfast*Znudgcof(i,j)
#  ifdef ZONAL_NUDGING
     &                                 *(sshzon(j)-zetazon(j))
#  else
     &                                 *(ssh(i,j)-zeta_new(i,j))
#  endif /* ZONAL_NUDGING */
        enddo
      enddo
# endif /* ZNUDGING */
!
!-----------------------------------------------------------------------
! Mask land cells & ensure that total depth in masked cells is > Dcrit
!-----------------------------------------------------------------------
!
# ifdef MASKING
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
#  ifdef WET_DRY
          cff=0.5+SIGN(0.5,Dcrit(i,j)-h(i,j))
          zeta_new(i,j)=zeta_new(i,j)+
     &                 cff*(Dcrit(i,j)-h(i,j))*(1.-rmask(i,j))
#  endif
        enddo
      enddo
# endif
!
!-----------------------------------------------------------------------
! Computes zetarhs to use in mometum equations
!-----------------------------------------------------------------------
!
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          zwrk(i,j)=cff0*zeta_new(i,j) +cff1*zeta(i,j,kstp)
     &             +cff2*zeta(i,j,kbak)+cff3*zeta(i,j,kold)
# if defined VAR_RHO_2D && defined SOLVE3D
          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
!
!-----------------------------------------------------------------------
! Load new free-surface values into shared array
!-----------------------------------------------------------------------
!
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          Dnew(i,j)=zeta_new(i,j)+h(i,j)
          zeta(i,j,knew)=zeta_new(i,j)
        enddo
      enddo
!
!-----------------------------------------------------------------------
! Load rhs values into additional AGRIF shared array for nesting
!-----------------------------------------------------------------------
!
# ifdef AGRIF
      if (FIRST_2D_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 (index 0) step2d','r',
     &   ondevice=.TRUE.)
       call check_tab2d(DU_avg3(:,:,0),'DU_avg3 (index 0) step2d','u',
     &   ondevice=.TRUE.)
       call check_tab2d(DV_avg3(:,:,0),'DV_avg3 (index 0) step2d','v',
     &   ondevice=.TRUE.)
C$OMP END MASTER
       endif
#  endif
# endif /* AGRIF */
!
!
!-----------------------------------------------------------------------
! Compute wet/dry masks
!-----------------------------------------------------------------------
!
!$acc end kernels

# ifdef WET_DRY
      call wetdry_tile (Istr,Iend,Jstr,Jend)
!RB this barrier seems to have a benefic effect
! in some cases with wetdry but not sure it is
! the proper fix
C$OMP BARRIER
# endif

!-----------------------------------------------------------------------
! Apply mass point sources (volume vertical influx) 
!-----------------------------------------------------------------------
!
# ifdef PSOURCE_MASS
      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
          zeta(i,j,knew)=zeta(i,j,knew)+ 
     &    Qbar(is)*pm(i,j)*pn(i,j)*dtfast 
        endif
      enddo
# endif

!-----------------------------------------------------------------------
! Apply Rain flux 
!-----------------------------------------------------------------------
# ifdef RAIN_FLUX
      do j=JstrR,JendR
       do i=IstrR,IendR
        zeta(i,j,knew)=zeta(i,j,knew)-EmP(i,j)*dtfast 
       enddo
      enddo
# endif
!
!-----------------------------------------------------------------------
! Debug zeta
!-----------------------------------------------------------------------
!
# ifdef CVTK_DEBUG_ADVANCED
C$OMP BARRIER
C$OMP MASTER
!$acc wait
       call check_tab2d(zeta(:,:,knew),'zeta step2d #2','rint',
     &   ondevice=.TRUE.)
C$OMP END MASTER
# endif

!
!-----------------------------------------------------------------------
! Set boundary conditions for the free-surface
!-----------------------------------------------------------------------
!
      call zetabc_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
! Compute time averaged fields over all short timesteps.
!
! Reset/initialise arrays for averaged fields during the first
! barotropic time step; Accumulate averages after that. Include
! physical boundary points, but not periodic ghost points or
! computation  MPI computational margins.
!----------------------------------------------------------------------
!
# ifdef SOLVE3D
      cff1=weight(1,iif)
      cff2=weight(2,iif)
!$acc kernels if(compute_on_device) default(present) async(1)
      if (FIRST_2D_STEP) then
        do j=JstrR,JendR
          do i=IstrR,IendR
            Zt_avg1(i,j)=cff1*zeta(i,j,knew)
            DU_avg1(i,j,nnew)=0.
            DV_avg1(i,j,nnew)=0.
            DU_avg2(i,j)=cff2*DUon(i,j)
            DV_avg2(i,j)=cff2*DVom(i,j)
          enddo
        enddo
      else
        do j=JstrR,JendR
          do i=IstrR,IendR
            Zt_avg1(i,j)=Zt_avg1(i,j)+cff1*zeta(i,j,knew)
            DU_avg2(i,j)=DU_avg2(i,j)+cff2*DUon(i,j)
            DV_avg2(i,j)=DV_avg2(i,j)+cff2*DVom(i,j)
          enddo
        enddo
      endif
# else      
!$acc kernels if(compute_on_device) default(present) async(1)
# endif

!
!
!======================================================================
! AM4 Backward Step: Computes ubar,vbar at time n+1
!                    first computes rubar,rvbar
!======================================================================
!
!
! Compute pressure-gradient terms  NOTE that "rubar" and "rvbar"
!-------- -------- -------- -----  are computed within the same
! fused loop despite the fact that their normal index ranges are
! different. Fusing loops causes redundant computation of one
! column of "rubar" on the western physical boundary and one row
! of "rvbar" on the southern, but, at the same time it allows to
! share references to array elements (i,j) which results in an
! increase of computational density by almost a factor of 1.5
! resulting in overall more efficient code pipelined in 26 cycles
! (61% of peak speed) on R10000 vs. 16+16 cycles of separate loop
! version for the case when both CPP switches below are defined.
!
      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)
# if defined VAR_RHO_2D && defined SOLVE3D
     &              +(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
# if defined POT_TIDES && !defined SOLVE3D
     &                  + ( h(i-1,j)+h(i,j)+rzeta(i-1,j)+rzeta(i,j) )
     &                                   *( Ptide(i,j)-Ptide(i-1,j) )
# endif
# if defined READ_PATM && !defined SOLVE3D
     &                  + ( h(i-1,j)+h(i,j)+rzeta(i-1,j)+rzeta(i,j) )
     &                                   *( patm2d(i,j)-patm2d(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)
# if defined VAR_RHO_2D && defined SOLVE3D
     &              +(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
# if defined POT_TIDES && !defined SOLVE3D
     &                  + ( h(i,j-1)+h(i,j)+rzeta(i,j-1)+rzeta(i,j) )
     &                                   *( Ptide(i,j)-Ptide(i,j-1) )
# endif
# if defined READ_PATM && !defined SOLVE3D
     &                  + ( h(i,j-1)+h(i,j)+rzeta(i,j-1)+rzeta(i,j) )
     &                                   *( patm2d(i,j)-patm2d(i,j-1) )
# endif
     &                                                              )
        enddo
      enddo            !--> discard  zwrk, rzeta, rzeta2, rzetaSA

# undef rzetaSA
# undef rzeta2
# undef rzeta
# undef zwrk
!
!========================================================================
! Compute horizontal advection terms for momentum equations (2D only)
!-------- ---------- --------- ----- --- -------- --------- --- -----
! 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 the purpose 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.
!========================================================================
!
# ifdef UV_ADV
!
#  if defined SOLVE3D || !defined M2_HADV_UP3
!-----------------------------------------------------------------------
! 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)
!-----------------------------------------------------------------------
      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
#   ifdef WET_DRY0
     &                                 *pmask_wet(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
#   ifdef WET_DRY0
     &                                 *pmask_wet(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

#  else
!-----------------------------------------------------------------------
!  UP3 (QUICK) spatial advection scheme as in 3D part
!
!  in this case, there is no numerical requirement for explicit viscosity
!-----------------------------------------------------------------------
#   define uxx wrk1
#   define Huxx wrk2
!
#   ifdef EW_PERIODIC
#    define IU_EXT_RANGE IstrU-1,Iend+1
#   else
#    ifdef MPI
        if (WEST_INTER) then
          imin=IstrU-1
        else
          imin=max(IstrU-1,2)
        endif
        if (EAST_INTER) then
          imax=Iend+1
        else
          imax=min(Iend+1,Lmmpi)
        endif
#     define IU_EXT_RANGE imin,imax
#    else
#     define IU_EXT_RANGE max(IstrU-1,2),min(Iend+1,Lm)
#    endif
#   endif

        do j=Jstr,Jend
          do i=IU_EXT_RANGE
            uxx(i,j)=urhs(i-1,j)-2.*urhs(i,j)+urhs(i+1,j)
            Huxx(i,j)=Duon(i-1,j)-2.*Duon(i,j)+Duon(i+1,j)
          enddo
        enddo
#   undef IU_EXT_RANGE
#   ifndef EW_PERIODIC
        if (WESTERN_EDGE) then
          do j=Jstr,Jend
            uxx(1,j)=uxx(2,j)
            Huxx(1,j)=Huxx(2,j)
          enddo
        endif
        if (EASTERN_EDGE) then
#    ifdef MPI
          do j=Jstr,Jend
            uxx(Lmmpi+1,j)=uxx(Lmmpi,j)
            Huxx(Lmmpi+1,j)=Huxx(Lmmpi,j)
          enddo
#    else
          do j=Jstr,Jend
            uxx(Lm+1,j)=uxx(Lm,j)
            Huxx(Lm+1,j)=Huxx(Lm,j)
          enddo
#    endif
        endif
#   endif

        do j=Jstr,Jend
          do i=IstrU-1,Iend
            cffX=urhs(i,j)+urhs(i+1,j)
            if (cffX.gt.0.) then
              curvX=uxx(i,j)
            else
              curvX=uxx(i+1,j)
            endif
            UFx(i,j)=0.25*(cffX+gamma*curvX)*( Duon(i,j)+
     &               Duon(i+1,j)-0.125*(Huxx(i,j)+Huxx(i+1,j)))
          enddo
        enddo
#   undef Huxx
#   undef uxx
#   define vee wrk1
#   define Hvee wrk2

#   ifdef NS_PERIODIC
#    define JV_EXT_RANGE JstrV-1,Jend+1
#   else
#    ifdef MPI
        if (SOUTH_INTER) then
          jmin=JstrV-1
        else
          jmin=max(JstrV-1,2)
        endif
        if (NORTH_INTER) then
          jmax=Jend+1
        else
          jmax=min(Jend+1,Mmmpi)
        endif
#     define JV_EXT_RANGE jmin,jmax
#    else
#     define JV_EXT_RANGE max(JstrV-1,2),min(Jend+1,Mm)
#    endif
#   endif

        do j=JV_EXT_RANGE
          do i=Istr,Iend
            vee(i,j) =vrhs(i,j-1)-2.*vrhs(i,j)+vrhs(i,j+1)
            Hvee(i,j)=Dvom(i,j-1)-2.*Dvom(i,j)+Dvom(i,j+1)
          enddo
        enddo
#   undef JV_EXT_RANGE
#   ifndef NS_PERIODIC
        if (SOUTHERN_EDGE) then
          do i=Istr,Iend
            vee(i,1)=vee(i,2)
            Hvee(i,1)=Hvee(i,2)
          enddo
        endif
        if (NORTHERN_EDGE) then
#    ifdef MPI
          do i=Istr,Iend
            vee(i,Mmmpi+1)=vee(i,Mmmpi)
            Hvee(i,Mmmpi+1)=Hvee(i,Mmmpi)
          enddo
#    else
          do i=Istr,Iend
            vee(i,Mm+1)=vee(i,Mm)
            Hvee(i,Mm+1)=Hvee(i,Mm)
          enddo
#    endif
        endif
#   endif
        do j=JstrV-1,Jend
          do i=Istr,Iend
            cffE=vrhs(i,j)+vrhs(i,j+1)
            if (cffE.gt.0.) then
              curvE=vee(i,j)
            else
              curvE=vee(i,j+1)
            endif
            VFe(i,j)=0.25*(cffE+gamma*curvE)*( Dvom(i,j)+
     &               Dvom(i,j+1)-0.125*(Hvee(i,j)+Hvee(i,j+1)))
          enddo
        enddo
#   undef Hvee
#   undef vee
#   define uee wrk1
#   define Hvxx wrk2

#   ifdef NS_PERIODIC
#    define JU_EXT_RANGE Jstr-1,Jend+1
#   else
#    ifdef MPI
        if (SOUTH_INTER) then
          jmin=Jstr-1
        else
          jmin=max(Jstr-1,1)
        endif
        if (NORTH_INTER) then
          jmax=Jend+1
        else
          jmax=min(Jend+1,Mmmpi)
        endif
#     define JU_EXT_RANGE jmin,jmax
#    else
#     define JU_EXT_RANGE max(Jstr-1,1),min(Jend+1,Mm)
#    endif
#   endif

        do j=JU_EXT_RANGE
          do i=IstrU,Iend
            uee(i,j)=urhs(i,j-1)-2.*urhs(i,j)+urhs(i,j+1)
          enddo
        enddo
#   undef JU_EXT_RANGE
#   ifndef NS_PERIODIC
        if (SOUTHERN_EDGE) then
          do i=IstrU,Iend
            uee(i,0)=uee(i,1)
          enddo
        endif
        if (NORTHERN_EDGE) then
#    ifdef MPI
          do i=IstrU,Iend
            uee(i,Mmmpi+1)=uee(i,Mmmpi)
          enddo
#    else
          do i=IstrU,Iend
            uee(i,Mm+1)=uee(i,Mm)
          enddo
#    endif
        endif
#   endif
        do j=Jstr,Jend+1
          do i=IstrU-1,Iend
           Hvxx(i,j)=Dvom(i-1,j)-2.*Dvom(i,j)+Dvom(i+1,j)
          enddo
        enddo
        do j=Jstr,Jend+1
          do i=IstrU,Iend
            cffX=urhs(i,j)+urhs(i,j-1)
            cffE=Dvom(i,j)+Dvom(i-1,j)
            if (cffE.gt.0.) then
              curvX=uee(i,j-1)
            else
              curvX=uee(i,j)
            endif
            UFe(i,j)=0.25*(cffX+gamma*curvX)*(cffE-0.125*(
     &                             Hvxx(i,j)+Hvxx(i-1,j) ))
          enddo
        enddo
#   undef Hvxx
#   undef uee
#   define vxx wrk1
#   define Huee wrk2

#   ifdef EW_PERIODIC
#    define IV_EXT_RANGE Istr-1,Iend+1
#   else
#    ifdef MPI
        if (WEST_INTER) then
          imin=Istr-1
        else
          imin=max(Istr-1,1)
        endif
        if (EAST_INTER) then
          imax=Iend+1
        else
          imax=min(Iend+1,Lmmpi)
        endif
#     define IV_EXT_RANGE imin,imax
#    else
#     define IV_EXT_RANGE max(Istr-1,1),min(Iend+1,Lm)
#    endif
#   endif

        do j=JstrV,Jend
          do i=IV_EXT_RANGE
            vxx(i,j)=vrhs(i-1,j)-2.*vrhs(i,j)+vrhs(i+1,j)
          enddo
        enddo
#   undef IV_EXT_RANGE
#   ifndef EW_PERIODIC
        if (WESTERN_EDGE) then
          do j=JstrV,Jend
            vxx(0,j)=vxx(1,j)
          enddo
        endif
        if (EASTERN_EDGE) then
#    ifdef MPI
          do j=JstrV,Jend
            vxx(Lmmpi+1,j)=vxx(Lmmpi,j)
          enddo
#    else
          do j=JstrV,Jend
            vxx(Lm+1,j)=vxx(Lm,j)
          enddo
#    endif
        endif

#   endif
        do j=JstrV-1,Jend
          do i=Istr,Iend+1
           Huee(i,j)=Duon(i,j-1)-2.*Duon(i,j)+Duon(i,j+1)
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend+1
            cffE=vrhs(i,j)+vrhs(i-1,j)
            cffX=Duon(i,j)+Duon(i,j-1)
            if (cffX.gt.0.) then
              curvE=vxx(i-1,j)
            else
              curvE=vxx(i,j)
            endif
            VFx(i,j)=0.25*(cffE+gamma*curvE)*(cffX-0.125*(
     &                             Huee(i,j)+Huee(i,j-1) ))
          enddo
        enddo
#   undef Huee
#   undef vxx
        do j=Jstr,Jend
          do i=IstrU,Iend
            rubar(i,j)=rubar(i,j)-UFx(i,j  )+UFx(i-1,j)
     &                           -UFe(i,j+1)+UFe(i  ,j)
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            rvbar(i,j)=rvbar(i,j)-VFx(i+1,j)+VFx(i,j  )
     &                           -VFe(i  ,j)+VFe(i,j-1)
          enddo
        enddo
#  endif /* !SOLVE3D && M2_HADV_UP3 */
# 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
!
!-----------------------------------------------------------------------
! Compute horizontal viscous stress terms (2D only).
!-----------------------------------------------------------------------
!
# ifdef SOLVE3D
#  undef UV_VIS2
# endif
# ifdef UV_VIS2
      do j=Jstr-1,Jend
        do i=Istr-1,Iend
          cff=2.*Drhs(i,j)*visc2_r(i,j)
          UFx(i,j)=cff*(ubar(i+1,j,kstp)-ubar(i,j,kstp))
     &                                 *pm(i,j)*on_r(i,j)
          VFe(i,j)=cff*(vbar(i,j+1,kstp)-vbar(i,j,kstp))
     &                                 *pn(i,j)*om_r(i,j)

          cff1=0.0625*visc2_p(i+1,j+1)*( Drhs(i,j)
     &       +Drhs(i+1,j)+Drhs(i,j+1)+Drhs(i+1,j+1) )*(
     &          (pn(i+1,j+1)+pn(i,j+1)+pn(i+1,j)+pn(i,j))
     &             *(ubar(i+1,j+1,kstp)-ubar(i+1,j,kstp))
     &         +(pm(i+1,j+1)+pm(i,j+1)+pm(i+1,j)+pm(i,j))
     &             *(vbar(i+1,j+1,kstp)-vbar(i,j+1,kstp))
     &                                                  )
#  ifdef MASKING
     &                     *pmask(i+1,j+1)
#  endif
          UFe(i+1,j+1)=cff1*om_p(i+1,j+1)
          VFx(i+1,j+1)=cff1*on_p(i+1,j+1)
        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
# endif /* UV_VIS2 */
!
!-----------------------------------------------------------------------
! Linear and/or quadratic bottom stress.
!-----------------------------------------------------------------------
!
# 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
!
!-----------------------------------------------------------------------
! 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
!$acc end kernels
!
!-----------------------------------------------------------------------
! Coupling between 2D and 3D parts.
!--------- ------- -- --- -- ------
! Before the predictor step of the first barotropic time step
! arrays "rufrc" and "rvfrc" contain vertically integrals of the
! 3D right-hand-side terms for the momentum equations (including
! surface and bottom stresses, if so prescribed).
!
! During the first barotropic time step connvert them 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 added to "rubar" and "rvbar"
! during all subsequent barotropic time steps.
!-----------------------------------------------------------------------
!
# ifdef SOLVE3D
      if (FIRST_2D_STEP) then

#  ifdef PRED_COUPLED_MODE
        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)
!$acc kernels if(compute_on_device) default(present) async(1)
        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
#  else
        do j=Jstr,Jend                       ! This version is
          do i=Istr,Iend                     ! designed for coupling
            rufrc(i,j)=rufrc(i,j)-rubar(i,j) ! during 3D corrector
            rvfrc(i,j)=rvfrc(i,j)-rvbar(i,j) ! sub-step: no forward
          enddo                              ! extrapolation needs
        enddo                                ! to be performed.
#  endif
!
!-----------------------------------------------------------------------
! Since coupling requires that pressure gradient term is computed
! using zeta(:,:,kstp) instead of zeta_new(:,:) needed to achieve
! numerical stability, apply compensation to shift pressure gradient
! terms from "kstp" to "knew": in essense, convert the fist 2D step
! from Forward Euler to Forward-Backward].
!-----------------------------------------------------------------------
!
#  define zwrk UFx
#  define rzeta  UFe
#  define rzeta2  VFe
#  define rzetaSA VFx

        do j=JstrV-1,Jend
          do i=IstrU-1,Iend
            zwrk(i,j)=zeta_new(i,j)-zeta(i,j,kstp)
#  if defined VAR_RHO_2D && defined SOLVE3D
            rzeta(i,j)=(1.+rhoS(i,j))*zwrk(i,j)
            rzeta2(i,j)=rzeta(i,j)*(zeta_new(i,j)+zeta(i,j,kstp))
            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)*(zeta_new(i,j)+zeta(i,j,kstp))
#  endif
          enddo
        enddo

        cff=0.5*g
        do j=Jstr,Jend
          do i=Istr,Iend
            rubar(i,j)=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)
#  if defined VAR_RHO_2D && defined SOLVE3D
     &              +(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
     &                                                              )
!
            rvbar(i,j)=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)

#  if defined VAR_RHO_2D && defined SOLVE3D
     &              +(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
     &                                                              )
          enddo
        enddo            !--> discard  zwrk, rzeta, rzeta2, rzetaSA

#  undef rzetaSA
#  undef rzeta2
#  undef rzeta
#  undef zwrk
!$acc end kernels
      endif   !<-- FIRST_2D_STEP

# endif /* SOLVE3D */
!
!-----------------------------------------------------------------------
! Perform time step for the 2D momentum equations.
!
! Also compute fast-time averaged barotropic mass fluxes.
! Doing so on the fly yields a more computationally dense code and
! eliminates repeated multiplication by Dnew (since mass fluxes are
! actually available as volatile variables DUnew,DVnew at this moment.
! However, as a result of this arrangement, a special code is needed
! to compute fast-time averages along the physical boundaries, which is
! done below.
!-----------------------------------------------------------------------
!

# define Dstp DUon
!$acc kernels if(compute_on_device) default(present) async(1)
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          Dstp(i,j)=zeta(i,j,kstp)+h(i,j)
        enddo
      enddo

      cff=0.5*dtfast
# ifdef SOLVE3D
      cff1=0.5*weight(1,iif)
# else
      cff2=2.*dtfast
# endif
      do j=Jstr,Jend
        do i=IstrU,Iend
          DUnew=( (Dstp(i,j)+Dstp(i-1,j))*ubar(i,j,kstp)
     &        +cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
# ifdef SOLVE3D
     &                             *(rubar(i,j)+rufrc(i,j))
# else
     &                          *rubar(i,j)+cff2*sustr(i,j)
#  ifdef MRL_WCI
     &                                    +cff2*brk2dx(i,j)
#   ifdef WAVE_STREAMING
     &                                    +cff2*frc2dx(i,j)
#   endif
#  endif
# endif
     &                                                    )
# ifdef MASKING
     &                                         *umask(i,j)
# endif
# ifdef WET_DRY
          cff1_WD=ABS(ABS(umask_wet(i,j))-1.)
          cff2_WD=0.5+SIGN(0.5,DUnew)*umask_wet(i,j)
          umask_wet(i,j)=0.5*umask_wet(i,j)*cff1_WD
     &                         +cff2_WD*(1.-cff1_WD)
          DUnew=DUnew*umask_wet(i,j)
#  ifdef MRL_WCI
          ust2d(i,j)=ust2d(i,j)*umask_wet(i,j)
#  endif
# endif
          ubar(i,j,knew)=DUnew/(Dnew(i,j)+Dnew(i-1,j))
# ifdef SOLVE3D
          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
     &                                                  )
# endif
        enddo
      enddo

      do j=JstrV,Jend
        do i=Istr,Iend
          DVnew=( (Dstp(i,j)+Dstp(i,j-1))*vbar(i,j,kstp)
     &        +cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
# ifdef SOLVE3D
     &                             *(rvbar(i,j)+rvfrc(i,j))
# else
     &                          *rvbar(i,j)+cff2*svstr(i,j)
#  ifdef MRL_WCI
     &                                    +cff2*brk2de(i,j)
#   ifdef WAVE_STREAMING
     &                                    +cff2*frc2de(i,j)
#   endif
#  endif
# endif
     &                                                    )
# ifdef MASKING
     &                                         *vmask(i,j)
# endif
# ifdef WET_DRY
          cff1_WD=ABS(ABS(vmask_wet(i,j))-1.)
          cff2_WD=0.5+SIGN(0.5,DVnew)*vmask_wet(i,j)
          vmask_wet(i,j)=0.5*vmask_wet(i,j)*cff1_WD
     &                        +cff2_WD*(1.-cff1_WD)
          DVnew=DVnew*vmask_wet(i,j)
#  ifdef MRL_WCI
          vst2d(i,j)=vst2d(i,j)*vmask_wet(i,j)
#  endif
# endif
          vbar(i,j,knew)=DVnew/(Dnew(i,j)+Dnew(i,j-1))
# ifdef SOLVE3D
          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
     &                                                   )
# endif
        enddo
      enddo
!
!-----------------------------------------------------------------------
!  Set 2D Momemtum nudging
!-----------------------------------------------------------------------
!
# if defined M2NUDGING && defined M2CLIMATOLOGY

#  ifdef ZONAL_NUDGING
!$acc end kernels
      if (iic.eq.ntstart .or. mod(iic,10).eq.0) then
        if (FIRST_2D_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 (iic.eq.ntstart) then
        if (FIRST_2D_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
!$acc kernels if(compute_on_device) default(present) async(1)
#  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
#  ifdef SOLVE3D
          DU_avg1(i,j,nnew)=DU_avg1(i,j,nnew) +cff1*DUnew*
     &                         (Dnew(i,j)+Dnew(i-1,j))*on_u(i,j)
#  endif
        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
#  ifdef SOLVE3D
          DV_avg1(i,j,nnew)=DV_avg1(i,j,nnew) +cff1*DVnew*
     &                         (Dnew(i,j)+Dnew(i,j-1))*om_v(i,j)
#  endif

        enddo
      enddo
# endif
!
!-----------------------------------------------------------------------
!  Body force for the Internal Tide test case
!-----------------------------------------------------------------------
!
# if defined INTERNAL || defined BODYTIDE
      omega=2.*pi/(12.4*3600.)
      U0=0.02
      do j=Jstr,Jend
        do i=IstrU,Iend
          DUnew = dtfast*omega*U0*cos(omega*time)
#  ifdef MASKING
     &                 * umask(i,j)
#  endif
#  ifdef WET_DRY
     &                 * umask_wet(i,j)
#  endif
          ubar(i,j,knew)=ubar(i,j,knew) + DUnew
#  ifdef SOLVE3D
          DU_avg1(i,j,nnew)=DU_avg1(i,j,nnew) +cff1*DUnew*
     &                         (Dnew(i,j)+Dnew(i-1,j))*on_u(i,j)
#  endif
        enddo
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
          DVnew = dtfast*0.5*(f(i,j)+f(i,j-1))*
     &                   U0*sin(omega*time)
#  ifdef MASKING
     &                 * vmask(i,j)
#  endif
#  ifdef WET_DRY
     &                 * vmask_wet(i,j)
#  endif
          vbar(i,j,knew)=vbar(i,j,knew) + DVnew
#  ifdef SOLVE3D
          DV_avg1(i,j,nnew)=DV_avg1(i,j,nnew) +cff1*DVnew*
     &                         (Dnew(i,j)+Dnew(i,j-1))*om_v(i,j)
#  endif

        enddo
      enddo
# endif
!$acc end kernels
!
!-----------------------------------------------------------------------
! Set boundary conditions and compute integral mass flux accross
! all open boundaries, if any.
!-----------------------------------------------------------------------
!
      call u2dbc_tile (Istr,Iend,Jstr,Jend, UFx)
      call v2dbc_tile (Istr,Iend,Jstr,Jend, UFx)

# ifdef OBC_VOLCONS
      call obc_flux_tile (Istr,Iend,Jstr,Jend)
# endif
!
!-----------------------------------------------------------------------
! Compute fast-time averaged barotropic mass fluxes along physical
! boundaries.
!-----------------------------------------------------------------------
!
!$acc kernels if(compute_on_device) default(present) async(1)
# ifdef SOLVE3D
#  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)
        enddo
      endif
      if (EASTERN_EDGE) then
        do j=Jstr-1,JendR
          Dnew(Iend+1,j)=h(Iend+1,j)+zeta(Iend+1,j,knew)
        enddo
      endif
#  endif
#  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)
        enddo
      endif
      if (NORTHERN_EDGE) then
        do i=Istr-1,IendR
          Dnew(i,Jend+1)=h(i,Jend+1)+zeta(i,Jend+1,knew)
        enddo
      endif
#  endif
!$acc end kernels

      cff1=0.5*weight(1,iif)  
!$acc kernels if(compute_on_device) default(present) async(1)
#  ifndef EW_PERIODIC
      if (WESTERN_EDGE) then
        do j=JstrR,JendR
          DU_avg1(IstrU-1,j,nnew)=DU_avg1(IstrU-1,j,nnew)
     &         +cff1*(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)
        enddo
        do j=JstrV,Jend
          DV_avg1(Istr-1,j,nnew)=DV_avg1(Istr-1,j,nnew)
     &       +cff1*(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)
        enddo
      endif

      if (EASTERN_EDGE) then
        do j=JstrR,JendR
          DU_avg1(Iend+1,j,nnew)=DU_avg1(Iend+1,j,nnew)
     &            +cff1*( 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)
        enddo
        do j=JstrV,Jend
          DV_avg1(Iend+1,j,nnew)=DV_avg1(Iend+1,j,nnew)
     &        +cff1*( 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)
        enddo
      endif
#  endif
#  ifndef NS_PERIODIC
      if (SOUTHERN_EDGE) then
        do i=IstrU,Iend
          DU_avg1(i,Jstr-1,nnew)=DU_avg1(i,Jstr-1,nnew)
     &        +cff1*( 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)
        enddo
        do i=IstrR,IendR
          DV_avg1(i,JstrV-1,nnew)=DV_avg1(i,JstrV-1,nnew)
     &         +cff1*(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)
        enddo
      endif
      if (NORTHERN_EDGE) then
        do i=IstrU,Iend
          DU_avg1(i,Jend+1,nnew)=DU_avg1(i,Jend+1,nnew)
     &        +cff1*( 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)
        enddo
        do i=IstrR,IendR
          DV_avg1(i,Jend+1,nnew)=DV_avg1(i,Jend+1,nnew)
     &            +cff1*( 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)
        enddo
      endif
#  endif
# endif
!
!-----------------------------------------------------------------------
! 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
            ubar(i,j,knew)=2.*Qbar(is)/( on_u(i,j)
     &                       *(Dnew(i-1,j)+Dnew(i,j)) )
#  ifdef SOLVE3D
            DU_avg1(i,j,nnew)=Qbar(is)
#  endif
          elseif (Dsrc(is).eq.1) then
            vbar(i,j,knew)=2.*Qbar(is)/( om_v(i,j)
     &                       *(Dnew(i,j-1)+Dnew(i,j)) )
#  ifdef SOLVE3D
            DV_avg1(i,j,nnew)=Qbar(is)
#  endif
          endif
        endif
      enddo
# endif
!
!-----------------------------------------------------------------------
!  Diagnostics
!-----------------------------------------------------------------------
!
!$acc end kernels
# ifndef SOLVE3D
      call diag_tile (Istr,Iend,Jstr,Jend, UFx,UFe)
# endif
!
!-----------------------------------------------------------------------
!  Exchange boundary information.
!-----------------------------------------------------------------------
!
# if defined BD_EXCHANGE
!# ifdef BASIN      
!      call exchange_r2d_1pts_tile (Istr,Iend,Jstr,Jend,
!     &                   zeta(START_2D_ARRAY,knew))
!# else
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                   zeta(START_2D_ARRAY,knew))
!#endif      
      call exchange_u2d_1pts_tile (Istr,Iend,Jstr,Jend,
     &                   ubar(START_2D_ARRAY,knew))
      call exchange_v2d_1pts_tile (Istr,Iend,Jstr,Jend,
     &                   vbar(START_2D_ARRAY,knew))
#  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
!-----------------------------------------------------------------------
!
!$acc kernels if(compute_on_device) default(present) async(1)
# 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 step2d','r',
     &   ondevice=.TRUE.)
       call check_tab2d(DU_avg3(:,:,iif),'DU_avg3 step2d','u',
     &   ondevice=.TRUE.)
       call check_tab2d(DV_avg3(:,:,iif),'DV_avg3 step2d','v',
     &   ondevice=.TRUE.)
C$OMP END MASTER
#  endif
      endif
#  ifdef AGRIF_CONSERV_VOL
       if (iif == 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 */
!
!-----------------------------------------------------------------------
!  TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP
!-----------------------------------------------------------------------
!
      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
!
!-----------------------------------------------------------------------
!  Debugging ubar,vbar
!-----------------------------------------------------------------------
!
# ifdef CVTK_DEBUG_ADVANCED
C$OMP BARRIER
C$OMP MASTER
       call check_tab2d(ubar(:,:,knew),'ubar step2d #2','u',
     &   ondevice=.TRUE.)
       call check_tab2d(vbar(:,:,knew),'vbar step2d #2','v',
     &   ondevice=.TRUE.)
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/))')
     &     '                                         ',
     &     '                                         ',
     &     ' ======================================= ',
     &     ' =                                     = ',
     &     ' =   STEP2D:   ABNORMAL JOB END        = ',
     &     ' =                 BLOW UP             = ',
     &     ' =                                     = ',
     &     ' ======================================= ',
     &     '                                         '
        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
# if defined MPI
        write(stdout,'(A,I6)')
     &                                       '  NODE  =',mynode
        write(stdout,'(A,2I6)')
     &       '  IMAX JMAX  =',imax-iminmpi+1,jmax-jminmpi+1
# endif
# ifdef SOLVE3D
        write(stdout,'(A,2I6/)')
     &                                         '  IINT IEXT  =',iic,iif
# else
        write(stdout,'(A,I6/)')    '  IIC        =',iic
# endif
        may_day_flag=1
# ifdef MPI
        call mpi_abort (MPI_COMM_WORLD, err)
# else
        stop                            !-->  EXIT
# endif
      ENDIF

# ifdef SOLVE3D
      if (LAST_FAST_STEP) then
!$acc wait
      endif
# endif       
      return
      end
#else
      subroutine step2d_empty
      end
#endif /* M3FAST */
