ERF_TurbStruct.H

Go to the documentation of this file.

#ifndef ERF_TURB_STRUCT_H_
#define ERF_TURB_STRUCT_H_
 
#include <AMReX_ParallelDescriptor.H>
#include <ERF_MYNNStruct.H>
#include <ERF_RANSClosure.H>
 
#include <iomanip>
#include <sstream>
 
/**
 * @brief Large-eddy simulation closure type.
 */
AMREX_ENUM(LESType, None, Smagorinsky, Smagorinsky2D, Deardorff);
 
/**
 * @brief Reynolds-averaged turbulence closure type.
 */
AMREX_ENUM(RANSType, None, kEqn);
 
/**
 * @brief Planetary boundary-layer closure type.
 */
AMREX_ENUM(PBLType, None, MYJ, MYNN25, MYNNEDMF, YSU, YSUNew, MRF, SHOC, EAMXX_SHOC, NATIVE_SHOC);
 
/**
 * @brief Thermodynamic variable used for stability stratification.
 */
AMREX_ENUM(StratType, theta, thetav, thetal);
 
/**
 * @brief Query a scalar or per-level input value.
 * @tparam T Value type accepted by ParmParse::query.
 * @param pp ParmParse object to query.
 * @param query_string Input key to query.
 * @param query_var Destination value for the requested level.
 * @param lev AMR level index.
 * @param maxlev Maximum AMR level configured for the run.
 * @return ParmParse query result, or zero if the key is absent.
 */
template <typename T>
int
query_one_or_per_level (
  const amrex::ParmParse& pp,
  const char* query_string,
  T& query_var,
  const int lev,
  const int maxlev)
{
  int count = pp.countval(query_string);
  if (count == 0) {
    return 0; // nothing to do
  } else if (count == 1) {
    // In this case, we assume the same value is used at every level
    return pp.query(query_string, query_var);
  } else if (count >= maxlev + 1) {
    // In this case, there may be more values than levels so we
    // will just read the number of values we need and ignore the
    // extra levels
    return pp.query(query_string, query_var, lev);
  } else {
    // In this case, there is more than one value AND there
    // are more levels than values; as an example, if max_level = 2
    // but count = 2 so we have three levels but only two values --
    // it is not clear how to interpret this so we abort
    amrex::Error(
      "For parmparse variable " + pp.prefixedName(query_string) +
      ": if specified, specify once total or at least once for each level");
    return 0; // avoid compiler warning
  }
}
 
/**
 * @brief Query a scalar or per-level enum input value using case-insensitive matching.
 * @tparam T Enum type accepted by ParmParse::query_enum_case_insensitive.
 * @param pp ParmParse object to query.
 * @param query_string Input key to query.
 * @param query_var Destination enum value for the requested level.
 * @param lev AMR level index.
 * @param maxlev Maximum AMR level configured for the run.
 * @return ParmParse query result, or zero if the key is absent.
 */
template <typename T>
int
query_one_or_per_level_enum_case_insensitive (
  const amrex::ParmParse& pp,
  const char* query_string,
  T& query_var,
  const int lev,
  const int maxlev)
{
  int count = pp.countval(query_string);
  if (count == 0) {
    return 0; // nothing to do
  } else if (count == 1) {
    // In this case, we assume the same value is used at every level
    return pp.query_enum_case_insensitive(query_string, query_var);
  } else if (count >= maxlev + 1) {
    // In this case, there may be more values than levels so we
    // will just read the number of values we need and ignore the
    // extra levels
    return pp.query_enum_case_insensitive(query_string, query_var, lev);
  } else {
    // In this case, there is more than one value AND there
    // are more levels than values; as an example, if max_level = 2
    // but count = 2 so we have three levels but only two values --
    // it is not clear how to interpret this so we abort
    amrex::Error(
      "For parmparse variable " + pp.prefixedName(query_string) +
      ": if specified, specify once total or at least once for each level");
    return 0; // avoid compiler warning
  }
}
 
/**
 * Container holding quantities related to turbulence parametrizations
 */
