!======================================================================
! 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 SOLVE3D
# ifndef VADV_ADAPT_PRED
#  undef VADV_ADAPT_IMP
# endif

      subroutine pre_step3d (tile)
!
      implicit none
      integer tile,  trd,omp_get_thread_num
# include "param.h"
# include "private_scratch.h"
# include "ocean3d.h"
# include "compute_tile_bounds.h"
      trd=omp_get_thread_num()
      call pre_step3d_tile (Istr,Iend,Jstr,Jend,
     &                A3d(1,1,trd), A3d(1,2,trd), A3d(1,3,trd),
     &                A2d(1,1,trd), A2d(1,2,trd), A2d(1,3,trd),
# ifdef VADV_ADAPT_IMP
     &                                            A2d(1,4,trd),
# endif
     &                A2d(1,1,trd), A2d(1,2,trd), A2d(1,3,trd)
# ifdef OPENACC
     &   ,A3d(1,4,trd),A3d(1,5,trd),A3d(1,6,trd),A3d(1,7,trd)
# endif
     &                                          , A3d(1,8,trd)         
     &                                                        )
      return
      end

      subroutine pre_step3d_tile (Istr,Iend,Jstr,Jend, ru,rv,rw,
     &                                                 FC,CF,DC,
# ifdef VADV_ADAPT_IMP
     &                                                       WC,
# endif
     &                                               FX,FE,WORK
# ifdef OPENACC
     &                           ,FX_3D,FE_3D,WORK_3D,WORK2_3D
# endif
     &                                                 ,Hz_half
     &                                                         )
!
!--------------------------------------------------------------------
! Preliminary step: initialize computations of the new time step
! 3D primitive variables.
!
! Since r.h.s. arrays ru,rv,rt(:,:,???[,:]), which at this moment
! contain r.h.s at time step n-2 will be overwritten by the
! subsequent routines within rhs3d driver, both [n-1 and n-2] old-
! time-step r.h.s. term in Adams-Bashforth stepping scheme are
! added at this time to the time step [n] fields and the result is
! stored as u,v,t(:,:,???[,:]).
!
! The actual time step will be completed in step3d, after the
! time step [n] r.h.s. terms and new-time step Hz will be available
! after the completion rhs3d computations and the 2D (barotropic
! mode) computations.
!--------------------------------------------------------------------
!
#ifdef CVTK_DEBUG
      use debug
#endif              
      implicit none
# include "param.h"
      integer Istr,Iend,Jstr,Jend, itrc, i,j,k, indx
     &       ,imin,imax,jmin,jmax,nadv,iAkt
# if defined PSOURCE || defined PSOURCE_MASS
     &       ,is
