anchor.diff:                  U2=u(:,k,j)
                   l =q(:,k,j)
                   call flux1(q1(1:n-1),q1(2:n),U2,f)
  -
  +                
                   include "adv_mpdata_FCT_core.f90"
                   u(:,k,j)=U2
               end do
  -
  +            
               n=nz
               q1=>q1z
               l =>lz
  
positive.code:                 U2=v(j,k,:)
                  l =q(j,k,:)
                  call flux1(q1(1:n-1),q1(2:n),U2,f)
  
                  include "adv_mpdata_FCT_core.f90"
                  v(j,k,:)=U2
              end do
          end do
          !$omp end do
          !$omp end parallel
  
      end subroutine flux_limiter
  
      subroutine advect3d(q,u,v,w,rho,dz,dx,nx,nz,ny,jaco,options,err)
          implicit none
          real,dimension(1:nx,1:nz,1:ny), intent(inout) :: q
          real,dimension(1:nx-1,1:nz,1:ny),intent(in) :: u
          real,dimension(1:nx,1:nz,1:ny-1),intent(in) :: v
          real,dimension(1:nx,1:nz,1:ny), intent(in) :: w
          real,dimension(1:nx,1:nz,1:ny), intent(in) :: rho
          real,dimension(1:nx,1:nz,1:ny), intent(in) :: dz
          integer, intent(in) :: ny,nz,nx
          real,dimension(1:nx,1:nz,1:ny), intent(in) :: jaco
          type(options_t), intent(in)::options
          integer, intent(inout) :: err
          real, intent(in) :: dx
  
          ! used for intermediate values in the mpdata calculation
          real,dimension(1:nx,1:nz,1:ny)   :: q2
          real,dimension(1:nx-1,1:nz,1:ny) :: u2
          real,dimension(1:nx,1:nz,1:ny-1) :: v2
          real,dimension(1:nx,1:nz,1:ny)   :: w2
  
          integer :: iord, i
  
          do iord=1,options%adv_options%mpdata_order
              if (iord==1) then
                  call upwind_advection(q, u, v, w, q2, dx,dz,nx,nz,ny,jaco)
              else
                  call mpdata_fluxes(q2, u, v, w, u2,v2,w2, nx,nz,ny)
                  if (this_image()==100) write(*,*) maxval(u2)
                  if ( (sum(abs(u2))+sum(abs(v2))+sum(abs(w2)) < 0.01)) write(*,*) "no ADV corr--1"
  
                  if (options%adv_options%flux_corrected_transport) then
                      call flux_limiter(q, q2, u2,v2,w2, nx,nz,ny)
                  endif
                  if ( (sum(abs(u2))+sum(abs(v2))+sum(abs(w2)) < 0.01)) write(*,*) "no ADV corr--2"
                  call upwind_advection(q2, u2, v2, w2, q, dx,dz,nx,nz,ny,jaco)
              endif
  
              !
              if (iord/=options%adv_options%mpdata_order) then
                  if (iord>1) then
                      !$omp parallel shared(q,q2) firstprivate(ny) private(i)
                      !$omp do schedule(static)
                      do i=1,ny
                          q2(:,:,i)=q(:,:,i)
                      enddo
                      !$omp end do
                      !$omp end parallel
                  endif
              else
                  if (iord==1) then
                      !$omp parallel shared(q,q2) firstprivate(ny) private(i)
                      !$omp do schedule(static)
                      do i=1,ny
                          q(:,:,i)=q2(:,:,i)
                      enddo
                      !$omp end do
                      !$omp end parallel
                  endif
  
              endif
          end do
  
  
      end subroutine advect3d
  
      subroutine mpdata_init(domain,options)
          type(domain_t), intent(in) :: domain
          type(options_t), intent(in) :: options
  
          ! originally used to permit the order of dimensions in advection to be rotated
          order    = 0
      end subroutine mpdata_init
  
  !   primary entry point, advect all scalars in domain
      subroutine mpdata(domain,options,dt)
          implicit none
          type(domain_t),intent(inout)::domain
          type(options_t), intent(in)::options
          real,intent(in)::dt
  
          real::dx
          integer::nx,nz,ny,i, error
          integer :: ims, ime, jms, jme, kms, kme
  
          ims = domain%grid%ims
          ime = domain%grid%ime
  
negative.code:         real, dimension(nz),   target :: q1z,lz
          real, dimension(nz-1), target :: fz, U2z
          logical :: flux_is_w
          
          real :: qmax_i,qmin_i,qmax_i2,qmin_i2
          real :: beta_in_i, beta_out_i, beta_in_i2, beta_out_i2
          real :: fin_i, fout_i, fin_i2, fout_i2
          
          ! NOTE: before inclusion of FCT_core the following variables must be setup: 
          ! q1 and l (l=q0)
          !$omp parallel shared(q2,q,u,v,w) firstprivate(nx,ny,nz) default(private)
          !$omp do schedule(static)
  