#include "options.inc"
    
c=======================================================================
c      Different advection schemes for buoyancy and other tracers
c=======================================================================


      subroutine adv_flux(adv_fe,adv_fn,adv_ft,var)
c=======================================================================
c       advective tracer flux:  call different schemes
c=======================================================================
      use cpflame_module
      implicit none
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), var(imt,jmt,km,0:2)

      if (enable_quicker_advection) then
       call setcyclic3D_j2(var(:,:,:,taum1) )
       call adv_flux_quicker(adv_fe,adv_fn,adv_ft,var)
      else if (enable_4th_advection) then
       call setcyclic3D_j2(var(:,:,:,tau) )
       call adv_flux_4th(adv_fe,adv_fn,adv_ft,var)
      else if (enable_upwind_advection) then
       call adv_flux_upw(adv_fe,adv_fn,adv_ft,var)
      else if (enable_lax_advection) then
       call adv_flux_lax_wendroff(adv_fe,adv_fn,adv_ft,var)
      else
       call adv_flux_2th(adv_fe,adv_fn,adv_ft,var)
      endif
      end subroutine adv_flux




      subroutine adv_flux_2th(adv_fe,adv_fn,adv_ft,var)
c=======================================================================
c      2th order advective tracer flux
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), var(imt,jmt,km,0:2)
      real :: totvel

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         totvel = u(i,j,k,1,tau)
         adv_fe(i,j,k)=
     &  0.5*(var(i,j,k,tau) + var(i+1,j,k,tau) )*totvel
     &    *maskU(i,j,k)
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         totvel = u(i,j,k,2,tau)
         adv_fn(i,j,k)=
     &  0.5*( var(i,j,k,tau) + var(i,j+1,k,tau) )*totvel
     &    *maskV(i,j,k)
        enddo
       enddo
      enddo
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         totvel = u(i,j,k,3,tau)
         adv_ft(i,j,k)=
     &  0.5*( var(i,j,k,tau) + var(i,j,k+1,tau) )*totvel
     &    *maskW(i,j,k)
        enddo
       enddo
      enddo
      end subroutine adv_flux_2th





      subroutine adv_flux_4th(adv_fe,adv_fn,adv_ft,var)