# endif
      real   ru(PRIVATE_2D_SCRATCH_ARRAY,N),    cff,
     &       rv(PRIVATE_2D_SCRATCH_ARRAY,N),    cff1,
     &       rw(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &       FC(PRIVATE_1D_SCRATCH_ARRAY,0:N),  cff2,
     &       CF(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &       DC(PRIVATE_1D_SCRATCH_ARRAY,0:N),  gamma,
     &       FX(PRIVATE_2D_SCRATCH_ARRAY),      epsil,
     &       FE(PRIVATE_2D_SCRATCH_ARRAY),        cdt,
     &     WORK(PRIVATE_2D_SCRATCH_ARRAY)
# ifdef OPENACC
      real FX_3D(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &     FE_3D(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &     WORK_3D(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &     WORK2_3D(PRIVATE_2D_SCRATCH_ARRAY,0:N)
# endif
# ifdef VADV_ADAPT_IMP      
      real   WC(PRIVATE_1D_SCRATCH_ARRAY,0:N) 
# endif 
      real Hz_half(PRIVATE_2D_SCRATCH_ARRAY,N)
# ifdef WET_DRY
      real cff3,cff4,cff5
# endif
# ifdef M3FAST
      real qdm_nstp,qdm_indx
# endif
      parameter (gamma=1./6., epsil=1.E-16)
# include "grid.h"
# include "ocean3d.h"
# include "coupling.h"
# include "ocean2d.h"
# include "forces.h"
# include "mixing.h"
# include "scalars.h"
# include "weno_interfaces.h"
# if defined PSOURCE || defined PSOURCE_MASS
#  include "sources.h"
# endif
# ifdef M3FAST
#  include "nbq.h"
# endif
# include "openacc_work.h"
# ifdef MASKING
#  define SWITCH *
# else
#  define SWITCH !
# endif

# ifdef M3FAST
#  define ru_nbq_avg1 ru_nbq
#  define rv_nbq_avg1 rv_nbq
#  ifdef NBQ
#   define rw_nbq_avg1 rw_nbq
#  endif
# endif
!
!--------------------------------------------------------------------
! Definition of flux operators: 1st, 2nd, 3rd, 4th, 5th or 6th order,
! used in UP5 and C6 advection schemes (and order degradation near
! land masks). cdiff is part of laplacian diffusion in flux1 (used
! near mask):
!    0 --> flux1=flux2 (second order C2 advection scheme)
!    1 --> flux1 gives 1st order monotonic UP1 advection scheme
!--------------------------------------------------------------------
!
# if defined TS_HADV_UP5 || defined TS_HADV_C6 \
      || defined TS_HADV_WENO5 || defined BIO_HADV_WENO5 \
      || defined TS_VADV_WENO5

      REAL    :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2
      REAL    :: ua, vel, cdiff, cdif
      REAL    :: flux1, flux2, flux3, flux4, flux5, flux6
      REAL    :: flx2, flx3, flx4, flx5
      REAL    :: mask0, mask1, mask2, mask3

      flux2(q_im1, q_i, ua, cdiff) = 0.5*( q_i + q_im1 )

      flux1(q_im1, q_i, ua, cdiff) = flux2(q_im1, q_i, ua, cdiff)
     &                          - 0.5*cdiff*sign(1.,ua)*(q_i - q_im1)

      flux4(q_im2, q_im1, q_i, q_ip1, ua) =
     &                     ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.

      flux3(q_im2, q_im1, q_i, q_ip1, ua) =
     &                            flux4(q_im2, q_im1, q_i, q_ip1, ua)
     &           + sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i - q_im1))/12.

      flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =
     &                       ( 37.*(q_i + q_im1) - 8.*(q_ip1 + q_im2)
     &                                        + (q_ip2 + q_im3) )/60.

      flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =
     &              flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)
     &             - sign(1.,ua)*( (q_ip2 - q_im3)-5.*(q_ip1 - q_im2)
     &                                      + 10.*(q_i - q_im1) )/60.

# endif
!
!--------------------------------------------------------------------
!
# include "compute_auxiliary_bounds.h"
!
      indx=3-nstp
      nadv=  nstp   !<-- do not remove: used in compute_tracer_fluxes.h

      if (FIRST_TIME_STEP) then
        cff=0.5*dt
        cff1=1.
        cff2=0.
      else
        cff=(1.-gamma)*dt
        cff1=0.5+gamma
        cff2=0.5-gamma
      endif

# if defined CVTK_DEBUG_PERFRST
       call check_tab3d(Hz(:,:,:),'Hz in prestep3d','r')
       call check_tab3d(Hz_bak(:,:,:),'Hz_bak in prestep3d','r')
       call check_tab3d(Huon(:,:,:),'Huon in prestep3d','r')
# endif
!$acc kernels if(compute_on_device) default(present)

      do k=1,N
        do j=JstrV-1,Jend
          do i=IstrU-1,Iend
            Hz_half(i,j,k)=cff1*Hz(i,j,k)+cff2*Hz_bak(i,j,k)
     &        -cff*pm(i,j)*pn(i,j)*( Huon(i+1,j,k)-Huon(i,j,k)
     &                              +Hvom(i,j+1,k)-Hvom(i,j,k)
     &                                  +We(i,j,k)-We(i,j,k-1)
# ifdef VADV_ADAPT_IMP
     &                                  +Wi(i,j,k)-Wi(i,j,k-1)
# endif
     &                                                       )
          enddo
        enddo
      enddo
!$acc end kernels
!
!  Set extended range
!
# ifdef EW_PERIODIC
      imin=Istr-2
      imax=Iend+2
# else
      if (WESTERN_EDGE) then
        imin=Istr-1
      else
        imin=Istr-2
      endif
      if (EASTERN_EDGE) then
        imax=Iend+1
      else
        imax=Iend+2
      endif
# endif
# ifdef NS_PERIODIC
      jmin=Jstr-2
      jmax=Jend+2
# else
      if (SOUTHERN_EDGE) then
        jmin=Jstr-1
      else
        jmin=Jstr-2
      endif
      if (NORTHERN_EDGE) then
        jmax=Jend+1
      else
        jmax=Jend+2
      endif
# endif
!
!--------------------------------------------------------------------
!  Start computation of the auxiliary tracer field
!  -----------------------------------------------
!  Once it will be completed, t(:,:,???,:) is effectively halfway
!  between time steps n and n+1. A high spatial order, centered,
!  non-conservative [but constancy preserving!] scheme is used for
!  this auxiliary step. This is done by introducing an artificial
!  continuity equation [Easter, 1993].
!
!  Since this field will be used exclussively to compute the high-
!  order fluxes during subsequent step3d_t operation, the final
!  values of t(i,j,k,??,itrc) alfer step3d_t will be computed in
!  a flux-conservative manner. The overall time step will be both
!  conservative and constancy preserving.
!
!  This preliminary step shall be done before field t(:,:,:,???,:)
!  loses its meaningful values during the pre-step operation.
!
!======================================================================
!
!  Compute horizontal advection
!
!======================================================================
!

# ifdef TRACERS
      do itrc=1,NT
#  if !defined OPENACC
      do k=1,N
#  endif
#  ifdef BIO_HADV_WENO5
#   if defined TEMPERATURE && defined SALINITY

          if (itrc .le. 2) then
#   elif defined TEMPERATURE || defined SALINITY
          if (itrc .le. 1) then
#   endif
#  endif
!
!----------------------------------------------------------
! Advection for active tracers: Compute fluxes FX and FE
!----------------------------------------------------------
!
#  define PREDICTOR
#  if !defined OPENACC
#   include "compute_horiz_tracer_fluxes.h"
#  else
#   include "compute_horiz_tracer_fluxes_3D.h"
#  endif
#  undef PREDICTOR

!
!----------------------------------------------------------
! WENO5 advection for passive tracers (biology ...)
!----------------------------------------------------------
!
#  ifdef BIO_HADV_WENO5
          else  ! if (itrc .gt. 2)
#   define FLUX5 flux5_weno
#   define FLUX3 flux3_weno
#   define FLUX2 flux1
#   define UP5_MASKING
!
            cdif=1.
#   include "t3dadv_order5.h"
!
#   undef FLUX5
#   undef FLUX3
#   undef FLUX2
#   undef UP5_MASKING
          endif
#  endif
!
!----------------------------------------------------------
! Apply point sources for river runoff simulations
!----------------------------------------------------------
!
#  ifdef PSOURCE
!$acc parallel loop if(compute_on_device) default(present)
!$acc& collapse(2)
          do is=1,Nsrc
#   if defined OPENACC
            DOEXTEND(k,1,N,FX,FE,WORK)
#   endif
#   ifdef MPI
           i=Isrc_mpi(is,mynode)
           j=Jsrc_mpi(is,mynode)
#   else
           i=Isrc(is)
           j=Jsrc(is)
#   endif
            if (Istr.le.i .and. i.le.Iend+1
     &                   .and. Jstr.le.j .and. j.le.Jend+1) then
              if (Dsrc(is).eq.0) then
                if (Lsrc(is,itrc)) then
                  FX(i,j)=Huon(i,j,k)*Tsrc(is,k,itrc)
#   ifdef MASKING
                else
                  if (rmask(i,j).eq.0 .and. rmask(i-1,j).eq.1) then
                    FX(i,j)=Huon(i,j,k)*t(i-1,j,k,nstp,itrc)
                  elseif (rmask(i,j).eq.1. .and. rmask(i-1,j).eq.0) then
                    FX(i,j)=Huon(i,j,k)*t(i  ,j,k,nstp,itrc)
                  endif
#   endif
                endif
              elseif(Dsrc(is).eq.1) then
                if (Lsrc(is,itrc)) then
                  FE(i,j)=Hvom(i,j,k)*Tsrc(is,k,itrc)
#   ifdef MASKING
                else
                  if (rmask(i,j).eq.0 .and. rmask(i,j-1).eq.1) then
                    FE(i,j)=Hvom(i,j,k)*t(i,j-1,k,nstp,itrc)
                  elseif (rmask(i,j).eq.1 .and. rmask(i,j-1).eq.0) then
                    FE(i,j)=Hvom(i,j,k)*t(i,j  ,k,nstp,itrc)
                  endif
#   endif
                endif
              endif
            endif
#   if defined OPENACC
            ENDDOEXTEND
#   endif
          enddo
#  endif /* PSOURCE */
!
!----------------------------------------------------------
! Finalize horizontal advection: compute flux divergences
!----------------------------------------------------------
!
          if (FIRST_TIME_STEP) then
            cff=0.5*dt
!$acc parallel loop if(compute_on_device) default(present)
#  if defined OPENACC
            DOEXTEND(k,1,N,FX,FE,WORK)
#  endif
            do j=Jstr,Jend
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=Hz(i,j,k)*t(i,j,k,nstp,itrc)
     &                 -cff*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j)
     &                                      +FE(i,j+1)-FE(i,j))
              enddo
            enddo
