anchor.diff:  			allocate(U_m(nx-1,nz,ny))
   			allocate(V_m(nx,nz,ny-1))
   			allocate(W_m(nx,nz,ny))
  +			allocate(rho_m(nx,nz,ny))
  +			allocate(pii_m(nx,nz,ny))
   		endif
  +! 		rho_m(:,:,:)=1
  +        !$omp parallel firstprivate(ny) &
  +		!$omp private(i) shared(pii_m,rho_m,domain)
  +        !$omp do schedule(static)
  +		do i=1,ny
  +			pii_m(:,:,i)=(domain%p(:,:,i)/100000.0)**(R/cp)
  +	        rho_m(:,:,i)=domain%p(:,:,i)/(R*domain%th(:,:,i)*pii_m(:,:,i)) ! kg/m^3
  +		enddo
  +        !$omp end do
  +        !$omp end parallel
   		
   ! 		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
  +		U_m=domain%u(1:nx-1,:,:)*dt/dx! * (rho_m(1:nx-1,:,:)+rho_m(2:nx,:,:))/2.0
  +		V_m=domain%v(:,:,1:ny-1)*dt/dx! * (rho_m(:,:,1:ny-1)+rho_m(:,:,2:ny))/2.0
   ! 		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
  +! 		W_m(:,:nz-1,:) = W_m(:,:nz-1,:) * (rho_m(:,1:nz-1,:)+rho_m(:,2:nz,:))/2.0
  +! 		W_m(:,nz,:) = W_m(:,nz,:) * rho_m(:,nz,:)
   ! 		should probably be converting to mass (q*rho) before advecting, then back again... but testing showed minimal difference
  +		domain%qv=domain%qv*rho_m
  +		call advect3d(domain%qv,   U_m,V_m,W_m,nx,nz,ny,0)
  +		domain%qv=domain%qv/rho_m
  +		
  +		domain%cloud=domain%cloud*rho_m
  +		call advect3d(domain%cloud,U_m,V_m,W_m,nx,nz,ny,0)
  +		domain%cloud=domain%cloud/rho_m
  +		
  +		domain%ice=domain%ice*rho_m
  +		call advect3d(domain%ice,  U_m,V_m,W_m,nx,nz,ny,0)
  +		domain%ice=domain%ice/rho_m
  +		
  +		domain%qrain=domain%qrain*rho_m
  +		call advect3d(domain%qrain,U_m,V_m,W_m,nx,nz,ny,0)
  +		domain%qrain=domain%qrain/rho_m
  +		
  +		domain%qsnow=domain%qsnow*rho_m
  +		call advect3d(domain%qsnow,U_m,V_m,W_m,nx,nz,ny,0)
  +		domain%qsnow=domain%qsnow/rho_m
  +		
  +		domain%qgrau=domain%qgrau*rho_m
  +		call advect3d(domain%qgrau,U_m,V_m,W_m,nx,nz,ny,0)
  +		domain%qgrau=domain%qgrau/rho_m
  +		
   		call advect3d(domain%th,   U_m,V_m,W_m,nx,nz,ny,0)
  -		call advect3d(domain%qv,   U_m,V_m,W_m,nx,nz,ny,0)
  -		call advect3d(domain%cloud,U_m,V_m,W_m,nx,nz,ny,0)
  -		call advect3d(domain%ice,  U_m,V_m,W_m,nx,nz,ny,0)
   		call advect3d(domain%nice, U_m,V_m,W_m,nx,nz,ny,0)
  -		call advect3d(domain%qrain,U_m,V_m,W_m,nx,nz,ny,0)
   		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)
  -		
   	end subroutine advect
   
   
  
positive.code: 		real,intent(in)::dt
  		
  		real::dx
  		integer::nx,nz,ny
  			
  		dx=domain%dx
  		nx=size(domain%dz,1)
  
negative.code: !            call flux2(qin(2:nx-1,1:nz-1,i),qin(2:nx-1,2:nz,i),w(2:nx-1,1:nz-1,i),nx-2,nz-1,1,f5)
  		   
             ! perform horizontal advection
             q(2:nx-1,:,i)=q(2:nx-1,:,i) - ((f1(2:nx-1,:)-f1(1:nx-2,:)) + (f3(:,:)-f4(:,:)))!/rho(2:nx-1,:,i)
             ! then vertical (order doesn't matter because fluxes f1-6 are calculated before applying them)
             ! add fluxes to middle layers
             q(2:nx-1,2:nz-1,i)=q(2:nx-1,2:nz-1,i)-(f5(:,2:nz-1)-f5(:,1:nz-2))!/rho(2:nx-1,2:nz-1,i)
             ! add fluxes to bottom layer
             q(2:nx-1,1,i)=q(2:nx-1,1,i)-f5(:,1)!/rho(2:nx-1,1,i)
             ! add fluxes to top layer
             q(2:nx-1,nz,i)=q(2:nx-1,nz,i)-(qin(2:nx-1,nz,i)*w(2:nx-1,nz,i)-f5(:,nz-1))!/rho(2:nx-1,nz,i)
          enddo
          !$omp end do
          !$omp end parallel
  