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

      subroutine step3d_t (tile)
!
      implicit none
      integer tile
# ifdef TRACERS
      integer 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 step3d_t_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),
     &                                                A3d(1,1,trd)
#  ifdef OPENACC
     &   ,A3d(1,2,trd),A3d(1,3,trd),A3d(1,4,trd),A3d(1,5,trd)
#  endif
     &                            )
# endif /* TRACERS */      
      return
      end

# ifdef TRACERS
      subroutine step3d_t_tile (Istr,Iend,Jstr,Jend,
     &                          FX,FE, WORK, FC,CF,BC,DC,EC,GC, swdk
#  ifdef OPENACC
     &                           ,FX_3D,FE_3D,WORK_3D,WORK2_3D
#  endif
     &  )

#  ifdef SUBSTANCE
      USE comsubstance, ONLY : nv_grav, nvp, ws_part, nv_sand
#   ifdef MUSTANG
      USE comMUSTANG,  ONLY : flx_s2w_CROCO, flx_w2s_CROCO
      USE comMUSTANG,  ONLY : flx_w2s_sum_CROCO
      USE comMUSTANG,  ONLY : corflux, corfluy
      USE comsubstance, ONLY : irkm_var_assoc, nvpc, nvp
#    ifdef key_sand2D
      USE comsubstance, ONLY : l_subs2D
#    endif
#   endif
#   ifdef SUBSTANCE_SUBMASSBALANCE
      USE submassbalance, ONLY :  submassbalance_flxcum_kiv
      USE submassbalance, ONLY :  submassbalance_psource
      USE submassbalance, ONLY :  submassbalance_flxcum_ws
#   endif
#  endif
      implicit none
#  include "param.h"
#  include "ncscrum.h"
#  include "grid.h"
#  include "ocean3d.h"
#  include "mixing.h"
#  include "climat.h"
#  include "scalars.h"
#  include "sources.h"
#  include "forces.h"
#  ifdef DIAGNOSTICS_TS
#   include "diagnostics.h"
#  endif
#  ifdef DIAGNOSTICS_PV
#   include "diags_pv.h"
#  endif
      integer Istr,Iend,Jstr,Jend, itrc, i,j,k, indx, kmld
     &       ,imin,imax,jmin,jmax,iAkt,nadv
#  if defined PSOURCE || defined PSOURCE_MASS
     &       ,is
#  endif
#  ifdef PSOURCE_MASS
     &       ,iis,jjs
#  endif
#  ifdef SUBSTANCE
      integer isubstep,nbsubstep
      real :: cflm,cfl_loc,dz_cflm,ws_cflm
#  endif
#  ifdef AGRIF
#   include "zoom.h"
#  endif
      real FX(PRIVATE_2D_SCRATCH_ARRAY),
     &     FE(PRIVATE_2D_SCRATCH_ARRAY),   cff,
     &     WORK(PRIVATE_2D_SCRATCH_ARRAY), epsil,
     &     FC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     CF(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     BC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     DC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     EC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     GC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &     swdk(PRIVATE_2D_SCRATCH_ARRAY,0:N)
#  ifdef DIAGNOSTICS_TS_MLD_DENS
      integer kblt(PRIVATE_2D_SCRATCH_ARRAY)
#  endif
#  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
      real cff1,cff2,gama,dRz,hbltmp,sig,dXmax,dEmax,
     &     dpth,smax,amax,amaxx
#  ifdef TS_VADV_FCT
      real tmax,tmin,Ppos,Qpos,Pneg,Qneg,lmdmax
#  endif
#  include "openacc_work.h"
      parameter (epsil=1.E-16)
!
#  ifdef MASKING
#   define SWITCH *
#  else
#   define SWITCH !
#  endif
!
#  if defined TS_HADV_UP5 || defined TS_HADV_C6 || \
     defined TS_HADV_WENO5 || defined BIO_HADV_WENO5 || \
     defined TS_VADV_WENO5
!
!--------------------------------------------------------------------
! 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
!--------------------------------------------------------------------
!
      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.0

      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.0

      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.0

      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.0
#  endif
#  if defined TS_HADV_WENO5 || defined BIO_HADV_WENO5 || defined TS_VADV_WENO5
      REAL    :: flux3_weno, flux5_weno
#  endif
!
#  include "compute_auxiliary_bounds.h"

      nadv = 3   !<-- do not remove: used in compute_tracer_fluxes.h
!
!======================================================================
!
!  Compute horizontal advection
!
!======================================================================
!

      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 .eq. 1) then
#   endif
#  endif
!
!----------------------------------------------------------
! Compute fluxes FX and FE
!----------------------------------------------------------
!
#  undef PREDICTOR
#  if !defined OPENACC
#   include "compute_horiz_tracer_fluxes.h"
#  else
#   include "compute_horiz_tracer_fluxes_3D.h"
#  endif
!
!----------------------------------------------------------
! 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
!
#  if defined MUSTANG && defined MUSTANG_CORFLUX
          if (k.eq.1 .and.
     &        (itrc.ge.itsubs1+nv_grav .and.
     &        itrc.le.itsubs1+nv_grav+nv_sand-1) ! isand1, isand2
     &        .or.
     &        (itrc.ge.itsubs1+nvpc .and.
     &        itrc.le.itsubs1+nvp-1) ! nvpc+1, nvp
     &        ) then
            do j = Jstr, Jend
              do i = Istr, Iend+1
                FX(i, j) = FX(i, j) *
     &                     corflux(itrc-itsubs1+1, i, j)
              enddo
            enddo
            do j = Jstr, Jend+1
              do i = Istr, Iend
                FE(i, j) = FE(i, j) *
     &                     corfluy(itrc-itsubs1+1, i, j)
              enddo
            enddo
          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)
#   if defined MUSTANG && defined key_sand2D
!     To do: sum FX for k=1
#   endif
#   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,3,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,3,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,3,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,3,itrc)
                  endif
#   endif
                endif
              endif
            endif
#   if defined OPENACC
            ENDDOEXTEND
#   endif
          enddo
