#include "options.inc"




      subroutine adv_vel()
      use spflame_module
      use freesurf_module
      implicit none
      integer is,ie,js,je,i,j,k
      real dys,dyr,dyn,asw,ase,anw,ane
c
c=======================================================================
c     calculate advection velocities for momentum and tracer equations     
c     SPFLAME version:   c.eden
c=======================================================================
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     advection velocity on northern face of "T" cells. Note the
c     imbedded cosine.
c     adv_vnt = WT_AVG_X(u(1,1,1,2,tau))
c-----------------------------------------------------------------------
c
      do j=js-1,je+1
        do k=1,km
          do i=is,ie+1
#ifdef partial_cell
            adv_vnt(i,k,j) = (u(i  ,k,j,2,tau)*dxu(i  )*dhu(i  ,k,j) +
     &                        u(i-1,k,j,2,tau)*dxu(i-1)*dhu(i-1,k,j)
     &                       )*csu(j)*dxt2r(i)
#else
            adv_vnt(i,k,j) = (u(i,k,j,2,tau)*dxu(i) +   
     &                     u(i-1,k,j,2,tau)*dxu(i-1))*csu(j)*dxt2r(i)
#endif
          enddo
        enddo
      enddo

      call set_cyclic(adv_vnt,km,1)
      if (my_blk_i == 1 .and. enable_obc_west) 
     &              adv_vnt(1,:,:) = adv_vnt(2,:,:)
      if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &              adv_vnt(imt,:,:) = adv_vnt(imt-1,:,:)
c
c-----------------------------------------------------------------------
c     advection velocity on the eastern face of "T" cells
c     adv_vnt = WT_AVG_Y(u(1,1,1,1,tau))
c-----------------------------------------------------------------------
c
      do j=js,je+1
        do k=1,km
          do i=is-1,ie+1
#ifdef partial_cell
            adv_vet(i,k,j) = (u(i,k,j  ,1,tau)*dyu(j  )*dhu(i,k,j  )
     &                      + u(i,k,j-1,1,tau)*dyu(j-1)*dhu(i,k,j-1)
     &                       )*dyt2r(j)
#else
            adv_vet(i,k,j) = (u(i,k,j,1,tau)*dyu(j) + 
     &                     u(i,k,j-1,1,tau)*dyu(j-1))*dyt2r(j)
#endif
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     construct vertical velocity on the bottom face of "T" cells
c-----------------------------------------------------------------------
c
      if (enable_freesurf) then

       call calc_convU

       do j=js,je+1
        do i=is,ie+1
          adv_vbt(i,0,j)   = convU(i,j)
        enddo     
       enddo     

      else

       do j=js,je+1
        do i=is,ie+1
          adv_vbt(i,0,j)   = 0.
        enddo     
       enddo     

      endif

      do j=js,je+1
        do k=1,km
          do i=is,ie+1
            adv_vbt(i,k,j) =
     &                   ((adv_vet(i,k,j) - adv_vet(i-1,k,j))*dxtr(i)
     &                   +(adv_vnt(i,k,j) - adv_vnt(i,k,j-1))*dytr(j)
#ifdef partial_cell
     &                   )*cstr(j)
#else
     &                   )*cstr(j)*dzt(k)
#endif
          enddo
        enddo
        do k=1,km
          do i=is,ie+1
            adv_vbt(i,k,j) = adv_vbt(i,k,j) + adv_vbt(i,k-1,j)
          enddo
        enddo
      end do

      call set_cyclic(adv_vbt,km+1,1)
      if (my_blk_i == 1 .and. enable_obc_west)        
     &                 adv_vbt(1,:,:) = adv_vbt(2,:,:)
      if (my_blk_i == n_pes_i .and. enable_obc_east)  
     &                 adv_vbt(imt,:,:) = adv_vbt(imt-1,:,:)
