anchor.diff:  	private
   	public :: simple_pbl, finalize_simple_pbl, init_simple_pbl
   	
  -! 	note _m indicates module level variables
  +! 	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
  -	real, allocatable, dimension(:,:,:) :: rig_m,shear_m
  -	real, allocatable, dimension(:,:,:) :: stability_m,prandtl_m,l_m,K_m,Kq_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
  
positive.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
  		do t=1,nsubsteps
  			! First water vapor
  			call diffuse_variable(domain%qv,rhomean,rho_dz,j)
  
negative.code: 			call diffuse_variable(domain%ice,rhomean,rho_dz,j)
  			! and snow
  			call diffuse_variable(domain%qsnow,rhomean,rho_dz,j)
  			! don't bother with rain or graupel assuming they are falling fast *enough* not entirely fair...
  		enddo
  
  	end subroutine pbl_diffusion
  	
  