#  endif /* PSOURCE */
!
!----------------------------------------------------------
! Finalize horizontal advection: compute flux divergences
!----------------------------------------------------------
!
!$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
#  if !(defined MUSTANG && defined key_sand2D)
              t(i,j,k,nnew,itrc)=Hz_bak(i,j,k)*t(i,j,k,nstp,itrc)
     &                    -dt*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j)
     &                                         +FE(i,j+1)-FE(i,j))
#  else
              if (l_subs2D(itrc-itsubs1+1)) then
                if (k.eq.1) then
                  t(i,j,k,nnew,itrc)=
     &                         ( Hz_bak(i,j,k)*t(i,j,k,nstp,itrc)
     &                     -dt*pm(i,j)*pn(i,j)*(FX(i+1,j)-FX(i,j)
     &                                         +FE(i,j+1)-FE(i,j))
     &                               +dt*flx_s2w_CROCO(i,j,itrc) )
     &                     /(Hz(i,j,k)+dt*flx_w2s_CROCO(i,j,itrc))
                  flx_w2s_sum_CROCO(i,j,itrc)=flx_w2s_CROCO(i,j,itrc)
     &                                         *dt*t(i,j,1,nnew,itrc)
                endif
              else
                t(i,j,k,nnew,itrc)=
     &                           Hz_bak(i,j,k)*t(i,j,k,nstp,itrc)
     &                    -dt*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j)
     &                                         +FE(i,j+1)-FE(i,j))
              endif
#  endif
            enddo
          enddo                 !-->  discard FX,FE

#  if defined AGRIF && defined AGRIF_CONSERV_TRA
          if (Agrif_Root()) then
            MYFX(IstrR:IendR,JstrR:JendR,k,itrc)=
     &        dt*FX(IstrR:IendR,JstrR:JendR)
            MYFY(IstrR:IendR,JstrR:JendR,k,itrc)=
     &        dt*FE(IstrR:IendR,JstrR:JendR)
          else
            MYFX(IstrR:IendR,JstrR:JendR,k,itrc)=
     &       MYFX(IstrR:IendR,JstrR:JendR,k,itrc)
     &       +dt*FX(IstrR:IendR,JstrR:JendR)
            MYFY(IstrR:IendR,JstrR:JendR,k,itrc)=
     &       MYFY(IstrR:IendR,JstrR:JendR,k,itrc)
     &      +dt*FE(IstrR:IendR,JstrR:JendR)
          endif
#  endif

!add submassbalance !! in loop k
! this has to be checked with OPENACC
! as well as OPENACC and MUSTANG in general          
#  if defined SUBSTANCE_SUBMASSBALANCE
        if (itrc .ge. itsubs1) then
            call submassbalance_flxcum_kiv(Istr, Iend, Jstr, Jend, 
     &        FX, FE, itrc-itsubs1+1)
        endif
#  endif
#  if defined OPENACC
            ENDDOEXTEND
# endif  
!
!----------------------------------------------------------
! Store diagnostic transport term
!----------------------------------------------------------
!
#  ifdef DIAGNOSTICS_TS
          do j=Jstr,Jend
            do i=Istr,Iend
              TXadv(i,j,k,itrc)=-(FX(i+1,j)-FX(i,j)) SWITCH rmask(i,j)
              TYadv(i,j,k,itrc)=-(FE(i,j+1)-FE(i,j)) SWITCH rmask(i,j)
          enddo
        enddo
!
#   if defined TS_HADV_WENO5 || defined TS_HADV_UP5 || defined TS_HADV_UP3
!
! Compute centered scheme fluxes FX and FE to separate
! advection and diffusion parts of the transport term
!
#    undef PREDICTOR

#    if defined TS_HADV_WENO5
#     define TS_HADV_C6
#     undef  TS_HADV_WENO5
#     include "compute_horiz_tracer_fluxes.h"
#     undef  TS_HADV_C6
#     define TS_HADV_WENO5

#    elif defined TS_HADV_UP5
#     define TS_HADV_C6
#     undef  TS_HADV_UP5
#     include "compute_horiz_tracer_fluxes.h"
#     undef  TS_HADV_C6
#     define TS_HADV_UP5

#    elif defined TS_HADV_UP3
#     define TS_HADV_C4
#     undef  TS_HADV_UP3
#     include "compute_horiz_tracer_fluxes.h"
#     undef  TS_HADV_C4
#     define TS_HADV_UP3

#    endif
        do j=Jstr,Jend
          do i=Istr,Iend
              THmix(i,j,k,itrc)=( TXadv(i,j,k,itrc)+TYadv(i,j,k,itrc)
     &                          +(FX(i+1,j)-FX(i,j))
     &                          +(FE(i,j+1)-FE(i,j)) ) SWITCH rmask(i,j)
              TXadv(i,j,k,itrc)=-(FX(i+1,j)-FX(i,j))   SWITCH rmask(i,j)
              TYadv(i,j,k,itrc)=-(FE(i,j+1)-FE(i,j))   SWITCH rmask(i,j)
          enddo
        enddo
#   else
        do j=Jstr,Jend
          do i=Istr,Iend
              THmix(i,j,k,itrc)=0.
          enddo
        enddo
#   endif
#  endif /* DIAGNOSTICS_TS */
#  if !defined OPENACC
        enddo    !<-- k
#  endif
      enddo    !<-- itrc

!
!----------------------------------------------------------
! Compute fraction of the solar shortwave flux "swdk"
! penetrating to grid level depth (at vertical w-points).
! (swdk is used later in this routine)
!----------------------------------------------------------
!
#  if defined LMD_SKPP || defined LMD_BKPP || defined GLS_MIXING
#   define wrk1 FX
#   define wrk2 FE
      do k=1,N-1
!$acc parallel loop if(compute_on_device) default(present) async(1)
        do j=Jstr,Jend
          do i=Istr,Iend
            wrk1(i,j)=z_w(i,j,k)-z_w(i,j,N)
          enddo
        enddo
        call lmd_swfrac_tile (Istr,Iend,Jstr,Jend,1.,wrk1,wrk2)