struct TurbChoice
{
public:
  /**
   * @brief Read turbulence options for one AMR level from the input parameter database.
   * @param lev AMR level index.
   * @param max_level Maximum AMR level configured for the run.
   * @param pp_prefix ParmParse prefix for the ERF input namespace.
   */
  void init_params (int lev, int max_level, std::string pp_prefix)
  {
    amrex::ParmParse pp(pp_prefix);
 
    // Which LES closure?
    query_one_or_per_level(pp, "les_type", les_type, lev, max_level);
    // Immersed-boundary awareness of the MRF and YSUNew schemes (off by default):
    // the surface of a column with immersed solid cells is the first fluid cell
    // above them, heights and the boundary-layer depth are measured from it,
    // the eddy diffusivities are zero inside the solid, and the surface
    // scales over such a column come from a neutral log law at its top with
    // the roughness pbl_ib_z0. Only MRF and YSUNew honour it.
    query_one_or_per_level(pp, "pbl_ib_aware", pbl_ib_aware, lev, max_level);
    query_one_or_per_level(pp, "pbl_ib_z0",    pbl_ib_z0,    lev, max_level);
 
    // Handle 2-D Smag
    if (les_type == LESType::Smagorinsky2D) {
        les_type = LESType::Smagorinsky;
        smag2d = true;
    }
 
    // Which RANS closure?
    query_one_or_per_level(pp, "rans_type", rans_type, lev, max_level);
 
    if ((rans_type != RANSType::None) && (les_type != LESType::None)) {
        amrex::Error("Hybrid RANS-LES not implemented");
    }
 
    // Which PBL Closure
    query_one_or_per_level_enum_case_insensitive(pp, "pbl_type", pbl_type, lev, max_level);
    if (pbl_type == PBLType::SHOC) {
      static bool warned_legacy_shoc = false;
      if (!warned_legacy_shoc) {
        amrex::Warning("erf.pbl_type = SHOC is deprecated; use erf.pbl_type = EAMXX_SHOC");
        warned_legacy_shoc = true;
      }
      pbl_type = PBLType::EAMXX_SHOC;
    }
 
    // Do some more stuff for PBL Modeling
    if (pbl_type != PBLType::None) {
      // Check for compatibility between PBL, LES, Molec Transport
      if (les_type != LESType::None) {
        amrex::Print() << "Selected a PBL model and an LES model: "
                       << "Using PBL for vertical transport, LES for horizontal"
                       << std::endl;
      }
      if (les_type == LESType::Smagorinsky) {
        if (!smag2d)
          amrex::Error("If using Smagorinsky with a PBL model, the 2-D "
                       "formulation should be used");
      } else if (les_type == LESType::Deardorff) {
        amrex::Error(
          "It is not recommended to use Deardorff LES and a PBL model");
      }
      // Both the k-eqn RANS closure and these PBL schemes transport RhoKE
      if ((rans_type == RANSType::kEqn) &&
          ((pbl_type == PBLType::MYJ) || (pbl_type == PBLType::MYNN25) ||
           (pbl_type == PBLType::MYNNEDMF) || uses_shoc_family())) {
        amrex::Error("erf.rans_type = kEqn cannot be combined with erf.pbl_type = " +
                     std::string(amrex::getEnumNameString(pbl_type)) +
                     ": both transport the turbulent kinetic energy");
      }
 
      if (pbl_type == PBLType::MYNN25 || pbl_type == PBLType::MYNNEDMF) {
        query_one_or_per_level(
          pp, "pbl_mynn_Lt_taper_exp", pbl_mynn.Lt_taper_exp, lev, max_level);
        if (pbl_mynn.Lt_taper_exp < 0) {
            amrex::Error("erf.pbl_mynn_Lt_taper_exp must be >= 0 "
                         "(0 disables the taper)");
        }
        // NN09 (the default) uses Lt_alpha = 0.23 and a 1/l sum for the master
        // length scale; Chen2021 uses Lt_alpha = 0.10 and a 1/l^2 sum.
        query_one_or_per_level_enum_case_insensitive(
          pp, "pbl_mynn_config", pbl_mynn.config, lev, max_level);
        // l_T = Lt_alpha*<zq>/<q> follows the config unless set explicitly.
        pbl_mynn.Lt_alpha = (pbl_mynn.config == MYNNConfigType::Chen2021)
                              ? amrex::Real(0.1) : amrex::Real(0.23);
        query_one_or_per_level(
          pp, "pbl_mynn_Lt_alpha", pbl_mynn.Lt_alpha, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_A1", pbl_mynn.A1, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_A2", pbl_mynn.A2, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_B1", pbl_mynn.B1, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_B2", pbl_mynn.B2, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_C1", pbl_mynn.C1, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_C2", pbl_mynn.C2, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_C3", pbl_mynn.C3, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_C4", pbl_mynn.C4, lev, max_level);
        query_one_or_per_level(pp, "pbl_mynn_C5", pbl_mynn.C5, lev, max_level);
        pbl_mynn_level2.init_coeffs(
          pbl_mynn.A1, pbl_mynn.A2, pbl_mynn.B1, pbl_mynn.B2, pbl_mynn.C1,
          pbl_mynn.C2, pbl_mynn.C3, pbl_mynn.C4, pbl_mynn.C5);
        // Moisture is now always diffused by MYNN, as it is by every other
        // PBL scheme; catch the retired flag rather than ignoring it silently.
        if (pp.contains("pbl_mynn_diffuse_moistvars")) {
            amrex::Abort("erf.pbl_mynn_diffuse_moistvars has been removed. "
                         "MYNN always diffuses the moisture variables now, as "
                         "MYJ, YSU, MRF and SHOC always have; there is no "
                         "longer a way to turn it off.");
        }
        query_one_or_per_level(
          pp, "pbl_mynn_SMmin", pbl_mynn.SMmin, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mynn_SMmax", pbl_mynn.SMmax, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mynn_SHmin", pbl_mynn.SHmin, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mynn_SHmax", pbl_mynn.SHmax, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mynn_SQfactor", pbl_mynn.SQfac, lev, max_level);
 
        // SQ = SQfac*SM, so its bound is tied to SHmax and SQfac; the default member
        // initializer was evaluated before either could be overridden here.
        pbl_mynn.SQmax = pbl_mynn.SQfac * pbl_mynn.SHmax;
 
      } else if (pbl_type == PBLType::YSU || pbl_type == PBLType::YSUNew) {
        query_one_or_per_level(
          pp, "pbl_ysu_coriolis_freq", pbl_ysu_coriolis_freq, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_ysu_use_consistent_coriolis",
          pbl_ysu_use_consistent_coriolis, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_ysu_force_over_water", pbl_ysu_force_over_water, lev,
          max_level);
        query_one_or_per_level(
          pp, "pbl_ysu_land_Ribcr", pbl_ysu_land_Ribcr, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_ysu_unst_Ribcr", pbl_ysu_unst_Ribcr, lev, max_level);
        // YSU now uses theta_v at WRF's thvx sites and theta_l at its thlix
        // sites; catch the retired switch between them rather than ignoring it.
        if (pp.contains("enable_ysu_liquid_theta")) {
            amrex::Abort("erf.enable_ysu_liquid_theta has been removed. YSU uses "
                         "theta_v (all condensate loading) for the bulk Richardson "
                         "number and clear-sky entrainment, and theta_l (cloud "
                         "liquid and ice) for the cloud-top extension scan and "
                         "cloudy entrainment, as WRF does with thvx and thlix.");
        }
        query_one_or_per_level(
          pp, "enable_ysu_countergradient", enable_ysu_countergradient, lev, max_level);
        query_one_or_per_level(
          pp, "enable_ysu_terrain_pblh_floor", enable_ysu_terrain_pblh_floor, lev, max_level);
        query_one_or_per_level(
          pp, "enable_ysu_sat_limiter", enable_ysu_sat_limiter, lev, max_level);
        query_one_or_per_level(
          pp, "enable_ysu_topdown", enable_ysu_topdown, lev, max_level);
        query_one_or_per_level(
          pp, "enable_ysu_entrainment", enable_ysu_entrainment, lev, max_level);
        query_one_or_per_level(
          pp, "enable_ysu_cloud_pblh", enable_ysu_cloud_pblh, lev, max_level);
        query_one_or_per_level(
          pp, "enable_mrf_unbounded_vpert", enable_mrf_unbounded_vpert, lev, max_level);
        query_one_or_per_level(
          pp, "ysu_qcloud_threshold", ysu_qcloud_threshold, lev, max_level);
        query_one_or_per_level(
          pp, "ysu_moistvars", ysu_moistvars, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_ysunew_highres_bounds", pbl_ysunew_highres_bounds, lev, max_level);
        // Radiative tendency limiter (YSUNew-only)
        query_one_or_per_level(
          pp, "enable_ysu_rad_tend_limiter", enable_ysu_rad_tend_limiter, lev, max_level);
        query_one_or_per_level(
          pp, "ysu_rad_tend_limiter_magnitude", ysu_rad_tend_limiter_magnitude, lev, max_level);
        if (!(ysu_rad_tend_limiter_magnitude > 0.0)) {
            if (amrex::ParallelDescriptor::IOProcessor()) {
                amrex::Print() << "WARNING: erf.ysu_rad_tend_limiter_magnitude = "
                               << ysu_rad_tend_limiter_magnitude
                               << " is not positive; using the default of 1.0 K/s\n";
            }
            ysu_rad_tend_limiter_magnitude = amrex::Real(1.0);
        }
        // Vogelezang & Holtslag (1996) shear correction
        query_one_or_per_level(
          pp, "enable_vh96_shear_correction", enable_vh96_shear_correction, lev, max_level);
        query_one_or_per_level(
          pp, "vh96_shear_const_b", vh96_shear_const_b, lev, max_level);
        // PBLH spatial smoothing
        query_one_or_per_level(
          pp, "enable_pblh_smoothing", enable_pblh_smoothing, lev, max_level);
        query_one_or_per_level(
          pp, "pblh_smoothing_passes", pblh_smoothing_passes, lev, max_level);
        query_one_or_per_level(
          pp, "pblh_smoothing_weight", pblh_smoothing_weight, lev, max_level);
        // Validate pblh_smoothing_weight is in [0,1]
        if (enable_pblh_smoothing && (pblh_smoothing_weight < 0 || pblh_smoothing_weight > 1)) {
          amrex::Error("pblh_smoothing_weight must be in [0, 1]");
        }
        // QNSE stable stability functions (Sukoriansky et al. 2005)
        query_one_or_per_level(
          pp, "enable_qnse_stable_functions", enable_qnse_stable_functions, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_am", qnse_am, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_bm", qnse_bm, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_ah", qnse_ah, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_bh", qnse_bh, lev, max_level);
      } else if (pbl_type == PBLType::MRF) {
        query_one_or_per_level(
          pp, "pbl_mrf_coriolis_freq", pbl_mrf_coriolis_freq, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mrf_Ribcr", pbl_mrf_Ribcr, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mrf_const_b", pbl_mrf_const_b, lev, max_level);
        query_one_or_per_level(pp, "pbl_mrf_sf", pbl_mrf_sf, lev, max_level);
        query_one_or_per_level(
          pp, "enable_mrf_countergradient", enable_mrf_countergradient, lev, max_level);
        query_one_or_per_level(
          pp, "enable_mrf_cloud_adjustment", enable_mrf_cloud_adjustment, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mrf_highres_bounds", pbl_mrf_highres_bounds, lev, max_level);
        query_one_or_per_level(
          pp, "enable_mrf_unbounded_vpert", enable_mrf_unbounded_vpert, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_mrf_use_zero_ri_extent", pbl_mrf_use_zero_ri_extent, lev, max_level);
        // QNSE stable stability functions (Sukoriansky et al. 2005)
        query_one_or_per_level(
          pp, "enable_qnse_stable_functions", enable_qnse_stable_functions, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_am", qnse_am, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_bm", qnse_bm, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_ah", qnse_ah, lev, max_level);
        query_one_or_per_level(
          pp, "qnse_bh", qnse_bh, lev, max_level);
        // Scale-aware PBL-LES blending parameters.
        query_one_or_per_level(
          pp, "pbl_blend_length", pbl_blend_length, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_blend_cs", pbl_blend_cs, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_blend_c_max", pbl_blend_c_max, lev, max_level);
        query_one_or_per_level(
          pp, "pbl_blend_use_smag", pbl_blend_use_smag, lev, max_level);
        // Vogelezang & Holtslag (1996) shear correction
        query_one_or_per_level(
          pp, "enable_vh96_shear_correction", enable_vh96_shear_correction, lev, max_level);
        query_one_or_per_level(
          pp, "vh96_shear_const_b", vh96_shear_const_b, lev, max_level);
        // PBLH spatial smoothing
        query_one_or_per_level(
          pp, "enable_pblh_smoothing", enable_pblh_smoothing, lev, max_level);
        query_one_or_per_level(
          pp, "pblh_smoothing_passes", pblh_smoothing_passes, lev, max_level);
        query_one_or_per_level(
          pp, "pblh_smoothing_weight", pblh_smoothing_weight, lev, max_level);
        // Validate pblh_smoothing_weight is in [0,1]
        if (enable_pblh_smoothing && (pblh_smoothing_weight < 0 || pblh_smoothing_weight > 1)) {
          amrex::Error("pblh_smoothing_weight must be in [0, 1]");
        }
      }
      if (uses_eamxx_shoc()) {
#ifndef ERF_USE_EAMXX_SHOC
          amrex::Abort("PBLType::EAMXX_SHOC requested, but ERF was not built with ERF_ENABLE_EAMXX_SHOC=ON");
#endif
      }
 
      if (pbl_requires_surface_layer()) {
          std::string zlo_bc = "none";
          amrex::ParmParse pp_bc("zlo");
          pp_bc.get("type",zlo_bc);
          if (amrex::toLower(zlo_bc) != "surface_layer") {
               amrex::Abort("You must use the surface_layer BC at zlo with the selected PBL.");
          }
      }
    }
 
    // Flags for QKE/TKE equation
    if (pbl_type == PBLType::MYJ || pbl_type == PBLType::MYNN25 || pbl_type == PBLType::MYNNEDMF) {
        // Add sources/sinks to QKE/TKE? (MYJ does this inline)
        if (pbl_type == PBLType::MYNN25 || pbl_type == PBLType::MYNNEDMF) {
            use_pbl_tke = true;
        }
        // Advect QKE/TKE?
        query_one_or_per_level(pp, "advect_tke"    , advect_tke    , lev, max_level);
        // Apply numerical diffusion to QKE/TKE?
        query_one_or_per_level(pp, "diffuse_tke_3D", diffuse_tke_3D, lev, max_level);
    }
 
    // There is a default value of 1e-6 but the user can override the value here,
    //   and in addition can add perturbational values.  and in addition can add perturbational values.
    pp.queryAdd("tke_min",tke_min);
 
    // LES constants...
    query_one_or_per_level(pp, "Cs", Cs, lev, max_level);
 
    query_one_or_per_level(pp, "Pr_t", Pr_t, lev, max_level);
    query_one_or_per_level(pp, "Sc_t", Sc_t, lev, max_level);
 
    // Compute relevant forms of diffusion parameters
    Pr_t_inv = one / Pr_t;
    Sc_t_inv = one / Sc_t;
 
    if (les_type == LESType::Deardorff) {
        query_one_or_per_level(pp, "Ck", Ck, lev, max_level);
        query_one_or_per_level(pp, "Ce", Ce, lev, max_level);
        query_one_or_per_level(pp, "Ce_wall", Ce_wall, lev, max_level);
    }
 
    // To quantify atmospheric stability for subgrid modeling
    const bool strat_type_set = (pp.countval("thermal_stratification") > 0);
    query_one_or_per_level(pp, "thermal_stratification", strat_type, lev, max_level);
    // The PBL schemes have always used theta_v for their buoyancy gradient,
    // which is the physically right choice there.  Honour an explicit
    // erf.thermal_stratification everywhere, but when the user has not set it
    // keep the historical split: theta for the LES closures, theta_v for the
    // PBL ones.  That way the knob reaches the PBL schemes without silently
    // changing any existing moist answer.
    pbl_strat_type = strat_type_set ? strat_type : StratType::thetav;
    amrex::Print() << "PBL thermal stratification variable: "
                   << amrex::getEnumNameString(pbl_strat_type) << std::endl;
    if (strat_type == StratType::theta) {
        amrex::Print() << "Thermal stratification based on gradient of potential temperature" << std::endl;
    } else if (strat_type == StratType::thetav) {
        amrex::Print() << "Thermal stratification based on gradient of virtual potential temperature" << std::endl;
    } else if (strat_type == StratType::thetal) {
        amrex::Print() << "Thermal stratification based on gradient of linearized liquid-water potential temperature" << std::endl;
    }
 
    // k-eqn constants
    query_one_or_per_level(pp, "Cmu0", Cmu0, lev, max_level);
    query_one_or_per_level(pp, "Cb", Cb, lev, max_level);
    query_one_or_per_level(pp, "Rt_crit", Rt_crit, lev, max_level);
    query_one_or_per_level(pp, "Rt_min", Rt_min, lev, max_level);
    query_one_or_per_level(pp, "max_geom_lscale", l_g_max, lev, max_level);
    query_one_or_per_level(pp, "dirichlet_k", dirichlet_k, lev, max_level);
    query_one_or_per_level(pp, "tke_floor", tke_floor, lev, max_level);
    query_one_or_per_level(pp, "implicit_tke_dissipation", implicit_tke_dissipation, lev, max_level);
    query_one_or_per_level(pp, "rans_consistent_diffusivities", rans_consistent_diffusivities, lev, max_level);
    query_one_or_per_level(pp, "rans_lscale_from_pblh", rans_lscale_from_pblh, lev, max_level);
    query_one_or_per_level(pp, "rans_lscale_min", rans_lscale_min, lev, max_level);
    query_one_or_per_level(pp, "rans_realizability", rans_realizability, lev, max_level);
    if (rans_realizability != zero) {
        if (!(rans_realizability >= amrex::Real(0.01) && rans_realizability <= amrex::Real(10))) {
            amrex::Error("erf.rans_realizability must be 0 (off) or in [0.01, 10]; Cmu0^2 = 0.309 is the neutral equilibrium value");
        }
        if (rans_type != RANSType::kEqn) {
            // One value for the whole hierarchy (a single entry, or a list whose entries are all
            // equal) is meant for the levels that run the k-eqn closure and is reset here; a list
            // with different entries names this level on purpose and must not switch the cap on
            // where there is no closure.
            amrex::Vector<amrex::Real> all_levels;
            pp.getarr("rans_realizability", all_levels);
            bool uniform = true;
            for (auto v : all_levels) { uniform = uniform && (v == all_levels[0]); }
            if (uniform) {
                amrex::Print() << "erf.rans_realizability reset to 0 on level " << lev
                               << " (no k-eqn closure there)" << std::endl;
                rans_realizability = zero;
            } else {
                amrex::Error("erf.rans_realizability needs erf.rans_type = kEqn (it caps the k-eqn eddy viscosity; level "
                             + std::to_string(lev) + ")");
            }
        }
    }
    // After the reset above, so a hierarchy with a PBL scheme on an LES level and the cap on
    // its k-eqn level is accepted
    if ((rans_realizability != zero) && (pbl_type != PBLType::None)) {
        amrex::Error("erf.rans_realizability needs erf.pbl_type = None: a PBL scheme overwrites the vertical "
                     "diffusivities the cap acts on (level " + std::to_string(lev) + ")");
    }
    if (rans_lscale_from_pblh && !(rans_lscale_min > 0)) {
        amrex::Error("erf.rans_lscale_min must be > 0 with erf.rans_lscale_from_pblh");
    }
 
    if (rans_type == RANSType::kEqn) {
        // The stability functions (AL01 Eqs. 31, 32) have poles at Rt = -3.59869
        // (cmu) and -3.61011 (cmu'); the Burchard & Petersen smoothing maps
        // Rt < Rt_crit onto (Rt_min, Rt_crit) and needs Rt_min < Rt_crit <= 0,
        // with Rt_min above the pole of cmu (AL01::Rt_min_lower_bound) and above
        // the zero of its numerator (AL01::cmu_numerator_root), so cmu > 0.
        if (!(Cmu0 > 0)) { amrex::Error("erf.Cmu0 must be > 0 for the k-eqn RANS"); }
        if (!(Cb > 0)) { amrex::Error("erf.Cb must be > 0 for the k-eqn RANS"); }
        if (!(l_g_max > 0)) { amrex::Error("erf.max_geom_lscale must be > 0 for the k-eqn RANS"); }
        if (!(Rt_crit <= 0)) { amrex::Error("erf.Rt_crit must be <= 0 for the k-eqn RANS"); }
        if (!(Rt_min < Rt_crit)) { amrex::Error("erf.Rt_min must be < erf.Rt_crit for the k-eqn RANS"); }
        if (!(Rt_min > AL01::Rt_min_lower_bound)) {
            std::ostringstream msg;
            msg << std::setprecision(8) << "erf.Rt_min must be > " << AL01::Rt_min_lower_bound
                << " for the k-eqn RANS: the momentum stability function cmu has a pole at Rt = -3.5986889850,"
                << " and just below it cmu and the eddy viscosity are negative (level " << lev << ")";
            amrex::Error(msg.str());
        }
        if (!(Rt_min > AL01::cmu_numerator_root(Cmu0))) {
            std::ostringstream msg;
            msg << std::setprecision(8) << "erf.Rt_min must be > -erf.Cmu0/0.108 = " << AL01::cmu_numerator_root(Cmu0)
                << " for the k-eqn RANS: below it the numerator Cmu0 + 0.108 Rt of cmu is negative;"
                << " raise erf.Cmu0 or erf.Rt_min (level " << lev << ")";
            amrex::Error(msg.str());
        }
    }
    if (tke_floor < 0) { amrex::Error("erf.tke_floor must be >= 0"); }
 
    // Common inputs (LES or RANS)
    if (!query_one_or_per_level(pp, "sigma_k", sigma_k, lev, max_level) && rans_type == RANSType::kEqn) {
        amrex::Print() << "Overriding default sigma_k for k-eqn RANS" << std::endl;
        sigma_k = one;
    };
    query_one_or_per_level(pp, "theta_ref", theta_ref, lev, max_level);
 
    query_one_or_per_level(pp, "mix_isotropic", mix_isotropic, lev, max_level);
    query_one_or_per_level(pp, "use_Ri_correction", use_Ri_correction, lev, max_level);
    query_one_or_per_level(pp, "Ri_crit", Ri_crit, lev, max_level);
 
    // Set common flags
    use_kturb =
      ((les_type != LESType::None) || (rans_type != RANSType::None) ||
       (pbl_type != PBLType::None));
    use_keqn =
      ((les_type == LESType::Deardorff) || (rans_type == RANSType::kEqn));
    use_tke =
      ((les_type == LESType::Deardorff) || (rans_type == RANSType::kEqn) ||
       (pbl_type == PBLType::MYJ)       || (pbl_type == PBLType::MYNN25) ||
       (pbl_type == PBLType::MYNNEDMF)  || uses_shoc_family());
 
    if (use_tke) {
        query_one_or_per_level(pp, "init_tke_from_ustar", init_tke_from_ustar, lev, max_level);
    }
 
    // Validate inputs
    if (les_type == LESType::Smagorinsky) {
      if (Cs == 0) {
        amrex::Error("Need to specify Cs for Smagorsinky LES");
      }
      if (smag2d && mix_isotropic) {
        amrex::Print() << "Turning off mix_isotropic for 2-D Smagorinsky" << std::endl;
        mix_isotropic = false;
      }
    }
 
    // Opt-in limits of WRF's Smagorinsky-2D closure (dyn_em/module_diffusion_em.F, smag2d_km):
    // the slope limiter, which divides K_h by alpha or alpha^2 on terrain-fitted meshes, and the
    // cap K_h <= smag2d_kh_cap * Delta_h (WRF: 10 m/s).  Both are off by default.  The mesh check
    // for the slope limiter is in SolverChoice::init_params, where the terrain type is known.
    query_one_or_per_level(pp, "smag2d_slope_limiter", smag2d_slope_limiter, lev, max_level);
    query_one_or_per_level(pp, "smag2d_kh_cap", smag2d_kh_cap, lev, max_level);
    if (!(smag2d_kh_cap >= zero)) {
      amrex::Error("erf.smag2d_kh_cap must be >= 0 (0 = no cap; WRF uses 10 m/s)");
    }
    smag2d_limits_requested = (smag2d_slope_limiter || smag2d_kh_cap > zero);
    if (!smag2d) {
      // A single value applies to every level; it is meant for the levels that run
      // Smagorinsky2D (SolverChoice aborts if none does).  A per-level list must not
      // switch the limits on at a level that does not.
      if (pp.countval("smag2d_slope_limiter") == 1) { smag2d_slope_limiter = false; }
      if (pp.countval("smag2d_kh_cap") == 1)        { smag2d_kh_cap = zero; }
      if (smag2d_slope_limiter || smag2d_kh_cap > zero) {
        amrex::Error("erf.smag2d_slope_limiter and erf.smag2d_kh_cap apply only to "
                     "erf.les_type = Smagorinsky2D (level " + std::to_string(lev) + ")");
      }
    }
  }
 
  /**
   * @brief Validate turbulence options against physical boundary conditions.
   * @param phys_bc_type Physical boundary-condition types.
   */
  void check_params (amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2>& phys_bc_type)
  {
      // BC compatibility
      if (pbl_requires_surface_layer() &&
          phys_bc_type[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)] != ERF_BC::surface_layer ) {
          amrex::Abort("The selected PBL model requires MOST at lower boundary");
      }
      if ( (les_type == LESType::Deardorff) && (Ce_wall > 0) &&
           (phys_bc_type[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)] != ERF_BC::surface_layer) &&
           (phys_bc_type[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)] != ERF_BC::slip_wall) &&
           (phys_bc_type[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)] != ERF_BC::no_slip_wall) )
      {
          amrex::Warning("Deardorff LES assumes wall at zlo when applying Ce_wall");
      }
  }
 
  /**
   * @brief Print turbulence settings for one AMR level.
   * @param lev AMR level index.
   */
  void display (int lev)
  {
    amrex::Print() << "Turbulence Settings at level " << lev << std::endl;
 
    if (
      les_type == LESType::None && rans_type == RANSType::None &&
      pbl_type == PBLType::None) {
      amrex::Print() << "    Using DNS model at level " << lev << std::endl;
    } else if (les_type == LESType::Smagorinsky) {
      if (smag2d) {
        amrex::Print() << "    Using 2D Smagorinsky LES model at level " << lev << std::endl;
      } else {
        amrex::Print() << "    Using Smagorinsky LES model at level " << lev << std::endl;
      }
      if (use_Ri_correction) {
        amrex::Print() << "    Smagorinsky uses Richardson number correction with Ri_crit = "
                       << Ri_crit << std::endl;
      }
      if (smag2d) {
        amrex::Print() << "    WRF slope limiter       : " << smag2d_slope_limiter << std::endl;
        if (smag2d_kh_cap > zero) {
          amrex::Print() << "    K_h cap                 : " << smag2d_kh_cap << " * Delta_h" << std::endl;
        } else {
          amrex::Print() << "    K_h cap                 : none" << std::endl;
        }
      }
    } else if (les_type == LESType::Deardorff) {
      amrex::Print() << "    Using Deardorff LES model at level " << lev << std::endl;
    } else if (rans_type == RANSType::kEqn) {
      amrex::Print()
        << "    Using Axell & Liungman one-equation RANS k model at level " << lev
        << std::endl;
    } else if (pbl_type == PBLType::MYJ) {
      amrex::Print() << "    Using MYJ PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::MYNN25) {
      amrex::Print() << "    Using MYNN2.5 PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::MYNNEDMF) {
      amrex::Print() << "    Using MYNNEDMF PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::YSU) {
      amrex::Print() << "    Using YSU PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::YSUNew) {
      amrex::Print() << "    Using YSU PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::MRF) {
      amrex::Print() << "    Using MRF PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::EAMXX_SHOC) {
      amrex::Print() << "    Using EAMxx SHOC PBL model at level " << lev << std::endl;
    } else if (pbl_type == PBLType::NATIVE_SHOC) {
      amrex::Print() << "    Using native SHOC PBL model at level " << lev << std::endl;
    } else {
      amrex::Error("Unknown turbulence model");
    }
 
    if (les_type != LESType::None) {
      if (les_type == LESType::Smagorinsky) {
        amrex::Print() << "    Cs                      : " << Cs << std::endl;
      }
      if (les_type == LESType::Deardorff) {
        amrex::Print() << "    Ce                      : " << Ce      << std::endl;
        amrex::Print() << "    Ce at wall              : " << Ce_wall << std::endl;
        amrex::Print() << "    Ck                      : " << Ck      << std::endl;
        amrex::Print() << "    sigma_k                 : " << sigma_k << std::endl;
 
        // Sullivan et al 1994, Eqn 14
        amrex::Real Cs_equiv = std::sqrt(Ck * std::sqrt(Ck / Ce));
        amrex::Print() << "    equivalent Cs           : " << Cs_equiv
                       << std::endl;
      }
      amrex::Print() << "    isotropic mixing        : " << mix_isotropic
                     << std::endl;
    }
 
    if (rans_type != RANSType::None) {
      if (rans_type == RANSType::kEqn) {
        amrex::Print() << "Cmu0                        : " << Cmu0 << std::endl;
        amrex::Print() << "sigma_k                     : " << sigma_k << std::endl;
        amrex::Print() << "Cb                          : " << Cb << std::endl;
        amrex::Print() << "Rt_crit                     : " << Rt_crit << std::endl;
        amrex::Print() << "Rt_min                      : " << Rt_min << std::endl;
        amrex::Print() << "max_geom_lscale             : " << l_g_max << std::endl;
        amrex::Print() << "rans_realizability          : " << rans_realizability << std::endl;
        amrex::Print() << "dirichlet_k                 : " << dirichlet_k << std::endl;
        amrex::Print() << "tke_floor                   : " << tke_floor << std::endl;
        amrex::Print() << "implicit_tke_dissipation    : " << implicit_tke_dissipation << std::endl;
        amrex::Print() << "consistent diffusivities    : " << rans_consistent_diffusivities << std::endl;
        amrex::Print() << "lscale cap from PBL height  : " << rans_lscale_from_pblh
                       << " (floor " << rans_lscale_min << " m)" << std::endl;
      }
    }
 
    if ((les_type  == LESType::Deardorff) ||
        (rans_type == RANSType::kEqn)) {
      if (theta_ref > 0) {
        amrex::Print() << "    reference theta         : " << theta_ref << std::endl;
      } else {
        amrex::Print() << "    reference theta         : n/a" << std::endl;
      }
    }
 
    if ((les_type != LESType::None) || (rans_type != RANSType::None)) {
      amrex::Print() << "    Pr_t                    : " << Pr_t << std::endl;
      amrex::Print() << "    Sc_t                    : " << Sc_t << std::endl;
    }
 
    if (pbl_type == PBLType::MYNN25 || pbl_type == PBLType::MYNNEDMF) {
      amrex::Print() << "    pbl_mynn_Lt_taper_exp : "
                     << pbl_mynn.Lt_taper_exp << std::endl;
      amrex::Print() << "    pbl_mynn_config : "
                     << amrex::getEnumNameString(pbl_mynn.config) << std::endl;
      amrex::Print() << "    pbl_mynn_Lt_alpha : " << pbl_mynn.Lt_alpha << std::endl;
      amrex::Print() << "    pbl_mynn_A1 : " << pbl_mynn.A1 << std::endl;
      amrex::Print() << "    pbl_mynn_A2 : " << pbl_mynn.A2 << std::endl;
      amrex::Print() << "    pbl_mynn_B1 : " << pbl_mynn.B1 << std::endl;
      amrex::Print() << "    pbl_mynn_B2 : " << pbl_mynn.B2 << std::endl;
      amrex::Print() << "    pbl_mynn_C1 : " << pbl_mynn.C1 << std::endl;
      amrex::Print() << "    pbl_mynn_C2 : " << pbl_mynn.C2 << std::endl;
      amrex::Print() << "    pbl_mynn_C3 : " << pbl_mynn.C3 << std::endl;
      amrex::Print() << "    pbl_mynn_C4 : " << pbl_mynn.C4 << std::endl;
      amrex::Print() << "    pbl_mynn_C5 : " << pbl_mynn.C5 << std::endl;
    } else if (pbl_type == PBLType::YSU || pbl_type == PBLType::YSUNew) {
      amrex::Print() << "    pbl_ysu_coriolis_freq              : "
                     << pbl_ysu_coriolis_freq << std::endl;
      amrex::Print() << "    pbl_ysu_use_consistent_coriolis    : "
                     << pbl_ysu_use_consistent_coriolis << std::endl;
      amrex::Print() << "    pbl_ysu_force_over_water           : "
                     << pbl_ysu_force_over_water << std::endl;
      amrex::Print() << "    pbl_ysu_land_Ribcr                 : "
                     << pbl_ysu_land_Ribcr << std::endl;
      amrex::Print() << "    pbl_ysu_unst_Ribcr                 : "
                     << pbl_ysu_unst_Ribcr << std::endl;
      amrex::Print() << "    enable_ysu_countergradient         : "
                     << enable_ysu_countergradient << std::endl;
      amrex::Print() << "    enable_ysu_terrain_pblh_floor      : "
                     << enable_ysu_terrain_pblh_floor << std::endl;
      amrex::Print() << "    enable_ysu_sat_limiter             : "
                     << enable_ysu_sat_limiter << std::endl;
      amrex::Print() << "    enable_ysu_topdown                 : "
                     << enable_ysu_topdown << std::endl;
      amrex::Print() << "    enable_ysu_entrainment             : "
                     << enable_ysu_entrainment << std::endl;
      amrex::Print() << "    enable_ysu_cloud_pblh              : "
                     << enable_ysu_cloud_pblh << std::endl;
      amrex::Print() << "    ysu_qcloud_threshold               : "
                     << ysu_qcloud_threshold << std::endl;
      amrex::Print() << "    ysu_moistvars                      : "
                     << ysu_moistvars << std::endl;
      amrex::Print() << "    pbl_ysunew_highres_bounds          : "
                     << pbl_ysunew_highres_bounds << std::endl;
      // Radiative tendency limiter/smoothing (YSUNew-only)
      amrex::Print() << "    enable_ysu_rad_tend_limiter        : "
                     << enable_ysu_rad_tend_limiter << std::endl;
      amrex::Print() << "    ysu_rad_tend_limiter_magnitude     : "
                     << ysu_rad_tend_limiter_magnitude << std::endl;
      amrex::Print() << "    enable_vh96_shear_correction       : "
                     << enable_vh96_shear_correction << std::endl;
      amrex::Print() << "    vh96_shear_const_b                 : "
                     << vh96_shear_const_b << std::endl;
      amrex::Print() << "    enable_pblh_smoothing              : "
                     << enable_pblh_smoothing << std::endl;
      amrex::Print() << "    pblh_smoothing_passes              : "
                     << pblh_smoothing_passes << std::endl;
      amrex::Print() << "    pblh_smoothing_weight              : "
                     << pblh_smoothing_weight << std::endl;
      amrex::Print() << "    enable_mrf_unbounded_vpert         : "
                     << enable_mrf_unbounded_vpert << std::endl;
      amrex::Print() << "    enable_qnse_stable_functions       : "
                     << enable_qnse_stable_functions << std::endl;
      if (enable_qnse_stable_functions) {
        amrex::Print() << "      qnse_am (momentum coeff)         : " << qnse_am << std::endl;
        amrex::Print() << "      qnse_bm (momentum coeff)         : " << qnse_bm << std::endl;
        amrex::Print() << "      qnse_ah (heat coeff)             : " << qnse_ah << std::endl;
        amrex::Print() << "      qnse_bh (heat coeff)             : " << qnse_bh << std::endl;
      }
    } else if (pbl_type == PBLType::MRF) {
      amrex::Print() << "    pbl_mrf_coriolis_freq : " << pbl_mrf_coriolis_freq
                     << std::endl;
      amrex::Print() << "    pbl_mrf_Ribcr         : " << pbl_mrf_Ribcr
                     << std::endl;
      amrex::Print() << "    pbl_mrf_const_b       : " << pbl_mrf_const_b
                     << std::endl;
      amrex::Print() << "    pbl_mrf_sf            : " << pbl_mrf_sf
                     << std::endl;
      amrex::Print() << "    enable_mrf_countergradient : " << enable_mrf_countergradient
                     << std::endl;
      amrex::Print() << "    enable_mrf_cloud_adjustment : " << enable_mrf_cloud_adjustment
                     << std::endl;
      amrex::Print() << "    pbl_mrf_highres_bounds     : " << pbl_mrf_highres_bounds
                     << std::endl;
      amrex::Print() << "    enable_mrf_unbounded_vpert : " << enable_mrf_unbounded_vpert
                     << std::endl;
      amrex::Print() << "    pbl_mrf_use_zero_ri_extent : " << pbl_mrf_use_zero_ri_extent
                     << std::endl;
      amrex::Print() << "    pbl_blend_length   : " << pbl_blend_length
                     << std::endl;
      amrex::Print() << "    pbl_blend_cs       : " << pbl_blend_cs
                     << std::endl;
      amrex::Print() << "    pbl_blend_c_max    : " << pbl_blend_c_max
                     << std::endl;
      amrex::Print() << "    pbl_blend_use_smag : " << pbl_blend_use_smag
                     << std::endl;
      amrex::Print() << "    enable_vh96_shear_correction : " << enable_vh96_shear_correction
                     << std::endl;
      amrex::Print() << "    vh96_shear_const_b           : " << vh96_shear_const_b
                     << std::endl;
      amrex::Print() << "    enable_pblh_smoothing        : " << enable_pblh_smoothing
                     << std::endl;
      amrex::Print() << "    pblh_smoothing_passes        : " << pblh_smoothing_passes
                     << std::endl;
      amrex::Print() << "    pblh_smoothing_weight        : " << pblh_smoothing_weight
                     << std::endl;
      amrex::Print() << "    enable_qnse_stable_functions       : "
                     << enable_qnse_stable_functions << std::endl;
      if (enable_qnse_stable_functions) {
        amrex::Print() << "      qnse_am (momentum coeff)         : " << qnse_am << std::endl;
        amrex::Print() << "      qnse_bm (momentum coeff)         : " << qnse_bm << std::endl;
        amrex::Print() << "      qnse_ah (heat coeff)             : " << qnse_ah << std::endl;
        amrex::Print() << "      qnse_bh (heat coeff)             : " << qnse_bh << std::endl;
      }
    }
  }
 
  // LES model
  LESType les_type = LESType::None;
 
  // Turbulent Prandtl number
  amrex::Real Pr_t = one / three;
  amrex::Real Pr_t_inv = three;
 
  // Turbulent Schmidt number
  amrex::Real Sc_t = one;
  amrex::Real Sc_t_inv = one;
 
  // Smagorinsky
  amrex::Real Cs = zero;
  bool smag2d = false;
  // WRF smag2d_km limits (opt-in): slope limiter K_h /= alpha or alpha^2 on terrain-fitted
  // meshes, and the cap K_h <= smag2d_kh_cap * Delta_h (m/s; 0 = no cap)
  bool smag2d_slope_limiter = false;
  amrex::Real smag2d_kh_cap = zero;
  // whether either limit was given as on for this level, before a single value is
  // dropped at a level without Smagorinsky2D
  bool smag2d_limits_requested = false;
 
  // Deardorff
  amrex::Real Ce = amrex::Real(0.93);
  amrex::Real Ce_wall = zero; // if > 0, then set Ce to this at k=0
  amrex::Real Ck = amrex::Real(0.1);
 
  // k-eqn RANS coefficients (Axell & Liungman 2001)
  amrex::Real Cmu0 = amrex::Real(0.5562);
  amrex::Real Cb = amrex::Real(0.35);
  amrex::Real Rt_crit = -one;
  amrex::Real Rt_min = -three;
  amrex::Real l_g_max = amrex::Real(30.0); // ~ kappa * (amrex::Real(0.1) * zi)
 
  // Deardorff or k-eqn RANS
  // - diffusivity of tke is 1/sigma_k times the eddy viscosity
  amrex::Real sigma_k = amrex::Real(0.5);
  // - reference potential temperature used to quantify the stratification in
  //   a stable region
  amrex::Real theta_ref = zero;
  // - how to quantify stratification effects
  StratType strat_type = StratType::theta;
  // Stratification variable used by the PBL schemes; see init_params.
  StratType pbl_strat_type = StratType::thetav;
 
  // Anisotropic length scales
  bool mix_isotropic = true;
 
  bool use_Ri_correction = true;
  amrex::Real Ri_crit = fourth;
 
  // RANS type
  RANSType rans_type = RANSType::None;
 
  bool dirichlet_k = false;
 
  // Runtime floor on the turbulent kinetic energy [m2/s2] in the k-equation
  // update and the RANS closure; 0 means machine epsilon on rho*k (the
  // historical floor). Distinct from tke_min, the initial TKE value.
  amrex::Real tke_floor = zero;
 
  // Treat the TKE dissipation implicitly: eps = (Cmu0^3 sqrt(k_old)/L) * k_new
  // folded into the update (Deardorff and k-eqn RANS)
  bool implicit_tke_dissipation = false;
 
  // k-eqn RANS: horizontal heat, scalar and moisture diffusivities follow the
  // AL01 scalar stability function (cmu') like the vertical heat diffusivity,
  // instead of Pr_t and Sc_t times the eddy viscosity
  bool rans_consistent_diffusivities = false;
 
  // k-eqn RANS: cap the geometric length at kappa * 0.1 * PBL height (from
  // erf.most.pblh_calc), clamped to [rans_lscale_min, max_geom_lscale],
  // instead of the fixed max_geom_lscale
  bool rans_lscale_from_pblh = false;
  amrex::Real rans_lscale_min = one;
 
  // k-eqn RANS: realizability cap nu_t <= C_r k / S on the eddy viscosity
  // (and every eddy diffusivity through the same ratio), 0 = off.  C_r = Cmu0^2
  // leaves a neutral equilibrium surface layer untouched and limits the
  // viscosity in the shear layer that leaves a crest (AL01 Eqs. 11 and 19;
  // the SST limiter form).
  amrex::Real rans_realizability = zero;
 
  // PBL model
  PBLType pbl_type = PBLType::None;
 
  /**
   * @brief Query whether this level uses the EAMxx SHOC PBL scheme.
   * @return True if PBLType::EAMXX_SHOC is selected.
   */
  bool uses_eamxx_shoc () const noexcept
  {
      return pbl_type == PBLType::EAMXX_SHOC;
  }
 
  /**
   * @brief Query whether this level uses the native SHOC PBL scheme.
   * @return True if PBLType::NATIVE_SHOC is selected.
   */
  bool uses_native_shoc () const noexcept
  {
      return pbl_type == PBLType::NATIVE_SHOC;
  }
 
  /**
   * @brief Query whether this level uses any SHOC-family PBL scheme.
   * @return True if native or EAMxx SHOC is selected.
   */
  bool uses_shoc_family () const noexcept
  {
      return uses_eamxx_shoc() || uses_native_shoc();
  }
 
  /**
   * @brief Query whether the horizontal and vertical eddy viscosities can differ on this level.
   * @return True for Smagorinsky-2D, anisotropic Smagorinsky or Deardorff, and any PBL scheme
   *         (which sets the vertical diffusivities on its own).
   */
  bool kh_kv_can_differ () const noexcept
  {
      return ( (les_type == LESType::Smagorinsky && (smag2d || !mix_isotropic)) ||
               (les_type == LESType::Deardorff   && !mix_isotropic) ||
               (pbl_type != PBLType::None) );
  }
 
  /**
   * @brief Query whether the selected PBL scheme requires surface-layer boundary conditions.
   * @return True if the PBL scheme requires a lower surface-layer boundary.
   */
  bool pbl_requires_surface_layer () const noexcept
  {
      return (pbl_type == PBLType::MYNN25)   ||
             (pbl_type == PBLType::MYNNEDMF) ||
             (pbl_type == PBLType::YSU)      ||
             (pbl_type == PBLType::YSUNew)      ||
             (pbl_type == PBLType::MRF)       ||
             uses_shoc_family();
  }
 
  /**
   * @brief Query whether the selected PBL scheme suppresses microphysics condensation.
   * @return True when the PBL scheme owns condensation handling.
   */
  bool pbl_suppresses_microphysics_condensation () const noexcept
  {
      return uses_shoc_family();
  }
 
  MYNNLevel25 pbl_mynn;
  MYNNLevel2 pbl_mynn_level2; // limiting in the decaying turbulence regime
 
  // Common Flags
  bool use_kturb = false; // Any turbulence modeling?
  bool use_keqn =
    false; // Any microscale turbulence modeling (LES, RANS) with TKE closure?
           // Then need to populate SmnSmn_lev for production term.
  bool use_pbl_tke =
    false; // Any mesoscale turbulence modeling (PBL) with TKE closure?
  bool use_tke = false; // Any TKE closure (meso or microscale)?
 
  // Initialize TKE/QKE with linear profiles whose surface values are a
  // function of the friction velocity calculated by the surface layer scheme
  // (e.g., as done in MYNN-EDMF)
  bool init_tke_from_ustar = false;
 
  // This is the value of tke_min in WRF 4.5
  amrex::Real tke_min = amrex::Real(1.e-6);
 
  // Model coefficients - YSU
  // TODO: Add parmparse for all of these above
  amrex::Real pbl_ysu_coriolis_freq =
    amrex::Real(1.0e-4); // 1e-4 is hardcoded in WRF, we let the user specify or take the
            // value from ERF coriolis forcing
  bool pbl_ysu_use_consistent_coriolis =
    false; // ignore input pbl_ysu_coriolis_freq, take value from ERF coriolis
           // forcing instead
  bool pbl_ysu_force_over_water =
    false; // Force YSU to act as if it is over water regardless of other inputs
           // (for testing)
  amrex::Real pbl_ysu_land_Ribcr =
    fourth; // Critical Bulk Richardson number of Land for stable conditions
  amrex::Real pbl_ysu_unst_Ribcr =
    zero; // Critical Bulk Richardson number for unstable conditions
  bool enable_ysu_countergradient =
    true; // if true, apply YSU-style HGAMT/HGAMQ/HGAMU/HGAMV countergradient correction
  bool enable_ysu_terrain_pblh_floor =
    true; // if true, apply terrain-following floor on PBL height
  bool enable_ysu_sat_limiter =
    false; // if true, apply saturation limiter to moisture countergradient
  bool enable_ysu_topdown =
    true; // if true, enable top-down mixing for cloud-topped boundary layers
  bool enable_ysu_entrainment =
    true; // if true, enable entrainment layer parameterization
  bool enable_ysu_cloud_pblh =
    true; // if true, enable cloud-based PBL height detection
  amrex::Real ysu_qcloud_threshold =
    amrex::Real(1.0e-4); // Cloud liquid water threshold for YSUNew (kg/kg)
  bool ysu_moistvars = false; // if true, adds turbulence to moisture
  bool pbl_ysunew_highres_bounds = false; // if true, apply high-resolution grid-dependent diffusivity bounds
  // Radiative tendency limiter (YSUNew-only)
  bool enable_ysu_rad_tend_limiter = false; // if true, guard and bound the qheating_rates top-down forcing
  amrex::Real ysu_rad_tend_limiter_magnitude = amrex::Real(1.0); // Bound on the limited radiative tendency (K/s)
  // Vogelezang & Holtslag (1996) shear correction and PBLH smoothing for YSU/YSUNew
  bool enable_vh96_shear_correction = false; // if true, use VH96 shear correction term instead of ad-hoc floor
  amrex::Real vh96_shear_const_b = amrex::Real(100.0); // VH96 shear correction constant b
  bool enable_pblh_smoothing = false; // if true, apply spatial smoothing to diagnosed PBLH
  int pblh_smoothing_passes = 1; // number of smoothing iterations
  amrex::Real pblh_smoothing_weight = amrex::Real(0.5); // center-cell weight in smoothing stencil
  // Model coefficients - MRF
  amrex::Real pbl_mrf_coriolis_freq = amrex::Real(1.0e-4);
  amrex::Real pbl_mrf_Ribcr   = amrex::Real(0.5); // Critical Bulk Richardson number for MRF PBL model
  amrex::Real pbl_mrf_const_b = amrex::Real(7.8); // Constant b in MRF PBL model, used to compute the PBL height
  amrex::Real pbl_mrf_sf      = amrex::Real(0.1); // MRF surface flux, used to compute the PBL height
  // Scale-aware PBL-LES blending parameters.
  // Boutle et al. (2014): https://doi.org/10.1175/MWR-D-13-00229.1
  // Set pbl_blend_length = 0 to disable (default). Recovers original behaviour.
  amrex::Real pbl_blend_length   = amrex::Real(0.0);   ///< Boutle blending length L [m]. 0 = off.
  amrex::Real pbl_blend_cs       = amrex::Real(0.17);  ///< Smagorinsky coeff for K_h ceiling.
  amrex::Real pbl_blend_c_max    = amrex::Real(0.1);   ///< Power-law ceiling coeff [m^(2/3)/s].
  bool        pbl_blend_use_smag = true;               ///< Use Smagorinsky ceiling (else power-law).
  bool enable_mrf_countergradient = false; // if true, apply WRF-style HGAMT/HGAMQ correction to PBL height
  bool        pbl_ib_aware = false;              // MRF / YSUNew: measure heights from the top of immersed solid columns
  amrex::Real pbl_ib_z0    = amrex::Real(0.01);  // roughness [m] of the neutral log law at the top of an immersed column
  bool enable_mrf_cloud_adjustment = false; // if true, apply cloud-aware stability function adjustments in MRF scheme
  bool pbl_mrf_highres_bounds = false; // if true, apply high-resolution grid-dependent diffusivity bounds instead of global forecast bounds
  bool enable_mrf_unbounded_vpert = false; // if true, does not limit VPERT to GAMCRT
  bool pbl_mrf_use_zero_ri_extent = false; // if true, use Ri=0 (pbli_zero_arr) for K-profile extent; default false uses Ri=0.5 corrector (pbli_arr) for WRF compatibility
  // QNSE (Quasi-Normal Scale Elimination) spectral stable stability functions
  // Sukoriansky, S., B. Galperin, and V. Perov, 2005: Application of a New Spectral
  // Theory of Stably Stratified Turbulence to the Atmospheric Boundary Layer over
  // Sea Ice. Boundary-Layer Meteorology, 117, 231–257.
  // https://doi.org/10.1007/s10546-004-6848-4
  bool enable_qnse_stable_functions = false; // if true, use QNSE spectral stability functions
                                              // instead of linear Businger-Dyer form for stable conditions
  amrex::Real qnse_am = amrex::Real(2.5);    // QNSE momentum stability function coefficient (dimensionless)
  amrex::Real qnse_bm = amrex::Real(0.2);    // QNSE momentum stability function coefficient (dimensionless)
  amrex::Real qnse_ah = amrex::Real(2.5);    // QNSE heat stability function coefficient (dimensionless)
  amrex::Real qnse_bh = amrex::Real(0.2);    // QNSE heat stability function coefficient (dimensionless)
  // MYNN2.5 PBL model
  // TKE/QKE stuff
  bool advect_tke = true; // if MYNN2.5 PBL is used default is turb transport in
                          // Z-direction only
  bool diffuse_tke_3D = true; // if numerical diffusion is turned on
};
 
/**
 * @var TurbChoice::les_type
 * @brief Selected LES closure.
 * @var TurbChoice::Pr_t
 * @brief Turbulent Prandtl number.
 * @var TurbChoice::Pr_t_inv
 * @brief Inverse turbulent Prandtl number.
 * @var TurbChoice::Sc_t
 * @brief Turbulent Schmidt number.
 * @var TurbChoice::Sc_t_inv
 * @brief Inverse turbulent Schmidt number.
 * @var TurbChoice::Cs
 * @brief Smagorinsky model coefficient.
 * @var TurbChoice::smag2d
 * @brief Whether the 2-D Smagorinsky formulation is used.
 * @var TurbChoice::smag2d_slope_limiter
 * @brief Apply WRF's Smagorinsky-2D slope limiter (K_h / alpha or / alpha^2) on terrain-fitted meshes.
 * @var TurbChoice::smag2d_kh_cap
 * @brief Cap K_h <= smag2d_kh_cap * Delta_h for Smagorinsky-2D (m/s); 0 means no cap.
 * @var TurbChoice::smag2d_limits_requested
 * @brief Whether either Smagorinsky-2D limit was given as on for this level.
 * @var TurbChoice::Ce
 * @brief Deardorff dissipation coefficient.
 * @var TurbChoice::Ce_wall
 * @brief Wall value for the Deardorff dissipation coefficient.
 * @var TurbChoice::Ck
 * @brief Deardorff eddy-viscosity coefficient.
 * @var TurbChoice::Cmu0
 * @brief One-equation RANS Cmu0 coefficient.
 * @var TurbChoice::Cb
 * @brief One-equation RANS buoyancy coefficient.
 * @var TurbChoice::Rt_crit
 * @brief Critical turbulent Reynolds number.
 * @var TurbChoice::Rt_min
 * @brief Minimum turbulent Reynolds number.
 * @var TurbChoice::l_g_max
 * @brief Maximum geometric length scale.
 * @var TurbChoice::sigma_k
 * @brief TKE diffusivity coefficient denominator.
 * @var TurbChoice::theta_ref
 * @brief Reference potential temperature for stable stratification.
 * @var TurbChoice::strat_type
 * @brief Thermodynamic variable used for stability stratification.
 * @var TurbChoice::mix_isotropic
 * @brief Whether subgrid mixing uses isotropic length scales.
 * @var TurbChoice::use_Ri_correction
 * @brief Whether Richardson-number correction is applied.
 * @var TurbChoice::Ri_crit
 * @brief Critical Richardson number for stability correction.
 * @var TurbChoice::rans_type
 * @brief Selected RANS closure.
 * @var TurbChoice::dirichlet_k
 * @brief Whether TKE uses Dirichlet boundary treatment.
 * @var TurbChoice::pbl_type
 * @brief Selected PBL closure.
 * @var TurbChoice::pbl_mynn
 * @brief MYNN level-2.5 closure coefficients.
 * @var TurbChoice::pbl_mynn_level2
 * @brief MYNN level-2 closure coefficients for limiting.
 * @var TurbChoice::use_kturb
 * @brief Whether any turbulence model is active.
 * @var TurbChoice::use_keqn
 * @brief Whether a microscale TKE closure is active.
 * @var TurbChoice::use_pbl_tke
 * @brief Whether a mesoscale PBL TKE closure is active.
 * @var TurbChoice::use_tke
 * @brief Whether any TKE or QKE closure is active.
 * @var TurbChoice::init_tke_from_ustar
 * @brief Whether initial TKE/QKE profiles are based on surface friction velocity.
 * @var TurbChoice::tke_min
 * @brief Minimum TKE/QKE value.
 * @var TurbChoice::pbl_ysu_coriolis_freq
 * @brief Coriolis frequency used by YSU-family PBL schemes.
 * @var TurbChoice::pbl_ysu_use_consistent_coriolis
 * @brief Whether YSU uses the ERF Coriolis frequency.
 * @var TurbChoice::pbl_ysu_force_over_water
 * @brief Whether YSU is forced to use over-water behavior.
 * @var TurbChoice::pbl_ysu_land_Ribcr
 * @brief Critical bulk Richardson number over land for stable YSU conditions.
 * @var TurbChoice::pbl_ysu_unst_Ribcr
 * @brief Critical bulk Richardson number for unstable YSU conditions.
 * @var TurbChoice::enable_ysu_countergradient
 * @brief Whether YSU countergradient corrections are enabled.
 * @var TurbChoice::enable_ysu_terrain_pblh_floor
 * @brief Whether YSU applies a terrain-following PBL-height floor.
 * @var TurbChoice::enable_ysu_sat_limiter
 * @brief Whether YSU applies a saturation limiter to moisture countergradient terms.
 * @var TurbChoice::enable_ysu_topdown
 * @brief Whether YSU top-down mixing is enabled.
 * @var TurbChoice::enable_ysu_entrainment
 * @brief Whether YSU entrainment-layer parameterization is enabled.
 * @var TurbChoice::enable_ysu_cloud_pblh
 * @brief Whether YSU cloud-based PBL-height detection is enabled.
 * @var TurbChoice::ysu_qcloud_threshold
 * @brief Cloud liquid water threshold for YSUNew [kg/kg].
 * @var TurbChoice::ysu_moistvars
 * @brief Whether YSU applies turbulence to moisture variables.
 * @var TurbChoice::pbl_ysunew_highres_bounds
 * @brief Whether YSUNew applies high-resolution grid-dependent diffusivity bounds.
 * @var TurbChoice::enable_vh96_shear_correction
 * @brief Whether Vogelezang & Holtslag (1996) shear-correction term is enabled.
 * @var TurbChoice::vh96_shear_const_b
 * @brief Vogelezang & Holtslag (1996) shear-correction constant b.
 * @var TurbChoice::enable_pblh_smoothing
 * @brief Whether spatial smoothing of diagnosed PBLH is enabled.
 * @var TurbChoice::pblh_smoothing_passes
 * @brief Number of PBLH smoothing iterations to apply.
 * @var TurbChoice::pblh_smoothing_weight
 * @brief Center-cell weight in PBLH smoothing stencil (must be in [0,1]).
 * @var TurbChoice::pbl_mrf_coriolis_freq
 * @brief Coriolis frequency used by the MRF PBL scheme.
 * @var TurbChoice::pbl_mrf_Ribcr
 * @brief Critical bulk Richardson number for the MRF PBL scheme.
 * @var TurbChoice::pbl_mrf_const_b
 * @brief MRF constant used to compute PBL height.
 * @var TurbChoice::pbl_mrf_sf
 * @brief MRF surface flux value used to compute PBL height.
 * @var TurbChoice::enable_mrf_countergradient
 * @brief Whether MRF countergradient corrections are enabled.
 * @var TurbChoice::pbl_ib_aware
 * When true, the MRF and YSUNew schemes treat the first fluid cell above the
 * immersed solid of a column (buildings or terrain by immersed forcing) as
 * that column's surface: the bulk Richardson heights, the boundary-layer
 * depth and the K profile are measured from it, the diffusivities vanish
 * inside the solid, and the surface scales of such a column are a neutral
 * log law at its top (u* from the wind there and pbl_ib_z0, theta* = 0). The
 * boundary-layer height both schemes store (pblh, Lturb) stays the absolute
 * height above the domain bottom; with enable_pblh_smoothing the stencil runs
 * on that absolute height, and the smoothed value is kept at or above each
 * column's own floor (half a cell above its surface, at least 10 m). Off by
 * default; without immersed cells the results are identical. Not supported
 * with terrain-fitted coordinates.
 * @var TurbChoice::pbl_ib_z0
 * Roughness length [m] of that log law.
 * @var TurbChoice::enable_mrf_cloud_adjustment
 * @brief Whether MRF cloud-aware stability adjustments are enabled.
 * @var TurbChoice::pbl_mrf_highres_bounds
 * @brief Whether MRF applies high-resolution grid-dependent diffusivity bounds.
 * @var TurbChoice::enable_mrf_unbounded_vpert
 * @brief Whether MRF leaves VPERT unlimited by GAMCRT.
 * @var TurbChoice::pbl_mrf_use_zero_ri_extent
 * @brief Whether MRF uses the Ri=0 K-profile extent.
 * @var TurbChoice::advect_tke
 * @brief Whether TKE/QKE is advected.
 * @var TurbChoice::diffuse_tke_3D
 * @brief Whether three-dimensional numerical diffusion is applied to TKE/QKE.
 */
#endif