anchor.diff:  !            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(:,:))
  +           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))
  +           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)
  +           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))
  +           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
  
positive.code: 	    real,dimension(1:nx,1:nz,1:ny), intent(in) :: w
          real,dimension(1:nx-1,1:nz,1:ny),intent(in) :: u
          real,dimension(1:nx,1:nz,1:ny-1),intent(in) :: v
  		integer, intent(in) :: ny,nz,nx,debug
          ! interal parameters
          integer :: err,i
  
negative.code: 			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! * (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%nice, U_m,V_m,W_m,nx,nz,ny,0)
  		call advect3d(domain%nrain,U_m,V_m,W_m,nx,nz,ny,0)
  	end subroutine advect
  
  
  