anchor.diff:                  ! diffusion for scalars
                   Kq_m(its:ite,k,j) = K_m(its:ite,k,j) / prandtl_m(its:ite,k,j)
   
  -                ! rescale diffusion to cut down on excessive mixing
  -                Kq_m(its:ite,k,j) = Kq_m(its:ite, k,j) / diffusion_reduction
                   Kq_m(its:ite,k,j) = Kq_m(its:ite, k,j) * dt / ((dz(its:ite,k,j) + dz(its:ite,k+1,j))/2)
   
                   ! enforce limits specified in HP96
  
positive.code:     ! note, they actually use 30m because they only use this for free-atmosphere mixing
      ! but they note that 250m is used in the operational model for the full PBL mixing
      real, parameter :: N_substeps=10. ! number of substeps to allow (puts a cap on K to match CFL)
      real, parameter :: diffusion_reduction=10.0 ! used to reduce diffusion rates
  
  
  contains
      subroutine simple_pbl(th, qv, cloud, ice, qrain, qsnow, um, vm, pii, rho, z, dz, terrain, its, ite, jts, jte, kts, kte_in, dt, tend_qv_pbl)
          real,   intent(inout), dimension(ims:ime, kms:kme, jms:jme) :: th            ! potential temperature [K]
          real,   intent(inout), dimension(ims:ime, kms:kme, jms:jme) :: qv            ! water vapor mixing ratio [kg/kg]
          real,   intent(inout), dimension(ims:ime, kms:kme, jms:jme) :: cloud         ! cloud water mixing ratio [kg/kg]
  
negative.code:         real,   intent(in),    dimension(ims:ime, kms:kme, jms:jme) :: z             ! model level heights [m]
          real,   intent(in),    dimension(ims:ime, kms:kme, jms:jme) :: dz            ! model level thickness [m]
          real,   intent(in),    dimension(ims:ime, jms:jme)          :: terrain       ! terrain height above sea level [m]
          integer,intent(in),    dimension(ims:ime, jms:jme)          :: land_mask      ! water = 2 (kLC_WATER) land = 1 (kLC_LAND)
          integer,intent(in) :: its, ite, jts, jte, kts, kte_in
          real,   intent(in) :: dt                                 ! time step [s]
          real,   intent(inout), dimension(ims:ime,kms:kme,jms:jme), optional :: tend_qv_pbl   ! output water vapor tendency [kg/kg/s] (for use in other physics)
  