      subroutine adv_flux (joff, js, je, is, ie, n)

#if defined O_mom
# if defined O_linearized_advection
!=======================================================================
!     Linearized advective tracer flux

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!=======================================================================

      implicit none

      integer istrt, iend, j, js, je, k, i, n, joff, ie, is
      parameter (istrt=2, iend=imt-1)

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "mw.h"

!-----------------------------------------------------------------------
!     limit the longitude indices based on those from the argument list
!     Note: this is currently bypassed. istrt and iend are set as
!           parameters to optimize performance
!-----------------------------------------------------------------------

!     istrt = max(2,is)
!     iend  = min(imt-1,ie)

!-----------------------------------------------------------------------
!     advective flux across eastern and northern face of "T" cells.
!     is zero due to linearization about state of no motion
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!     calculate 2*advective flux across bottom face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      do j=js,je
        do k=1,km-1
          do i=istrt,iend
            adv_fb(i,k,j)  = adv_vbt(i,k,j)*(tbarz(k,n) + tbarz(k+1,n))
          enddo
        enddo
      enddo

# elif defined O_quicker
!=======================================================================
!     3rd order advective tracer flux

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!        n   = tracer
!=======================================================================

      implicit none

      integer istrt, iend, i, k, j, ip, kr, jq, lag, js, je, ip2, n
      integer joff, jrow, jp2, ie, is
      parameter (istrt=2, iend=imt-1)

      real totvel, upos, uneg, vpos, vneg, wpos, wneg

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "grdvar.h"
      include "mw.h"

#  if defined O_gent_mcwilliams
      include "isopyc.h"
#  endif

#  if defined O_ncar_upwind3
      lag = tau
#  else
      lag = taum1
#  endif

!-----------------------------------------------------------------------
!     calculate 2*advective flux across eastern face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      do j=js,je
        do k=1,km
          do i=istrt,iend-1
            ip2 = i + 2
            totvel = adv_vet(i,k,j)
#  if defined O_gent_mcwilliams
     &              + adv_vetiso(i,k,j)
#  endif
            upos = p5*(totvel + abs(totvel))
     &             *tmask(i-1,k,j)*tmask(i,k,j)*tmask(i+1,k,j)
            uneg = p5*(totvel - abs(totvel))
     &             *tmask(ip2,k,j)*tmask(i+1,k,j)*tmask(i,k,j)

            adv_fe(i,k,j) = totvel*(
     &                          quick_x(i,1)*t(i,  k,j,n,tau)
     &                        + quick_x(i,2)*t(i+1,k,j,n,tau))
     &                  - upos*(curv_xp(i,1)*t(i+1,k,j,n,lag)
     &                         +curv_xp(i,2)*t(i  ,k,j,n,lag)
     &                         +curv_xp(i,3)*t(i-1,k,j,n,lag))
     &                  - uneg*(curv_xn(i,1)*t(ip2,k,j,n,lag)
     &                         +curv_xn(i,2)*t(i+1,k,j,n,lag)
     &                         +curv_xn(i,3)*t(i  ,k,j,n,lag))
          enddo
        enddo

        do k=1,km
          i=iend
          ip2 = 3
          totvel = adv_vet(i,k,j)
#  if defined O_gent_mcwilliams
     &              + adv_vetiso(i,k,j)
