anchor.diff:          kme = domain%grid%kme
   
           dx=domain%dx
  -        nx = domain%grid%nx
  -        nz = domain%grid%nz
  -        ny = domain%grid%ny
   
           if (.not.allocated(domain%advection_dz)) then
               allocate(domain%advection_dz(ims:ime,kms:kme,jms:jme))
  
positive.code:         real,dimension(1:nx,1:nz,1:ny-1),intent(inout) :: v
          real,dimension(1:nx,1:nz,1:ny),  intent(inout) :: w
          integer, intent(in) :: nx,nz,ny
  
          integer :: i,j,k,n
          real, dimension(:), pointer :: q1, U2, l, f
          ! q1 = q after applying previous iteration advection
  
negative.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,rho,ims,ime,kms,kme,jms,jme,jaco,dz,options)
          implicit none
          real,dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: q
          real,dimension(ims:ime,kms:kme,jms:jme), intent(in) :: rho
          integer, intent(in) :: ims,ime,kms,kme,jms,jme
          real,dimension(ims:ime,kms:kme,jms:jme), intent(in) :: jaco,dz
          type(options_t), intent(in)::options
  
          ! used for intermediate values in the mpdata calculation
          real,dimension(ims:ime,kms:kme,jms:jme)   :: q2
          real,dimension(ims+1:ime,kms:kme,jms:jme) :: u2
          real,dimension(ims:ime,kms:kme,jms+1:jme) :: v2
          real,dimension(ims:ime,kms:kme,jms:jme)   :: w2
          
          integer :: iord, i
                  
          do iord=1,options%adv_options%mpdata_order
              if (iord==1) then
                  call upwind_advection(q, U_m, V_m, W_m, rho, q2,dz,ims,ime,kms,kme,jms,jme,jaco)
              else
                  ! Due to an uneven vertical grid spacing, W_m is not normalized by dz, because this causes problems in
                  ! advection code. However, pseudo-velocity expects, all velocities normalized by time and dx/z
                  ! so do that here before passing to pseudo-velocity calculations
                  call mpdata_fluxes(q2, U_m, V_m, (W_m/dz), u2,v2,w2, ims,ime,kms,kme,jms,jme,(jaco*rho))
                  ! and un-normalize, since upwind advection scheme includes dz
                  ! Since pseudo-velocities cannot be gaurenteed to be non-divergent, we assume worst-case and multiply by 0.5 to
                  ! ensure stability (from Smolarkiewicz 1984, after Eq. 24)
                  call upwind_advection(q2, 0.5*u2, 0.5*v2, 0.5*(w2*dz), rho, q,dz,ims,ime,kms,kme,jms,jme,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=jms,jme
                          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=jms,jme
                          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
      
      
      subroutine test_divergence(dz, ims, ime, kms, kme, jms, jme)
          implicit none
          real, intent(in) :: dz(ims:ime,kms:kme,jms:jme)
          integer, intent(in) :: ims, ime, jms, jme, kms, kme
  
          real, allocatable :: du(:,:), dv(:,:), dw(:,:)
          integer :: i,j,k
  
          allocate(du(ims+1:ime-1,jms+1:jme-1))
          allocate(dv(ims+1:ime-1,jms+1:jme-1))
          allocate(dw(ims+1:ime-1,jms+1:jme-1))
  
          do i=ims+1,ime-1
              do j=jms+1,jme-1
                  do k=kms,kme
                      du(i,j) = (U_m(i+1,k,j)-U_m(i,k,j))
                      dv(i,j) = (V_m(i,k,j+1)-V_m(i,k,j))
                      if (k==kms) then
                          dw(i,j) = (W_m(i,k,j))/dz(i,k,j)
                      else
                          dw(i,j) = (W_m(i,k,j)-W_m(i,k-1,j))/dz(i,k,j)
                      endif
                      if (abs(du(i,j) + dv(i,j) + dw(i,j)) > 1e-3) then
                          print*, this_image(), i,k,j , abs(du(i,j) + dv(i,j) + dw(i,j))
                          print*, "Winds are not balanced on entry to advect"
                          !error stop
                      endif
                  enddo
              enddo
          enddo
  
      end subroutine test_divergence
      
  !   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, allocatable, dimension(:,:,:)   :: rho
          real::dx
          integer :: i, ims, ime, jms, jme, kms, kme
  
          ims = domain%grid%ims
          ime = domain%grid%ime
  