anchor.diff:  			call diffuse_variable(domain%ice,rhomean,rho_dz,j)
   			! and snow
   			call diffuse_variable(domain%qsnow,rhomean,rho_dz,j)
  -		enddo
  -		! don't bother with rain or graupel assuming they are falling fast *enough* not entirely fair...
  +			! don't bother with rain or graupel assuming they are falling fast *enough* not entirely fair...
  +		enddo
   
   	end subroutine pbl_diffusion
   	
  
positive.code: 		integer,intent(in) :: j
  		real,dimension(nx,nz)::fluxes
  		
  		! if gradient is downward (qv[z+1]>qv[z]) then flux is negative
  		fluxes=Kq_m(:,:,j)*rhomean*(q(2:nx+1,:nz,j+1)-q(2:nx+1,2:,j+1))
  		! first layer assumes no flow through the surface, that comes from the LSM
  		q(2:nx+1,1,j+1) = q(2:nx+1,1,j+1) - &
  								fluxes(:,1) / rho_dz(:,1)
  		! middle layers (no change for top layer assuming flux in = flux out)
  		q(2:nx+1,2:nz,j+1) = q(2:nx+1,2:nz,j+1) - &
  		   					   (fluxes(:,2:nz)-fluxes(:,:nz-1)) / rho_dz
  		
  	end subroutine diffuse_variable
  	
  
negative.code: 		! for most regions it is < 0.5, for the small regions it isn't just substep for now. 
  		nsubsteps=ceiling(2*maxval(Kq_m(:,:,j)/domain%dz(2:nx+1,:nz,j)))
  		Kq_m(:,:,j)=Kq_m(:,:,j)/nsubsteps
  ! 		nsubsteps will typically be 1
  		do t=1,nsubsteps
  			! First water vapor
  			call diffuse_variable(domain%qv,rhomean,rho_dz,j)
  