!$acc parallel loop if(compute_on_device) default(present) async(1)
        do j=Jstr,Jend
          do i=Istr,Iend
            swdk(i,j,k)=wrk2(i,j)
          enddo
        enddo
      enddo
!$acc wait
#   undef wrk1
#   undef wrk2
#  endif /* LMD_SKPP || LMD_BKPP || GLS_MIXING */
!
!----------------------------------------------------------
! Add tracer divergence due to cell-centered point sources
!----------------------------------------------------------
!
#  ifdef PSOURCE_MASS
        do itrc=1,NT
          do is=1,Nsrc
#   ifdef MPI
           iis=Isrc_mpi(is,mynode)
           jjs=Jsrc_mpi(is,mynode)
#   else
           iis=Isrc(is)
           jjs=Jsrc(is)
#   endif
           if (Istr.le.iis .and. iis.le.Iend+1
     &        .and. Jstr.le.jjs .and. jjs.le.Jend+1) then
              cff1=dt*pm(iis,jjs)*pn(iis,jjs)
              do k=1,N
                if (Lsrc(is,itrc)) then
                  cff2=Tsrc(is,k,itrc)
                else
                  cff2=t(iis,jjs,k,3,itrc)
                endif
                 t(iis,jjs,k,nnew,itrc)=t(iis,jjs,k,nnew,itrc)+
     &           cff1*Qsrc(is,k)*cff2
              enddo
           endif
          enddo ! source
        enddo ! itrc
#  endif /* PSOURCE_MASS */

!
!======================================================================
!
!  Compute vertical advection
!
!======================================================================
!
!$acc kernels if(compute_on_device) default(present)
      DOLOOP2D(Istr,Iend,Jstr,Jend)
        do itrc=1,NT
!
!----------------------------------------------------------
! Compute vertical fluxes FC
!----------------------------------------------------------
!
#  undef PREDICTOR
#  include "compute_vert_tracer_fluxes.h"

#  ifdef TS_HADV_TEST
        do i=Istr,Iend
          FC(i,0:N)=0.
        enddo
#  endif

!----------------------------------------------------------
! Apply flux correction for monotonicity
!----------------------------------------------------------
!
#  ifdef TS_VADV_FCT
#   include "FCT.h"
#  endif
!
!----------------------------------------------------------
! Apply vertical advective flux divergence
!----------------------------------------------------------
!
#  if defined MUSTANG && defined key_sand2D
          if (.not.l_subs2D(itrc-itsubs1+1)) then
#  endif
            do k=1,N
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)
     &                             -CF(i,0)*(FC(i,k)-FC(i,k-1))
              enddo
            enddo               !--> discard FC
#  ifdef DIAGNOSTICS_TS
            do k=1,N
              do i=Istr,Iend
                TVadv(i,j,k,itrc) = -(FC(i,k)-FC(i,k-1))
     &                                 SWITCH rmask(i,j)
              enddo
            enddo
#  endif /* DIAGNOSTICS_TS */
#  if defined MUSTANG && defined key_sand2D
          endif
#  endif
!
!
!======================================================================
!
!  Compute surface and Bottom forcing
!
!======================================================================
!
!----------------------------------------------------------
! Add top and bottom fluxes
!----------------------------------------------------------
!
          do i=Istr,Iend
            FC(i,N)=dt*stflx(i,j,itrc)
            FC(i,0)=-dt*btflx(i,j,itrc)
          enddo
!
!----------------------------------------------------------
! Add solar radiation flux in temperature equation
! Also compute the nonlocal transport flux for unstable
! (convective) forcing conditions into matrix DC when using
! the Large et al. 1994 KPP scheme.
!----------------------------------------------------------
!
#  ifdef TEMPERATURE
          if (itrc.eq.itemp) then
            do k=1,N-1
              do i=Istr,Iend
                FC(i,k)=0.
#   if defined LMD_SKPP || defined LMD_BKPP \
                        || defined GLS_MIXING
     &    +dt*srflx(i,j)*swdk(i,j,k)
#    ifdef LMD_NONLOCAL
     &    -dt*ghats(i,j,k)*(stflx(i,j,itemp)-srflx(i,j))
#    endif
#   endif
              enddo
            enddo
          endif
#  endif /* TEMPERATURE */
#  ifdef SALINITY
          if (itrc.eq.isalt) then
            do k=1,N-1
              do i=Istr,Iend
                FC(i,k)=0.
#   if defined LMD_SKPP || defined LMD_BKPP
#    ifdef LMD_NONLOCAL
     &    -dt*ghats(i,j,k)*stflx(i,j,isalt)
#    endif
#   endif
              enddo
            enddo
          endif
#  endif /* SALINITY */
!
#  if defined SALINITY || defined TEMPERATURE
#   if defined SALINITY && defined TEMPERATURE
          if (itrc.eq.itemp .or. itrc.eq.isalt) then
#   elif defined SALINITY
          if (itrc.eq.isalt) then
#   elif defined TEMPERATURE
          if (itrc.eq.itemp) then
#   endif
            do k=1,N
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+FC(i,k )
     &                                               -FC(i,k-1)
#   if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV
                TForc(i,j,k,itrc)=(FC(i,k)-FC(i,k-1))
     &                                    /(dt*pm(i,j)*pn(i,j))
     &                                        SWITCH rmask(i,j)
#   endif  /* DIAGNOSTICS_TS */
              enddo
            enddo
#   if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV
          else
            do k=1,N
              do i=Istr,Iend
                TForc(i,j,k,itrc)=0.
              enddo
            enddo
#   endif  /* DIAGNOSTICS_TS */
          endif
#  endif  /* defined SALINITY || TEMPERATURE  */
!
!
!======================================================================
!
! Compute vertical mixing (implicit step)
!
! If lateral rotated diffusion is chosen, this step is performed
! after lateral diffusion to include implicit correction to lateral
! diffusion
!
! If MUSTANG is used, mixing is done here for all passive tracers
! but sediment tracers include a settling process with substeps
!
!======================================================================
!
#  ifdef VADV_ADAPT_IMP
          do i=Istr,Iend
            DC(i,0)=dt*pn(i,j)*pm(i,j)
          enddo