#  if defined OPENACC
            ENDDOEXTEND
#  endif
          else
            cff=(1.-gamma)*dt
            cff1=0.5+gamma
            cff2=0.5-gamma
!$acc parallel loop if(compute_on_device) default(present)
!$acc& collapse(3)
#  if defined OPENACC
            DOEXTEND(k,1,N,FX,FE,WORK)
#  endif
            do j=Jstr,Jend
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=cff1*Hz(i,j,k)*t(i,j,k,nstp,itrc)
     &                            +cff2*Hz_bak(i,j,k)*t(i,j,k,indx,itrc)
     &                -cff*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j)
     &                                     +FE(i,j+1)-FE(i,j))
              enddo
            enddo
#  if defined OPENACC
            ENDDOEXTEND
#  endif
          endif
#  if !defined OPENACC
        enddo   ! <-- k
#  endif
      enddo       ! <-- itrc
!
!----------------------------------------------------------
! Add tracer divergence due to cell-centered point sources
!----------------------------------------------------------
!
#  if defined PSOURCE_MASS
        do itrc=1,NT
          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 (Istr.le.i .and. i.le.Iend+1
     &        .and. Jstr.le.j .and. j.le.Jend+1) then
             cff1=dt*pm(i,j)*pn(i,j)
             do k=1,N
               if (Lsrc(is,itrc)) then
                 cff2=Tsrc(is,k,itrc)
               else
                 cff2=t(i,j,k,nstp,itrc)
               endif
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+
     &          cff1*Qsrc(is,k)*cff2
             enddo
           endif
          enddo ! source
        enddo ! itrc
