anchor.diff:  !           note, even though dz!=dx, W is computed from the divergence in U/V so it is scaled by dx/dz already
               W_m=domain%w*(dt/dx)
           endif
  -
  +        
  +        error=0
  +        
   !         print*, "qv"
  -        call advect3d(domain%qv,   U_m,V_m,W_m, domain%rho,domain%dz,nx,nz,ny,0,options)
  +        call advect3d(domain%qv,   U_m,V_m,W_m, domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +        if (error/=0) then
  +            print*, "qv", error
  +            if (error>0) then
  +                stop "MPDATA Error"
  +            endif
  +        endif
   !         print*, "qc"
  -        call advect3d(domain%cloud,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +        call advect3d(domain%cloud,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +        if (error/=0) then
  +            print*, "qc", error
  +            if (error>0) then
  +                stop "MPDATA Error"
  +            endif
  +        endif
   !         print*, "qr"
  -        call advect3d(domain%qrain,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +        call advect3d(domain%qrain,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +        if (error/=0) then
  +            print*, "qr", error
  +            if (error>0) then
  +                stop "MPDATA Error"
  +            endif
  +        endif
   !         print*, "qs"
  -        call advect3d(domain%qsnow,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +        call advect3d(domain%qsnow,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +        if (error/=0) then
  +            print*, "qs", error
  +            if (error>0) then
  +                stop "MPDATA Error"
  +            endif
  +        endif
   !         print*, "th"
  -        call advect3d(domain%th,   U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +        call advect3d(domain%th,   U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +        if (error/=0) then
  +            print*, "th", error
  +            if (error>0) then
  +                stop "MPDATA Error"
  +            endif
  +        endif
           if (options%physics%microphysics==kMP_THOMPSON) then
   !             print*, "qi"
  -            call advect3d(domain%ice,  U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +            call advect3d(domain%ice,  U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +            if (error/=0) then
  +                print*, "qi", error
  +                if (error>0) then
  +                    stop "MPDATA Error"
  +                endif
  +            endif
   !             print*, "qg"
  -            call advect3d(domain%qgrau,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +            call advect3d(domain%qgrau,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +            if (error/=0) then
  +                print*, "qg", error
  +                if (error>0) then
  +                    stop "MPDATA Error"
  +                endif
  +            endif
   !             print*, "ni"
  -            call advect3d(domain%nice, U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  +            call advect3d(domain%nice, U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +            if (error/=0) then
  +                print*, "ni", error
  +                if (error>0) then
  +                    stop "MPDATA Error"
  +                endif
  +            endif
   !             print*, "nr"
  -            call advect3d(domain%nrain,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  -        endif
  -!         print*, "CLOUD", minval(domain%cloud), maxval(domain%cloud)
  +            call advect3d(domain%nrain,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options, error)
  +            if (error/=0) then
  +                print*, "nr", error
  +                if (error>0) then
  +                    stop "MPDATA Error"
  +                endif
  +            endif
  +        endif
           order=mod(order+1,3)
       end subroutine mpdata
   end module adv_mpdata
  
positive.code:     end subroutine advect_u
  
  
      subroutine advect3d(q,u,v,w,rho,dz,nx,nz,ny,debug,options)
          implicit none
          real,dimension(1:nx,1:nz,1:ny), intent(inout) :: q
          real,dimension(1:nx,1:nz,1:ny), intent(in) :: w
  
negative.code:     integer :: order
      public:: mpdata, mpdata_init
      public:: advect_u, advect_w, advect_v ! for test_mpdata testing only!
      
  
  contains
      subroutine flux2(l,r,U,nx,nz,ny,f)
  