#  endif

#  ifdef TS_MIX_IMP
#   if defined TEMPERATURE && defined SALINITY
          if (itrc .le. 2) then                  ! active tracers
#   elif defined TEMPERATURE || defined SALINITY
          if (itrc .le. 1) then                  ! active tracers
#   endif
            do k=1,N
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc) / Hz(i,j,k)
              enddo
            enddo
          else                                   ! passive tracers
!                                                 (sed. and others)
#   ifdef MUSTANG
!
            if (itrc.ge.itsubs1+nv_grav .and.
     &          itrc.le.itsubs1+nvp) then        ! sediment tracers
#    ifdef key_sand2D
              if (.not.l_subs2D(itrc-itsubs1+1)) then
#    endif
#    undef TS_MIX_IMP
#    define CD EC
#    define NBSUBSTEP GC
#    include "t3dmix_tridiagonal_settling.h"
#    undef NBSUBSTEP
#    undef CD
#    define TS_MIX_IMP
#    ifdef key_sand2D
              endif
#    endif
            else                                 ! other tracers
#    undef TS_MIX_IMP
#    include "t3dmix_tridiagonal.h"
#    define TS_MIX_IMP
            endif
!
#   else  /* MUSTANG */
!                                                ! passive tracers
#    ifdef BIO_HADV_WENO5
#     undef TS_MIX_IMP
#     include "t3dmix_tridiagonal.h"
#     define TS_MIX_IMP
#    else
            do k=1,N
              do i=Istr,Iend
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc) / Hz(i,j,k)
              enddo
            enddo
#    endif /* BIO_HADV_WENO5 */
#   endif  /* MUSTANG */
!
          endif ! itrc <= 2
!
#  else /* TS_MIX_IMP */

#   ifndef MUSTANG
!                                                ! all tracers
#    include "t3dmix_tridiagonal.h"
!
#   else
          if (itrc.ge.itsubs1+nv_grav .and.
     &        itrc.le.itsubs1+nvp) then          ! sediment tracers
#    ifdef key_sand2D
            if (.not.l_subs2D(itrc-itsubs1+1)) then
#    endif
#    define CD EC
#    define NBSUBSTEP GC
#    include "t3dmix_tridiagonal_settling.h"
#    undef NBSUBSTEP
#    undef CD
#    ifdef key_sand2D
            endif
#    endif
          else                                   ! all other tracers
#    include "t3dmix_tridiagonal.h"
          endif
#   endif /* MUSTANG */

#  endif /* TS_MIX_IMP */

!
!  Test of constant tracers
!
#  ifdef CONST_TRACERS
          do k=1,N
            do i=Istr,Iend
              t(i,j,k,nnew,itrc)=t(i,j,k,nstp,itrc)
            enddo
          enddo
#   if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV
          do k=1,N
            do i=Istr,Iend
               TVmix(i,j,k,itrc)=0.0
#    ifdef MASKING
     &              * rmask(i,j)
#    endif
            enddo
          enddo
#   endif /* DIAGNOSTICS_TS */
#  endif /* CONST_TRACERS */


        enddo      ! <-- itrc

!
!----------------------------------------------------------
! Compute vertical gradient in vertical T-diffusion
! coefficient for floats random walk.
!----------------------------------------------------------
!
#  if defined FLOATS && defined RANDOM_WALK
        do k=1,N
          do i=Istr,Iend
            dAktdz(i,j,k)=(Akt(i,j,k,1)-Akt(i,j,k-1,1))/Hz(i,j,k)
          enddo
        enddo
#  endif

!
      ENDDOLOOP2D ! <-- j


!add submassbalance for river!! 
!mustang and openacc not compatible         
# if defined SUBSTANCE_SUBMASSBALANCE
      call submassbalance_psource(Istr, Iend, Jstr, Jend)
#  ifdef MUSTANG
      call submassbalance_flxcum_ws(Istr, Iend, Jstr, Jend)
#  endif
# endif

!$acc end kernels
!======================================================================
!  Set lateral boundary conditions; nudge toward tracer climatology;
!  apply land-sea mask and exchange periodic boundary conditions.
!======================================================================
!
      do itrc=1,NT
        call t3dbc_tile (Istr,Iend,Jstr,Jend, nnew,itrc, WORK)
      enddo
!
!======================================================================
!  Compute biological fluxes
!======================================================================
!
#  ifdef BIOLOGY
      call biology_tile (Istr,Iend,Jstr,Jend)
#  endif


!
!======================================================================
!  Compute inverse vertical Rho gradients for use in t3dmix
!======================================================================
!
#  if !defined TS_MIX_S && \
      (defined TS_DIF4 || defined TS_DIF2 || defined SPONGE_DIF2)
!$acc kernels if(compute_on_device) default(present)
#   ifdef TS_MIX_GEO
!
! --- Geopotential diffusion: idRz contains idZ values ---
!
      do k=1,N-1
        do j=jstr,jend
          do i=istr,iend
            idRz(i,j,k) = 1./(z_r(i,j,k+1)-z_r(i,j,k))
          enddo
        enddo
      enddo
      do j=jstr,jend
        do i=istr,iend
          idRz(i,j,N) = 0.
          idRz(i,j,0) = 0.
        enddo
      enddo

#   elif defined TS_MIX_ISO
!
!  --- Isopycnal diffusion: idRz is inverse vertical rho gradient ---
!
      do k=1,N-1
        do j=jstr,jend
          do i=istr,iend
#    ifdef SPLIT_EOS
            dpth=z_w(i,j,N)-0.5*(z_r(i,j,k+1)+z_r(i,j,k))
            dRz =rho1(i,j,k+1)-rho1(i,j,k)
     &          +(qp1(i,j,k+1)- qp1(i,j,k))
     &               *dpth*(1.- qp2*dpth)
#    else
            dRz = rho1(i,j,k+1)-rho1(i,j,k)
#    endif
            cff  = 1./min(dRz, -1.E-14)    !  minimum stratification imposed
            cff1 = 1./(z_r(i,j,k+1)-z_r(i,j,k))