c
c-----------------------------------------------------------------------
c     construct advection velocity on the northern face of "u" cells by
c     averaging advection velocity on northern face of "t" cells 
c     note: je-1 is used instead of jemw to account for possible non
c           integral number of MW`s in jmt 
c     adv_vnu = LINEAR_INTRP_Y(WT_AVG_X(adv_vnt))
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
	    adv_vnu(i,k,j) = ((adv_vnt(i,k,j)*duw(i) 
     &                       + adv_vnt(i+1,k,j)*due(i)
     &                        )*dus(j+1) +
     &                        (adv_vnt(i,k,j+1)*duw(i) 
     &                       + adv_vnt(i+1,k,j+1)*due(i)
     &                        )*dun(j))*dytr(j+1)*dxur(i)
          enddo
        enddo
      enddo

      call set_cyclic(adv_vnu,km,1)
      if (my_blk_i == 1 .and. enable_obc_west)    
     &                   adv_vnu(1,:,:) = adv_vnu(2,:,:)
      if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &                   adv_vnu(imt,:,:) = adv_vnu(imt-1,:,:)
c
c-----------------------------------------------------------------------
c     construct advection velocity on the eastern face of "u" cells by
c     averaging advection velocity on eastern face of "t" cells 
c     note: take special care of zonal b.c. on this term.
c     adv_veu = LINEAR_INTRP_X(WT_AVG_Y(adv_vet))
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
	    adv_veu(i,k,j) = ((adv_vet(i,k,j)*dus(j)
     &                       + adv_vet(i,k,j+1)*dun(j)
     &                        )*duw(i+1) +    
     &                        (adv_vet(i+1,k,j)*dus(j)
     &                       + adv_vet(i+1,k,j+1)*dun(j)
     &                        )*due(i))*dyur(j)*dxtr(i+1)
          enddo
        enddo
      enddo

      call set_cyclic(adv_veu,km,1)
      if (my_blk_i == 1 .and. enable_obc_west) 
     &           adv_veu(1,:,:) = adv_veu(2,:,:)
      if (my_blk_i == n_pes_i .and. enable_obc_east) then
          adv_veu(imt,:,:) = adv_veu(imt-1,:,:)
      elseif (my_blk_i == n_pes_i) then
          adv_veu(imt,:,:) = 0.
      endif
c
c-----------------------------------------------------------------------
c     construct advection velocity on the bottom face of "u" cells by
c     averaging advection velocity on bottom face of "t" cells 
c-----------------------------------------------------------------------
c
c
#ifdef partial_cell
c
c      this is the way MOM 3 calculates adv_cbu
c      use it for partial cell only, it should be however the same, 
c      maybe faster
c
      do j=js,je
c
c       construct vertical velocity at top of first level
c
	dyn  = dun(j)*cst(j+1)
	dys  = dus(j)*cst(j)
	dyr  = dyur(j)*csur(j)
        k=0
        do i=is,ie
	  asw = duw(i)*dys
	  anw = duw(i)*dyn
	  ase = due(i)*dys
	  ane = due(i)*dyn
	  adv_vbu(i,k,j) = dyr*dxur(i)*(
     &                    adv_vbt(i,k,j)*asw + adv_vbt(i+1,k,j)*ase
     &                  + adv_vbt(i,k,j+1)*anw + adv_vbt(i+1,k,j+1)*ane)
        enddo
c
c       construct divergence of advection velocity * level thickness
c
        do k=1,km
          do i=is,ie
            adv_vbu(i,k,j) =
     &                   ((adv_veu(i,k,j) - adv_veu(i-1,k,j))*dxur(i)
     &                   +(adv_vnu(i,k,j) - adv_vnu(i,k,j-1))*dyur(j)
#ifdef partial_cell
     &                   )*csur(j)
#else
     &                   )*csur(j)*dzt(k)
#endif
          enddo
        enddo
c
c       integrate downward to define "adv_vbu" at the bottom of levels
c
        do k=1,km
          do i=is,ie
            adv_vbu(i,k,j) = adv_vbu(i,k,j) + adv_vbu(i,k-1,j)
          enddo
        enddo
      enddo
c
#else
c
c      this is the way MOM 2 calculates adv_cbu
c

      do j=js,je
        do k=0,km
          do i=is,ie
	    adv_vbu(i,k,j) = dyur(j)*csur(j)*dxur(i)*(
     &                    adv_vbt(i,k,j)*duw(i)*dus(j)*cst(j)
     &                  + adv_vbt(i+1,k,j)*due(i)*dus(j)*cst(j)
     &                  + adv_vbt(i,k,j+1)*duw(i)*dun(j)*cst(j+1)
     &                  + adv_vbt(i+1,k,j+1)*due(i)*dun(j)*cst(j+1) )
          enddo
        enddo
      enddo
#endif

      call set_cyclic(adv_vbu,km+1,1)
      if (my_blk_i == 1 .and. enable_obc_west) 
     &         adv_vbu(1,:,:) = adv_vbu(2,:,:)
      if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &         adv_vbu(imt,:,:) = adv_vbu(imt-1,:,:)


      if (enable_bbl) call adv_vel_bbl

      end subroutine adv_vel
      

c
c=======================================================================
c     Module Quicker Advection for SPFLAME
c     taken from MOM 3, originally written by
c       r.c.pacanowski       e-mail rcp@gfdl.gov
c
c     it is linked in the main code in setup
c     and tracer (see also topog for array tmask)
c
c       c.eden
c=======================================================================
c
      module quicker_module
      implicit none

      real, allocatable :: quick_x(:,:)
      real, allocatable :: curv_xp(:,:),curv_xn(:,:)

      real, allocatable :: quick_y(:,:)
      real, allocatable :: curv_yp(:,:),curv_yn(:,:)

      real, allocatable :: quick_z(:,:)
      real, allocatable :: curv_zp(:,:),curv_zn(:,:)

      end module



      subroutine init_quicker
c=======================================================================
c     Initialisation of Quicker Advection module
c=======================================================================
      use spflame_module
      use quicker_module
      implicit none
      integer i,j,k,ip2,jp1,jp2,jm1,kp1,kp2,km1

      if (my_pe==0) then
       print*,' Initialisation of quicker advection scheme'
       print*,''
       print*,' using Quicker advection for all tracers'
       print*,' instead of centered differences'
       print*,''
      endif

      if (enable_4th_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled 4th order scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_quiekcer')
      endif

      if (enable_fct_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled FCT scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_quiekcer')
      endif

      if (enable_upstream_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled upstream scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_quiekcer')
      endif

      allocate( quick_x(imt,2) ); quick_x=0.
      allocate( curv_xp(imt,3) ); curv_xp=0.
      allocate( curv_xn(imt,3) ); curv_xn=0.

      do i=2,imt-1
        ip2=min(i+2,imt)
        if (i == imt-1 .and. cyclic) ip2=3
        quick_x(i,1) = 2.*dxt(i+1)/(dxt(i+1)+dxt(i))
        quick_x(i,2) = 2.*dxt(i  )/(dxt(i+1)+dxt(i))
        curv_xp(i,1) = 2.*dxt(i)*dxt(i+1)/
     &               ((dxt(i-1)+2.0*dxt(i)+dxt(i+1))*(dxt(i)+dxt(i+1)))
        curv_xp(i,2) =-2.*dxt(i)*dxt(i+1)/((dxt(i)+dxt(i+1))
     &                                  *(dxt(i-1)+dxt(i)))
        curv_xp(i,3) = 2.*dxt(i)*dxt(i+1)/
     &               ((dxt(i-1)+2.0*dxt(i)+dxt(i+1))*(dxt(i-1)+dxt(i)))
        curv_xn(i,1) = 2.*dxt(i)*dxt(i+1)/
     &             ((dxt(i)+2.0*dxt(i+1)+dxt(ip2))*(dxt(i+1)+dxt(ip2)))
        curv_xn(i,2) =-2.*dxt(i)*dxt(i+1)/((dxt(i+1)+dxt(ip2))
     &                                  *(dxt(i)+dxt(i+1)))
        curv_xn(i,3) = 2.*dxt(i)*dxt(i+1)/
     &               ((dxt(i)+2.0*dxt(i+1)+dxt(ip2))*(dxt(i)+dxt(i+1)))
      enddo

      if (cyclic) then
        quick_x(1,:)   = quick_x(imt-1,:)
	quick_x(imt,:) = quick_x(2,:)
        curv_xp(1,:)   = curv_xp(imt-1,:)
	curv_xp(imt,:) = curv_xp(2,:)
        curv_xn(1,:)   = curv_xn(imt-1,:)
	curv_xn(imt,:) = curv_xn(2,:)
      endif

      allocate( quick_y(jmt,2) ); quick_y=0.
      allocate( curv_yp(jmt,3) ); curv_yp=0.
      allocate( curv_yn(jmt,3) ); curv_yn=0.

      do j=1,jmt
        jp2 = min(j+2,jmt); jp1 = min(j+1,jmt); jm1 = max(j-1,1)
        quick_y(j,1) = 2.*dyt(jp1)/(dyt(jp1)+dyt(j))
        quick_y(j,2) = 2.*dyt(j  )/(dyt(jp1)+dyt(j))
        curv_yp(j,1) = 2.*dyt(j)*dyt(jp1)/
     &               ((dyt(jm1)+2.0*dyt(j)+dyt(jp1))*(dyt(j)+dyt(jp1)))
        curv_yp(j,2) =-2.*dyt(j)*dyt(jp1)/((dyt(j)+dyt(jp1))
     &                                  *(dyt(jm1)+dyt(j)))
        curv_yp(j,3) = 2.*dyt(j)*dyt(jp1)/
     &               ((dyt(jm1)+2.0*dyt(j)+dyt(jp1))*(dyt(jm1)+dyt(j)))
        curv_yn(j,1) = 2.*dyt(j)*dyt(jp1)/
     &             ((dyt(j)+2.0*dyt(jp1)+dyt(jp2))*(dyt(jp1)+dyt(jp2)))
        curv_yn(j,2) =-2.*dyt(j)*dyt(jp1)/((dyt(jp1)+dyt(jp2))
     &                                  *(dyt(j)+dyt(jp1)))
        curv_yn(j,3) = 2.*dyt(j)*dyt(jp1)/
     &               ((dyt(j)+2.0*dyt(jp1)+dyt(jp2))*(dyt(j)+dyt(jp1)))
      enddo
c
      allocate( quick_z(km,2) ); quick_z=0.
      allocate( curv_zp(km,3) ); curv_zp=0.
      allocate( curv_zn(km,3) ); curv_zn=0.

      do k=1,km
        kp2 = min(k+2,km); kp1 = min(k+1,km); km1 = max(k-1,1)
        quick_z(k,1) = 2.*dzt(kp1)/(dzt(kp1)+dzt(k))
        quick_z(k,2) = 2.*dzt(k  )/(dzt(kp1)+dzt(k))
        curv_zp(k,1) = 2.*dzt(k)*dzt(kp1)/
     &               ((dzt(km1)+2.0*dzt(k)+dzt(kp1))*(dzt(k)+dzt(kp1)))
        curv_zp(k,2) =-2.*dzt(k)*dzt(kp1)/((dzt(k)+dzt(kp1))
     &                                  *(dzt(km1)+dzt(k)))
        curv_zp(k,3) = 2.*dzt(k)*dzt(kp1)/
     &               ((dzt(km1)+2.0*dzt(k)+dzt(kp1))*(dzt(km1)+dzt(k)))
c
        curv_zn(k,1) = 2.*dzt(k)*dzt(kp1)/
     &             ((dzt(k)+2.0*dzt(kp1)+dzt(kp2))*(dzt(kp1)+dzt(kp2)))
        curv_zn(k,2) =-2.*dzt(k)*dzt(kp1)/((dzt(kp1)+dzt(kp2))
     &                                  *(dzt(k)+dzt(kp1)))
        curv_zn(k,3) = 2.*dzt(k)*dzt(kp1)/
     &               ((dzt(k)+2.0*dzt(kp1)+dzt(kp2))*(dzt(k)+dzt(kp1)))
      enddo
      if (my_pe==0) print*,' done'
      end subroutine init_quicker


      subroutine adv_east_flux_quicker(n,adv_fe)
c=======================================================================
c     3rd order advective tracer flux
c=======================================================================
      use spflame_module
      use quicker_module
      implicit none
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      integer is,ie,js,je,n, i,j,k,ip2
      real totvel,upos,uneg,eastmsk,westmsk

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across eastern face of T cells.  (It`s
c     done this way for performance issues).  In grid points near land
c     points and near bottom and surface, project the curvature coeff of
c     the land point to the center point.
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
            ip2=min(i+2,imt)
            if (i == imt-1 .and. cyclic) ip2=i+2
            totvel = adv_vet(i,k,j)
            upos = .5*(totvel + abs(totvel))
     &                             *tmask(i,k,j)*tmask(i+1,k,j)
            uneg = .5*(totvel - abs(totvel))
     &                             *tmask(i+1,k,j)*tmask(i,k,j)
