#include "options.inc"

c=======================================================================
c      Advection schemes for momentum
c=======================================================================


      subroutine adv_flux_u_2nd(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of zonal momentum
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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) 
      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
         adv_fe(i,j,k)=
     &  0.5*( u(i  ,j,k,1,tau) + u(i+1,j,k,1,tau) )*
     &   (u(i,j,k,1,tau)+u(i+1,j,k,1,tau))
     &        *0.5*maskU(i+1,j,k)*maskU(i,j,k)
        enddo
       enddo
      enddo
      call setcyclic3D(adv_fe)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         adv_fn(i,j,k)=
     &  0.5*( u(i,j,k,1,tau) + u(i,j+1,k,1,tau) )*
     &   (u(i,j,k,2,tau)+u(i+1,j,k,2,tau))
     &        *0.5*maskU(i,j+1,k)*maskU(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         adv_ft(i,j,k)=
     &  0.5*(u(i,j,k,1,tau)+u(i,j,k+1,1,tau))*
     &   (u(i,j,k,3,tau)+u(i+1,j,k,3,tau))
     &        *0.5*maskU(i,j,k+1)*maskU(i,j,k)
        enddo
       enddo
      enddo
      end subroutine adv_flux_u_2nd


      subroutine adv_flux_v_2nd(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of meridional momentum
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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) 
      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
         adv_fe(i,j,k)=
     &  0.5*( u(i  ,j,k,2,tau) + u(i+1,j,k,2,tau) )*
     &    (u(i,j,k,1,tau)+u(i,j+1,k,1,tau))
     &        *0.5*maskV(i+1,j,k)*maskV(i,j,k)
        enddo
       enddo
      enddo
      call setcyclic3D(adv_fe)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         adv_fn(i,j,k)=
     &  0.5*( u(i,j,k,2,tau) + u(i,j+1,k,2,tau) )*
     &   (u(i,j,k,2,tau)+u(i,j+1,k,2,tau))
     &        *0.5*maskV(i,j+1,k)*maskV(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         adv_ft(i,j,k)=
     &  0.5*( u(i,j,k,2,tau) + u(i,j,k+1,2,tau) )*
     &   (u(i,j,k,3,tau)+u(i,j+1,k,3,tau))
     &        *0.5*maskV(i,j,k+1)*maskV(i,j,k)
        enddo
       enddo
      enddo
      end subroutine adv_flux_v_2nd




      subroutine adv_flux_w_2nd(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of vertical momentum
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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) 
      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
         adv_fe(i,j,k)=
     &     0.5*( u(i  ,j,k,3,tau) + u(i+1,j,k,3,tau) )*
     &         (u(i,j,k,1,tau)+u(i,j,k+1,1,tau))
     &        *0.5*maskW(i+1,j,k)*maskW(i,j,k)
        enddo
       enddo
      enddo
      call setcyclic3D(adv_fe)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         adv_fn(i,j,k)=
     &     0.5*( u(i,j,k,3,tau) + u(i,j+1,k,3,tau) )*
     &     (u(i,j,k,2,tau)+u(i,j,k+1,2,tau))
     &        *0.5*maskW(i,j+1,k)*maskW(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         adv_ft(i,j,k)=
     &     0.5*( u(i,j,k,3,tau) + u(i,j,k+1,3,tau) )*
     &     (u(i,j,k,3,tau)+u(i,j,k+1,3,tau))
     &        *0.5*maskW(i,j,k+1)*maskW(i,j,k)
        enddo
       enddo
      enddo
      end subroutine adv_flux_w_2nd



      subroutine adv_flux_u_4th(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of zonal momentum 4th order
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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)
      real,parameter  :: a2nd = 1.0,     b2nd = 0.0
      real,parameter  :: a4th = 7.0/6.0, b4th = -1.0/6.0
      real :: mask,totvel

      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 = maskU(im1,j,k)*maskU(ip2,j,k)
c         totvel = 0.5*(u(i+1,j,k,1,tau)+u(i,j,k,1,tau))
         totvel = (u(i+1,j,k,1,tau)*maskU(i+1,j,k)
     &             +u(i,j,k,1,tau)*maskU(i,j,k))/2.0
         adv_fe(i,j,k) = 0.5*totvel*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(u(i,  j,k,1,tau) + 
     &                                      u(i+1,j,k,1,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(u(im1,j,k,1,tau) + 
     &                                      u(ip2,j,k,1,tau)))
        enddo
       enddo
      enddo
      call setcyclic3D(adv_fe)

      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 = maskU(i,jm1,k)*maskU(i,jp2,k)
c         totvel = 0.5*(u(i+1,j,k,2,tau)+u(i,j,k,2,tau))
         totvel = (u(i+1,j,k,2,tau)*maskV(i+1,j,k)
     &             +u(i,j,k,2,tau)*maskV(i,j,k))/2.0
         adv_fn(i,j,k) = 0.5*totvel*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(u(i,  j,k,1,tau) + 
     &                                      u(i,j+1,k,1,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(u(i,jm1,k,1,tau) + 
     &                                      u(i,jp2,k,1,tau)))
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)
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 = maskU(i,j,km1)*maskU(i,j,kp2)
c         totvel = 0.5*(u(i+1,j,k,3,tau)+u(i,j,k,3,tau))
         totvel = (u(i+1,j,k,3,tau)*maskW(i+1,j,k)
     &             +u(i,j,k,3,tau)*maskW(i,j,k))/2.0
         adv_ft(i,j,k) = 0.5*totvel*(
     &         (a2nd*(1.0-mask) + a4th*mask)*(u(i,j,  k,1,tau) + 
     &                                        u(i,j,k+1,1,tau))+
     &         (b2nd*(1.0-mask) + b4th*mask)*(u(i,j,km1,1,tau) + 
     &                                        u(i,j,kp2,1,tau)))
        enddo
       enddo
      enddo
      end subroutine adv_flux_u_4th



      subroutine adv_flux_v_4th(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of meridional momentum 4th order
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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)
      real,parameter  :: a2nd = 1.0,     b2nd = 0.0
      real,parameter  :: a4th = 7.0/6.0, b4th = -1.0/6.0
      real :: mask,totvel

      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 = maskV(im1,j,k)*maskV(ip2,j,k)
c         totvel = 0.5*(u(i,j+1,k,1,tau)+u(i,j,k,1,tau))
         totvel = (u(i,j+1,k,1,tau)*maskU(i,j+1,k)
     &             +u(i,j,k,1,tau)*maskU(i,j,k))/2.0
         adv_fe(i,j,k) = 0.5*totvel*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(u(i,  j,k,2,tau) + 
     &                                      u(i+1,j,k,2,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(u(im1,j,k,2,tau) + 
     &                                      u(ip2,j,k,2,tau)))
        enddo
       enddo
      enddo
      call setcyclic3D(adv_fe)

      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 = maskV(i,jm1,k)*maskV(i,jp2,k)
c         totvel = 0.5*(u(i,j+1,k,2,tau)+u(i,j,k,2,tau))
         totvel = (u(i,j+1,k,2,tau)*maskV(i,j+1,k)
     &             +u(i,j,k,2,tau)*maskV(i,j,k))/2.0
         adv_fn(i,j,k) = 0.5*totvel*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(u(i,  j,k,2,tau) + 
     &                                      u(i,j+1,k,2,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(u(i,jm1,k,2,tau) + 
     &                                      u(i,jp2,k,2,tau)))
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)
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 = maskV(i,j,km1)*maskV(i,j,kp2)
c         totvel = 0.5*(u(i,j+1,k,3,tau)+u(i,j,k,3,tau))
         totvel = (u(i,j+1,k,3,tau)*maskW(i,j+1,k)
     &             +u(i,j,k,3,tau)*maskW(i,j,k))/2.0
         adv_ft(i,j,k) = 0.5*totvel*(
     &         (a2nd*(1.0-mask) + a4th*mask)*(u(i,j,  k,2,tau) + 
     &                                        u(i,j,k+1,2,tau))+
     &         (b2nd*(1.0-mask) + b4th*mask)*(u(i,j,km1,2,tau) + 
     &                                        u(i,j,kp2,2,tau)))
        enddo
       enddo
      enddo
      end subroutine adv_flux_v_4th




      subroutine adv_flux_w_4th(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of vertical momentum 4th order
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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)
      real,parameter  :: a2nd = 1.0,     b2nd = 0.0
      real,parameter  :: a4th = 7.0/6.0, b4th = -1.0/6.0
      real :: mask,totvel

      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 = maskW(im1,j,k)*maskW(ip2,j,k)
         totvel = (u(i,j,k+1,1,tau)*maskU(i,j,k+1)
     &             +u(i,j,k,1,tau)*maskU(i,j,k))/2.0
         adv_fe(i,j,k) = 0.5*totvel*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(u(i,  j,k,3,tau) + 
     &                                      u(i+1,j,k,3,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(u(im1,j,k,3,tau) + 
     &                                      u(ip2,j,k,3,tau)))
        enddo
       enddo
      enddo
      call setcyclic3D(adv_fe)

      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 = maskW(i,jm1,k)*maskW(i,jp2,k)
         totvel = (u(i,j,k+1,2,tau)*maskV(i,j,k+1)
     &             +u(i,j,k,2,tau)*maskV(i,j,k))/2.0
         adv_fn(i,j,k) = 0.5*totvel*(
     &       (a2nd*(1.0-mask) + a4th*mask)*(u(i,  j,k,3,tau) + 
     &                                      u(i,j+1,k,3,tau))+
     &       (b2nd*(1.0-mask) + b4th*mask)*(u(i,jm1,k,3,tau) + 
     &                                      u(i,jp2,k,3,tau)))
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)
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 = maskW(i,j,km1)*maskW(i,j,kp2)
         totvel = (u(i,j,k+1,3,tau)*maskW(i,j,k+1)
     &             +u(i,j,k,3,tau)*maskW(i,j,k))/2.0
         adv_ft(i,j,k) = 0.5*totvel*(
     &         (a2nd*(1.0-mask) + a4th*mask)*(u(i,j,  k,3,tau) + 
     &                                        u(i,j,k+1,3,tau))+
     &         (b2nd*(1.0-mask) + b4th*mask)*(u(i,j,km1,3,tau) + 
     &                                        u(i,j,kp2,3,tau)))
        enddo
       enddo
      enddo
      end subroutine adv_flux_w_4th






      subroutine adv_flux_u_quicker(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of zonal momentum,   3th order upwind
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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)
      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 = 0.5*(u(i+1,j,k,1,tau)+u(i,j,k,1,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskU(i,j,k)*maskU(i+1,j,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskU(i+1,j,k)*maskU(i,j,k)
         eastmsk=maskU(i,j,k)*(1-maskU(im1,j,k))
         westmsk=maskU(i+1,j,k)*(1-maskU(ip2,j,k))
	 adv_fe(i,j,k) =
     &        0.5*totvel*(u(i,j,k,1,tau)+u(i+1,j,k,1,tau))
     &        - 0.5*upos*(0.25*u(i+1,j,k,1,taum1)-0.5*u(i,j,k,1,taum1)
     &               +0.25*(u(im1,j,k,1,taum1)*(1-eastmsk)
     &                          +u(i,j,k,1,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(ip2,j,k,1,taum1)*(1-westmsk)+
     &          u(i+1,j,k,1,taum1)*westmsk)-0.5*u(i+1,j,k,1,taum1)
     &                         +0.25*u(i,j,k,1,taum1))
          enddo
        enddo
      enddo
      call setcyclic3D(adv_fe)

      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 = 0.5*(u(i+1,j,k,2,tau)+u(i,j,k,2,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskU(i,j,k)*maskU(i,j+1,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskU(i,j+1,k)*maskU(i,j,k)
         eastmsk=maskU(i,j,k)*(1-maskU(i,jm1,k))
         westmsk=maskU(i,j+1,k)*(1-maskU(i,jp2,k))
	 adv_fn(i,j,k) =
     &        0.5*totvel*(u(i,j,k,1,tau)+u(i,j+1,k,1,tau))
     &        - 0.5*upos*(0.25*u(i,j+1,k,1,taum1)-0.5*u(i,j,k,1,taum1)
     &               +0.25*(u(i,jm1,k,1,taum1)*(1-eastmsk)
     &                          +u(i,j,k,1,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(i,jp2,k,1,taum1)*(1-westmsk)+
     &          u(i,j+1,k,1,taum1)*westmsk)-0.5*u(i,j+1,k,1,taum1)
     &                         +0.25*u(i,j,k,1,taum1))
          enddo
        enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)

      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 = 0.5*(u(i+1,j,k,3,tau)+u(i,j,k,3,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskU(i,j,k)*maskU(i,j,k+1)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskU(i,j,k+1)*maskU(i,j,k)
         eastmsk=maskU(i,j,k)*(1-maskU(i,j,km1))
         westmsk=maskU(i,j,k+1)*(1-maskU(i,j,kp2))
	 adv_ft(i,j,k) =
     &        0.5*totvel*(u(i,j,k,1,tau)+u(i,j,k+1,1,tau))
     &        - 0.5*upos*(0.25*u(i,j,k+1,1,taum1)-0.5*u(i,j,k,1,taum1)
     &               +0.25*(u(i,j,km1,1,taum1)*(1-eastmsk)
     &                          +u(i,j,k,1,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(i,j,kp2,1,taum1)*(1-westmsk)+
     &          u(i,j,k+1,1,taum1)*westmsk)-0.5*u(i,j,k+1,1,taum1)
     &                         +0.25*u(i,j,k,1,taum1))
          enddo
        enddo
      enddo
      end subroutine adv_flux_u_quicker





      subroutine adv_flux_v_quicker(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of meridional momentum,   3th order upwind
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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)
      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 = 0.5*(u(i,j+1,k,1,tau)+u(i,j,k,1,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskV(i,j,k)*maskV(i+1,j,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskV(i+1,j,k)*maskV(i,j,k)
         eastmsk=maskV(i,j,k)*(1-maskV(im1,j,k))
         westmsk=maskV(i+1,j,k)*(1-maskV(ip2,j,k))
	 adv_fe(i,j,k) =
     &        0.5*totvel*(u(i,j,k,2,tau)+u(i+1,j,k,2,tau))
     &        - 0.5*upos*(0.25*u(i+1,j,k,2,taum1)-0.5*u(i,j,k,2,taum1)
     &               +0.25*(u(im1,j,k,2,taum1)*(1-eastmsk)
     &                          +u(i,j,k,2,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(ip2,j,k,2,taum1)*(1-westmsk)+
     &          u(i+1,j,k,2,taum1)*westmsk)-0.5*u(i+1,j,k,2,taum1)
     &                         +0.25*u(i,j,k,2,taum1))
          enddo
        enddo
      enddo
      call setcyclic3D(adv_fe)

      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 = 0.5*(u(i,j+1,k,2,tau)+u(i,j,k,2,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskV(i,j,k)*maskV(i,j+1,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskV(i,j+1,k)*maskV(i,j,k)
         eastmsk=maskV(i,j,k)*(1-maskV(i,jm1,k))
         westmsk=maskV(i,j+1,k)*(1-maskV(i,jp2,k))
	 adv_fn(i,j,k) =
     &        0.5*totvel*(u(i,j,k,2,tau)+u(i,j+1,k,2,tau))
     &        - 0.5*upos*(0.25*u(i,j+1,k,2,taum1)-0.5*u(i,j,k,2,taum1)
     &               +0.25*(u(i,jm1,k,2,taum1)*(1-eastmsk)
     &                          +u(i,j,k,2,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(i,jp2,k,2,taum1)*(1-westmsk)+
     &          u(i,j+1,k,2,taum1)*westmsk)-0.5*u(i,j+1,k,2,taum1)
     &                         +0.25*u(i,j,k,2,taum1))
          enddo
        enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)

      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 = 0.5*(u(i,j+1,k,3,tau)+u(i,j,k,3,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskV(i,j,k)*maskV(i,j,k+1)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskV(i,j,k+1)*maskV(i,j,k)
         eastmsk=maskV(i,j,k)*(1-maskV(i,j,km1))
         westmsk=maskV(i,j,k+1)*(1-maskV(i,j,kp2))
	 adv_ft(i,j,k) =
     &        0.5*totvel*(u(i,j,k,2,tau)+u(i,j,k+1,2,tau))
     &        - 0.5*upos*(0.25*u(i,j,k+1,2,taum1)-0.5*u(i,j,k,2,taum1)
     &               +0.25*(u(i,j,km1,2,taum1)*(1-eastmsk)
     &                          +u(i,j,k,2,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(i,j,kp2,2,taum1)*(1-westmsk)+
     &          u(i,j,k+1,2,taum1)*westmsk)-0.5*u(i,j,k+1,2,taum1)
     &                         +0.25*u(i,j,k,2,taum1))
          enddo
        enddo
      enddo
      end subroutine adv_flux_v_quicker



      subroutine adv_flux_w_quicker(adv_fe,adv_fn,adv_ft)
c=======================================================================
c       advection of vertical momentum,   3th order upwind
c       fluxes are stored in adv_fe,adv_fn and adv_ft
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)
      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 = 0.5*(u(i,j,k+1,1,tau)+u(i,j,k,1,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskW(i,j,k)*maskW(i+1,j,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskW(i+1,j,k)*maskW(i,j,k)
         eastmsk=maskW(i,j,k)*(1-maskW(im1,j,k))
         westmsk=maskW(i+1,j,k)*(1-maskW(ip2,j,k))
	 adv_fe(i,j,k) =
     &        0.5*totvel*(u(i,j,k,3,tau)+u(i+1,j,k,3,tau))
     &        - 0.5*upos*(0.25*u(i+1,j,k,3,taum1)-0.5*u(i,j,k,3,taum1)
     &               +0.25*(u(im1,j,k,3,taum1)*(1-eastmsk)
     &                          +u(i,j,k,3,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(ip2,j,k,3,taum1)*(1-westmsk)+
     &          u(i+1,j,k,3,taum1)*westmsk)-0.5*u(i+1,j,k,3,taum1)
     &                         +0.25*u(i,j,k,3,taum1))
          enddo
        enddo
      enddo
      call setcyclic3D(adv_fe)

      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 = 0.5*(u(i,j,k+1,2,tau)+u(i,j,k,2,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskW(i,j,k)*maskW(i,j+1,k)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskW(i,j+1,k)*maskW(i,j,k)
         eastmsk=maskW(i,j,k)*(1-maskW(i,jm1,k))
         westmsk=maskW(i,j+1,k)*(1-maskW(i,jp2,k))
	 adv_fn(i,j,k) =
     &        0.5*totvel*(u(i,j,k,3,tau)+u(i,j+1,k,3,tau))
     &        - 0.5*upos*(0.25*u(i,j+1,k,3,taum1)-0.5*u(i,j,k,3,taum1)
     &               +0.25*(u(i,jm1,k,3,taum1)*(1-eastmsk)
     &                          +u(i,j,k,3,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(i,jp2,k,3,taum1)*(1-westmsk)+
     &          u(i,j+1,k,3,taum1)*westmsk)-0.5*u(i,j+1,k,3,taum1)
     &                         +0.25*u(i,j,k,3,taum1))
          enddo
        enddo
      enddo
      call border_exchg3D(adv_fn,1); call setcyclic3D(adv_fn)

      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 = 0.5*(u(i,j,k+1,3,tau)+u(i,j,k,3,tau))
         upos = .5*(totvel + abs(totvel))
     &                  *maskW(i,j,k)*maskW(i,j,k+1)
         uneg = .5*(totvel - abs(totvel))
     &                  *maskW(i,j,k+1)*maskW(i,j,k)
         eastmsk=maskW(i,j,k)*(1-maskW(i,j,km1))
         westmsk=maskW(i,j,k+1)*(1-maskW(i,j,kp2))
	 adv_ft(i,j,k) =
     &        0.5*totvel*(u(i,j,k,3,tau)+u(i,j,k+1,3,tau))
     &        - 0.5*upos*(0.25*u(i,j,k+1,3,taum1)-0.5*u(i,j,k,3,taum1)
     &               +0.25*(u(i,j,km1,3,taum1)*(1-eastmsk)
     &                          +u(i,j,k,3,taum1)*eastmsk))
     &        - 0.5*uneg*(0.25*(u(i,j,kp2,3,taum1)*(1-westmsk)+
     &          u(i,j,k+1,3,taum1)*westmsk)-0.5*u(i,j,k+1,3,taum1)
     &                         +0.25*u(i,j,k,3,taum1))
          enddo
        enddo
      enddo
      end subroutine adv_flux_w_quicker