#  endif /* PSOURCE_MASS */
# endif /* TRACERS */

!$acc kernels if(compute_on_device) default(present)
!
! Continue computation of the auxiliary tracer field: auxiliary
! continuity equation (the same for all tracers): DC=1/Hz_half_new,
! where Hz_half_new is Hz at time step n+1/2 as it would be computed
! from three-dimensional divergence of volume fluxes Huon,Hvom and W.
!
!======================================================================
!
! Compute vertical advection
!
!======================================================================
!
! Finalize computation of the auxiliary tracer field: compute its
! change due to vertical advection. Computation of vertical advective
! fluxes requires interpolation of tracer values to the verical grid-
! box interfaces (W-points). This can be done by either using
! parabolic spline interpolation or, more simple local cubic
! polynomial [linear interpolation is considered obsolete].
!
      if (FIRST_TIME_STEP) then
            cdt=0.5*dt
      else
            cdt=(1.-gamma)*dt
      endif

# ifdef M3FAST
      do j=Jstr,Jend
        do i=IstrU,Iend
           ru_int_nbq_2d(i,j)=0.
        enddo
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
           rv_int_nbq_2d(i,j)=0.
        enddo
      enddo
# endif

# if defined TS_HADV_TEST && defined SOLID_BODY_PER
      do k=1,N
        do j=JstR,JendR
          do i=IstrR,IendR
             u(i,j,k,nstp) = u(i,j,k,nnew)
          enddo
        enddo
        do j=JstrR,JendR
          do i=IstR,IendR
             v(i,j,k,nstp) = v(i,j,k,nnew)
          enddo
        enddo
      enddo

#  if defined BD_EXCHANGE
#   ifdef THREE_GHOST_POINTS_UV
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             u(START_2D_ARRAY,1,nstp))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             v(START_2D_ARRAY,1,nstp))
#   else
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        u(START_2D_ARRAY,1,nstp))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        v(START_2D_ARRAY,1,nstp))
#   endif
#  endif
# endif

# ifdef TRACERS
        DOLOOP2D(Istr,Iend,Jstr,Jend)
#  ifndef VADV_ADAPT_IMP 
        do k=1,N
          do i=Istr,Iend
             DC(i,k)=1./Hz_half(i,j,k)
          enddo
        enddo
#  endif
        do i=Istr,Iend
           DC(i,0)=cdt*pn(i,j)*pm(i,j)
        enddo

        do itrc=1,NT
!
!----------------------------------------------------------------------
! Compute vertical fluxes FC
!----------------------------------------------------------------------
!
#  ifdef TS_VADV_WENO5
#   undef PREDICTOR
#  else
#   define PREDICTOR
#  endif
#  include "compute_vert_tracer_fluxes.h"
#  undef PREDICTOR
!
#  ifdef TS_HADV_TEST
      FC(istr:iend,0:N) = 0.d0
#  endif
!
!----------------------------------------------------------------------
!++ Apply vertical advection for tracers (DC corresponds to 1/Hz_half)
!----------------------------------------------------------------------
!
          do k=1,N
            do i=Istr,Iend
#  ifdef VADV_ADAPT_IMP
              t(i,j,k,nnew,itrc)= t(i,j,k,nnew,itrc)  !<-- in this case division by Hz
     &               -DC(i,0)*(FC(i,k)-FC(i,k-1))     !<-- in the tridiagonal matrix problem
#  else
              t(i,j,k,nnew,itrc)=DC(i,k)*( t(i,j,k,nnew,itrc)
     &               -DC(i,0)*(FC(i,k)-FC(i,k-1)))
#   ifdef CONST_TRACERS
              t(i,j,k,nnew,itrc)=t(i,j,k,nstp,itrc)
#   endif
#  endif
            enddo
          enddo            !--> discard FC
!
!----------------------------------------------------------------------
!  Compute implicit vertical advection/diffusion
!----------------------------------------------------------------------
!
#  ifdef VADV_ADAPT_IMP
#   ifdef SALINITY
        iAKt=min(itrc,isalt)
#   else
        iAKt=min(itrc,itemp)
#   endif
#   define TRIDIAG_TRA
#   undef  TRIDIAG_U
#   undef  TRIDIAG_V
#   include "tridiag_pred.h"
#   undef TRIDIAG_TRA
#  endif
!
        enddo   !<-- itrc  !--> discard DC
        ENDDOLOOP2D