c           eastmsk is one only in grid points just east of land:
            eastmsk=tmask(i,k,j)*(1-tmask(i-1,k,j))
c           westmsk is one only in grid points just west of land:
            westmsk=tmask(i+1,k,j)*(1-tmask(ip2,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,taum1)
     &                         +curv_xp(i,2)*t(i  ,k,j,n,taum1)
     &                         +curv_xp(i,3)
     &        *(t(i-1,k,j,n,taum1)*(1-eastmsk)
     &         +t(i,k,j,n,taum1)*eastmsk))
     &                  - uneg*(curv_xn(i,1)
     &        *(t(ip2,k,j,n,taum1)*(1-westmsk)+
     &          t(i+1,k,j,n,taum1)*westmsk)
     &                         +curv_xn(i,2)*t(i+1,k,j,n,taum1)
     &                         +curv_xn(i,3)*t(i  ,k,j,n,taum1))
          enddo
        enddo
      enddo
      end subroutine adv_east_flux_quicker




      subroutine adv_north_flux_quicker(n,adv_fn)
c=======================================================================
c     3rd order advective tracer flux
c=======================================================================
      use spflame_module
      use quicker_module
      implicit none
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      integer is,ie,js,je,n,i,j,k,jp2
      real totvel,vpos,vneg,rnormsk,soutmsk

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across northern face of T cells.
c     (It`s done this way for performance issues)
c-----------------------------------------------------------------------
c
      do j=js-1,je
        jp2 = min(jmt,j+2)
        do k=1,km
          do i=is,ie
            totvel = adv_vnt(i,k,j)
            vpos = .5*(totvel + abs(totvel))
     &                             *tmask(i,k,j)*tmask(i,k,j+1)
            vneg = .5*(totvel - abs(totvel))
     &                             *tmask(i,k,j+1)*tmask(i,k,j)
c           rnormsk is one only in grid points just north of land:
            rnormsk=tmask(i,k,j)*(1-tmask(i,k,j-1))
c           soutmsk is one only in grid points just south of land:
            soutmsk=tmask(i,k,j+1)*(1-tmask(i,k,jp2))
	    adv_fn(i,k,j) = totvel*(
     &                          quick_y(j,1)*t(i,k,j  ,n,tau)
     &                        + quick_y(j,2)*t(i,k,j+1,n,tau))
     &                  - vpos*(curv_yp(j,1)*t(i,k,j+1,n,taum1)
     &                         +curv_yp(j,2)*t(i,k,j  ,n,taum1)
     &                         +curv_yp(j,3)*
     &         (t(i,k,j-1,n,taum1)*(1-rnormsk)
     &         +t(i,k,j  ,n,taum1)*rnormsk))
     &                  - vneg*(curv_yn(j,1)*
     &         (t(i,k,jp2,n,taum1)*(1-soutmsk) +
     &          t(i,k,j+1,n,taum1)*soutmsk)
     &                         +curv_yn(j,2)*t(i,k,j+1,n,taum1)
     &                         +curv_yn(j,3)*t(i,k,j  ,n,taum1))
          enddo
        enddo
      enddo
      end subroutine adv_north_flux_quicker


      subroutine adv_bottom_flux_quicker(n,adv_fb)
c=======================================================================
c     3rd order advective tracer flux
c=======================================================================
      use spflame_module
      use quicker_module
      implicit none
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
      integer is,ie,js,je,n,i,j,k
      real totvel,wpos,wneg,upmsk

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across bottom face of T cells.
c     (It`s done this way for performance issues)
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=2,km-2
          do i=is,ie
            totvel = adv_vbt(i,k,j)
            wpos = .5*(totvel + abs(totvel))
     &             *tmask(i,k+2,j)*tmask(i,k+1,j)             
            wneg = .5*(totvel - abs(totvel))
     &                            *tmask(i,k,j)*tmask(i,k+1,j)
