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: !!----------------------------------------------------------
  module pbl_simple
      use data_structures
      use domain_interface,   only : domain_t
      use options_interface,  only : options_t
  
  
negative.code:         ! don't process the top model level regardless, there shouldn't be any real pbl diffusion occuring there
          kte = min(kme-1, kte_in)
  !       OpenMP parallelization small static chunk size because we typically get a small area that takes most of the time (because of substepping)
          !$omp parallel shared(land_mask, th, qv, cloud, ice, qrain, qsnow, um, vm, pii, rho, z, dz, terrain) & !, tend_qv_pbl)     &
          !$omp shared(l_m, K_m, Kq_m, stability_m, prandtl_m, virt_pot_temp_zgradient_m, rig_m, shear_m, lastqv_m, ims, ime, jms, jme, kms, kme) &
          !$omp firstprivate(its, ite, jts, jte, kts, kte, dt) private(i, k, j)
  
  