# endif /* TRACERS */
!
!======================================================================
!
! Momentum equations: time stepping to time n+1/2
!
!======================================================================
!
# undef DC
# undef FC
# undef CF

      do j=Jstr,Jend
        do i=IstrU,Iend
          WORK(i,j)=pm_u(i,j)*pn_u(i,j)
        enddo
      enddo

        if (FIRST_TIME_STEP) then
          do k=1,N
           do j=Jstr,Jend
            do i=IstrU,Iend
# ifdef VADV_ADAPT_IMP
              cff = 1.
# else
              cff = 2./(Hz_half(i,j,k)+Hz_half(i-1,j,k))
# endif
              u(i,j,k,nnew)=(u(i,j,k,nstp)*0.5*(Hz(i,j,k)+Hz(i-1,j,k))
     &                                    +cdt*WORK(i,j)*ru(i,j,k)
     &                      )*cff
# ifdef M3FAST
              ru_int_nbq(i,j,k) = 0.
!              ru_int_nbq_2d(i,j)=ru_int_nbq_2d(i,j)+ru_int_nbq(i,j,k)
# endif
              u(i,j,k,indx)=u(i,j,k,nstp)*0.5*(Hz(i,j,k)+
     &                                         Hz(i-1,j,k))
            enddo
          enddo
         enddo

!# ifdef M3FAST
!          do k=1,N
!            do i=IstrU,Iend
!               ru_int_nbq(i,j,k)=ru_int_nbq(i,j,k)
!     &            - ru_int_nbq_2d(i,j) *(Hz(i,j,k)+Hz(i-1,j,k))
!     &                               /( Zt_avg1(i  ,j)+h(i  ,j)
!     &                                 +Zt_avg1(i-1,j)+h(i-1,j))
!            enddo
!          enddo
!# endif

        else

          cff1=0.5+gamma
          cff2=0.5-gamma
          do k=1,N
            do j=Jstr,Jend
             do i=IstrU,Iend

# ifdef VADV_ADAPT_IMP
              cff = 1.
# else
              cff = 2./(Hz_half(i,j,k)+Hz_half(i-1,j,k))
# endif

# ifdef M3FAST
              qdm_nstp = u(i,j,k,nstp)*0.5*(Hz(i  ,j,k)+
     &                                      Hz(i-1,j,k))
              qdm_indx = u(i,j,k,indx)*0.5*(Hz_bak(i  ,j,k)+
     &                                      Hz_bak(i-1,j,k))
              u(i,j,k,nnew)=( cff1*qdm_nstp
     &                       +cff2*qdm_indx
     &                       +cdt*WORK(i,j)*(ru(i,j,k)+ru_nbq_avg1(i,j,k))
     &                                                     )*cff
              ru_int_nbq(i,j,k) = (qdm_nstp-qdm_indx) /dt
     &                                -WORK(i,j)*ru_nbq_avg2(i,j,k)
       
              ru_int_nbq_2d(i,j)=ru_int_nbq_2d(i,j)+ru_int_nbq(i,j,k)
# else
              u(i,j,k,nnew)=( cff1*u(i,j,k,nstp)*0.5*(Hz(i  ,j,k)+
     &                                                Hz(i-1,j,k))
     &                       +cff2*u(i,j,k,indx)*0.5*(Hz_bak(i  ,j,k)+
     &                                                Hz_bak(i-1,j,k))
     &                       +cdt*WORK(i,j)*ru(i,j,k)
     &                                                     )*cff
# endif
              u(i,j,k,indx)=u(i,j,k,nstp)*0.5*(Hz(i,j,k)+
     &                                         Hz(i-1,j,k))

            enddo
          enddo  
        enddo

# ifdef M3FAST
          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_int_nbq_2d(i,j) *(Hz(i,j,k)+Hz(i-1,j,k))
     &                               /( Zt_avg1(i  ,j)+h(i  ,j)
     &                                 +Zt_avg1(i-1,j)+h(i-1,j))
            enddo
          enddo
        enddo
# endif
         endif

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

        do j=JstrV,Jend
          do i=Istr,Iend
            WORK(i,j)=pm_v(i,j)*pn_v(i,j)
          enddo
        enddo
          if (FIRST_TIME_STEP) then
            do k=1,N
             do j=JstrV,Jend
              do i=Istr,Iend
# ifdef VADV_ADAPT_IMP
                cff = 1.
# else
                cff = 2./(Hz_half(i,j,k)+Hz_half(i,j-1,k))
# endif
                v(i,j,k,nnew)=(v(i,j,k,nstp)*0.5*(Hz(i,j,k)+Hz(i,j-1,k))
     &                                    +cdt*WORK(i,j)*rv(i,j,k)
     &                           )*cff

# ifdef M3FAST
                rv_int_nbq(i,j,k) = 0.
!                rv_int_nbq_2d(i,j)=rv_int_nbq_2d(i,j)+rv_int_nbq(i,j,k)
# endif
                v(i,j,k,indx)=v(i,j,k,nstp)*0.5*(Hz(i,j  ,k)+
     &                                           Hz(i,j-1,k))
              enddo
            enddo
          enddo