c=======================================================================
c      4th order advective tracer flux
c      author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je,ip2,im1,jp2,jm1,kp2,km1
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), var(imt,jmt,km,0:2)
      real,parameter  :: a2nd = 1.0,     b2nd = 0.0
      real,parameter  :: a4th = 7.0/6.0, b4th = -1.0/6.0
      real :: mask

      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
       do k=2,km-1
        do i=1,imt-1
         ip2=i+2
         if (ip2>imt.and..not.enable_cyclic_x) ip2=imt
         if (ip2>imt.and.     enable_cyclic_x) ip2=3
         im1=i-1
         if (im1<1.and..not. enable_cyclic_x) im1=1
         if (im1<1.and.      enable_cyclic_x) im1=imt-2
         mask = maskT(im1,j,k)*maskT(ip2,j,k)
         adv_fe(i,j,k) = 0.5*u(i,j,k,1,tau)*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(var(i,  j,k,tau) + 
     &                                      var(i+1,j,k,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(var(im1,j,k,tau) + 
     &                                      var(ip2,j,k,tau)))
        enddo
       enddo
      enddo

      do j=js,je
       do k=2,km-1
        do i=2,imt-1
         jp2=j+2
         if (jp2>jmt.and..not.enable_cyclic_y) jp2=jmt
         if (jp2>jmt.and.     enable_cyclic_y) jp2=3
         jm1=j-1
         if (jm1<1.and..not. enable_cyclic_y) jm1=1
         if (jm1<1.and.      enable_cyclic_y) jm1=jmt-2
         mask = maskT(i,jm1,k)*maskT(i,jp2,k)
         adv_fn(i,j,k) = 0.5*u(i,j,k,2,tau)*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(var(i,  j,k,tau) + 
     &                                      var(i,j+1,k,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(var(i,jm1,k,tau) + 
     &                                      var(i,jp2,k,tau)))
        enddo
       enddo
      enddo
c
      do j=js,je
       do k=1,km-1
        do i=2,imt-1
         kp2=min(k+2,km)
         km1=max(1,k-1)
         mask = maskT(i,j,km1)*maskT(i,j,kp2)
         adv_ft(i,j,k) = 0.5*u(i,j,k,3,tau)*(
     &         (a2nd*(1.0-mask) + a4th*mask)*(var(i,j,  k,tau) + 
     &                                        var(i,j,k+1,tau))+
     &         (b2nd*(1.0-mask) + b4th*mask)*(var(i,j,km1,tau) + 
     &                                        var(i,j,kp2,tau)))
        enddo
       enddo
      enddo

      end subroutine adv_flux_4th






      subroutine adv_flux_quicker(adv_fe,adv_fn,adv_ft,var)
c=======================================================================
c      3th order advective tracer flux
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,j,k,ip2,jp2,kp2,im1,jm1,km1,js,je
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), var(imt,jmt,km,0:2)
      real :: totvel,upos,uneg,eastmsk,westmsk
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         ip2=i+2
         if (ip2>imt.and..not.enable_cyclic_x) ip2=imt
         if (ip2>imt.and.     enable_cyclic_x) ip2=3
         im1=i-1
         if (im1<1.and..not. enable_cyclic_x) im1=1
         if (im1<1.and.      enable_cyclic_x) im1=imt-2

         totvel = u(i,j,k,1,tau)
         upos = .5*(totvel + abs(totvel))
     &                  *maskt(i,j,k)*maskt(i+1,j,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskt(i+1,j,k)*maskt(i,j,k)
         eastmsk=maskt(i,j,k)*(1-maskt(im1,j,k))
         westmsk=maskt(i+1,j,k)*(1-maskt(ip2,j,k))
	 adv_fe(i,j,k) =
     &        0.5*totvel*(var(i,j,k,tau)+var(i+1,j,k,tau))
     &        - 0.5*upos*(0.25*var(i+1,j,k,taum1)-0.5*var(i,j,k,taum1)
     &               +0.25*(var(im1,j,k,taum1)*(1-eastmsk)
     &                          +var(i,j,k,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(var(ip2,j,k,taum1)*(1-westmsk)+
     &          var(i+1,j,k,taum1)*westmsk)-0.5*var(i+1,j,k,taum1)
     &                         +0.25*var(i,j,k,taum1))
          enddo
        enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         jp2=j+2
         if (jp2>jmt.and..not.enable_cyclic_y) jp2=jmt
         if (jp2>jmt.and.     enable_cyclic_y) jp2=3
         jm1=j-1
         if (jm1<1.and..not. enable_cyclic_y) jm1=1
         if (jm1<1.and.      enable_cyclic_y) jm1=jmt-2

         totvel = u(i,j,k,2,tau)
         upos = .5*(totvel + abs(totvel))
     &                  *maskt(i,j,k)*maskt(i,j+1,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskt(i,j+1,k)*maskt(i,j,k)
         eastmsk=maskt(i,j,k)*(1-maskt(i,jm1,k))
         westmsk=maskt(i,j+1,k)*(1-maskt(i,jp2,k))
	 adv_fn(i,j,k) =
     &        0.5*totvel*(var(i,j,k,tau)+var(i,j+1,k,tau))
     &        - 0.5*upos*(0.25*var(i,j+1,k,taum1)-0.5*var(i,j,k,taum1)
     &               +0.25*(var(i,jm1,k,taum1)*(1-eastmsk)
     &                          +var(i,j,k,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(var(i,jp2,k,taum1)*(1-westmsk)+
     &          var(i,j+1,k,taum1)*westmsk)-0.5*var(i,j+1,k,taum1)
     &                         +0.25*var(i,j,k,taum1))
          enddo
        enddo
      enddo
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         kp2=min(k+2,km)
         km1=max(1,k-1)
         totvel = u(i,j,k,3,tau)
         upos = .5*(totvel + abs(totvel))
     &                  *maskt(i,j,k)*maskt(i,j,k+1)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskt(i,j,k+1)*maskt(i,j,k)
         eastmsk=maskt(i,j,k)*(1-maskt(i,j,km1))
         westmsk=maskt(i,j,k+1)*(1-maskt(i,j,kp2))
	 adv_ft(i,j,k) =
     &        0.5*totvel*(var(i,j,k,tau)+var(i,j,k+1,tau))
     &        - 0.5*upos*(0.25*var(i,j,k+1,taum1)-0.5*var(i,j,k,taum1)
     &               +0.25*(var(i,j,km1,taum1)*(1-eastmsk)
     &                          +var(i,j,k,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(var(i,j,kp2,taum1)*(1-westmsk)+
     &          var(i,j,k+1,taum1)*westmsk)-0.5*var(i,j,k+1,taum1)
     &                         +0.25*var(i,j,k,taum1))
          enddo
        enddo
      enddo
      end subroutine adv_flux_quicker






      subroutine adv_flux_lax_wendroff(adv_fe,adv_fn,adv_ft,var)
c=======================================================================
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), var(imt,jmt,km,0:2)
      real :: totvel,d0,d1,absc

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         totvel = u(i,j,k,1,tau)
         absc=abs(totvel*dt/dx)
         adv_fe(i,j,k)=maskU(i,j,k)*0.5*( 
     &        totvel*(var(i+1,j,k,tau)+var(i,j,k,tau))
     &     -abs(totvel)*absc*(var(i+1,j,k,taum1)-var(i,j,k,taum1)) )
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         totvel = u(i,j,k,2,tau)
         absc=abs(totvel*dt/dx)
         adv_fn(i,j,k)=maskV(i,j,k)*0.5*( 
     &        totvel*(var(i,j+1,k,tau)+var(i,j,k,tau))
     &     -abs(totvel)*absc*(var(i,j+1,k,taum1)-var(i,j,k,taum1)) )
        enddo
       enddo
      enddo
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         totvel = u(i,j,k,3,tau)
         absc=abs(totvel*dt/dz)
         adv_ft(i,j,k)=maskW(i,j,k)*0.5*( 
     &        totvel*(var(i,j,k+1,tau)+var(i,j,k,tau))
     &     -abs(totvel)*absc*(var(i,j,k+1,taum1)-var(i,j,k,taum1)) )
        enddo
       enddo
      enddo
      end subroutine adv_flux_lax_wendroff




      subroutine adv_flux_upw(adv_fe,adv_fn,adv_ft,var)
c=======================================================================
c        upwind tracer flux:   F = U b^i 
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), var(imt,jmt,km,0:2)
      real :: totvel

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         totvel = u(i,j,k,1,tau)
         adv_fe(i,j,k)=maskU(i,j,k)*( 
     &     0.5*totvel*(var(i,j,k,taum1)+var(i+1,j,k,taum1))
     &   - 0.5*abs(totvel)*(var(i+1,j,k,taum1)-var(i,j,k,taum1) ) )
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         totvel = u(i,j,k,2,tau)
         adv_fn(i,j,k)=maskV(i,j,k)*( 
     &     0.5*totvel*(var(i,j,k,taum1)+var(i,j+1,k,taum1))
     &   - 0.5*abs(totvel)*(var(i,j+1,k,taum1)-var(i,j,k,taum1) ) )
        enddo
       enddo
      enddo
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         totvel = u(i,j,k,3,tau)
         adv_ft(i,j,k)=maskW(i,j,k)*( 
     &     0.5*totvel*(var(i,j,k,taum1)+var(i,j,k+1,taum1))
     &   - 0.5*abs(totvel)*(var(i,j,k+1,taum1)-var(i,j,k,taum1) ) )
        enddo
       enddo
      enddo
      end subroutine adv_flux_upw