#    ifdef LMD_SKPP2005
            hbltmp=hbls(i,j,3-nstp)
#    elif defined LMD_SKPP1994 || defined GLS_MIXING
            hbltmp=hbl(i,j)
#    else
            hbltmp=10.
#    endif
            gama = 1.                                      ! mixing becomes
            sig = (z_w(i,j,N)-z_w(i,j,k))/max(hbltmp,10.)  ! isosigma within
            if (sig .lt. 1.) gama = sig*sig*(3.-2.*sig)    ! PBL depth (>10m)
            dXmax = max(abs(dRdx(i,j,k  )),abs(dRdx(i+1,j,k  )),
     &                  abs(dRdx(i,j,k+1)),abs(dRdx(i+1,j,k+1)),1E-14)
            dEmax = max(abs(dRde(i,j,k  )),abs(dRde(i,j+1,k  )),
     &                  abs(dRde(i,j,k+1)),abs(dRde(i,j+1,k+1)),1E-14)
            smax=min(Rslope_max,Gslope_max*min(pm(i,j),pn(i,j))/cff1)
!
            idRz(i,j,k) = max(cff,   -smax*gama*cff1 / dXmax,
     &                               -smax*gama*cff1 / dEmax )
          enddo
        enddo
      enddo
      do j=jstr,jend
        do i=istr,iend
          idRz(i,j,N) = 0.
          idRz(i,j,0) = 0.
        enddo
      enddo
#   endif /* TS_MIX_ISO */
!
!$acc end kernels
#   if defined BD_EXCHANGE
      call exchange_w3d_tile (  istr,iend,jstr,jend,
     &                          idRz(START_2D_ARRAY,0)  )
#   endif
!
#  endif /* TS_DIF2 || TS_DIF4 || SPONGE_DIF2 */
!
! 

!
!======================================================================
!  Compute Nudging terms
!======================================================================
!
#  if defined TNUDGING && defined TCLIMATOLOGY
#   ifdef ZONAL_NUDGING
      do itrc=1,NT   !<-- itrc
        if (iic.eq.ntstart .or. mod(iic,10).eq.0) then
          call zonavg_3d(istr,iend,jstr,jend,
     &               t(START_2D_ARRAY,1,nnew,itrc),tzon(:,:,itrc))
        endif
        if (iic.eq.ntstart) then
          call zonavg_3d(istr,iend,jstr,jend,
     &              tclm(START_2D_ARRAY,1,itrc),tclmzon(:,:,itrc))
        endif
      enddo   ! <-- itrc
#   endif /* ZONAL_NUDGING */
#  endif
!$acc kernels if(compute_on_device) default(present)
!$acc loop seq
      do itrc=1,NT   !<-- itrc
!
!======================================================================
!  Compute Nudging terms
!======================================================================
!
#  if defined TNUDGING && defined TCLIMATOLOGY
        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
              t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)
     &                                   +dt*Tnudgcof(i,j,k,itrc)*
#   ifdef ZONAL_NUDGING
     &                          (tclmzon(j,k,itrc)-tzon(j,k,itrc))
#   else
     &                       (tclm(i,j,k,itrc)-t(i,j,k,nnew,itrc))
#   endif /* ZONAL_NUDGING */
#   if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV
              TForc(i,j,k,itrc)=(TForc(i,j,k,itrc)
     &                                       +Tnudgcof(i,j,k,itrc)
     &                                          *(tclm(i,j,k,itrc)
#    ifdef ZONAL_NUDGING
     &                                            -tzon(j,k,itrc))
#    else
     &                                        -t(i,j,k,nnew,itrc))
#    endif /* ZONAL_NUDGING */
     &                             *(Hz(i,j,k)/(pm(i,j)*pn(i,j))))
     &                                           SWITCH rmask(i,j)
#   endif  /* DIAGNOSTICS_TS */
            enddo
          enddo
        enddo