c           upmsk is one only in grid points just above land:
            upmsk=tmask(i,k,j)*(1-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,taum1)
     &                         +curv_zp(k,2)*t(i,k  ,j,n,taum1)
     &                         +curv_zp(k,3)*t(i,k-1,j,n,taum1))
     &                  - wpos*(curv_zn(k,1)*t(i,k+2,j,n,taum1)
     &                         +curv_zn(k,2)*t(i,k+1,j,n,taum1)
     &                         +curv_zn(k,3)*
     &        (t(i,k  ,j,n,taum1)*(1-upmsk)+t(i,k+1,j,n,taum1)*upmsk))
          enddo
        enddo

	k=1
        do i=is,ie
          totvel = adv_vbt(i,k,j)
          wpos = .5*(totvel + abs(totvel))
     &         *tmask(i,k+2,j)*tmask(i,k+1,j)*tmask(i,k,j)

          wneg = .5*(totvel - abs(totvel))
     &         *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,taum1)
     &                         +curv_zp(k,2)*t(i,k  ,j,n,taum1)
     &                         +curv_zp(k,3)*t(i,k  ,j,n,taum1))
     &                - wpos*(curv_zn(k,1)*t(i,k+2,j,n,taum1)
     &                       +curv_zn(k,2)*t(i,k+1,j,n,taum1)
     &                       +curv_zn(k,3)*t(i,k  ,j,n,taum1))
        enddo

	k=km-1
        do i=is,ie
          totvel = adv_vbt(i,k,j)
          wneg = .5*(totvel - abs(totvel))
     &             *tmask(i,k-1,j)*tmask(i,k,j)*tmask(i,k+1,j)
          wpos = .5*(totvel + abs(totvel))
     &                            *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))
     &                - wneg*(curv_zp(k,1)*t(i,k+1,j,n,taum1)
     &                       +curv_zp(k,2)*t(i,k  ,j,n,taum1)
     &                       +curv_zp(k,3)*t(i,k-1,j,n,taum1))
     &                - wpos*(curv_zn(k,1)*t(i,k+1,j,n,taum1)
     &                       +curv_zn(k,2)*t(i,k+1,j,n,taum1)
     &                       +curv_zn(k,3)*t(i,k  ,j,n,taum1))
        enddo
      enddo

      end subroutine adv_bottom_flux_quicker



      subroutine init_adv_flux_4th
