anchor.diff:  		integer,intent(in) :: j
   		real,dimension(nx,nz)::fluxes
   		
  +! 		Eventually this should be made into an implicit solution to avoid substepping
   		! 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)
  +		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
  +		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
   	
  
positive.code: 			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...
  
  	end subroutine pbl_diffusion
  	
  
negative.code: 	private
  	public :: simple_pbl, finalize_simple_pbl, init_simple_pbl
  	
  ! 	NOTE *_m indicates module level variables
  ! 	these variables are declared as module level variables so that they do not need to be allocated
  !   deallocated, and re-allocated all the time. 
  
  ! 	gradient in the virtual potential temperature
  	real, allocatable, dimension(:,:,:) :: virt_pot_temp_zgradient_m
  ! 	gradient in the richardson number
  	real, allocatable, dimension(:,:,:) :: rig_m
  ! 	vertical wind shear (dwind / dz)
  	real, allocatable, dimension(:,:,:) :: shear_m
  ! 	atmospheric stability
  	real, allocatable, dimension(:,:,:) :: stability_m
  ! 	length scale that asymptotes from kappa*z to l0 (250m below)
  	real, allocatable, dimension(:,:,:) :: l_m
  ! 	diffusion term for momentum
  	real, allocatable, dimension(:,:,:) :: K_m
  ! 	diffusion term for scalars (K/prandtl)
  	real, allocatable, dimension(:,:,:) :: Kq_m
  ! 	prandtl number to convert K for momentum to K for scalars
  	real, allocatable, dimension(:,:,:) :: prandtl_m
  	integer :: nx,nz,ny !NOTE these are subset from full domain e.g. nx-2,ny-2,nz-1
  	
  !	limits on Pr noted in HP96 page 2325 below eqn 13
  