!# ifdef M3FAST
!          do k=1,N
!            do i=Istr,Iend
!               rv_int_nbq(i,j,k)=rv_int_nbq(i,j,k)
!     &            - rv_int_nbq_2d(i,j) *(Hz(i,j,k)+Hz(i,j-1,k))
!     &                               /( Zt_avg1(i,j  )+h(i,j  )
!     &                                 +Zt_avg1(i,j-1)+h(i,j-1))
!            enddo
!          enddo
!# endif

          else

            cff1=0.5+gamma
            cff2=0.5-gamma
            do k=1,N
              do j=JstrV,Jend
               do i=Istr,Iend              
# ifdef VADV_ADAPT_IMP              
                cff = 1.
# else
                cff = 2./(Hz_half(i,j,k)+Hz_half(i,j-1,k))
# endif

# ifdef M3FAST
                qdm_nstp = v(i,j,k,nstp)*0.5*(Hz(i,j  ,k)+
     &                                        Hz(i,j-1,k))
                qdm_indx = v(i,j,k,indx)*0.5*(Hz_bak(i  ,j,k)+
     &                                        Hz_bak(i,j-1,k))
                v(i,j,k,nnew)=( cff1*qdm_nstp
     &                         +cff2*qdm_indx
     &                     +cdt*WORK(i,j)*(rv(i,j,k)+rv_nbq_avg1(i,j,k))
     &                                                           )*cff

                rv_int_nbq(i,j,k) = (qdm_nstp-qdm_indx)/dt
     &                                 -WORK(i,j)*rv_nbq_avg2(i,j,k)

                rv_int_nbq_2d(i,j)=rv_int_nbq_2d(i,j)+rv_int_nbq(i,j,k)
# else
                v(i,j,k,nnew)=( cff1*v(i,j,k,nstp)*0.5*(Hz(i,j  ,k)+
     &                                                  Hz(i,j-1,k))
     &                         +cff2*v(i,j,k,indx)*0.5*(Hz_bak(i,j  ,k)+
     &                                                  Hz_bak(i,j-1,k))
     &                         +cdt*WORK(i,j)*rv(i,j,k)
     &                                                           )*cff
# endif

                v(i,j,k,indx)=v(i,j,k,nstp)*0.5*(Hz(i,j,k)+
     &                                           Hz(i,j-1,k))
              enddo
            enddo               !--> discard DC(:,:)
          enddo

# ifdef M3FAST
            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_int_nbq_2d(i,j) *(Hz(i,j,k)+Hz(i,j-1,k))
     &                               /( Zt_avg1(i,j  )+h(i,j  )
     &                                 +Zt_avg1(i,j-1)+h(i,j-1))
            enddo
          enddo
        enddo
# endif

          endif
!
!======================================================================
!

# ifdef VADV_ADAPT_IMP
!$acc loop private(DC)
      do j=Jstr,Jend
        do i=IstrU,Iend
          DC(i,0)=pm_u(i,j)*pn_u(i,j)
        enddo
!
!======================================================================
!
!  Compute implicit vertical advection/viscosity
!

#  undef  TRIDIAG_TRA
#  define TRIDIAG_U
#  undef  TRIDIAG_V
#  include "tridiag_pred.h"

!
!======================================================================
!
        if (j.ge.JstrV) then
          do i=Istr,Iend
            DC(i,0)=pm_v(i,j)*pn_v(i,j)
          enddo
!
!======================================================================
!
!  Compute implicit vertical advection/viscosity
!

#  undef  TRIDIAG_TRA
#  undef  TRIDIAG_U
#  define TRIDIAG_V

#  include "tridiag_pred.h"
!
!======================================================================
!
      enddo     !<-- j
# endif

# ifdef NBQ
!----------------------------------------------------------
!<-- at this point
!<-- rw            contains internal 3D advection + Coriolis
!<-- rw_nbq_avg1   contains NBQ pressure gradient (+ gravity)
!                                        + second viscosity
!----------------------------------------------------------
      do j=Jstr,Jend
!**************************************
      if (FIRST_TIME_STEP) then
!**************************************
        do k=1,N-1
          do i=Istr,Iend
            wz(i,j,k,nnew)=( wz(i,j,k,nstp)*(Hz(i,j,k)+Hz(i,j,k+1))
     &                                          +dt*pn(i,j)*pm(i,j)
     &                              *(rw(i,j,k)+rw_nbq_avg1(i,j,k))
     &                     )   /  (Hz_half(i,j,k)+Hz_half(i,j,k+1))

            rw_int_nbq(i,j,k)= 0.
            wz(i,j,k,indx)=wz(i,j,k,nstp)*0.5*(Hz(i,j,k  )+
     &                                         Hz(i,j,k+1))
          enddo
        enddo
!==
        do i=Istr,Iend
