anchor.diff:  		call advect3d(domain%nrain,U_m,V_m,W_m,nx,nz,ny,0)
   		call advect3d(domain%qsnow,U_m,V_m,W_m,nx,nz,ny,0)
   		call advect3d(domain%qgrau,U_m,V_m,W_m,nx,nz,ny,0)
  -		
  -! 		deallocate(U,V,W)
   		
   	end subroutine advect
   
  
positive.code:         ! interal parameters
          integer :: err,i
  	    real,dimension(1:nx,1:nz,1:ny) :: qin
          real, dimension(1:nx-1,1:nz,1) :: f1 ! there used to be an f2 to store f[x+1]
          real, dimension(1:nx-2,1:nz,1) :: f3,f4
          real, dimension(1:nx-2,1:nz-1,1) ::f5
          !$omp parallel shared(q,u,v,w) firstprivate(nx,ny,nz) private(i,f1,f3,f4,f5)
          !$omp do
          do i=1,ny
  
negative.code: 		nz=size(domain%dz,2)
  		ny=size(domain%dz,3)
  		
  ! 		if this if the first time we are called, we need to allocate the module level arrays
  		if (.not.allocated(U_m)) then
  			allocate(U_m(nx-1,nz,ny))
  			allocate(V_m(nx,nz,ny-1))
  			allocate(W_m(nx,nz,ny))
  		endif
  		
  ! 		calculate U,V,W normalized for dt/dx
  		U_m=domain%u(1:nx-1,:,:)*dt/dx
  		V_m=domain%v(:,:,1:ny-1)*dt/dx
  ! 		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
  ! 		should probably be converting to mass (q*rho) before advecting, then back again... but testing showed minimal difference
  