anchor.diff:  	    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
  +! 	    real,dimension(1:nx,1:nz,1:ny), intent(in) :: rho
   		integer, intent(in) :: ny,nz,nx,debug
           ! interal parameters
           integer :: err,i
  
positive.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(:,:))
             ! 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))
             ! add fluxes to bottom layer
             q(2:nx-1,1,i)=q(2:nx-1,1,i)-f5(:,1)
             ! 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))
          enddo
          !$omp end do
          !$omp end parallel
  
negative.code: 		real,intent(in)::dt
  		
  		real::dx
  		integer::nx,nz,ny,i
  			
  		dx=domain%dx
  		nx=size(domain%dz,1)
  