anchor.diff:          real,intent(in)::dt
           
           real::dx
  -        integer::nx,nz,ny,i
  -            
  +        integer::nx,nz,ny,i, error
  +        
           dx=domain%dx
           nx=size(domain%dz,1)
           nz=size(domain%dz,2)
  
positive.code: !           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
  
  !         print*, "qv"
          call advect3d(domain%qv,   U_m,V_m,W_m, domain%rho,domain%dz,nx,nz,ny,0,options)
  !         print*, "qc"
          call advect3d(domain%cloud,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  !         print*, "qr"
          call advect3d(domain%qrain,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  !         print*, "qs"
          call advect3d(domain%qsnow,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  !         print*, "th"
          call advect3d(domain%th,   U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
          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)
  !             print*, "qg"
              call advect3d(domain%qgrau,U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  !             print*, "ni"
              call advect3d(domain%nice, U_m,V_m,W_m,domain%rho,domain%dz,nx,nz,ny,0,options)
  !             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)
          order=mod(order+1,3)
      end subroutine mpdata
  end module adv_mpdata
  
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)
  