c=======================================================================
c     Initialisation of 4.th order advection module
c=======================================================================
      use spflame_module
      implicit none

      if (my_pe==0) then
       print*,' Initialisation of 4.th order advection scheme'
       print*,''
       print*,' using 4.th order advection for all tracers'
       print*,' instead of centered differences'
       print*,''
      endif

      if (enable_quicker_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled Quicker scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_4th')
      endif

      if (enable_fct_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled FCT scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_4th')
      endif

      if (enable_upstream_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled upstream scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_4th')
      endif

      if (my_pe==0) print*,' done'
      end subroutine init_adv_flux_4th



      subroutine adv_flux_4th (n,adv_fe,adv_fn,adv_fb)
      use spflame_module
      implicit none
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
c
c=======================================================================
c     4th order advective tracer flux
c     author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c=======================================================================
c
      integer :: n,i,k,j,is,ie,js,je
      real    :: mask
      real,parameter  :: a2nd = 1.0,     b2nd = 0.0
      real,parameter  :: a4th = 7.0/6.0, b4th = -1.0/6.0

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across eastern face of T cells.
c     (It`s done this way for performance issues)
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
            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
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     2*advective flux across northern face of T cells 
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
            mask = tmask(i,k,j-1)*tmask(i,k,j+2)
            adv_fn(i,k,j) = 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,j+2,n,tau)))
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across bottom face of T cells.
c     (It`s done this way for performance issues)
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=2,km-2
            do i=is,ie
              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
        do i=is,ie
         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,km-1,j)=adv_vbt(i,km-1,j)*(t(i,km-1,j,n,tau) + 
     &                                      t(i,km  ,j,n,tau))
        enddo
      enddo

      end subroutine adv_flux_4th



      subroutine init_adv_flux_upstream
c=======================================================================
c     Initialisation of upstream advection module
c=======================================================================
      use spflame_module
      implicit none

      if (my_pe==0) then
       print*,' Initialisation of upstream advection scheme'
       print*,''
       print*,' using upstream advection for all tracers'
       print*,' instead of centered differences'
       print*,''
      endif

      if (enable_quicker_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled Quicker scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_4th')
      endif

      if (enable_fct_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled FCT scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_4th')
      endif

      if (enable_4th_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled 4th order scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_upstream')
      endif

      if (my_pe==0) print*,' done'
      end subroutine init_adv_flux_upstream



      subroutine adv_flux_upstream (n,adv_fe,adv_fn,adv_fb)
      use spflame_module
      implicit none
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
c
c=======================================================================
c     upstream advective tracer flux
c     author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c=======================================================================
c
      integer :: n,i,k,j,is,ie,js,je
      real    :: totadv
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across northern, 
c     eastern and bottom faces of T cells
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
            totadv = adv_vnt(i,k,j)
            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

      do j=js,je
        do k=1,km
          do i=is-1,ie
            totadv = adv_vet(i,k,j)
            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
      enddo

      do j=js,je
        do k=1,km-1
          do i=is,ie
            totadv = adv_vbt(i,k,j)
            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=is,ie
          adv_fb(i,0,j) = adv_vbt(i,0,j)*2.0*t(i,1,j,n,taum1)
          adv_fb(i,km,j) = 0.
        enddo
      enddo

      end subroutine adv_flux_upstream




      subroutine adv_flux_2nd (n,adv_fe,adv_fn,adv_fb)
      use spflame_module
      implicit none
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
c
c=======================================================================
c     second order advective tracer flux
c     author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c=======================================================================
c
      integer :: n,i,k,j,is,ie,js,je

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across northern, 
c     eastern and bottom faces of T cells
c-----------------------------------------------------------------------
c
       do j=js,je
         do k=1,km
            do i=is-1,ie
	      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

       do j=js-1,je
         do k=1,km
           do i=is,ie
	    adv_fn(i,k,j) = adv_vnt(i,k,j)*(t(i,k,j  ,n,tau) + 
     &                                      t(i,k,j+1,n,tau))
           enddo
         enddo
       enddo

       do j=js,je
         do k=1,km-1
            do i=is,ie
	      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
         do i=is,ie
	    adv_fb(i,0,j)   = adv_vbt(i,0,j)*2.*t(i,1,j,n,tau) 
	    adv_fb(i,km,j)  = 0.
         enddo
       enddo

      end subroutine adv_flux_2nd



      subroutine init_adv_flux_fct
c=======================================================================
c     Initialisation of FCT advection module
c=======================================================================
      use spflame_module
      implicit none
c
c     some compiler switches for FCT module, 
c     could become namelist parameter as well
c
c#define dlimiter_1
#define dlimiter_2
#define fct_3d
#define fct_4th

      if (my_pe==0) then
       print*,' Initialisation of FCT advection scheme'
       print*,''
       print*,' using FCT advection for all tracers'
       print*,' instead of centered differences'
       print*,''
#ifdef dlimiter_1
       print*,' using first form of delimiter '
#endif
#ifdef dlimiter_2
       print*,' using second form of delimiter '
#endif
#ifdef fct_3d
       print*,' delimiting fluxes 3 dimensional as well '
#endif
#ifdef fct_4th
       print*,' using 4.th order scheme for antidiff. fluxes '
#else
       print*,' using 2.th order scheme for antidiff. fluxes '
#endif
      endif

      if (enable_quicker_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled Quicker scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_fct')
      endif

      if (enable_4th_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled 4th order scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_fct')
      endif
      if (enable_upstream_advection) then
       if (my_pe==0) then
        print*,' ERROR:'
        print*,' You have enabled upstream scheme as well '
        print*,' this will not work '
       endif
       call halt_stop(' in init_adv_flux_fct')
      endif

      if (my_pe==0) print*,' done'

      end subroutine init_adv_flux_fct





      subroutine adv_flux_fct (n,adv_fe,adv_fn,adv_fb)
      use spflame_module
      implicit none
      integer :: n
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
c
c=======================================================================
c     computes advective fluxes using a flux-corrected transport scheme
c
c        for reference see
c        Gerdes, R., C. Koeberle and J. Willebrand, 1991
c        the influence of numerical advection schemes on the results of 
c        ocean general circulation models. Clim Dynamics 5, 211-226 
c        and
c        Zalesak, S. T., 1979: Fully multidimensional flux-corrected 
c        transport algorithms for fluids. J. Comp. Phys. 31, 335-362. 
c
c     input: n
c
c     output:
c       adv_fn = 2*advective flux across northern face of T-cell
c       adv_fe = 2*advective flux across eastern face of T-cell
c       adv_fb = 2*advective flux across bottom face of T-cell
c
c     author:   C. Koeberle       e-mail ckoeberl@awi-bremerhaven.de
c
c                                   SPFLAME version   C. Eden April 2002
c=======================================================================
c               
      integer :: i,k,j,is,ie,js,je
      real :: fxa,fxb,totadv
      real,parameter :: epsln   = 1.0e-15  !  for single precision
      real, dimension(is_pe-2:ie_pe+2) :: dcf, Trmin, Trmax,
     &                            Cpos, Cneg, flxlft, flxrgt
      real :: Rpl(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) 
      real :: Rmn(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) 
      real :: t_lo(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) 
      real :: anti_fn(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) 
      real :: anti_fe(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) 
      real :: anti_fb(is_pe-1:ie_pe+1,0:km,js_pe-1:je_pe+1) 

      real :: adv_fn_hi(is_pe:ie_pe,km,js_pe-1:je_pe)
      real :: adv_fe_hi(is_pe-1:ie_pe,km,js_pe:je_pe)
      real :: adv_fb_hi(is_pe:ie_pe,0:km,js_pe:je_pe)
# ifdef fct_3d
      real :: tmpext(is_pe:ie_pe,km,js_pe:je_pe,2)
# endif
c
c-----------------------------------------------------------------------
c      FCT statement functions
c
c      calculate some variables needed as one-dimensional delimiters to 
c      raw anti-diffusive fluxes 
c
c      1) calculate P+/P-: maximal/minimal possible growth/decay of T
c      2) calculate Q+/Q-: maximal/minimal allowed growth/decay of T 
c           the appropriate choice of this "allowed" is the key
c           for the success of this scheme...
c      3) calculate ration R+- of Q+- to P+-, that is maximal/minimal 
c         possible change of T if no limit would be active, 
c         must be at least 1
c
c-----------------------------------------------------------------------
c
c     here the calculations of the flux delimiters are done as statement
c     functions for simplicity always loop along one row of POINTS
c 
c     maximal/minimal possible change at grid point
c     this variable is related to a point
c
      real :: pplus,pminus
      Pplus(i)  = c2dt*dcf(i)*(max(0.,flxlft(i))-min(0.,flxrgt(i)))
      Pminus(i) = c2dt*dcf(i)*(max(0.,flxrgt(i))-min(0.,flxlft(i)))
c
c     maximal/minimal feasible change at grid point
c     this variable is related to a point
c
      real :: qplus,qminus
      Qplus(i)  = Trmax(i) - t_lo(i,k,j) 
      Qminus(i) = t_lo(i,k,j) - Trmin(i)
c
c     ratio of Q and P to get delimiter
c     this variable is related to a point
c
      real :: Rpls,Rmns
      Rpls(i)  = min(1.,tmask(i,k,j)*Qplus(i)/(Pplus(i)+epsln)) 
      Rmns(i)  = min(1.,tmask(i,k,j)*Qminus(i)/(Pminus(i)+epsln))
c
c-----------------------------------------------------------------------
c     advective terms
c-----------------------------------------------------------------------
c
      real :: adv_tx
      ADV_Tx(i,k,j) = (adv_fe(i,k,j) - adv_fe(i-1,k,j))*cstdxt2r(i,j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      real :: adv_ty
      ADV_Ty(i,k,j) = (adv_fn(i,k,j) - adv_fn(i,k,j-1))*cstdyt2r(j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      real :: adv_tz
      ADV_Tz(i,k,j) = (adv_fb(i,k-1,j) - adv_fb(i,k,j))
#ifdef partial_cell
     &               *.5/dht(i,k,j)
#else
     &               *dzt2r(k)
#endif

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective (low order scheme) flux across northern, 
c     eastern and bottom faces of T cells
c-----------------------------------------------------------------------
c
      call adv_flux_upstream (n,adv_fe,adv_fn,adv_fb)
c
c-----------------------------------------------------------------------
c       solve for "tau+1" tracer at center of T cells 
c       - low order solution -
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          do i=is,ie
	    t_lo(i,k,j) = (t(i,k,j,n,taum1) - c2dt
     &                      *(ADV_Tx(i,k,j) + ADV_Ty(i,k,j) + 
     &                        ADV_Tz(i,k,j))*tmask(i,k,j))
          enddo
        enddo
      enddo

#ifdef dlimiter_2
      call set_cyclic(t_lo,km,1)
      call border_exchg(t_lo,km,1)
      if (my_blk_j == n_pes_j ) t_lo(:,:,jmt)=0. 
      if (my_blk_j ==1 )        t_lo(:,:,1)=0.
#endif
c
c-----------------------------------------------------------------------
c       next calculate raw antidiffusive fluxes, that is high order 
c       scheme flux (leap frog) minus the low order (upstream) 
c-----------------------------------------------------------------------
c
#ifdef fct_4th
      call adv_flux_4th (n,adv_fe_hi,adv_fn_hi,adv_fb_hi)
#else
      call adv_flux_2nd (n,adv_fe_hi,adv_fn_hi,adv_fb_hi)
#endif

      do j=js,je
        do k=1,km
          do i=is-1,ie
            anti_fe(i,k,j) = adv_fe_hi(i,k,j)-adv_fe(i,k,j)
          enddo
        enddo
      enddo

      do j=js-1,je
        do k=1,km
          do i=is,ie
            anti_fn(i,k,j) = adv_fn_hi(i,k,j)-adv_fn(i,k,j)
          enddo
        enddo
      enddo
c
      do j=js,je
        do k=0,km
          do i=is,ie
            anti_fb(i,k,j) = adv_fb_hi(i,k,j)-adv_fb(i,k,j)
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c       now calculate and apply one-dimensional delimiters to these 
c       raw antidiffusive fluxes 
c
c       1) calculate T*, that are all halfway neighbors of T
c       2) calculate ratio R+- of Q+- to P+-, that is maximal/minimal 
c          possible change of T if no limit would be active, 
c          must be at least 1
c       3) choose correct ratio depending on direction of flow as a 
c          delimiter
c       4) apply this delimiter to raw antidiffusive flux
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c       delimit x-direction
c-----------------------------------------------------------------------
c
      do j=js,je 
        do k=1,km
c
c         prepare some data for use in statement function   
# ifdef dlimiter_1
c
c         running mean of two adjacent points
c
          do i=is,ie+1
            Trmax(i) = 0.5*(t(i-1,k,j,n,tau) + t(i,k,j,n,tau))
          enddo
# endif
c
c         extremum of low order solution central point and adjacent 
c         halfway neighbours; check for land
c
          do i=is,ie
# ifdef dlimiter_1
            fxa = tmask(i-1,k,j)*Trmax(i) +
     &            (1.-tmask(i-1,k,j))*t_lo(i,k,j)
            fxb = tmask(i+1,k,j)*Trmax(i+1) +
     &            (1.-tmask(i+1,k,j))*t_lo(i,k,j)
# endif
# ifdef dlimiter_2
            fxa = tmask(i-1,k,j)*t_lo(i-1,k,j) +
     &            (1.-tmask(i-1,k,j))*t_lo(i,k,j)
            fxb = tmask(i+1,k,j)*t_lo(i+1,k,j) +
     &            (1.-tmask(i+1,k,j))*t_lo(i,k,j)
# endif
            Trmax(i) = max(fxa,fxb,t_lo(i,k,j))
            Trmin(i) = min(fxa,fxb,t_lo(i,k,j))
# ifdef fct_3d
            tmpext(i,k,j,1) = Trmax(i)
            tmpext(i,k,j,2) = Trmin(i)
# endif
            dcf(i) = cstdxt2r(i,j)
            flxlft(i) = anti_fe(i-1,k,j)
            flxrgt(i) = anti_fe(i,k,j)
          enddo
c
c         calculate ratio R           
c
          do i=is,ie
            Rpl(i,k,j) = Rpls(i)
            Rmn(i,k,j) = Rmns(i)
          enddo
        enddo
      enddo

      call set_cyclic(Rmn,km,1)
      call border_exchg(Rmn,km,1)
      call set_cyclic(Rpl,km,1)
      call border_exchg(Rpl,km,1)

      do j=js,je 
        do k=1,km
c
c         calculate delimiter using ratio at adjacent points         
c
          do i=is-1,ie
            Cpos(i) = min(Rpl(i+1,k,j),Rmn(i,k,j))
            Cneg(i) = min(Rpl(i,k,j),Rmn(i+1,k,j))
          enddo
c
c         finally apply appropriate delimiter to flux
c
          do i=is-1,ie
            anti_fe(i,k,j) = 0.5*((Cpos(i) + Cneg(i))
     &                               *anti_fe(i,k,j) +
     &                               (Cpos(i) - Cneg(i))
     &                               *abs(anti_fe(i,k,j)))
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c       delimit y-direction
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
c
c         prepare some data for use in statement function       
c
          do i=is,ie
# ifdef dlimiter_1
            fxa = 0.5*tmask(i,k,j-1)*(t(i,k,j-1,n,tau) + 
     &            t(i,k,j,n,tau)) +
     &            (1.-tmask(i,k,j-1))*t_lo(i,k,j)
            fxb = 0.5*tmask(i,k,j+1)*(t(i,k,j,n,tau) +
     &            t(i,k,j+1,n,tau)) +
     &            (1.-tmask(i,k,j+1))*t_lo(i,k,j)
# endif
# ifdef dlimiter_2
            fxa = tmask(i,k,j-1)*t_lo(i,k,j-1) +
     &            (1.-tmask(i,k,j-1))*t_lo(i,k,j)
            fxb = 0.5*tmask(i,k,j+1)*(t(i,k,j,n,tau) +
     &            t(i,k,j+1,n,tau)) +
     &            (1.-tmask(i,k,j+1))*t_lo(i,k,j)
# endif
            Trmax(i) = max(fxa,fxb,t_lo(i,k,j))
            Trmin(i) = min(fxa,fxb,t_lo(i,k,j))
# ifdef fct_3d
            tmpext(i,k,j,1) = max(Trmax(i),tmpext(i,k,j,1))
            tmpext(i,k,j,2) = min(Trmin(i),tmpext(i,k,j,2))
# endif
            dcf(i) = cstdyt2r(j)
            flxlft(i) = anti_fn(i,k,j-1)
            flxrgt(i) = anti_fn(i,k,j)
          enddo
c
c         calculate ratio R, related to a point          
c
          do i=is,ie
            Rpl(i,k,j) = Rpls(i)
            Rmn(i,k,j) = Rmns(i)
          enddo
        enddo
      enddo

      call border_exchg(Rmn,km,1)
      call border_exchg(Rpl,km,1)

      if (my_blk_j == n_pes_j ) then
       Rmn(:,:,jmt)=0. !! really zero ??
       Rpl(:,:,jmt)=0.
      endif
      if (my_blk_j ==1 ) then
       Rmn(:,:,1)=0.
       Rpl(:,:,1)=0.
      endif

      do j=js-1,je 
        do k=1,km
c
c         calculate delimiter using ratio at adjacent points
c
          do i=is,ie
            Cpos(i) = min(Rpl(i,k,j+1),Rmn(i,k,j))
            Cneg(i) = min(Rpl(i,k,j),Rmn(i,k,j+1))
          enddo
c
c         finally get delimiter c dependent on direction of flux and 
c         apply it to raw antidiffusive flux
c
          do i=is,ie
            anti_fn(i,k,j) = 0.5*((Cpos(i) + Cneg(i))
     &                             *anti_fn(i,k,j) +
     &                             (Cpos(i) - Cneg(i))
     &                             *abs(anti_fn(i,k,j)))
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c       delimit z-direction
c-----------------------------------------------------------------------
c
      do j=js,je 
        do k=1,km
c
c         prepare some data for use in statement function
c
          do i=is,ie
            dcf(i) = 
# ifdef partial_cell
     &               0.5/dht(i,k,j) ! unchanged j, probably wrong in MOM3
# else
     &               dzt2r(k)
# endif
            flxlft(i) = anti_fb(i,k,j)
            flxrgt(i) = anti_fb(i,k-1,j)
            if (k .gt. 1)then
# ifdef dlimiter_1
              fxa = 0.5*tmask(i,k-1,j)*
     &              (t(i,k-1,j,n,tau) + t(i,k,j,n,tau)) + 
     &              (1.-tmask(i,k-1,j))*t_lo(i,k,j)
# endif
# ifdef dlimiter_2
              fxa = tmask(i,k-1,j)*t_lo(i,k-1,j) +
     &              (1.-tmask(i,k-1,j))*t_lo(i,k,j)
# endif
            else
              fxa = t_lo(i,k,j)
            endif
# ifdef dlimiter_1
            if (k .lt. km) then
              fxb = 0.5*tmask(i,k+1,j)*
     &              (t(i,k,j,n,tau)+t(i,k+1,j,n,tau)) +
     &              (1.-tmask(i,k+1,j))*t_lo(i,k,j)
# endif
# ifdef dlimiter_2
            if (k .lt. km) then
              fxb = tmask(i,k+1,j)*t_lo(i,k+1,j) +
     &              (1.-tmask(i,k+1,j))*t_lo(i,k,j)
# endif
            else
              fxb = t_lo(i,k,j)
            endif
            Trmax(i) = max(fxa,fxb,t_lo(i,k,j))
            Trmin(i) = min(fxa,fxb,t_lo(i,k,j))
# ifdef fct_3d
            tmpext(i,k,j,1) = max(Trmax(i),tmpext(i,k,j,1))
            tmpext(i,k,j,2) = min(Trmin(i),tmpext(i,k,j,2))
# endif
          enddo
c
c         calculate delimiter using ratio at adjacent points
c         this variable is related to an arc (between two points,
c         the same way as fluxes are defined.)
c
          do i=is,ie
            Rpl(i,k,j) = Rpls(i)
            Rmn(i,k,j) = Rmns(i)
          enddo
c
        enddo
      enddo

      do j=js,je 
        do k=1,km-1
c
c         calculate delimiter using ratio at adjacent points
c         this variable is related to an arc (between two points,
c         the same way as fluxes are defined.)
c
          do i=is,ie
            Cneg(i) = min(Rpl(i,k+1,j),Rmn(i,k,j))
            Cpos(i) = min(Rpl(i,k,j),Rmn(i,k+1,j))
          enddo
c
c         finally get delimiter c dependent on direction of flux and 
c         apply it to raw antidiffusive flux
c
          do i=is,ie
            anti_fb(i,k,j) = 0.5*((Cpos(i)+Cneg(i))
     &                               *anti_fb(i,k,j) +
     &                               (Cpos(i)-Cneg(i))
     &                               *abs(anti_fb(i,k,j)))
          enddo
        enddo
        do i=is,ie
          anti_fb(i,0,j)  = 0.
          anti_fb(i,km,j) = 0.
        enddo

      enddo ! j loop


# ifdef fct_3d
c
c-----------------------------------------------------------------------
c     then calculate and apply 3-d delimiter to just delimited 
c     antidiffusive fluxes   
c-----------------------------------------------------------------------
c
      do j=js,je 
        do k=1,km
c
c         prepare some data for use in statement function
c
          do i=is,ie
            Trmax(i) = tmpext(i,k,j,1)
            Trmin(i) = tmpext(i,k,j,2)
          enddo
c
          do i=is,ie
            Rpl(i,k,j) = min(1.,tmask(i,k,j)*Qplus(i)/
     &                           (epsln+c2dt*(
     &                           cstdxt2r(i,j)
     &                             *(max(0.,anti_fe(i-1,k,j)) - 
     &                               min(0.,anti_fe(i,k,j))) +
     &                           cstdyt2r(j)
     &                             *(max(0.,anti_fn(i,k,j-1)) -
     &                               min(0.,anti_fn(i,k,j))) +
#  ifdef partial_cell
     &                          (0.5/dht(i,k,j)) ! unchanged j
#  else
     &                           dzt2r(k)
#  endif
     &                             *(max(0.,anti_fb(i,k,j)) -
     &                               min(0.,anti_fb(i,k-1,j)))
     &                           )))
c
            Rmn(i,k,j) = min(1.,tmask(i,k,j)*Qminus(i)/
     &                            (epsln+c2dt*(
     &                            cstdxt2r(i,j)
     &                              *(max(0.,anti_fe(i,k,j)) -
     &                                min(0.,anti_fe(i-1,k,j))) +
     &                            cstdyt2r(j)
     &                              *(max(0.,anti_fn(i,k,j)) -
     &                                min(0.,anti_fn(i,k,j-1))) +
#  ifdef partial_cell
     &                          (0.5/dht(i,k,j)) ! unchanged j
#  else
     &                           dzt2r(k)
#  endif
     &                              *(max(0.,anti_fb(i,k-1,j)) -
     &                                min(0.,anti_fb(i,k,j)))
     &                            )))
c
          enddo
        enddo
      enddo

      call set_cyclic(Rmn,km,1)
      call border_exchg(Rmn,km,1)
      call set_cyclic(Rpl,km,1)
      call border_exchg(Rpl,km,1)

c
c       finally apply 3-d delimiters to precorrected fluxes
c
      do j=js,je 
        do k=1,km
          do i=is-1,ie
            Cpos(i) = min(Rpl(i+1,k,j),Rmn(i,k,j))
            Cneg(i) = min(Rpl(i,k,j),  Rmn(i+1,k,j))
          enddo
          do i=is-1,ie
            anti_fe(i,k,j) = 0.5*((Cpos(i) + Cneg(i))
     &                               *anti_fe(i,k,j) +
     &                               (Cpos(i) - Cneg(i))
     &                               *abs(anti_fe(i,k,j)))
          enddo
        enddo
      enddo
c
      do j=js-1,je 
        do k=1,km
          do i=is,ie
            Cpos(i) = min(Rpl(i,k,j+1),Rmn(i,k,j))
            Cneg(i) = min(Rpl(i,k,j),Rmn(i,k,j+1))
          enddo
          do i=is,ie
            anti_fn(i,k,j) = 0.5*((Cpos(i) + Cneg(i))
     &                             *anti_fn(i,k,j) +
     &                             (Cpos(i) - Cneg(i))
     &                             *abs(anti_fn(i,k,j)))
          enddo
        enddo
      enddo
c
      do j=js,je 
        do k=1,km-1
          do i=is,ie
            Cneg(i) = min(Rpl(i,k+1,j),Rmn(i,k,j))
            Cpos(i) = min(Rpl(i,k,j),Rmn(i,k+1,j))
          enddo
          do i=is,ie
            anti_fb(i,k,j) = 0.5*((Cpos(i) + Cneg(i))
     &                               *anti_fb(i,k,j) +
     &                               (Cpos(i) - Cneg(i))
     &                               *abs(anti_fb(i,k,j)))
          enddo
        enddo
      enddo
c
# endif
c
c-----------------------------------------------------------------------
c     complete advective fluxes by adding low order fluxes to 
c     delimited antidiffusive fluxes    
c-----------------------------------------------------------------------
c
      do j=js,je 
        do k=1,km
          do i=is-1,ie
            adv_fe(i,k,j) = (anti_fe(i,k,j) + adv_fe(i,k,j))
     &                         *tmask(i,k,j)
          enddo
        enddo
      enddo

      do j=js-1,je 
        do k=1,km
          do i=is,ie
            adv_fn(i,k,j) = (anti_fn(i,k,j) + adv_fn(i,k,j))
     &                         *tmask(i,k,j)
          enddo
        enddo
      enddo 

      do j=js,je 
        do k=1,km
          do i=is,ie
            adv_fb(i,k,j) = (anti_fb(i,k,j) + adv_fb(i,k,j))
     &                           *tmask(i,k,j)
          enddo
        enddo
      enddo 

      end subroutine adv_flux_fct