!== Special treatment for k=N because the control volume at the top is Hz(N)/2
            wz(i,j,N,nnew)=( wz(i,j,N,nstp)*Hz(i,j,N)
     &                            +dt*pn(i,j)*pm(i,j)
     &                *(rw(i,j,N)+rw_nbq_avg1(i,j,N))
     &                     )        /  Hz_half(i,j,N) !<-- wz has units m.s-1
            rw_int_nbq(i,j,N)= 0.
            wz(i,j,N,indx)=wz(i,j,N,nstp)*0.5*Hz(i,j,N) !<-- wz(indx) has units kg.m-1.s-1

!== Special treatment for k=0 because the control volume at the top is Hz(1)/2
#  ifdef NBQ_FREESLIP
!           wz(i,j,0,nnew)=( wz(i,j,0,nstp)*Hz(i,j,1)
!     &                            +dt*pn(i,j)*pm(i,j)
!     &                *(rw(i,j,0)+ rw_nbq_avg1(i,j,0))
!     &                     )        /  Hz_half(i,j,1) !<-- wz has units m.s-1
!           rw_int_nbq(i,j,0)=pn(i,j)*pm(i,j)*rw(i,j,0) ! ! Francis: switch from k=0 to k=1 here

            rw_int_nbq(i,j,0)=0.

            wz(i,j,0,nnew)=  ( u(i,j,1,nnew) * pm_u(i,j)
     &                                    * ( z_w(i  ,j,0)-z_w(i-1,j,0) )
     &                        +v(i,j,1,nnew) * pm_v(i,j)
     &                                    * ( z_w(i,j  ,0)-z_w(i,j-1,0) )  )
            wz(i,j,0,indx)=wz(i,j,0,nnew)*0.5*Hz(i,j,1) !<-- wz(indx) has units kg.m-1.s-1
#  else
            wz(i,j,0,nnew)=0.
            wz(i,j,0,indx)=0.
#  endif
        enddo

!**************************************
      else
!**************************************
        cff1=0.5+gamma
        cff2=0.5-gamma
!==
        do k=1,N-1
          do i=Istr,Iend
            qdm_nstp = wz(i,j,k,nstp)*(Hz(i,j,k  )+
     &                                             Hz(i,j,k+1))
            qdm_indx = wz(i,j,k,indx)*(Hz_bak(i,j,k  )+
     &                                             Hz_bak(i,j,k+1))
            wz(i,j,k,nnew)=( cff1*qdm_nstp         !
     &                      +cff2*qdm_indx         !
     &                      +2.*cdt*pn(i,j)*pm(i,j)                 ! n - 1/2
     &                             *(rw(i,j,k)+rw_nbq_avg1(i,j,k))  !
     &                     ) / (Hz_half(i,j,k)+Hz_half(i,j,k+1))  !<-- wz has units m.s-1

            rw_int_nbq(i,j,k)=(qdm_nstp-qdm_indx)*0.5/dt
     &                       - pn(i,j)*pm(i,j)*rw_nbq_avg2(i,j,k)
            wz(i,j,k,indx)=wz(i,j,k,nstp)*0.5*(Hz(i,j,k  )+
     &                                          Hz(i,j,k+1)) !<-- wz(indx) has units kg.m-1.s-1
          enddo
        enddo

        do i=Istr,Iend

!== Special treatment for k=N because the control volume at the top is Hz(N)/2
          qdm_nstp = wz(i,j,N,nstp)*Hz    (i,j,N)
          qdm_indx = wz(i,j,N,indx)*Hz_bak(i,j,N)
          wz(i,j,N,nnew)=( cff1*qdm_nstp
     &                    +cff2*qdm_indx
     &                    +2.*cdt*pn(i,j)*pm(i,j)
     &                       *(rw(i,j,N)+rw_nbq_avg1(i,j,N))
     &                   )/ Hz_half(i,j,N)          !<-- wz has units m.s-1

          rw_int_nbq(i,j,N)=(qdm_nstp-qdm_indx)*0.5/dt
     &                     -  pn(i,j)*pm(i,j)*rw_nbq_avg2(i,j,N)
          wz(i,j,N,indx)=wz(i,j,N,nstp)*0.5*Hz(i,j,N)   !<-- wz(indx) has units kg.m-1.s-1

!== Special treatment for k=0 because the control volume at the bottom
!    is Hz(1)/2
#  ifdef NBQ_FREESLIP
!         wz(i,j,0,nnew)=( cff1*wz(i,j,0,nstp)*Hz    (i,j,1)
!     &                   +cff2*wz(i,j,0,indx)*Hz_bak(i,j,1)
!     &                   +2.*cdt*pn(i,j)*pm(i,j)
!     &                   *(rw(i,j,0)+rw_nbq_avg1(i,j,0))
!     &                  )/ Hz_half(i,j,1)         !<-- wz has units m.s-1
!        rw_int_nbq(i,j,0)=pn(i,j)*pm(i,j)*rw(i,j,0)  !! Francis: switch from k=0 to k=1 here

          rw_int_nbq(i,j,0)=0.

          wz(i,j,0,nnew)= (  u(i,j,1,nnew) * pm_u(i,j)
     &                                  * ( z_w(i  ,j,0)-z_w(i-1,j,0) )
     &                     + v(i,j,1,nnew) * pm_v(i,j)
     &                                  * ( z_w(i,j  ,0)-z_w(i,j-1,0) )  )
          wz(i,j,0,indx)=wz(i,j,0,nnew)*0.5*Hz(i,j,1)   !<-- wz(indx) has units kg.m-1.s-1