#  endif /* TNUDGING && TCLIMATOLOGY */
!
!======================================================================
! Compute the tendency term of tracer diagnostics
! Divide all diagnostic terms by the cell volume
! (Hz(i,j,k,itrc)/(pm(i,j).*pn(i,j)). There after the unit
! of diagnostic terms will be: (unit of tracers)* s-1.
!
! Note: the Horizontal mixing term is computed in t3dmix
! where Trate is updated accordingly
!======================================================================
!
#  if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV
        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
              Trate(i,j,k,itrc)=(Hz(i,j,k)*t(i,j,k,nnew,itrc)
     &                          -Hz_bak(i,j,k)*t(i,j,k,nstp,itrc))
     &                                       /(dt*pm(i,j)*pn(i,j))
     &                                           SWITCH rmask(i,j)
!
#   ifdef DIAGNOSTICS_TS_ADV
              TXadv(i,j,k,itrc)=TXadv(i,j,k,itrc)+
     &           t(i,j,k,3,itrc)*(Huon(i+1,j,k)-Huon(i,j,k))
              TYadv(i,j,k,itrc)=TYadv(i,j,k,itrc)+
     &           t(i,j,k,3,itrc)*(Hvom(i,j+1,k)-Hvom(i,j,k))
              TVadv(i,j,k,itrc)=TVadv(i,j,k,itrc)-
     &           t(i,j,k,3,itrc)*(Huon(i+1,j,k)-Huon(i,j,k)+
     &                            Hvom(i,j+1,k)-Hvom(i,j,k))
#   endif /* DIAGNOSTICS_TS_ADV */
!
              cff=pm(i,j)*pn(i,j)/Hz(i,j,k)
#   ifdef DIAGNOSTICS_TS
              TXadv(i,j,k,itrc)=TXadv(i,j,k,itrc)*cff
              TYadv(i,j,k,itrc)=TYadv(i,j,k,itrc)*cff
              TVadv(i,j,k,itrc)=TVadv(i,j,k,itrc)*cff
#   endif
              TForc(i,j,k,itrc)=TForc(i,j,k,itrc)*cff
              Trate(i,j,k,itrc)=Trate(i,j,k,itrc)*cff
              TVmix(i,j,k,itrc)=TVmix(i,j,k,itrc)*cff
              THmix(i,j,k,itrc)=THmix(i,j,k,itrc)*cff

#   if defined DIAGNOSTICS_TSVAR && !defined SPONGE_DIF2 \
             && !defined TS_DIF2 && !defined TS_DIF4
#    include "finalize_diagnostics_tsadv.h"
#   endif
            enddo
          enddo
        enddo
!
!----------------------------------------------------------
! Compute tracer diagnostics averaged over the MLD
!----------------------------------------------------------
!
#   ifdef DIAGNOSTICS_TS_MLD
#    define T_mld_nnew FX
#    define T_mld_nstp FE
#    ifdef DIAGNOSTICS_TS_MLD_DENS
         do j=Jstr,Jend
            do i=Istr,Iend
               kbl1(i,j)=1
               kblt(i,j)=1
               cff=0
               cff1=0
               do k=N,1,-1
                  if (cff.eq.0.and.rho(i,j,k).lt.rho(i,j,N)+mld_crit_D) then
                     kbl1(i,j)=k
                  endif
                  if (rho(i,j,k).ge.rho(i,j,N)+mld_crit_D)  cff=1

                  if (cff1.eq.0.and.t(i,j,k,nnew,1).gt.
     &               t(i,j,N,nnew,1)-mld_crit_T) then
                     kblt(i,j)=k
                  endif
                  if (t(i,j,k,nnew,1).le.t(i,j,N,nnew,1)-mld_crit_T)  cff1=1
               enddo
               if (itrc.eq.1) Tcrit_mld(i,j,itrc)=
     &                        z_w(i,j,N)-z_r(i,j,kbl1(i,j))
               Tcrit_mld(i,j,1)=z_r(i,j,N)-z_r(i,j,kbl1(i,j))
               if (itrc.eq.2) Tcrit_mld(i,j,itrc)=
     &            (z_r(i,j,N)-z_r(i,j,kblt(i,j)))-Tcrit_mld(i,j,1)
            enddo
         enddo
#    endif /* DIAGNOSTICS_TS_MLD_DENS */
        do j=Jstr,Jend
          do i=Istr,Iend
            TXadv_mld(i,j,itrc)=0.
            TYadv_mld(i,j,itrc)=0.
            TVadv_mld(i,j,itrc)=0.
            TVmix_mld(i,j,itrc)=0.
            THmix_mld(i,j,itrc)=0.
            TForc_mld(i,j,itrc)=0.
            Trate_mld(i,j,itrc)=0.
            Tentr_mld(i,j,itrc)=0.
            T_mld_nnew(i,j)=0.
            T_mld_nstp(i,j)=0.
          enddo
        enddo
        do j=Jstr,Jend
          do i=Istr,Iend
#    ifdef LMD_SKPP
            kmld=kbl(i,j)
#     ifdef DIAGNOSTICS_TS_MLD_DENS
           kmld=kbl1(i,j)
#     endif
#    else
            kmld=N-5
#    endif
            do k=N,kmld,-1
              cff=Hz(i,j,k)/(z_w(i,j,N)-z_w(i,j,kmld-1))
              TXadv_mld(i,j,itrc)=TXadv_mld(i,j,itrc)+
     &                            TXadv(i,j,k,itrc)*cff
#  ifdef MASKING
     &                   *rmask(i,j)
#  endif
              TYadv_mld(i,j,itrc)=TYadv_mld(i,j,itrc)+
     &                            TYadv(i,j,k,itrc)*cff
#  ifdef MASKING
     &                   *rmask(i,j)
#  endif
              TVadv_mld(i,j,itrc)=TVadv_mld(i,j,itrc)+
     &                            TVadv(i,j,k,itrc)*cff
#  ifdef MASKING
     &                   *rmask(i,j)
#  endif
              TVmix_mld(i,j,itrc)=TVmix_mld(i,j,itrc)+
     &                            TVmix(i,j,k,itrc)*cff
#  ifdef MASKING
     &                   *rmask(i,j)
#  endif
              TForc_mld(i,j,itrc)=TForc_mld(i,j,itrc)+
     &                            TForc(i,j,k,itrc)*cff
#  ifdef MASKING
     &                   *rmask(i,j)
#  endif
              T_mld_nnew(i,j)    =T_mld_nnew(i,j)+
     &                            t(i,j,k,nnew,itrc)*cff
            enddo
          enddo
        enddo

#    if (!defined TS_DIF2 && !defined TS_DIF4 && !defined SPONGE_DIF2) \
       || defined DIAGNOSTICS_DEBUG
!
! Compute entrainement/detrainement term.
! If diffusion terms are computed in t3dmix routines, entrainement
! terms are also computed in those routines
!
        do j=Jstr,Jend
          do i=Istr,Iend
#     if defined LMD_SKPP || defined GLS_MIXING
            if (kbl_nstp(i,j).eq.0) kbl_nstp(i,j)=kbl(i,j)
#      ifdef DIAGNOSTICS_TS_MLD_DENS
            if (kbl_nstp(i,j).eq.0) kbl_nstp(i,j)=kbl1(i,j)
#      endif
            kmld=kbl_nstp(i,j)
#     else
            kmld=N-5
#     endif
            do k=N,kmld,-1
              cff=Hz_bak(i,j,k)/(z_w(i,j,N)-z_w(i,j,kmld-1))
              T_mld_nstp(i,j)=T_mld_nstp(i,j)+
     &                        t(i,j,k,nstp,itrc)*cff
            enddo
            if (itrc .eq. NT) kbl_nstp(i,j)=kbl(i,j)
#     ifdef DIAGNOSTICS_TS_MLD_DENS
            if (itrc .eq. NT) kbl_nstp(i,j)=kbl1(i,j)
#     endif
          enddo
        enddo
        do j=Jstr,Jend
          do i=Istr,Iend
            Trate_mld(i,j,itrc)=(T_mld_nnew(i,j)-T_mld_nstp(i,j))/dt
            Tentr_mld(i,j,itrc)=Trate_mld(i,j,itrc)-
     &                          TXadv_mld(i,j,itrc)-
     &                          TYadv_mld(i,j,itrc)-
     &                          TVadv_mld(i,j,itrc)-
     &                          TVmix_mld(i,j,itrc)-
     &                          TForc_mld(i,j,itrc)
          enddo
        enddo
#    endif /* TS_DIF2 && TS_DIF4 && SPONGE_DIF2 */
#    undef T_mld_nnew
#    undef T_mld_nstp
#   endif /* DIAGNOSTICS_TS_MLD */

#   ifdef DIAGNOSTICS_DEBUG
        if (istr.eq.1 .and. jstr.eq.1
#    ifdef TEMPERATURE
     &       .and. itrc.eq.itemp
#    endif /* TEMPERATURE */
     &       ) then
          i=5
          j=5
#    if defined DIAGNOSTICS_TS_MLD
          cff=Trate_mld(i,j,itrc)-
     &        TXadv_mld(i,j,itrc)-
     &        TYadv_mld(i,j,itrc)-
     &        TVadv_mld(i,j,itrc)-
     &        TVmix_mld(i,j,itrc)-
     &        THmix_mld(i,j,itrc)-
     &        Tentr_mld(i,j,itrc)-
     &        TForc_mld(i,j,itrc)
             MPI_master_only write(stdout,'(A,2x,1pe10.3)')
     &       ' STEP3D_T: T budget closure MLD     :', cff
!
#    endif /* DIAGNOSTICS_TS_MLD */
          cff=Trate(i,j,N-5,itrc)-
     &        TXadv(i,j,N-5,itrc)-
     &        TYadv(i,j,N-5,itrc)-
     &        TVadv(i,j,N-5,itrc)-
     &        TVmix(i,j,N-5,itrc)-
     &        THmix(i,j,N-5,itrc)-
     &        TForc(i,j,N-5,itrc)
            MPI_master_only write(stdout,'(A,2x,1pe10.3)')
     &      ' STEP3D_T: T budget closure at k=N-5:', cff
        endif
#   endif /* DIAGNOSTICS_DEBUG */


#  endif /* DIAGNOSTICS_TS */
!
!
!======================================================================
!  Set PHYSICAL lateral boundary conditions.
!======================================================================
!
!----------------------------------------------------------
! Apply land mask
!----------------------------------------------------------
!
#  ifdef MASKING
#   ifdef EW_PERIODIC
#    define I_RANGE Istr,Iend
#   else
#    define I_RANGE IstrR,IendR
#   endif
#   ifdef NS_PERIODIC
#    define J_RANGE Jstr,Jend
#   else
#    define J_RANGE JstrR,JendR
#   endif
        do k=1,N
          do j=J_RANGE
            do i=I_RANGE
              t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*rmask(i,j)
            enddo
          enddo
        enddo

#   undef I_RANGE
#   undef J_RANGE
#  endif /* MASKING */

      enddo   ! <-- itrc
!$acc end kernels

      do itrc=1,NT
!
!----------------------------------------------------------
! Exchange periodic boundaries and computational margins.
!----------------------------------------------------------
!
#  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

      return
      end
# endif  /* TRACERS */
!
!======================================================================
!
# if defined BIO_HADV_WENO5 || \
     defined TS_HADV_WENO5  || defined TS_VADV_WENO5 || \
     defined UV_HADV_WENO5  || defined UV_VADV_WENO5 || \
     defined W_HADV_WENO5   || defined W_VADV_WENO5 || \
     defined BEDLOAD_WENO5
!
      function flux5_weno(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)
!$acc routine(flux5_weno) seq
!
!     WENO5: Weighted Essentially Non-Oscillatory scheme
!            with 5th-order accuracy
!
! This function computes tracer reconstruction at the grid cell's left
! edge (u-point i-1/2 or v-point j-1/2). WENO5 uses a convex combination
! of the polynomials reconstructed on the three ENO3 stencils in order
! to achieve higher accuracy on smooth profiles. Both left and right
! combinations of mirror symetric stencils around i-1/2 (j-1/2) points
! are computed. For stability, the upstream stencil combination,
! identified by the sign of ua (u or v), is selected.
!
! The scheme includes improvements from Borges et al., 2008 (WENO_Z), i.e.,
! smoothness indicators of higher order (with new non-oscillatory weights)
! that provide a scheme with less dissipation, higher resolution and better
! monotonicity preservation than the classical WENO5. In the meantime, WENO_Z
! removes the need to tune parameter Eps for added dissipation (needed in the
! original scheme due to suboptimal performance on critical points).
!
! References:
!  Guang-Shan JIANG and Chi-Wang SHU, 1996: Efficient Implementation of
!  Weighted ENO Schemes. JOURNAL OF COMPUTATIONAL PHYSICS 126, 202–228.
!
!  Rong WANG and Raymond J. SPITERI, 2007: Linear instability of the
!  fifth-order WENO method, SIAM J. Numer. Anal., 45, 1871-1901.
!
!  Borges R., M. Carmona, B. Costa, W.S. Don, 2008: An improved weighted
!  essentially non-oscillatory scheme for hyperbolic conservation laws.
!  Journal of Computational Physics 227 (2008) 3191–3211
!
!  Implementation: P. Marchesiello and J. Demange 2013
!
      implicit none
      REAL    :: flux5_weno
      REAL    :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua
      REAL    :: IS0, IS1, IS2
      REAl    :: d0, d1, d2
      REAl    :: a0, a1, a2
      REAL    :: w0, w1, w2
      REAL    :: p0, p1, p2
      REAL    :: Eps, cff1, cff2, T5

#  define WENO_Z

! Epsilon value for computing weigths
#  ifdef WENO_Z
      Eps = 1.e-40
#  else
!  --> from 1.e-7 (smoother) to 1.e-5 (sharper)
      Eps = 1.e-6
#  endif
! Weigths coefficients
      d0=1./10.
      d1=6./10.
      d2=3./10.
! other coefficients
      cff1=13./12.
      cff2=1./6.
!
      if (ua .ge. 0.) then        ! Take upstream stencils for stability
!
! === Reconstruction at u_i-1/2 using LEFT stencils ===
!
! Nonlinear Smoothness Indicators
        IS0 = cff1*(q_im3 - 2.*q_im2 +   q_im1)**2
     &      + 0.25*(q_im3 - 4.*q_im2 + 3*q_im1)**2
        IS1 = cff1*(q_im2 - 2.*q_im1 +   q_i  )**2
     &      + 0.25*(q_im2            -   q_i  )**2
        IS2 = cff1*(q_im1 - 2.*q_i   +   q_ip1)**2
     &   + 0.25*(3.*q_im1 - 4.*q_i   +   q_ip1)**2
#  ifdef WENO_Z
! Non-normalized stencil weights
! with parameter T5 of new smoothness indicator
        T5  = abs(IS2-IS0)
        a0  = d0*(1+T5/(Eps+IS0))
        a1  = d1*(1+T5/(Eps+IS1))
        a2  = d2*(1+T5/(Eps+IS2))
#  else
        a0  = d0/(Eps+IS0)**2
        a1  = d1/(Eps+IS1)**2
        a2  = d2/(Eps+IS2)**2
#  endif
! Normalized Weigths
        w0  = a0/(a0+a1+a2)
        w1  = a1/(a0+a1+a2)
        w2  = a2/(a0+a1+a2)
! Polynomials
        p0  = cff2*(2.*q_im3 - 7.*q_im2 + 11.*q_im1)
        p1  = cff2*(  -q_im2 + 5.*q_im1 +  2.*q_i  )
        p2  = cff2*(2.*q_im1 + 5.*q_i   -     q_ip1)
! Combination
        flux5_weno = w0*p0 + w1*p1 + w2*p2

      else
!
! === Reconstruction at u_i-1/2 using RIGHT stencils ===
!
! Nonlinear Smoothness Indicators
        IS0 = cff1*(q_ip2 - 2.*q_ip1 +   q_i  )**2
     &      + 0.25*(q_ip2 - 4.*q_ip1 + 3*q_i  )**2
        IS1 = cff1*(q_ip1 - 2.*q_i   +   q_im1)**2
     &      + 0.25*(q_ip1            -   q_im1)**2
        IS2 = cff1*(q_i   - 2.*q_im1 +   q_im2)**2
     &   + 0.25*(3.*q_i   - 4.*q_im1 +   q_im2)**2
#  ifdef WENO_Z
! Non-normalized stencil weights
! with parameter T5 of new smoothness indicator
        T5  = abs(IS2-IS0)
        a0  = d0*(1+T5/(Eps+IS0))
        a1  = d1*(1+T5/(Eps+IS1))
        a2  = d2*(1+T5/(Eps+IS2))
#  else
        a0  = d0/(Eps+IS0)**2
        a1  = d1/(Eps+IS1)**2
        a2  = d2/(Eps+IS2)**2
#  endif
! Normalized Weigths
        w0  = a0/(a0+a1+a2)
        w1  = a1/(a0+a1+a2)
        w2  = a2/(a0+a1+a2)
! Polynomials
        p0  = cff2*(2.*q_ip2 - 7.*q_ip1 + 11.*q_i  )
        p1  = cff2*(  -q_ip1 + 5.*q_i   +  2.*q_im1)
        p2  = cff2*(2.*q_i   + 5.*q_im1 -     q_im2)
! Combination
        flux5_weno = w0*p0 + w1*p1 + w2*p2

      endif

      return
      end
!
!======================================================================
!
      function flux3_weno( q_im2, q_im1, q_i, q_ip1, ua)
!$acc routine(flux3_weno) seq

!
!     WENO3: Weighted Essentially Non-Oscillatory of 3rd-order Accuracy
!
! This function computes tracer reconstruction at the grid cell's left
! edge (u-point i-1/2 or v-point j-1/2).
!
! Implementation: P. Marchesiello and J. Demange 2013
!
      implicit none
      REAL    :: flux3_weno
      REAL    :: q_im2, q_im1, q_i, q_ip1, ua
      REAL    :: IS0, IS1
      REAl    :: a0, a1
      REAL    :: w0, w1
      REAL    :: p0, p1
      REAL    :: Eps, d0,d1, T3
!
! Epsilon value for computing weigths
#  ifdef WENO_Z
      Eps = 1.e-40
#  else
!  --> from 1.e-7 (smoother) to 1.e-5 (sharper)
      Eps = 1.e-6
#  endif
      d0=1./3.
      d1=2./3.
!
      if (ua .ge. 0.) then        ! Take upstream stencils for stability

        IS0 = (q_im1-q_im2)**2
        IS1 = (q_im1-q_i)**2
#  ifdef WENO_Z
        T3 = abs(IS1-IS0)
        a0  = d0*(1+T3/(Eps+IS0))
        a1  = d1*(1+T3/(Eps+IS1))
#  else
        a0  = 1./(3.*(Eps+IS0)**2)
        a1  = 2./(3.*(Eps+IS1)**2)
#  endif
        w0  = a0/(a0+a1)
        w1  = a1/(a0+a1)
        p0  = 1./2.*(3.*q_im1-q_im2)
        p1  = 1./2.*(q_im1+q_i)

        flux3_weno = w0*p0 + w1*p1

      else

        IS0 = (q_i-q_ip1)**2
        IS1 = (q_im1-q_i)**2
#  ifdef WENO_Z
        T3 = abs(IS1-IS0)
        a0  = d0*(1+T3/(Eps+IS0))
        a1  = d1*(1+T3/(Eps+IS1))
#  else
        a0  = 1./(3.*(Eps+IS0)**2)
        a1  = 2./(3.*(Eps+IS1)**2)
#  endif
        w0  = a0/(a0+a1)
        w1  = a1/(a0+a1)
        p0  = 1./2.*(3.*q_i-q_ip1)
        p1  = 1./2.*(q_im1+q_i)

        flux3_weno = w0*p0 + w1*p1

      endif

      return
      end
# endif /* WENO */
!
!======================================================================
!
#else
      subroutine step3d_t_empty
      return
      end
#endif  /* SOLVE3D */