#  endif
          upos = p5*(totvel + abs(totvel))
     &             *tmask(i-1,k,j)*tmask(i,k,j)*tmask(i+1,k,j)
          uneg = p5*(totvel - abs(totvel))
     &             *tmask(ip2,k,j)*tmask(i+1,k,j)*tmask(i,k,j)

          adv_fe(i,k,j) = totvel*(
     &                          quick_x(i,1)*t(i,  k,j,n,tau)
     &                        + quick_x(i,2)*t(i+1,k,j,n,tau))
     &                  - upos*(curv_xp(i,1)*t(i+1,k,j,n,lag)
     &                         +curv_xp(i,2)*t(i  ,k,j,n,lag)
     &                         +curv_xp(i,3)*t(i-1,k,j,n,lag))
     &                  - uneg*(curv_xn(i,1)*t(ip2,k,j,n,lag)
     &                         +curv_xn(i,2)*t(i+1,k,j,n,lag)
     &                         +curv_xn(i,3)*t(i  ,k,j,n,lag))
        enddo
        call setbcx (adv_fe(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     calculate 2*advective flux across northern face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      if (joff +js .eq. 2) then
        do j=1,1
          do k=1,km
            do i=2,imt-1
              adv_f4n(i,k,j,n) = c0
            enddo
          enddo
        enddo
      endif
      do j=js,je
        jrow = j + joff
        jp2 = min(j+2+joff,jmt) - joff
        do k=1,km
          do i=istrt,iend
            totvel = adv_vnt(i,k,j)
#  if defined O_gent_mcwilliams
     &              + adv_vntiso(i,k,j)
#  endif
            vpos = p5*(totvel + abs(totvel))
     &             *tmask(i,k,j-1)*tmask(i,k,j)*tmask(i,k,j+1)
            vneg = p5*(totvel - abs(totvel))
     &             *tmask(i,k,jp2)*tmask(i,k,j+1)*tmask(i,k,j)

            adv_f4n(i,k,j,n) = totvel*(
     &                          quick_y(jrow,1)*t(i,k,j  ,n,tau)
     &                        + quick_y(jrow,2)*t(i,k,j+1,n,tau))
     &                  - vpos*(curv_yp(jrow,1)*t(i,k,j+1,n,lag)
     &                         +curv_yp(jrow,2)*t(i,k,j  ,n,lag)
     &                         +curv_yp(jrow,3)*t(i,k,j-1,n,lag))
     &                  - vneg*(curv_yn(jrow,1)*t(i,k,jp2,n,lag)
     &                         +curv_yn(jrow,2)*t(i,k,j+1,n,lag)
     &                         +curv_yn(jrow,3)*t(i,k,j  ,n,lag))
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     calculate 2*advective flux across bottom face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      do j=js,je
        do k=2,km-2
          do i=istrt,iend
            totvel = adv_vbt(i,k,j)
#  if defined O_gent_mcwilliams
     &              + adv_vbtiso(i,k,j)
#  endif
            wpos = p5*(totvel + abs(totvel))
     &             *tmask(i,k+2,j)*tmask(i,k+1,j)*tmask(i,k,j)
            wneg = p5*(totvel - abs(totvel))
     &             *tmask(i,k-1,j)*tmask(i,k,j)*tmask(i,k+1,j)

            adv_fb(i,k,j)  = totvel*(
     &                          quick_z(k,1)*t(i,k  ,j,n,tau)
     &                        + quick_z(k,2)*t(i,k+1,j,n,tau))
     &                  - wneg*(curv_zp(k,1)*t(i,k+1,j,n,lag)
     &                         +curv_zp(k,2)*t(i,k  ,j,n,lag)
     &                         +curv_zp(k,3)*t(i,k-1,j,n,lag))
     &                  - wpos*(curv_zn(k,1)*t(i,k+2,j,n,lag)
     &                         +curv_zn(k,2)*t(i,k+1,j,n,lag)
     &                         +curv_zn(k,3)*t(i,k  ,j,n,lag))
          enddo
        enddo
        k=1
        do i=istrt,iend
          totvel = adv_vbt(i,k,j)
#  if defined O_gent_mcwilliams
     &            + adv_vbtiso(i,k,j)
#  endif
          wpos = p5*(totvel + abs(totvel))
     &             *tmask(i,k+2,j)*tmask(i,k+1,j)*tmask(i,k,j)

          adv_fb(i,k,j)  = totvel*(
     &                        quick_z(k,1)*t(i,k  ,j,n,tau)
     &                      + quick_z(k,2)*t(i,k+1,j,n,tau))
     &                - wpos*(curv_zn(k,1)*t(i,k+2,j,n,lag)
     &                       +curv_zn(k,2)*t(i,k+1,j,n,lag)
     &                       +curv_zn(k,3)*t(i,k  ,j,n,lag))
        enddo
        k=km-1
        do i=istrt,iend
          totvel = adv_vbt(i,k,j)
#  if defined O_gent_mcwilliams
     &            + adv_vbtiso(i,k,j)
#  endif
          wneg = p5*(totvel - abs(totvel))
     &             *tmask(i,k-1,j)*tmask(i,k,j)*tmask(i,k+1,j)

          adv_fb(i,k,j)  = totvel*(
     &                        quick_z(k,1)*t(i,k  ,j,n,tau)
     &                      + quick_z(k,2)*t(i,k+1,j,n,tau))
     &                - wneg*(curv_zp(k,1)*t(i,k+1,j,n,lag)
     &                       +curv_zp(k,2)*t(i,k  ,j,n,lag)
     &                       +curv_zp(k,3)*t(i,k-1,j,n,lag))
        enddo
      enddo

# elif defined O_fourth_order_tracer_advection
!=======================================================================
!     4th order advective tracer flux

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!        n   = tracer
!=======================================================================

      implicit none

      integer istrt, iend, j, js, je, k, i, mask, n, ip2, joff, jp2, m
      integerie, is
      parameter (istrt=2, iend=imt-1)

      real a2nd, b2nd, a4th, b4th

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "mw.h"

!-----------------------------------------------------------------------
!     limit the longitude indices based on those from the argument list
!     Note: this is currently bypassed. istrt and iend are set as
!           parameters to optimize performance
!-----------------------------------------------------------------------

!     istrt = max(2,is)
!     iend  = min(imt-1,ie)

!-----------------------------------------------------------------------
!     calculate 2*advective flux across eastern face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      a2nd = 1.0
      b2nd = 0.0
      a4th = 7.0/6.0
      b4th = -1.0/6.0

      do j=js,je
        do k=1,km
          do i=istrt,iend-1
            mask = tmask(i-1,k,j)*tmask(i+2,k,j)
            adv_fe(i,k,j) = adv_vet(i,k,j)*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(t(i,  k,j,n,tau) +
     &                                      t(i+1,k,j,n,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(t(i-1,  k,j,n,tau) +
     &                                      t(i+2,k,j,n,tau)))
          enddo
          i = iend
#  if defined O_cyclic
          ip2 = 3
#  else
          ip2 = imt
#  endif
          mask = tmask(i-1,k,j)*tmask(ip2,k,j)
          adv_fe(i,k,j) = adv_vet(i,k,j)*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(t(i,  k,j,n,tau) +
     &                                      t(i+1,k,j,n,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(t(i-1,  k,j,n,tau) +
     &                                      t(ip2,k,j,n,tau)))
          adv_fe(1,k,j) = adv_fe(imt-1,k,j)
        enddo
      enddo

!-----------------------------------------------------------------------
!     2*advective flux across northern face of "T" cells is built
!     into ADV_Ty. (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      if (joff + js .eq. 2) then
        do j=1,1
          do k=1,km
            do i=2,imt-1
              adv_f4n(i,k,j,n) = 0.0
            enddo
          enddo
        enddo
      endif
      do j=js,je
        jp2 = min(j+2+joff,jmt) - joff
        do k=1,km
          do i=istrt,iend
            mask = tmask(i,k,j-1)*tmask(i,k,jp2)
            adv_f4n(i,k,j,n) = adv_vnt(i,k,j)*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(t(i,  k,j,n,tau) +
     &                                      t(i,k,j+1,n,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(t(i,  k,j-1,n,tau) +
     &                                      t(i,k,jp2,n,tau)))
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     calculate 2*advective flux across bottom face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

        do j=js,je
          do k=2,km-2
            do i=istrt,iend
              mask = tmask(i,k-1,j)*tmask(i,k+2,j)
              adv_fb(i,k,j) = adv_vbt(i,k,j)*(
     &         (a2nd*(1.0-mask) + a4th*mask)*(t(i,  k,j,n,tau) +
     &                                        t(i,k+1,j,n,tau))+
     &         (b2nd*(1.0-mask) + b4th*mask)*(t(i,k-1,j,n,tau) +
     &                                        t(i,k+2,j,n,tau)))
            enddo
          enddo
          k = 1
          m = km-1
          do i=istrt,iend
            adv_fb(i,k,j) = adv_vbt(i,k,j)*(t(i,k  ,j,n,tau) +
     &                                      t(i,k+1,j,n,tau))
            adv_fb(i,m,j) = adv_vbt(i,m,j)*(t(i,m  ,j,n,tau) +
     &                                      t(i,m+1,j,n,tau))
          enddo
        enddo

# elif defined O_fct
!=======================================================================
!     computes advective fluxes using a flux-corrected transport scheme

!        for reference see
!        Gerdes, R., C. Koeberle and J. Willebrandt, 1991
!        the influence of numerical advection schemes on the results of
!        ocean general circulation models. Clim Dynamics 5, 211-226
!        and
!        Zalesak, S. T., 1979: Fully multidimensional flux-corrected
!        transport algorithms for fluids. J. Comp. Phys. 31, 335-362.

!     input:
!       joff  = offset relating "j" in the MW to latitude "jrow"
!       js    = starting row in the MW
!       je    = ending row in the MW
!       jstrt = starting row in the MW for fct
!       jend  = ending row in the MW for fct
!       is    = starting longitude index in the MW
!       ie    = ending longitude index in the MW
!       istrt = max(2,starting longitude index in the MW)
!       iend  = min(imt-1,ending longitude index in the MW)
!       n     = tracer index

!     output: ( via common mw in "mw.h" )
!       adv_fn = 2*advective flux across northern face of T-cell
!       adv_fe = 2*advective flux across eastern face of T-cell
!       adv_fb = 2*advective flux across bottom face of T-cell
!=======================================================================

      implicit none

      integer i, k, j, ip, kr, jq, n, jp, jrow, istrt, is, iend, ie
      integer istrtm1, iendp1, jstrt, js, jend, je, joff, jlast, jp2
      integer jp1

      real t_i, t_j, dz_t2r, dz_tr, dz_wtr, dx_t2r, dx_tr, dy_t2r
      real dy_tr, adv_tx, adv_ty, adv_tz, adv_txiso, adv_tyiso
      real adv_tziso, diff_tx, diff_ty, diff_tz, aidif, totadv
      real fxa, fxb, pplus, pminus, qplus, qminus, reltim

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "accel.h"
      include "coord.h"
      include "grdvar.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "tmngr.h"
#  if defined O_gent_mcwilliams
      include "isopyc.h"
#  endif

!---------------------------------------------------------------------
!     dimension local data
!---------------------------------------------------------------------

      real twodt(km), dcf(imt), Trmin(imt), Trmax(imt)
      real Cpos(imt), Cneg(imt), flxlft(imt), flxrgt(imt)
      real Rpl(imt,km), Rmn(imt,km), tmaski(imt,km,jmw)
#  if defined O_fct_3d
      real tmpext(imt,km,2)
#  endif
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
     &,         t_lo(imt,km)
#  endif
      include "fdift.h"

!-----------------------------------------------------------------------
!     limit the indices based on those from the argument list
!-----------------------------------------------------------------------

      istrt   = max(2,is)
      iend    = min(imt-1,ie)
      istrtm1 = istrt - 1
      iendp1  = iend + 1
      jstrt   = js
      jend    = min(je,jmt-1-joff)

!-----------------------------------------------------------------------
!     initialization when calculating jrow 2
!-----------------------------------------------------------------------

      if (joff + js .eq. 2) then
        jstrt = js - 1
        do k=1,km
          do i=istrt-1,iend
            adv_fn(i,k,1) = c0
            anti_fn(i,k,1,n) = c0
            R_plusY(i,k,1,n) = c0
            R_minusY(i,k,1,n) = c0
#  if defined O_fct_3d
            R_plus3(i,k,1,n) = c0
            R_minus3(i,k,1,n) = c0
#  endif
          enddo
        enddo
      endif

!-----------------------------------------------------------------------
!     create an inverse land mask
!-----------------------------------------------------------------------

      do j=1,jmw
        do k=1,km
          do i=1,imt
            tmaski(i,k,j) = c1 - tmask(i,k,j)
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     calculate 2*advective (low order scheme) flux across northern,
!     eastern and bottom faces of "T" cells
!-----------------------------------------------------------------------

      jlast = min(jend+1+joff,jmt-1) - joff
      do j=js-1,jlast
        do k=1,km
          do i=istrt,iend
            totadv = adv_vnt(i,k,j)
#  if defined O_gent_mcwilliams
     &             + adv_vntiso(i,k,j)
#  endif
            adv_fn(i,k,j) = totadv*
     &                       (t(i,k,j,n,taum1) + t(i,k,j+1,n,taum1))
     &                      + abs(totadv)*
     &                       (t(i,k,j,n,taum1) - t(i,k,j+1,n,taum1))
          enddo
        enddo
      enddo

      jlast = min(jend+1+joff,jmt-1) - joff
      do j=js,jlast
        do k=1,km
          do i=istrtm1,iend
            totadv = adv_vet(i,k,j)
#  if defined O_gent_mcwilliams
     &             + adv_vetiso(i,k,j)
#  endif
            adv_fe(i,k,j) = totadv*
     &                       (t(i,k,j,n,taum1) + t(i+1,k,j,n,taum1))
     &                      + abs(totadv)*
     &                       (t(i,k,j,n,taum1) - t(i+1,k,j,n,taum1))
          enddo
        enddo

        do k=1,kmm1
          do i=istrt,iend
            totadv = adv_vbt(i,k,j)
#  if defined O_gent_mcwilliams
     &             + adv_vbtiso(i,k,j)
#  endif
            adv_fb(i,k,j) = totadv*
     &                      (t(i,k+1,j,n,taum1) + t(i,k,j,n,taum1))
     &                    + abs(totadv)*
     &                      (t(i,k+1,j,n,taum1) - t(i,k,j,n,taum1))
          enddo
        enddo
        do i=istrt,iend
          adv_fb(i,0,j) = adv_vbt(i,0,j)*c2*t(i,1,j,n,taum1)
          adv_fb(i,km,j) = c0
        enddo
      enddo

!-----------------------------------------------------------------------
!     main j loop
!-----------------------------------------------------------------------

      do j=jstrt,jend
        jrow = (j+1) + joff
        jp2  = min(j+2+joff,jmt) - joff
        jp1  = min(j+1+joff,jmt-1) - joff

!-----------------------------------------------------------------------
!       solve for "tau+1" tracer at center of "T" cells in row j+1
!       - low order solution -
!-----------------------------------------------------------------------

        do k=1,km
          twodt(k) = c2dtts*dtxcel(k)
          do i=istrt,iend
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            t_lo(i,k) = (t(i,k,j+1,n,taum1) - twodt(k)
#  else
            t_lo(i,k,j+1,n) = (t(i,k,j+1,n,taum1) - twodt(k)
#  endif
     &        *(ADV_Tx(i,k,jp1) + ADV_Ty(i,k,jp1,jrow,n) +
     &        ADV_Tz(i,k,jp1))*tmask(i,k,j+1))
          enddo
        enddo
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
        call setbcx (t_lo(1,1), imt, km)
#  else
        call setbcx (t_lo(1,1,j+1,n), imt, km)
#  endif

!-----------------------------------------------------------------------
!       next calculate raw antidiffusive fluxes, that is high order
!       scheme flux (leap frog) minus the low order (upstream)
!-----------------------------------------------------------------------

        do k=1,km
          do i=istrtm1,iend
            totadv = adv_vet(i,k,jp1)
#  if defined O_gent_mcwilliams
     &             + adv_vetiso(i,k,jp1)
#  endif
            anti_fe(i,k,j+1,n) = totadv*(t(i,k,j+1,n,tau) +
     &                           t(i+1,k,j+1,n,tau)) - adv_fe(i,k,jp1)
          enddo
          do i=istrt,iend
            totadv = adv_vnt(i,k,jp1)
#  if defined O_gent_mcwilliams
     &             + adv_vntiso(i,k,jp1)
#  endif
            anti_fn(i,k,j+1,n) = totadv*(t(i,k,j+1,n,tau) +
     &                           t(i,k,jp2,n,tau)) - adv_fn(i,k,jp1)
          enddo
        enddo

        do k=1,kmm1
          do i=istrt,iend
            totadv = adv_vbt(i,k,jp1)
#  if defined O_gent_mcwilliams
     &             + adv_vbtiso(i,k,jp1)
#  endif
            anti_fb(i,k,j+1,n) = totadv*(t(i,k,j+1,n,tau) +
     &                           t(i,k+1,j+1,n,tau)) - adv_fb(i,k,jp1)
     &                           *tmask(i,k,j+1)
          enddo
        enddo
        do i=istrt,iend
          anti_fb(i,0,j+1,n) = adv_vbt(i,0,j+1)*c2*t(i,1,j+1,n,taum1)
          anti_fb(i,km,j+1,n) = c0
        enddo

!-----------------------------------------------------------------------
!       now calculate and apply one-dimensional delimiters to these
!       raw antidiffusive fluxes

!       1) calculate T*, that are all halfway neighbors of T
!       2) calculate ratio R+- of Q+- to P+-, that is maximal/minimal
!          possible change of T if no limit would be active,
!          must be at least 1
!       3) choose correct ratio depending on direction of flow as a
!          delimiter
!       4) apply this delimiter to raw antidiffusive flux
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!       delimit x-direction
!-----------------------------------------------------------------------

        do k=1,km

!         prepare some data for use in statement function
#  if defined O_fct_dlm1 || !defined O_fct_dlm2

!         running mean of two adjacent points

          do i=istrt,iendp1
            Trmax(i) = p5*(t(i-1,k,j+1,n,tau) + t(i,k,j+1,n,tau))
          enddo
#  endif

!         extremum of low order solution central point and adjacent
!         halfway neighbours; check for land

          do i=istrt,iend
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            fxa = tmask(i-1,k,j+1)*Trmax(i) +
     &            tmaski(i-1,k,j+1)*t_lo(i,k)
            fxb = tmask(i+1,k,j+1)*Trmax(i+1) +
     &            tmaski(i+1,k,j+1)*t_lo(i,k)
            Trmax(i) = max(fxa,fxb,t_lo(i,k))
            Trmin(i) = min(fxa,fxb,t_lo(i,k))
#  else
            fxa = tmask(i-1,k,j+1)*t_lo(i-1,k,j+1,n) +
     &            tmaski(i-1,k,j+1)*t_lo(i,k,j+1,n)
            fxb = tmask(i+1,k,j+1)*t_lo(i+1,k,j+1,n) +
     &            tmaski(i+1,k,j+1)*t_lo(i,k,j+1,n)
            Trmax(i) = max(fxa,fxb,t_lo(i,k,j+1,n))
            Trmin(i) = min(fxa,fxb,t_lo(i,k,j+1,n))
#  endif
#  if defined O_fct_3d
            tmpext(i,k,1) = Trmax(i)
            tmpext(i,k,2) = Trmin(i)
#  endif
            dcf(i) = cstdxt2r(i,j+1)
            flxlft(i) = anti_fe(i-1,k,j+1,n)
            flxrgt(i) = anti_fe(i,k,j+1,n)
          enddo

!         calculate ratio R

          do i=istrt,iend
            Pplus  = c2dtts*dcf(i)*(max(c0,flxlft(i))-min(c0,flxrgt(i)))
            Pminus = c2dtts*dcf(i)*(max(c0,flxrgt(i))-min(c0,flxlft(i)))
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            Qplus  = Trmax(i) - t_lo(i,k)
            Qminus = t_lo(i,k) - Trmin(i)
#  else
            Qplus  = Trmax(i) - t_lo(i,k,j+1,n)
            Qminus = t_lo(i,k,j+1,n) - Trmin(i)
#  endif
            Rpl(i,k) = min(1.,tmask(i,k,j+1)*Qplus/(Pplus+epsln))
            Rmn(i,k) = min(1.,tmask(i,k,j+1)*Qminus/(Pminus+epsln))
          enddo
          call setbcx (Rpl, imt, km)
          call setbcx (Rmn, imt, km)

!         calculate delimiter using ratio at adjacent points

          do i=istrt,iendp1
            Cpos(i-1) = min(Rpl(i,k),Rmn(i-1,k))
            Cneg(i-1) = min(Rpl(i-1,k),Rmn(i,k))
          enddo

!         finally apply appropriate delimiter to flux

          do i=istrtm1,iend
            anti_fe(i,k,j+1,n) = p5*((Cpos(i) + Cneg(i))
     &                               *anti_fe(i,k,j+1,n) +
     &                               (Cpos(i) - Cneg(i))
     &                               *abs(anti_fe(i,k,j+1,n)))
          enddo
        enddo

!-----------------------------------------------------------------------
!       delimit y-direction
!-----------------------------------------------------------------------

        do k=1,km

!         prepare some data for use in statement function

          do i=istrt,iend
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            fxa = p5*tmask(i,k,j)*(t(i,k,j,n,tau) +
     &            t(i,k,j+1,n,tau)) +
     &            tmaski(i,k,j)*t_lo(i,k)
            fxb = p5*tmask(i,k,jp2)*(t(i,k,j+1,n,tau) +
     &            t(i,k,jp2,n,tau)) +
     &            tmaski(i,k,jp2)*t_lo(i,k)
            Trmax(i) = max(fxa,fxb,t_lo(i,k))
            Trmin(i) = min(fxa,fxb,t_lo(i,k))
#  else
            fxa = tmask(i,k,j)*t_lo(i,k,j,n) +
     &            tmaski(i,k,j)*t_lo(i,k,j+1,n)
            fxb = p5*tmask(i,k,jp2)*(t(i,k,j+1,n,tau) +
     &            t(i,k,jp2,n,tau)) +
     &            tmaski(i,k,jp2)*t_lo(i,k,j+1,n)
            Trmax(i) = max(fxa,fxb,t_lo(i,k,j+1,n))
            Trmin(i) = min(fxa,fxb,t_lo(i,k,j+1,n))
#  endif
#  if defined O_fct_3d
            tmpext(i,k,1) = max(Trmax(i),tmpext(i,k,1))
            tmpext(i,k,2) = min(Trmin(i),tmpext(i,k,2))
#  endif
            dcf(i) = cstdyt2r(jrow)
            flxlft(i) = anti_fn(i,k,j,n)
            flxrgt(i) = anti_fn(i,k,j+1,n)
          enddo

!         calculate ratio R, related to a point

          do i=istrt,iend
            Pplus  = c2dtts*dcf(i)*(max(c0,flxlft(i))-min(c0,flxrgt(i)))
            Pminus = c2dtts*dcf(i)*(max(c0,flxrgt(i))-min(c0,flxlft(i)))
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            Qplus  = Trmax(i) - t_lo(i,k)
            Qminus = t_lo(i,k) - Trmin(i)
#  else
            Qplus  = Trmax(i) - t_lo(i,k,j+1,n)
            Qminus = t_lo(i,k,j+1,n) - Trmin(i)
#  endif
            R_plusY(i,k,j+1,n) =
     &        min(1.,tmask(i,k,j+1)*Qplus/(Pplus+epsln))
            R_minusY(i,k,j+1,n) =
     &        min(1.,tmask(i,k,j+1)*Qminus/(Pminus+epsln))
          enddo

!         calculate delimiter using ratio at adjacent points

          do i=istrt,iend
            Cpos(i) = min(R_plusY(i,k,j+1,n),R_minusY(i,k,j,n))
            Cneg(i) = min(R_plusY(i,k,j,n),R_minusY(i,k,j+1,n))
          enddo

!         finally get delimiter c dependent on direction of flux and
!         apply it to raw antidiffusive flux

          do i=istrt,iend
            anti_fn(i,k,j,n) = p5*((Cpos(i) + Cneg(i))
     &                             *anti_fn(i,k,j,n) +
     &                             (Cpos(i) - Cneg(i))
     &                             *abs(anti_fn(i,k,j,n)))
          enddo
        enddo

!-----------------------------------------------------------------------
!       delimit z-direction
!-----------------------------------------------------------------------

        do k=1,km

!         prepare some data for use in statement function

          do i=istrt,iend
            dcf(i) = dzt2r(k)
            flxlft(i) = anti_fb(i,k,j+1,n)
            flxrgt(i) = anti_fb(i,k-1,j+1,n)
            if (k .gt. 1)then
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
              fxa = p5*tmask(i,k-1,j+1)*
     &              (t(i,k-1,j+1,n,tau) + t(i,k,j+1,n,tau)) +
     &              tmaski(i,k-1,j+1)*t_lo(i,k)
#  else
              fxa = tmask(i,k-1,j+1)*t_lo(i,k-1,j+1,n) +
     &              tmaski(i,k-1,j+1)*t_lo(i,k,j+1,n)
#  endif
            else
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
              fxa = t_lo(i,k)
#  else
              fxa = t_lo(i,k,j+1,n)
#  endif
            endif
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            if (k .lt. km) then
              fxb = p5*tmask(i,k+1,j+1)*
     &              (t(i,k,j+1,n,tau)+t(i,k+1,j+1,n,tau)) +
     &              tmaski(i,k+1,j+1)*t_lo(i,k)
#  else
            if (k .lt. km) then
              fxb = tmask(i,k+1,j+1)*t_lo(i,k+1,j+1,n) +
     &              tmaski(i,k+1,j+1)*t_lo(i,k,j+1,n)
#  endif
            else
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
              fxb = t_lo(i,k)
#  else
              fxb = t_lo(i,k,j+1,n)
#  endif
            endif
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            Trmax(i) = max(fxa,fxb,t_lo(i,k))
            Trmin(i) = min(fxa,fxb,t_lo(i,k))
#  else
            Trmax(i) = max(fxa,fxb,t_lo(i,k,j+1,n))
            Trmin(i) = min(fxa,fxb,t_lo(i,k,j+1,n))
#  endif
#  if defined O_fct_3d
            tmpext(i,k,1) = max(Trmax(i),tmpext(i,k,1))
            tmpext(i,k,2) = min(Trmin(i),tmpext(i,k,2))
#  endif
          enddo

!         calculate delimiter using ratio at adjacent points
!         this variable is related to an arc (between two points,
!         the same way as fluxes are defined.)

          do i=istrt,iend
            Pplus  = c2dtts*dcf(i)*(max(c0,flxlft(i))-min(c0,flxrgt(i)))
            Pminus = c2dtts*dcf(i)*(max(c0,flxrgt(i))-min(c0,flxlft(i)))
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            Qplus = Trmax(i) - t_lo(i,k)
            Qminus = t_lo(i,k) - Trmin(i)
#  else
            Qplus = Trmax(i) - t_lo(i,k,j+1,n)
            Qminus = t_lo(i,k,j+1,n) - Trmin(i)
#  endif
            Rpl(i,k) = min(1.,tmask(i,k,j+1)*Qplus/(Pplus+epsln))
            Rmn(i,k) = min(1.,tmask(i,k,j+1)*Qminus/(Pminus+epsln))
          enddo

        enddo

        do k=1,kmm1

!         calculate delimiter using ratio at adjacent points
!         this variable is related to an arc (between two points,
!         the same way as fluxes are defined.)

          do i=istrt,iend
            Cneg(i) = min(Rpl(i,k+1),Rmn(i,k))
            Cpos(i) = min(Rpl(i,k),Rmn(i,k+1))
          enddo

!         finally get delimiter c dependent on direction of flux and
!         apply it to raw antidiffusive flux

          do i=istrt,iend
            anti_fb(i,k,j+1,n) = p5*((Cpos(i)+Cneg(i))
     &                               *anti_fb(i,k,j+1,n) +
     &                               (Cpos(i)-Cneg(i))
     &                               *abs(anti_fb(i,k,j+1,n)))
          enddo
        enddo
        do i=istrt,iend
          anti_fb(i,0,j+1,n)  = c0
          anti_fb(i,km,j+1,n) = c0
        enddo
#  if defined O_fct_3d

!-----------------------------------------------------------------------
!     then calculate and apply 3-d delimiter to just delimited
!     antidiffusive fluxes
!-----------------------------------------------------------------------

        do k=1,km

!         prepare some data for use in statement function

          do i=istrt,iend
            Trmax(i) = tmpext(i,k,1)
            Trmin(i) = tmpext(i,k,2)
          enddo

          do i=istrt,iend
#  if defined O_fct_dlm1 || !defined O_fct_dlm2
            Qplus  = Trmax(i) - t_lo(i,k)
            Qminus = t_lo(i,k) - Trmin(i)
#  else
            Qplus  = Trmax(i) - t_lo(i,k,j+1,n)
            Qminus = t_lo(i,k,j+1,n) - Trmin(i)
#  endif
            R_plus3(i,k,j+1,n) = min(1.,tmask(i,k,j+1)*Qplus/
     &                           (epsln+c2dtts*(
     &                           cstdxt2r(i,j+1)
     &                             *(max(c0,anti_fe(i-1,k,j+1,n)) -
     &                               min(c0,anti_fe(i,k,j+1,n))) +
     &                           cstdyt2r(jrow)
     &                             *(max(c0,anti_fn(i,k,j,n)) -
     &                               min(c0,anti_fn(i,k,j+1,n))) +
     &                           dzt2r(k)
     &                             *(max(c0,anti_fb(i,k,j+1,n)) -
     &                               min(c0,anti_fb(i,k-1,j+1,n)))
     &                           )))

            R_minus3(i,k,j+1,n) = min(1.,tmask(i,k,j+1)*Qminus/
     &                            (epsln+c2dtts*(
     &                            cstdxt2r(i,j+1)
     &                              *(max(c0,anti_fe(i,k,j+1,n)) -
     &                                min(c0,anti_fe(i-1,k,j+1,n))) +
     &                            cstdyt2r(jrow)
     &                              *(max(c0,anti_fn(i,k,j+1,n)) -
     &                                min(c0,anti_fn(i,k,j,n))) +
     &                            dzt2r(k)
     &                              *(max(c0,anti_fb(i,k-1,j+1,n)) -
     &                                min(c0,anti_fb(i,k,j+1,n)))
     &                            )))

          enddo
        enddo
        call setbcx (R_plus3(1,1,j+1,n), imt, km)
        call setbcx (R_minus3(1,1,j+1,n), imt, km)

!       finally apply 3-d delimiters to precorrected fluxes

        do k=1,km
          do i=istrt,iendp1
            Cpos(i-1) = min(R_plus3(i,k,j+1,n),R_minus3(i-1,k,j+1,n))
            Cneg(i-1) = min(R_plus3(i-1,k,j+1,n),R_minus3(i,k,j+1,n))
          enddo
          do i=istrtm1,iend
            anti_fe(i,k,j+1,n) = p5*((Cpos(i) + Cneg(i))
     &                               *anti_fe(i,k,j+1,n) +
     &                               (Cpos(i) - Cneg(i))
     &                               *abs(anti_fe(i,k,j+1,n)))
          enddo

          do i=istrt,iend
            Cpos(i) = min(R_plus3(i,k,j+1,n),R_minus3(i,k,j,n))
            Cneg(i) = min(R_plus3(i,k,j,n),R_minus3(i,k,j+1,n))
          enddo
          do i=istrt,iend
            anti_fn(i,k,j,n) = p5*((Cpos(i) + Cneg(i))
     &                             *anti_fn(i,k,j,n) +
     &                             (Cpos(i) - Cneg(i))
     &                             *abs(anti_fn(i,k,j,n)))
          enddo
        enddo

        do k=1,kmm1
          do i=istrt,iend
            Cneg(i) = min(R_plus3(i,k+1,j+1,n),R_minus3(i,k,j+1,n))
            Cpos(i) = min(R_plus3(i,k,j+1,n),R_minus3(i,k+1,j+1,n))
          enddo
          do i=istrt,iend
            anti_fb(i,k,j+1,n) = p5*((Cpos(i) + Cneg(i))
     &                               *anti_fb(i,k,j+1,n) +
     &                               (Cpos(i) - Cneg(i))
     &                               *abs(anti_fb(i,k,j+1,n)))
          enddo
        enddo

#  endif

!-----------------------------------------------------------------------
!       complete advective fluxes by adding low order fluxes to
!       delimited antidiffusive fluxes
!-----------------------------------------------------------------------

        do k=1,km
          do i=istrtm1,iend
            anti_fe(i,k,j+1,n) = anti_fe(i,k,j+1,n) + adv_fe(i,k,jp1)
          enddo
          do i=istrt,iend
            anti_fn(i,k,j,n) = (anti_fn(i,k,j,n) + adv_fn(i,k,j))
     &                         *tmask(i,k,j)
            anti_fb(i,k,j+1,n) = (anti_fb(i,k,j+1,n) + adv_fb(i,k,jp1))
     &                           *tmask(i,k,j+1)
          enddo
        enddo

      enddo

!-----------------------------------------------------------------------
!     set 2*corrected advective fluxes across northern, eastern and
!     bottom faces of "T" cells
!-----------------------------------------------------------------------

      do j=js-1,jend
        do k=1,km
          do i=istrt,iend
            adv_fn(i,k,j) = anti_fn(i,k,j,n)
          enddo
        enddo
      enddo

      do j=js,jend
        do k=1,km
          do i=istrtm1,iend
            adv_fe(i,k,j) = anti_fe(i,k,j,n)
          enddo
        enddo

        do k=1,kmm1
          do i=istrt,iend
            adv_fb(i,k,j) = anti_fb(i,k,j,n)
          enddo
        enddo
      enddo

# else
!=======================================================================
!     2nd order advective tracer flux

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!        n   = tracer
!=======================================================================

      implicit none

      integer istrt, iend, j, js, je, k, i, n, joff, ie, is, ip, kr, jq

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "mw.h"

      parameter (istrt=2, iend=imt-1)

!-----------------------------------------------------------------------
!     calculate 2*advective flux across eastern face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      do j=js,je
        do k=1,km
          do i=istrt-1,iend
            adv_fe(i,k,j) = adv_vet(i,k,j)*(t(i,  k,j,n,tau) +
     &                                      t(i+1,k,j,n,tau))
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     calculate 2*advective flux across bottom face of "T" cells.
!     (It`s done this way for performance issues)
!-----------------------------------------------------------------------

      do j=js,je
        do k=1,km-1
          do i=istrt,iend
            adv_fb(i,k,j)  = adv_vbt(i,k,j)*(t(i,k,  j,n,tau) +
     &                                       t(i,k+1,j,n,tau))
          enddo
        enddo
      enddo

# endif
#endif
      return
      end