#  else
          wz(i,j,0,nnew)=0.
          wz(i,j,0,indx)=0.
#  endif
        enddo
!==
      endif

      enddo
!==

!==
# endif /* NBQ */

# if defined M3FAST && defined M2FILTER_NONE
#  undef ru_nbq_avg1
#  undef rv_nbq_avg1
#  ifdef NBQ
#   undef rw_nbq_avg1
#  endif
# endif
!
!--------------------------------------------------------------------
!

# if defined TS_HADV_TEST && ( defined DIAGONAL_ADV || defined SOLID_BODY_ROT )
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
            u(i,j,k,nnew) = u(i,j,k,nstp)
          enddo
        enddo
        do j=JstrR,JendR
          do i=IstrR,IendR
            v(i,j,k,nnew) = v(i,j,k,nstp)
          enddo
        enddo
      enddo
# endif
# if defined TS_HADV_TEST && defined SOLID_BODY_PER
      cff = cos(2.*pi*(time+dt/2.)/(float(ntimes)*dt))
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
            u(i,j,k,nnew) = u(i,j,k,nstp) * cff
          enddo
        enddo
        do j=JstrR,JendR
          do i=IstrR,IendR
            v(i,j,k,nnew) = v(i,j,k,nstp) * cff
          enddo
        enddo
      enddo
# endif
!$acc end kernels
!
!--------------------------------------------------------------------
!
!======================================================================
! Set PHYSICAL lateral boundary conditions for tracers.
!======================================================================
!
# if defined TRACERS
      do itrc=1,NT
        call t3dbc_tile (Istr,Iend,Jstr,Jend, nnew,itrc, WORK)
#  if defined BD_EXCHANGE
#   ifdef THREE_GHOST_POINTS_TS
        call exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                               t(START_2D_ARRAY,1,nnew,itrc))
#   else
        call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                          t(START_2D_ARRAY,1,nnew,itrc))
#   endif
#  endif
      enddo
# endif

      call u3dbc_tile (Istr,Iend,Jstr,Jend, WORK)
      call v3dbc_tile (Istr,Iend,Jstr,Jend, WORK)
# ifdef NBQ
      call w3dbc_tile (Istr,Iend,Jstr,Jend, WORK)
# endif
!
!======================================================================
! Coupling, include ghost points associated with PHYSICAL
! boundaries ONLY. Do not touch periodic ghost points or
! internal computational margins (MPI code).
!======================================================================
!
# 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 PHYSICAL lateral boundary conditions for momentum.
!======================================================================
!
# if defined BD_EXCHANGE
#  ifdef THREE_GHOST_POINTS_UV
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             u(START_2D_ARRAY,1,nnew))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             v(START_2D_ARRAY,1,nnew))
#   ifdef NBQ
      call exchange_w3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             wz(START_2D_ARRAY,0,nnew))
#   endif
#  else
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        u(START_2D_ARRAY,1,nnew))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        v(START_2D_ARRAY,1,nnew))
#   ifdef NBQ
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
     &                        wz(START_2D_ARRAY,0,nnew))
#   endif
#  endif
# endif

!
!======================================================================
! Prepare viscosity/diffusivity array for GLS_MIXING
! (needed for OPENMP parallelization)
!======================================================================
!
!$acc kernels if(compute_on_device) default(present)

# ifdef GLS_MIXING
      do k=0,N
        do j=Jstr-1,Jend+1
          do i=Istr-1,Iend+1
            Akv_old(i,j,k)=Akv(i,j,k)
#  if defined TEMPERATURE
            Akt_old(i,j,k)=Akt(i,j,k,itemp)
#  elif defined SALINITY
            Akt_old(i,j,k)=Akt(i,j,k,isalt)
#  endif
          enddo
        enddo
      enddo
# endif
!
!======================================================================
! Prepare to start two-dimensional time stepping:
! set the initial values of the fast-time-step free-surface
! field to its fast-time-averaged values corresponding
! to the time step n (nstp).
!======================================================================
!
# ifndef M2FILTER_NONE 
#  if defined BD_EXCHANGE
      do j=JU_RANGE
        do i=IV_RANGE
          zeta(i,j,knew)=Zt_avg1(i,j)
        enddo
      enddo
!$acc end kernels
       call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          zeta(START_2D_ARRAY,knew))
#  else
      do j=JstrR,JendR
        do i=IstrR,IendR
          zeta(i,j,knew)=Zt_avg1(i,j)
        enddo
      enddo     
!$acc end kernels
#  endif
# else
!$acc end kernels
# endif
# undef IU_RANGE
# undef JU_RANGE
# undef IV_RANGE
# undef JV_RANGE
#else
      subroutine pre_step3d_empty
#endif
      return
      end
