ERF_DataStruct.H

Go to the documentation of this file.

#ifndef ERF_DATA_STRUCT_H_
#define ERF_DATA_STRUCT_H_
 
#include <string>
#include <iostream>
#include <cmath>
 
#include <AMReX_ParmParse.H>
#include <AMReX_Print.H>
#include <AMReX_Gpu.H>
#include <AMReX_Geometry.H>
 
#include "ERF_Constants.H"
#include "ERF_IndexDefines.H"
#include "ERF_AdvStruct.H"
#include "ERF_DampingStruct.H"
#include "ERF_DiffStruct.H"
#include "ERF_EBStruct.H"
#include "ERF_SpongeStruct.H"
#include "ERF_TurbStruct.H"
#include "ERF_TurbPertStruct.H"
#include "ERF_RadStruct.H"
#include "ERF_ShocTransportStruct.H"
#include "ERF_CanopyTurbStruct.H"
 
/**
 * Enum of possible plotfile types.
 */
AMREX_ENUM(PlotFileType,
    None, Amrex, Netcdf
);
 
/**
 * @brief Indices for map scale factors on mass and face grids.
 */
enum MapFacType {
//  This version assumes isotropic
    m_x, u_x, v_x, num,
    m_y = 0, u_y = 1, v_y = 2
//  This version allows for non-isotropic
//  m_x, u_x, v_x,
//  m_y, u_y, v_y, num
};
 
/**
 * @brief Stress tensor component indices.
 */
enum TauType {
    tau11, tau22, tau33, tau12, tau13, tau23, tau21, tau31, tau32
};
 
/**
 * @brief Initial-condition source used to populate the ERF state.
 */
AMREX_ENUM(InitType,
    None, Input_Sounding, NCFile, WRFInput, Metgrid, Uniform, ConstantDensity, ConstantDensityLinearTheta,
    Isentropic, MoistBaseState, HindCast
);
 
/**
 * @brief Interpretation of input sounding profiles.
 */
AMREX_ENUM(SoundingType,
    ConstantDensity, Ideal, Isentropic, DryIsentropic
);
 
/**
 * @brief Large-scale atmospheric boundary-layer driver type.
 */
AMREX_ENUM(ABLDriverType,
    None, PressureGradient, GeostrophicWind
);
 
/**
 * @brief Multilevel coupling strategy.
 */
AMREX_ENUM(CouplingType,
    OneWay, TwoWay
);
 
/**
 * @brief Acoustic substepping strategy.
 */
AMREX_ENUM(SubsteppingType,
    None, Implicit
);
 
/**
 * @brief Two-stage scheme used by the anelastic time integrator.
 *
 * RK2 is the SSP (Heun) scheme: both stages advance a full timestep from the old
 * state and the second averages the two slow sources.  MidPoint advances only a
 * half timestep in the first stage, which makes the vertical diffusion second
 * order in time with an implicit solve in the first stage alone.
 */
AMREX_ENUM(AnelasticType,
    RK2, MidPoint
);
 
/**
 * @brief Vertical mesh representation.
 */
AMREX_ENUM(MeshType,
    ConstantDz, StretchedDz, VariableDz
);
 
/**
 * @brief Linear solver for the anelastic projection on a terrain-fitted mesh.
 *
 * gmres_fft is GMRES with ERF's TerrainPoisson operator and an FFT preconditioner (needs the
 * FFT build).  mlmg and gmres_mlmg use amrex::MLTerrainPoisson (AMReX-Codes/amrex#6030), the
 * same stencil with zero flux through the Neumann faces, solved by geometric multigrid
 * alone or as the preconditioner of GMRES; they need no FFT.
 */
AMREX_ENUM(TerrainPoissonSolver,
    gmres_fft, mlmg, gmres_mlmg
);
 
/**
 * @brief Terrain and immersed-boundary representation.
 */
AMREX_ENUM(TerrainType,
    None, StaticFittedMesh, MovingFittedMesh, EB, ImmersedForcing
);
 
/**
 * @brief Building representation.
 */
AMREX_ENUM(BuildingsType,
    None, ImmersedForcing
);
 
/**
 * @brief Moisture transport representation.
 */
AMREX_ENUM(MoistureModelType,
    Eulerian, Lagrangian, Undefined
);
 
/**
 * @brief Moisture and microphysics model.
 */
AMREX_ENUM(MoistureType,
    SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, Morrison, Morrison_NoIce, WSM6, WDM6, SuperDroplets, MoistNoCondensation, SBM, None
);
 
/**
 * @brief Whether an anelastic SuperDroplets configuration is unsupported.
 *
 * The current SuperDroplets implementation diagnoses thermodynamics through
 * the compressible equation of state and has no BaseState pressure input.
 */
inline bool anelastic_superdroplets_configuration_invalid (
    const MoistureType moisture_type,
    const bool any_anelastic) noexcept
{
    return any_anelastic && moisture_type == MoistureType::SuperDroplets;
}
 
/**
 * @brief Wind-farm model.
 */
AMREX_ENUM(WindFarmType,
    Fitch, EWP, SimpleAD, GeneralAD, None
);
 
/**
 * @brief Coordinate system used for wind-farm locations.
 */
AMREX_ENUM(WindFarmLocType,
    lat_lon, x_y, None
);
 
/**
 * @brief Land-surface model.
 */
AMREX_ENUM(LandSurfaceType,
    SLM, NOAHMP, None
);
 
/**
 * @brief Urban model.
 */
AMREX_ENUM(UrbanType,
    None
);
 
/**
 * @brief Radiation model.
 */
AMREX_ENUM(RadiationType,
    None, RRTMGP, Simple, TwoStream
);
 
/**
 * @brief Coordinate-axis selector.
 */
enum struct Coord {
    x, y, z
};
 
// These are used as integers so must be enum not enum struct
/**
 * @brief Rayleigh damping profile component indices.
 */
enum Rayleigh {
    ubar, vbar, wbar, thetabar, nvars
};
 
// These are used as integers so must be enum not enum struct
/**
 * @brief Sponge reference-state component indices.
 */
enum Sponge {
    ubar_sponge, vbar_sponge, nvars_sponge
};
 
static const std::string BoundaryFaceName[] = {"xlo", "ylo", "zlo", "xhi", "yhi", "zhi"};
 
/**
 * @brief The moisture data carried by the active microphysics scheme.
 *
 * This struct is the single source of truth for the question "which moisture
 * quantities does this run actually have?".  Everything downstream of the
 * microphysics -- source terms, turbulence closures, and every output path
 * (3D plotfiles, 2D plotfiles, subvolumes) -- answers that question by querying
 * this struct instead of enumerating MoistureType values, so a newly added
 * scheme only has to declare its layout once, in from_moisture_model() below.
 *
 * Two kinds of index live here:
 *
 *   1. Conserved-state components (qv ... nn).  Absolute component indices into
 *      the conserved state array, i.e. one of RhoQ1_comp ... RhoQ11_comp.  They
 *      name the *rho-weighted* species: state(i,j,k,qv) is rho*qv, so an output
 *      path must divide by Rho_comp to recover the mixing ratio.
 *
 *   2. qmoist diagnostic slots (rain_accum ... cond_rate).  Indices into the
 *      per-level qmoist vector that is filled from
 *      Microphysics::Get_Qmoist_Ptr(lev,idx).  Each scheme chooses that
 *      layout independently, so the slot numbers differ between schemes and are
 *      recorded per scheme below.
 *
 * In both cases the value is `absent` when the active scheme does not carry the
 * quantity, so `index >= 0` is the criterion for "this exists and is meaningful
 * to write".  That is deliberately *not* the same test as "the storage is
 * allocated": several schemes allocate a wider state than they integrate (the
 * Morrison class always allocates eleven moist components, so Morrison_NoIce
 * owns ice slots it never fills), and writing those slots to a plotfile would
 * publish untouched memory as if it were data.
 */
struct MoistureComponentIndices {
 
    // Sentinel values shared by both kinds of index above.
    static constexpr int absent = -1;              //!< not carried by the active scheme
    static constexpr int computed_from_state = -2; //!< carried, but derived at output
                                                   //!< time from the conserved state
                                                   //!< rather than read from a qmoist slot
 
    // ---------------------------------------------------------------------
    // Conserved-state components: mass mixing ratios (rho-weighted)
    // ---------------------------------------------------------------------
    int qv = absent;  //!< water vapor
    int qc = absent;  //!< cloud liquid water
    int qi = absent;  //!< cloud ice
    int qr = absent;  //!< rain
    int qs = absent;  //!< snow
    int qg = absent;  //!< graupel
 
    // ---------------------------------------------------------------------
    // Conserved-state components: number concentrations (rho-weighted)
    // ---------------------------------------------------------------------
    int nc = absent;  //!< cloud liquid water number
    int ni = absent;  //!< cloud ice number
    int nr = absent;  //!< rain number
    int ns = absent;  //!< snow number
    int ng = absent;  //!< graupel number
    // CCN / total aerosol number. Unlike every other number concentration above,
    // this has no companion mass species: it is an aerosol reservoir that
    // exchanges number with nc and nr through activation and evaporation rather
    // than a hydrometeor count. Kept distinct from ni, which is cloud ice number,
    // even though some schemes place the two in the same conserved-state slot.
    int nn = absent;  //!< CCN / total aerosol number
 
    // ---------------------------------------------------------------------
    // qmoist diagnostic slots
    //
    // These are *not* conserved-state components; they index the qmoist vector
    // handed out by the microphysics interface.  Slot numbers are scheme
    // specific, which is exactly why they belong in this per-scheme map.
    // ---------------------------------------------------------------------
    //
    // The accumulations are what reached the ground, in mm of liquid water:
    // rain, snow and graupel each hold that species alone, and precip_accum
    // holds everything that fell out -- those three plus whatever cloud water
    // and cloud ice the scheme sediments -- so it can exceed their sum.
    int rain_accum    = absent;  //!< accumulated surface rain
    int snow_accum    = absent;  //!< accumulated surface snow
    int graupel_accum = absent;  //!< accumulated surface graupel
    int precip_accum  = absent;  //!< accumulated surface precipitation, all species
    int rel_hum       = absent;  //!< relative humidity
    int cond_rate     = absent;  //!< condensation rate
 
    /**
     * @brief Construct the conserved-state component map.
     *
     * Prefer from_moisture_model() for production code: it fills in the qmoist
     * diagnostic slots as well.  This constructor exists for the callers -- unit
     * tests and kernels exercised in isolation -- that only need the state map.
     *
     * @param qv_comp Water vapor component index.
     * @param qc_comp Cloud liquid water component index.
     * @param qi_comp Cloud ice component index.
     * @param qr_comp Rain component index.
     * @param qs_comp Snow component index.
     * @param qg_comp Graupel component index.
     * @param nc_comp Cloud liquid water number concentration component index.
     * @param ni_comp Cloud ice number concentration component index.
     * @param nr_comp Rain number concentration component index.
     * @param ns_comp Snow number concentration component index.
     * @param ng_comp Graupel number concentration component index.
     * @param nn_comp CCN / total aerosol number concentration component index.
     */
    MoistureComponentIndices (int qv_comp, int qc_comp,
                              int qi_comp=absent,
                              int qr_comp=absent,
                              int qs_comp=absent,
                              int qg_comp=absent,
                              int nc_comp=absent,
                              int ni_comp=absent,
                              int nr_comp=absent,
                              int ns_comp=absent,
                              int ng_comp=absent,
                              int nn_comp=absent)
        : qv(qv_comp), qc(qc_comp), qi(qi_comp), qr(qr_comp), qs(qs_comp), qg(qg_comp),
                       nc(nc_comp), ni(ni_comp), nr(nr_comp), ns(ns_comp), ng(ng_comp),
                       nn(nn_comp) {}
 
    /**
     * @brief Construct an empty map, i.e. a dry run.
     */
    MoistureComponentIndices () = default;
 
    /**
     * @brief Build the complete index map -- conserved-state components and
     *        qmoist diagnostic slots -- for a moisture model.
     *
     * This is the one place a scheme declares what it carries.  The conserved
     * layout must match the component the scheme reads and writes in its
     * Copy_State_to_Micro / Copy_Micro_to_State pair, and the qmoist slots must
     * match the order the scheme installs in its MicVarMap.  A slot that a
     * scheme allocates but never fills is left `absent` on purpose: see the
     * class comment for why allocation is not the criterion here.
     *
     * @param moisture_type Active moisture/microphysics model.
     */
    static MoistureComponentIndices
    from_moisture_model (const MoistureType moisture_type)
    {
        MoistureComponentIndices mi;
 
        switch (moisture_type) {
 
        case MoistureType::None:
            // Dry: no moisture data of any kind.
            break;
 
        case MoistureType::MoistNoCondensation:
        case MoistureType::SAM_NoPrecip_NoIce:
        case MoistureType::Kessler_NoRain:
            // Vapor plus suspended cloud water, no precipitation and no
            // accumulation diagnostics.
            mi.qv = RhoQ1_comp;
            mi.qc = RhoQ2_comp;
            break;
 
        case MoistureType::SatAdj:
            // As above, but relative humidity is recovered from the state at
            // output time rather than stored: SatAdj publishes no qmoist arrays.
            mi.qv      = RhoQ1_comp;
            mi.qc      = RhoQ2_comp;
            mi.rel_hum = computed_from_state;
            break;
 
        case MoistureType::Kessler:
            // Warm rain in three moist components; qmoist is {rain_accum}.
            // Rain is the only species that falls, so the precipitation total
            // is the rain accumulation itself.
            mi.qv           = RhoQ1_comp;
            mi.qc           = RhoQ2_comp;
            mi.qr           = RhoQ3_comp;
            mi.rain_accum   = 0;
            mi.precip_accum = 0;
            break;
 
        case MoistureType::SBM:
            // SBM stores the distribution in an auxiliary state.  These two
            // compact liquid lanes are projections for host compatibility.
            mi.qv = RhoQ1_comp;
            mi.qc = RhoQ2_comp;
            mi.qr = RhoQ3_comp;
            break;
 
        case MoistureType::SAM_NoIce:
            // The SAM class always allocates six moist components, so rain sits
            // in the fourth slot even with the ice species switched off.
            mi.qv           = RhoQ1_comp;
            mi.qc           = RhoQ2_comp;
            mi.qr           = RhoQ4_comp;
            mi.rain_accum   = 0;
            mi.precip_accum = 3;
            break;
 
        case MoistureType::SAM:
        case MoistureType::WSM6:
            // Single-moment mixed phase; qmoist is
            // {rain_accum, snow_accum, graupel_accum, precip_accum}.
            mi.qv            = RhoQ1_comp;
            mi.qc            = RhoQ2_comp;
            mi.qi            = RhoQ3_comp;
            mi.qr            = RhoQ4_comp;
            mi.qs            = RhoQ5_comp;
            mi.qg            = RhoQ6_comp;
            mi.rain_accum    = 0;
            mi.snow_accum    = 1;
            mi.graupel_accum = 2;
            mi.precip_accum  = 3;
            break;
 
        case MoistureType::Morrison:
            // Two-moment mixed phase: six mass species followed by five number
            // concentrations.
            mi.qv          = RhoQ1_comp;
            mi.qc          = RhoQ2_comp;
            mi.qi          = RhoQ3_comp;
            mi.qr          = RhoQ4_comp;
            mi.qs          = RhoQ5_comp;
            mi.qg          = RhoQ6_comp;
            mi.nc          = RhoQ7_comp;
            mi.ni          = RhoQ8_comp;
            mi.nr          = RhoQ9_comp;
            mi.ns          = RhoQ10_comp;
            mi.ng          = RhoQ11_comp;
            mi.rain_accum    = 0;
            mi.snow_accum    = 1;
            mi.graupel_accum = 2;
            mi.precip_accum  = 3;
            break;
 
        case MoistureType::Morrison_NoIce:
            // Liquid-only Morrison: the class still allocates all eleven moist
            // components, but only vapor, cloud water and rain -- and the cloud
            // and rain numbers -- are integrated, so the frozen species and
            // their numbers stay absent here.
            mi.qv         = RhoQ1_comp;
            mi.qc         = RhoQ2_comp;
            mi.qr         = RhoQ4_comp;
            mi.nc           = RhoQ7_comp;
            mi.nr           = RhoQ9_comp;
            mi.rain_accum   = 0;
            mi.precip_accum = 3;
            break;
 
        case MoistureType::WDM6:
            // Double-moment warm phase on top of single-moment ice: nn is the
            // aerosol reservoir and takes the slot Morrison uses for ice number,
            // and there are no snow or graupel numbers.
            mi.qv          = RhoQ1_comp;
            mi.qc          = RhoQ2_comp;
            mi.qi          = RhoQ3_comp;
            mi.qr          = RhoQ4_comp;
            mi.qs          = RhoQ5_comp;
            mi.qg          = RhoQ6_comp;
            mi.nc          = RhoQ7_comp;
            mi.nn          = RhoQ8_comp;
            mi.nr          = RhoQ9_comp;
            mi.rain_accum    = 0;
            mi.snow_accum    = 1;
            mi.graupel_accum = 2;
            mi.precip_accum  = 3;
            break;
 
        case MoistureType::SuperDroplets:
            // Lagrangian microphysics that still carries the six Eulerian mass
            // species.  Its qmoist layout is its own (see the MicVarMap in
            // SuperDropletsMoist::Init): the graupel accumulation slot exists
            // but is never filled, so it stays absent.
            mi.qv         = RhoQ1_comp;
            mi.qc         = RhoQ2_comp;
            mi.qi         = RhoQ3_comp;
            mi.qr         = RhoQ4_comp;
            mi.qs         = RhoQ5_comp;
            mi.qg         = RhoQ6_comp;
            mi.cond_rate    = 3;
            mi.rel_hum      = 7;
            mi.rain_accum   = 8;
            mi.snow_accum   = 10;
            mi.precip_accum = 11;
            break;
 
        default:
            amrex::Abort("Unknown MoistureType in MoistureComponentIndices::from_moisture_model");
        }
 
        return mi;
    }
 
    /**
     * @brief An ordered list of conserved-state components forming an aggregate
     *        moisture output variable (qt, qn, qp, moist_density).
     *
     * Only the species the active scheme carries are appended, so a consumer can
     * sum the list without first asking which ones exist.
     */
    struct CompList {
        static constexpr int max_size = 6; //!< qv, qc, qi, qr, qs, qg
        int comp[max_size] = {absent, absent, absent, absent, absent, absent};
        int size = 0;
 
        //! Append a component, ignoring the ones the scheme does not carry.
        void append (int c) {
            if ( (c >= 0) && (size < max_size) ) { comp[size++] = c; }
        }
        [[nodiscard]] bool empty () const { return (size == 0); }
    };
 
    /**
     * @brief Components summed for total water, "qt": every mass mixing ratio.
     *
     * Number concentrations are excluded -- they are counts, not masses.
     */
    [[nodiscard]] CompList total_water_comps () const {
        CompList list;
        list.append(qv); list.append(qc); list.append(qi);
        list.append(qr); list.append(qs); list.append(qg);
        return list;
    }
 
    /**
     * @brief Components summed for non-precipitating water, "qn": vapor plus the
     *        suspended condensate.  Also the moist part of "moist_density".
     */
    [[nodiscard]] CompList nonprecipitating_comps () const {
        CompList list;
        list.append(qv); list.append(qc); list.append(qi);
        return list;
    }
 
    /**
     * @brief Components summed for precipitating water, "qp": the falling species.
     */
    [[nodiscard]] CompList precipitating_comps () const {
        CompList list;
        list.append(qr); list.append(qs); list.append(qg);
        return list;
    }
 
    /**
     * @brief Test whether a conserved-state component holds one of the moisture
     *        variables this scheme carries.
     *
     * Use this rather than enumerating moisture models when deciding whether a
     * "rhoQn" state component is meaningful for the active scheme.
     *
     * @param comp Absolute component index into the conserved state array.
     */
    [[nodiscard]] bool has_comp (int comp) const {
        return ( (comp >= 0) &&
                 ( (comp == qv) || (comp == qc) || (comp == qi) ||
                   (comp == qr) || (comp == qs) || (comp == qg) ||
                   (comp == nc) || (comp == ni) || (comp == nr) ||
                   (comp == ns) || (comp == ng) || (comp == nn) ) );
    }
 
    /**
     * @brief Test whether the active scheme carries rain, snow and graupel in the
     *        RhoQ4/RhoQ5/RhoQ6 slots assumed by erf_derreflectivity and
     *        erf_dermaxreflectivity.
     *
     * Those kernels read the precipitating species by hardcoded component, so this
     * checks the layout and not just the presence of the species.
     */
    [[nodiscard]] bool has_reflectivity_species () const {
        return ( (qr == RhoQ4_comp) && (qs == RhoQ5_comp) && (qg == RhoQ6_comp) );
    }
 
    /**
     * @brief Test whether the active scheme carries any moisture at all.
     *
     * Every moist scheme allocates water vapor, so this is equivalent to testing
     * moisture_type != MoistureType::None without needing the enum.
     */
    [[nodiscard]] bool has_moisture () const { return (qv >= 0); }
 
    /**
     * @brief The conserved-state component behind a single-species output name,
     *        or `absent` if the name is not a single species or the scheme does
     *        not carry it.
     *
     * An output path can use this to copy a species without repeating the
     * name-to-component mapping.
     *
     * @param name Derived or plot variable name.
     */
    [[nodiscard]] int comp_for_var (const std::string& name) const {
        if (name == "qv")     { return qv; }
        if (name == "qc")     { return qc; }
        if (name == "qi")     { return qi; }
        if (name == "qrain")  { return qr; }
        if (name == "qsnow")  { return qs; }
        if (name == "qgraup") { return qg; }
        if (name == "nc")     { return nc; }
        if (name == "ni")     { return ni; }
        if (name == "nr")     { return nr; }
        if (name == "ns")     { return ns; }
        if (name == "ng")     { return ng; }
        if (name == "nn")     { return nn; }
        return absent;
    }
 
    /**
     * @brief The qmoist slot behind a moist diagnostic output name.
     *
     * Returns `absent` when the scheme does not provide the diagnostic, and
     * `computed_from_state` when it does provide it but not through a qmoist
     * array -- in that case the caller must derive it from the conserved state.
     *
     * @param name Derived or plot variable name.
     */
    [[nodiscard]] int qmoist_index_for_var (const std::string& name) const {
        if (name == "rain_accum")        { return rain_accum; }
        if (name == "snow_accum")        { return snow_accum; }
        if (name == "graupel_accum")     { return graupel_accum; }
        if (name == "precip_accum")      { return precip_accum; }
        if (name == "rel_humidity")      { return rel_hum; }
        if (name == "condensation_rate") { return cond_rate; }
        return absent;
    }
 
    /**
     * @brief Whether this map decides a variable's availability, and if so
     *        whether the variable is available.
     */
    struct VarAvailability {
        bool governed  = false; //!< the name is a moisture output variable
        bool available = false; //!< the active scheme carries the data behind it
    };
 
    /**
     * @brief Classify an output variable name against this index map.
     *
     * This is the one table pairing output names with the data behind them.  It
     * is driven off the scheme's index map rather than an enumerated list of
     * moisture models, so it stays correct as schemes are added.  Names with no
     * moisture dependence come back not-governed, which lets a caller walk a
     * mixed list of variable names and only defer to this map where it applies.
     *
     * @param name Derived or plot variable name.
     */
    [[nodiscard]] VarAvailability query_var (const std::string& name) const {
 
        // Single species, taken straight from the conserved state
        if ( (name == "qv") || (name == "qc") || (name == "qi") ||
             (name == "qrain") || (name == "qsnow") || (name == "qgraup") ||
             (name == "nc") || (name == "ni") || (name == "nr") ||
             (name == "ns") || (name == "ng") || (name == "nn") ) {
            return {true, comp_for_var(name) >= 0};
        }
 
        // Moist diagnostics, either from a qmoist array or derived at output time
        if ( (name == "rain_accum") || (name == "snow_accum") ||
             (name == "graupel_accum") || (name == "precip_accum") ||
             (name == "rel_humidity") ||
             (name == "condensation_rate") ) {
            return {true, qmoist_index_for_var(name) != absent};
        }
 
        // Aggregates over several species
        if (name == "qt") { return {true, !total_water_comps().empty()}; }
        if (name == "qp") { return {true, !precipitating_comps().empty()}; }
 
        // Non-precipitating water is vapor plus the suspended condensate, and
        // moist density adds that same sum to the dry density
        if ( (name == "qn") || (name == "moist_density") ) {
            return {true, (qv >= 0) && (qc >= 0)};
        }
 
        // Column and thermodynamic diagnostics that read vapor only
        if ( (name == "qsat") || (name == "precipitable") ) {
            return {true, has_moisture()};
        }
 
        return {false, false};
    }
 
    /**
     * @brief Test whether the moist species behind a derived/plot variable name is
     *        carried by the active scheme.
     *
     * Names this map does not govern return true, so this can be used directly as
     * a filter over a mixed list of variable names.  Use query_var() instead when
     * the caller needs to distinguish "not a moisture variable" from "available".
     *
     * @param name Derived or plot variable name.
     */
    [[nodiscard]] bool has_derived_var (const std::string& name) const {
        const VarAvailability status = query_var(name);
        return (!status.governed) || status.available;
    }
};
 
/**
 * Container holding many of the algorithmic options and parameters
 */
 
struct SolverChoice {
  public:
    /**
     * @brief Read solver-wide algorithmic options from the input parameter database.
     * @param max_level Maximum AMR level configured for the run.
     * @param pp_prefix ParmParse prefix for the ERF input namespace.
     */
    void init_params (int max_level, std::string pp_prefix)
    {
        amrex::ParmParse pp(pp_prefix);
 
        bool bogus_bool;
        if (pp.query("use_terrain",bogus_bool) > 0)  {
            amrex::Error("The input use_terrain is deprecated.  Set terrain_type instead.");
        }
 
        if (pp.query("use_moist_background",bogus_bool) > 0)  {
            amrex::Error("The input use_moist_background is deprecated.  Set init_type = MoistBaseState instead.");
        }
 
        // Do we set map scale factors to myhalf instead of 1 for testing?
        pp.queryAdd("test_mapfactor", test_mapfactor);
 
        // For testing: if nonzero, the map scale factors on level 0 vary smoothly in space,
        // m = 1 + test_mapfactor_amp sin(2 pi x / Lx) cos(2 pi y / Ly), as if read from an init file
        // there; finer levels take theirs from the coarser level, as a level without a file does
        pp.queryAdd("test_mapfactor_amp", test_mapfactor_amp);
        if (test_mapfactor && test_mapfactor_amp != amrex::Real(0)) {
            amrex::Abort("erf.test_mapfactor and erf.test_mapfactor_amp cannot be used together");
        }
 
        // Which horizontal pressure gradient formulation to use with terrain fitted coords?
        // 0: dp/dx with dp/dz correction (default)
        // 1: gradient of vertically interpolated p, see Klemp 2011
        pp.queryAdd("gradp_type", gradp_type);
        AMREX_ALWAYS_ASSERT(gradp_type == 0 || gradp_type == 1);
 
        // For the lateral pressure gradient to be used in the momentum equation, should we
        // take the x- and y-derivatives of the perturbational pressure or the full pressure?
        pp.queryAdd("use_pert_pres_gradient", use_pert_pres_gradient);
 
        // What type of moisture model to use?
        moisture_type = MoistureType::None; // Default
        if (pp.query("moisture_type",moisture_type) > 0) {
            amrex::Error("The input moisture_type is deprecated.  Set moisture_model instead.");
        }
 
        pp.query_enum_case_insensitive("moisture_model",moisture_type);
 
        // Test-only SBM failure injections. They remain disabled unless an
        // integration test explicitly enables the corresponding option.
        pp.queryAdd("sbm_test_native_qc_write_fault", sbm_test_native_qc_write_fault);
        pp.queryAdd("sbm_test_donor_cfl_fault", sbm_test_donor_cfl_fault);
        if (moisture_type != MoistureType::SBM &&
            (sbm_test_native_qc_write_fault || sbm_test_donor_cfl_fault)) {
            amrex::Error("SBM test fault-injection options require moisture_model=SBM");
        }
 
        // Register the moisture data this scheme carries: conserved-state
        // components and qmoist diagnostic slots alike.  MoistureComponentIndices
        // is the one place that mapping lives, and every consumer -- source terms,
        // turbulence closures, and the plotfile/subvolume writers -- queries it
        // instead of testing moisture_type.
        moisture_indices = MoistureComponentIndices::from_moisture_model(moisture_type);
 
        if (moisture_type == MoistureType::SBM) {
            if (pp.contains("sbm_zero_transport_fixture")) {
                bool legacy_zero_transport = false;
                pp.get("sbm_zero_transport_fixture", legacy_zero_transport);
                if (legacy_zero_transport) {
                    amrex::Error("erf.sbm_zero_transport_fixture=true is obsolete; M3 advances the mapped spectral state, so remove this option");
                }
            }
            pp.queryAdd("sbm_nbins", sbm_nbins);
            pp.queryAdd("sbm_moment_mode", sbm_moment_mode);
            pp.queryAdd("sbm_cloud_rain_split", sbm_cloud_rain_split);
            pp.queryarr("sbm_edges", sbm_edges);
            pp.queryarr("sbm_pivots", sbm_pivots);
            pp.queryarr("sbm_fixture_initial_state", sbm_fixture_initial_state);
 
            if (sbm_nbins < 2 || (sbm_moment_mode != 1 && sbm_moment_mode != 2)) {
                amrex::Error("SBM requires sbm_nbins >= 2 and sbm_moment_mode equal to 1 or 2");
            }
            if (sbm_edges.empty()) {
                sbm_edges.resize(static_cast<std::size_t>(sbm_nbins + 1));
                sbm_pivots.resize(static_cast<std::size_t>(sbm_nbins));
                for (int b = 0; b <= sbm_nbins; ++b) {
                    const amrex::Real fraction = static_cast<amrex::Real>(b) /
                                                 static_cast<amrex::Real>(sbm_nbins);
                    sbm_edges[static_cast<std::size_t>(b)] =
                        std::pow(amrex::Real(10.0), amrex::Real(-18.0) + amrex::Real(6.0) * fraction);
                }
                for (int b = 0; b < sbm_nbins; ++b) {
                    sbm_pivots[static_cast<std::size_t>(b)] = std::sqrt(
                        sbm_edges[static_cast<std::size_t>(b)] * sbm_edges[static_cast<std::size_t>(b + 1)]);
                }
            } else {
                sbm_nbins = static_cast<int>(sbm_edges.size()) - 1;
                if (sbm_nbins < 2) amrex::Error("sbm_edges must define at least two bins");
                if (sbm_pivots.empty()) {
                    sbm_pivots.resize(static_cast<std::size_t>(sbm_nbins));
                    for (int b = 0; b < sbm_nbins; ++b) {
                        sbm_pivots[static_cast<std::size_t>(b)] = std::sqrt(
                            sbm_edges[static_cast<std::size_t>(b)] * sbm_edges[static_cast<std::size_t>(b + 1)]);
                    }
                }
            }
            if (sbm_pivots.size() != static_cast<amrex::Long>(sbm_nbins)) {
                amrex::Error("sbm_pivots must contain one value per SBM bin");
            }
            if (sbm_cloud_rain_split < 0) sbm_cloud_rain_split = sbm_nbins / 2;
            if (sbm_cloud_rain_split <= 0 || sbm_cloud_rain_split >= sbm_nbins) {
                amrex::Error("sbm_cloud_rain_split must be strictly inside the liquid bin range");
            }
            const int expected_sbm_components = sbm_nbins * (sbm_moment_mode == 2 ? 2 : 1);
            if (!sbm_fixture_initial_state.empty() &&
                sbm_fixture_initial_state.size() != static_cast<amrex::Long>(expected_sbm_components)) {
                amrex::Error("sbm_fixture_initial_state must match the full runtime SBM component count");
            }
            for (const amrex::Real value : sbm_fixture_initial_state) {
                if (!std::isfinite(value) || value < amrex::Real(0.0)) {
                    amrex::Error("sbm_fixture_initial_state values must be finite and nonnegative");
                }
            }
        }
 
        // Set a default for both dry and moist
        buoyancy_type.resize(max_level+1);
        for (int i = 0; i <= max_level; ++i) {
            buoyancy_type[i] = 1; // uses Rhoprime
        }
 
        if (moisture_type != MoistureType::None) {
            pp.queryAdd("moisture_tight_coupling",moisture_tight_coupling);
        }
 
        // Which expression (1,2/3 or 4) to use for buoyancy
        int default_buoyancy_type = buoyancy_type[0];
        read_int_string(max_level, "buoyancy_type", buoyancy_type, default_buoyancy_type);
 
        // What type of land surface model to use
        lsm_type = LandSurfaceType::None; // Default
        pp.query_enum_case_insensitive("land_surface_model",lsm_type);
 
        // What type of urban model to use and whether it is active at the given level
        urban_type = UrbanType::None;
        pp.query_enum_case_insensitive("urban_model", urban_type);
        read_int_string(max_level, "urban_lev_active", urban_enabled_lev, 1);
        for (int lev = 0; lev <= max_level; ++lev) {
            if (urban_enabled_lev[lev] != 0 && urban_enabled_lev[lev] != 1) {
                amrex::Error("urban_lev_active should be 0 or 1");
            }
        }
 
        // Will an external ocean coupler supply SST for part or all of the water
        // cells this run? The surface layer needs to know this it in time to select
        // ThetaCalcType::SURFACE_TEMPERATURE, which is long before the first
        // ApplyOceanSurfaceState call.
        pp.queryAdd("use_coupled_sst", use_coupled_sst);
 
        read_int_string(max_level, "is_land", is_land, 1);
        for (int lev = 0; lev <= max_level; ++lev) {
            if (is_land[lev] == 1) {
                amrex::Print() << "Level " << lev << " is land" << std::endl;
            } else if (is_land[lev] == 0) {
                amrex::Print() << "Level " << lev << " is water" << std::endl;
            } else {
                amrex::Error("is_land should be 0 or 1");
            }
        }
 
        // What type of radiation model to use
        rad_type = RadiationType::None; // Default
        pp.query_enum_case_insensitive("radiation_model", rad_type);
 
        // Verify that radiation model cannot be RRTMGP if ERF was not compiled with RRTMGP
#ifndef ERF_USE_RRTMGP
        if (rad_type == RadiationType::RRTMGP)
        {
            amrex::Error("ERF was not compiled with RRTMGP enabled!");
        }
#endif
 
        // Is the terrain none, static or moving?
        std::string terrain_type_temp = "";
        pp.queryAdd("terrain_type", terrain_type_temp);
        if (terrain_type_temp == "Moving") {
            amrex::Error("erf.terrain_type = Moving is no longer supported; please replace Moving by MovingFittedMesh");
        } else if (terrain_type_temp == "Static") {
            amrex::Error("erf.terrain_type = Static is no longer supported; please replace Static by StaticFittedMesh");
        } else {
            pp.query_enum_case_insensitive("terrain_type",terrain_type);
        }
 
        // Whether the terrain is flat - this allows planar averages (lsf, nudging) with StaticFittedMesh
        pp.queryAdd("flat_terrain", flat_terrain);
 
        // Get buildings type
        std::string buildings_type_temp = "";
        pp.queryAdd("buildings_type", buildings_type_temp);
        if (buildings_type_temp == "ImmersedForcing") {
            buildings_type = BuildingsType::ImmersedForcing;
        }
        pp.queryAdd("buildings_fc_factories", buildings_fc_factories);
 
        //
        // Read the init_type here to make sure we correctly set the mesh and terrain types
        //
        std::string init_type_temp_string;
 
        int found = pp.queryAdd("init_type",init_type_temp_string);
 
        if ( (init_type_temp_string == "Real") || (init_type_temp_string == "real") ) {
            amrex::Error("erf.init_type = Real is deprecated; please replace Real by WRFInput");
        } else if ( (init_type_temp_string == "Ideal") || (init_type_temp_string == "ideal") ) {
            amrex::Error("erf.init_type = Ideal is deprecated; please replace Ideal by WRFInput");
        } else {
            pp.query_enum_case_insensitive("init_type",init_type);
            use_real_bcs = ( (init_type == InitType::WRFInput) || (init_type == InitType::Metgrid) );
        }
 
        if ( (init_type == InitType::WRFInput) || (init_type == InitType::Metgrid) ) {
            if (terrain_type != TerrainType::StaticFittedMesh) {
                amrex::Error("Only terrain_type = StaticFittedMesh are allowed with init_type = WRFInput or Metgrid");
            }
        }
 
        if (init_type == InitType::MoistBaseState) {
            if (moisture_type == MoistureType::None) {
                amrex::Error("Makes no sense to have moist base state with no moisture model");
            }
        }
 
        if (init_type == InitType::WRFInput) {
            if (moisture_type == MoistureType::None) {
                amrex::Error("Can't have moisture_type = None with init_type = WRFInput");
            }
 
            // NetCDF wrfbdy lateral boundary file
            std::string nc_bdy_file_temp_string;
            bool has_bdy = pp.queryAdd("nc_bdy_file", nc_bdy_file_temp_string);
            if (!has_bdy) use_real_bcs = false;
 
            bool use_real_bcs_temp = use_real_bcs;
            pp.queryAdd("use_real_bcs", use_real_bcs_temp);
            if (use_real_bcs && !use_real_bcs_temp) {
                use_real_bcs = false;
            }
        }
 
        if (found == 0 || init_type == InitType::None) {
            amrex::Print() << "init_type must now be set.  The options are " << std::endl;
            amrex::Print() << " Input_Sounding, NCFile, WRFInput, Metgrid, Uniform, " << std::endl;
            amrex::Print() << " ConstantDensity, Isentropic, MoistBaseState, or HindCast " << std::endl;
            amrex::Error("Please add a string for init_type to your inputs file or command line");
        }
 
        // How to interpret input_sounding
        if (init_type == InitType::Input_Sounding) {
            pp.query_enum_case_insensitive("sounding_type",sounding_type);
        }
 
        if (terrain_type == TerrainType::StaticFittedMesh ||
            terrain_type == TerrainType::MovingFittedMesh) {
            mesh_type = MeshType::VariableDz;
        }
 
        if (flat_terrain && terrain_type != TerrainType::StaticFittedMesh) {
            amrex::Error("erf.flat_terrain requires terrain_type = StaticFittedMesh");
        }
 
        // Refined levels over moving terrain are not supported: erf_substep_MT does not take the
        // fluxes through coarse/fine faces from the coarse level (cf_flux_history_active)
        if (terrain_type == TerrainType::MovingFittedMesh && max_level > 0) {
            amrex::Error("erf.terrain_type = MovingFittedMesh requires amr.max_level = 0");
        }
 
        // Immersed-forcing buildings need the terrain blanking that init_zphys builds, and
        // init_zphys returns before building it when the terrain is EB
        if (terrain_type == TerrainType::EB && buildings_type == BuildingsType::ImmersedForcing) {
            amrex::Error("erf.buildings_type = ImmersedForcing is not supported with erf.terrain_type = EB");
        }
 
        if ( (moisture_type == MoistureType::SAM                ||
              moisture_type == MoistureType::SAM_NoIce          ||
              moisture_type == MoistureType::SAM_NoPrecip_NoIce   ) &&
             (mesh_type == MeshType::VariableDz) )
        {
            amrex::Abort("SAM is not correct with variable dz -- choose another moisture model");
        }
 
        // queryAdd puts the default in the table, so ask whether the ratio was given first
        const bool have_stretching_ratio = pp.contains("grid_stretching_ratio");
        pp.queryAdd("grid_stretching_ratio", grid_stretching_ratio);
        if (grid_stretching_ratio != 0) {
            AMREX_ASSERT_WITH_MESSAGE((grid_stretching_ratio >= one),
                                      "The grid stretching ratio must be greater than 1");
        }
        // erf.initial_dz on its own asks for the stretching ratio to be solved for, so that
        //    a first layer that thick stretches geometrically to close on prob_hi[2].  An
        //    explicit erf.grid_stretching_ratio = 0 still means no stretching at all.
        const bool solve_stretching_ratio = !have_stretching_ratio && pp.contains("initial_dz");
        if (grid_stretching_ratio >= 1 || solve_stretching_ratio) {
            if (terrain_type == TerrainType::None) {
                terrain_type = TerrainType::StaticFittedMesh;
            }
            if (mesh_type == MeshType::ConstantDz) {
                mesh_type = MeshType::StretchedDz;
            }
            pp.get("initial_dz", dz0);
        }
 
        // The stretched levels start at geometry.prob_lo[2]; erf.zsurface used to set the
        //    bottom separately, and a flat surface did not follow it (#4442)
        if (pp.contains("zsurface")) {
            amrex::Abort("erf.zsurface is no longer supported: a stretched vertical grid now starts "
                         "at geometry.prob_lo[2], so set the bottom of the domain there instead");
        }
 
        int n_zlevels = pp.countval("terrain_z_levels");
        if (n_zlevels > 0)
        {
            if (terrain_type == TerrainType::None) {
                terrain_type = TerrainType::StaticFittedMesh;
            }
            if (mesh_type == MeshType::ConstantDz) {
                mesh_type = MeshType::StretchedDz;
            }
        }
 
        // Use lagged_delta_rt in the fast integrator?
        pp.queryAdd("use_lagged_delta_rt", use_lagged_delta_rt);
 
        // Average the input z-face heights onto the nodes, or reconstruct our own
        pp.queryAdd("avg_grid_faces_to_nodes",avg_grid_faces_to_nodes);
 
        // Build the ERF vertical grid from the file's own mean layer-thickness profile,
        // or from a single geometric stretch anchored on the file's first layer?
        pp.queryAdd("wrfinput_zlevels_from_file",wrfinput_zlevels_from_file);
 
        // Catch the old name for the flag above so the error says what to do
        // rather than silently ignoring the setting.
        if (pp.contains("use_wrf_height_grid")) {
            amrex::Abort("erf.use_wrf_height_grid has been renamed. "
                         "Set erf.avg_grid_faces_to_nodes instead; it carries "
                         "the same meaning (average the input z-face heights "
                         "onto the nodes rather than reconstructing them).");
        }
 
        // Interpolate atmospheric state from coarse level instead of reading from wrfinput
        // when creating a finer level from coarse at a later time?
        pp.queryAdd("interp_atmos_from_coarse", interp_atmos_from_coarse);
 
        // Rebalance wrf state?
        pp.queryAdd("rebalance_wrf_input", rebalance_wrf_input);
 
        // Must rebalance if the nodal heights were reconstructed rather than averaged
        if (!avg_grid_faces_to_nodes && !rebalance_wrf_input) {
            amrex::Print() << "Must rebalance when not averaging the input height grid. Setting rebalance_wrf_input to true \n";
            rebalance_wrf_input = true;
        }
 
        // These default to true but are used for unit testing
        pp.queryAdd("use_gravity", use_gravity);
 
        // Initializing from Metgrid or WRFInput without gravity makes no sense
        if ( !use_gravity && ( (init_type == InitType::WRFInput) || (init_type == InitType::Metgrid) ) ) {
            amrex::Print() << "Setting use_gravity to true due to init_type = WRFInput / Metgrid " << std::endl;
            use_gravity = true;
        }
 
        gravity = use_gravity? CONST_GRAV: zero;
 
        pp.queryAdd("c_p", c_p);
        rdOcp = R_d / c_p;
 
        // *******************************************************************************
        // Read anelastic etc and over-ride if necessary
        // *******************************************************************************
 
        read_int_string(max_level, "anelastic", anelastic, 0);
        const bool fixed_density_given =
            read_int_string(max_level, "fixed_density", fixed_density, 0);
        read_int_string(max_level, "project_initial_velocity", project_initial_velocity, 0);
 
        // The anelastic equations hold the density at its base-state value, so
        // fixed_density is forced to 1 at every anelastic level; a user value of 0
        // there is ignored, and we say so rather than drop it silently.
        bool any_anelastic = false;
        bool any_compress  = false;
        for (int i = 0; i <= max_level; ++i) {
            if (anelastic[i] == 1) {
                if (fixed_density_given && fixed_density[i] == 0) {
                    amrex::Print() << "WARNING: erf.fixed_density = 0 is ignored at level " << i
                                   << " because erf.anelastic = 1 there; the anelastic equations"
                                   << " always hold the density fixed in time." << std::endl;
                }
                project_initial_velocity[i] = 1;
                fixed_density[i]            = 1; // Always fixed when anelastic
                buoyancy_type[i]            = 3; // (This isn't actually used when anelastic is set)
                any_anelastic               = true;
            } else {
                any_compress                = true;
            }
        }
 
        if (anelastic_superdroplets_configuration_invalid(moisture_type, any_anelastic)) {
            amrex::Error("erf.moisture_model = SuperDroplets is not supported with "
                         "erf.anelastic = 1: the Super-Droplet path currently uses "
                         "the compressible equation of state. Use compressible "
                         "dynamics or a supported Eulerian moisture model.");
        }
 
        // Want to have different immersed forcing defaults depending on anelastic or fully compressible.
        // We should have different starting values because the stability of the method is dependent on dt.
        // dt is coarser for anelastic --> need to loosen stiffness of immersed forcing.
        if (any_anelastic) {
            immersed_forcing_substep  = false;
            if_Cd_momentum            = amrex::Real(50.0);
            if_Cd_scalar              = amrex::Real(5.0);
        } else {
            immersed_forcing_substep  = true;
            if_Cd_momentum            = amrex::Real(500.0);
            if_Cd_scalar              = amrex::Real(50.0);
        }
 
        // *******************************************************************************
        // Read anelastic_type and allow for different values at each level
        // *******************************************************************************
        anelastic_type.resize(max_level+1);
 
        for (int i = 0; i <= max_level; i++) {
            anelastic_type[i] = AnelasticType::RK2;
        }
 
        for (int i = 0; i <= max_level; i++) {
            query_one_or_per_level_enum_case_insensitive(
                pp, "anelastic_type", anelastic_type[i], i, max_level);
        }
 
        if (pp.contains("anelastic_type") && !any_anelastic) {
            amrex::Print() << "Ignoring erf.anelastic_type because no level is anelastic" << std::endl;
            for (int i = 0; i <= max_level; i++) { anelastic_type[i] = AnelasticType::RK2; }
        }
 
        // *******************************************************************************
        // Read substepping_type and allow for different values at each level
        // *******************************************************************************
        substepping_type.resize(max_level+1);
 
        for (int i = 0; i <= max_level; i++) {
            substepping_type[i] = SubsteppingType::Implicit;
        }
 
        for (int i = 0; i <= max_level; i++) {
            query_one_or_per_level_enum_case_insensitive(
                pp, "substepping_type", substepping_type[i], i, max_level);
        }
 
        pp.queryAdd("substepping_diag", substepping_diag);
 
        pp.queryAdd("beta_s", beta_s);
 
        // Divergence damping coefficient for the forward projection of (rho theta)
        pp.queryAdd("beta_d", beta_d);
        if (!(beta_d >= amrex::Real(0.0) && beta_d <= amrex::Real(1.0))) {
            amrex::Abort("erf.beta_d must lie in [0, 1]; got " + std::to_string(beta_d));
        }
 
        // Spatially varying beta_s on the w-faces (implicit substepping only):
        //    beta_s = max(beta_s_min, min(beta_s_max, h_xi^2 + h_eta^2))
        pp.queryAdd("use_variable_beta_s", use_variable_beta_s);
        pp.queryAdd("beta_s_min", beta_s_min);
        pp.queryAdd("beta_s_max", beta_s_max);
        if (use_variable_beta_s) {
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(beta_s_min <= beta_s_max,
                "erf.beta_s_min must be less than or equal to erf.beta_s_max");
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(beta_s_min >= -one && beta_s_max <= one,
                "erf.beta_s_min and erf.beta_s_max must lie in [-1, 1]");
        }
 
        // *******************************************************************************
        // Error check on deprecated input
        // *******************************************************************************
        int nvals_old = pp.countval("no_substepping");
        if (nvals_old > 0) {
            amrex::Error("The no_substepping flag is deprecated -- set substepping_type instead");
        }
 
        pp.queryAdd("ncorr", ncorr);
        // Wall distance on a terrain-fitted mesh (RANS): "poisson" (Tucker 2003
        // differential-equation distance) or "terrain_height" (height above the
        // local surface projected on its normal, no linear solve)
        pp.queryAdd("wall_dist_type", wall_dist_type);
        if (wall_dist_type != "poisson" && wall_dist_type != "terrain_height") {
            amrex::Error("erf.wall_dist_type must be poisson or terrain_height");
        }
        pp.queryAdd("poisson_abstol", poisson_abstol);
        pp.queryAdd("poisson_reltol", poisson_reltol);
        // Solver for the projection on a terrain-fitted mesh: FFT-preconditioned GMRES (default),
        // or amrex::MLTerrainPoisson with multigrid, alone or under GMRES
        pp.query_enum_case_insensitive("terrain_poisson_solver", terrain_poisson_solver);
        if (pp.contains("terrain_poisson_solver") &&
            (mesh_type != MeshType::VariableDz || terrain_type == TerrainType::EB)) {
            amrex::Print() << "WARNING: erf.terrain_poisson_solver is ignored: it applies only to the "
                           << "projection on a terrain-fitted mesh (erf.terrain_type = StaticFittedMesh "
                           << "or MovingFittedMesh)" << std::endl;
        }
#ifdef AMREX_USE_FLOAT
        poisson_abstol = amrex::max(poisson_abstol,amrex::Real(1e-6));
        poisson_reltol = amrex::max(poisson_reltol,amrex::Real(1e-6));
#endif
        for (int lev = 0; lev <= max_level; lev++) {
            if (anelastic[lev] != 0)
            {
                substepping_type[lev] = SubsteppingType::None;
            }
        }
 
        pp.queryAdd("force_stage1_single_substep", force_stage1_single_substep);
 
        // Include Coriolis forcing?
        pp.queryAdd("use_coriolis", use_coriolis);
        pp.queryAdd("variable_coriolis", variable_coriolis);
 
        // Include four stream radiation approximation
        pp.queryAdd("four_stream_radiation", four_stream_radiation);
 
        // flags for whether to apply other source terms in substep only
        pp.queryAdd("immersed_forcing_substep", immersed_forcing_substep); // apply immersed forcing source terms in substep only
        pp.queryAdd("forest_substep", forest_substep);                     // apply canopy-related source terms in substep only
 
        // Fixed-leaf-temperature canopy heat exchange
        pp.query("forest_biophysics",       forest_biophysics);
        pp.query("forest_biophysics_heat",  forest_biophysics_heat);
        pp.query("forest_leaf_theta_fixed", forest_leaf_theta_fixed);
 
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            !forest_biophysics_heat || forest_biophysics,
            "forest_biophysics_heat requires forest_biophysics = true");
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            !(forest_biophysics && forest_biophysics_heat) ||
                forest_leaf_theta_fixed > zero,
            "forest_biophysics_heat requires a positive forest_leaf_theta_fixed");
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            forest_leaf_theta_fixed <= zero || forest_biophysics_heat,
            "forest_leaf_theta_fixed requires forest_biophysics_heat = true");
 
        // Canopy terms of the k-eqn RANS model (erf.forest_canopy_*); the
        // checks that need the forest and the turbulence model are in
        // ERF::ReadParameters, once both are known
        canopyTurb.init_params(pp);
 
        // immersed forcing parameters
        pp.queryAdd("if_Cd_momentum", if_Cd_momentum);
        pp.queryAdd("if_Cd_scalar", if_Cd_scalar);
        pp.queryAdd("if_implicit_drag", if_implicit_drag); // flag for implicit vs. explicit drag formulation
        pp.queryAdd("if_z0", if_z0);
        pp.queryAdd("if_surf_temp_flux", if_surf_temp_flux);
        pp.queryAdd("if_init_surf_temp", if_init_surf_temp);
 
        pp.queryAdd("if_surf_heating_rate", if_surf_heating_rate);
 
        // Modify rate to be in units of K / s rather than K / hr
        if_surf_heating_rate /= amrex::Real(3600.0);        // [K/s]
 
        pp.queryAdd("if_Olen", if_Olen_in);
        pp.queryAdd("if_use_most",if_use_most);
        pp.queryAdd("if_snap_partial_cells",if_snap_partial_cells);
        // The snapped staircase puts full-strength forcing on every wall face,
        // whose explicit relaxation rate (Cd / dx per substep) sits at the
        // stability limit and fails where faces stack at rims and corners; the
        // point-implicit form holds at any rate, so the snap uses it.
        if (if_snap_partial_cells && !if_implicit_drag) {
            amrex::Print() << "erf.if_snap_partial_cells: using the point-implicit drag (erf.if_implicit_drag = true)\n";
            if_implicit_drag = true;
        }
        pp.queryAdd("if_stability_correction",if_stability_correction);
 
        // Wall law of the terrain immersed forcing (erf-model/ERF#4134): "legacy" (the default)
        // relaxes the top partial cell toward the log law at fixed 0.5/1.5 dz;
        // "fraction_stress" applies the wall stress of ERF_ImmersedWallCell.H in the wall cell
        // and closes the wall face to diffusion (form 7 of the kynema-sgf immersed terrain)
        std::string if_wall_form_in = "legacy";
        pp.queryAdd("if_wall_form", if_wall_form_in);
        if (if_wall_form_in == "fraction_stress") {
            if_fraction_stress = true;
        } else if (if_wall_form_in != "legacy") {
            amrex::Abort("erf.if_wall_form must be legacy or fraction_stress, not " + if_wall_form_in);
        }
        // Spacing of the face above the wall cell (fraction_stress only): "log" (default) gives
        // that face the relative error of flat ground's first interior face for a log profile,
        // "physical" the centroid spacing, "grid" the cell height (ib_wall::face_above_factor)
        // (the keys below are accepted at their default values with the legacy law, so one deck
        // can run both laws; any other value belongs to fraction_stress)
        std::string if_wall_face_spacing_in = "log";
        pp.queryAdd("if_wall_face_spacing", if_wall_face_spacing_in);
        if (if_wall_face_spacing_in != "log" && !if_fraction_stress) {
            amrex::Abort("erf.if_wall_face_spacing belongs to erf.if_wall_form = fraction_stress");
        }
        if      (if_wall_face_spacing_in == "grid")     { if_wall_face_spacing = 0; }
        else if (if_wall_face_spacing_in == "physical") { if_wall_face_spacing = 1; }
        else if (if_wall_face_spacing_in == "log")      { if_wall_face_spacing = 2; }
        else { amrex::Abort("erf.if_wall_face_spacing must be log, physical or grid"); }
        // Wall value of the k-eqn TKE in the fraction-stress wall cell, (u*^3 + kappa B d1)^(2/3)
        // / Cmu0^2 (the value erf.dirichlet_k holds over flat ground): with a time factor of 1
        // (the default) the wall cell is held at it after every stage and in the implicit TKE
        // solve, as flat ground is; a larger factor relaxes toward it as a source over that many
        // steps (the kynema-sgf form)
        pp.queryAdd("if_wall_tke", if_wall_tke);
        pp.queryAdd("if_wall_tke_time_factor", if_wall_tke_time_factor);
        if ((if_wall_tke || if_wall_tke_time_factor != one) && !if_fraction_stress) {
            amrex::Abort("erf.if_wall_tke = true and erf.if_wall_tke_time_factor != 1 belong to erf.if_wall_form = fraction_stress");
        }
        if (!(if_wall_tke_time_factor >= one)) {
            amrex::Abort("erf.if_wall_tke_time_factor must be >= 1 (steps over which the wall TKE is relaxed)");
        }
        // Partial wall cells of the fraction-stress law dissipate TKE with the length scale of flat
        // ground's first cell (kynema-sgf KLAxell.flat_wall_dissipation); on by default with it.
        // query, not queryAdd: the default depends on erf.if_wall_form, so storing it would make a
        // later init_params in the same process (the unit tests) see it as given (cf. #4940)
        if_flat_wall_dissipation = if_fraction_stress;
        pp.query("if_flat_wall_dissipation", if_flat_wall_dissipation);
        if (if_flat_wall_dissipation && !if_fraction_stress) {
            amrex::Abort("erf.if_flat_wall_dissipation = true belongs to erf.if_wall_form = fraction_stress");
        }
        if (if_fraction_stress) {
            if (terrain_type != TerrainType::ImmersedForcing) {
                amrex::Abort("erf.if_wall_form = fraction_stress needs erf.terrain_type = ImmersedForcing");
            }
            if (buildings_type == BuildingsType::ImmersedForcing) {
                amrex::Abort("erf.if_wall_form = fraction_stress is not available with erf.buildings_type = ImmersedForcing");
            }
            if (!if_use_most) {
                amrex::Abort("erf.if_wall_form = fraction_stress is a wall law: it needs erf.if_use_most = true");
            }
            if (mesh_type != MeshType::ConstantDz) {
                amrex::Abort("erf.if_wall_form = fraction_stress needs a constant-dz mesh");
            }
            // The wall heat flux and the wall TKE are written by the slow scalar source; the
            // compressible fast substep (erf.immersed_forcing_substep) does not carry them
            if (immersed_forcing_substep) {
                amrex::Abort("erf.if_wall_form = fraction_stress runs with the anelastic model (erf.anelastic = 1 on "
                             "every level); it does not support erf.immersed_forcing_substep = true");
            }
            // One surface condition (heat flux, surface temperature - optionally changing at
            // erf.if_surf_heating_rate - or Obukhov length; none is neutral, checked below); the
            // MOST stability functions always apply, as on the flat surface layer, so
            // erf.if_stability_correction (a legacy switch) must not be set
            if (if_surf_heating_rate != zero && !(if_init_surf_temp > zero)) {
                amrex::Abort("erf.if_surf_heating_rate needs erf.if_init_surf_temp");
            }
            if (if_stability_correction) {
                amrex::Abort("erf.if_wall_form = fraction_stress always applies the stability functions; "
                             "erf.if_stability_correction is for the legacy wall law");
            }
        }
        pp.queryAdd("if_ws_floor",if_ws_floor);
        pp.queryAdd("if_damp_alpha",if_damp_alpha);
        // Opt-in WRF-style bounds on the wall law (erf-model/ERF#4016): the wind floor behind the
        // friction velocity of its stability estimate, and the cap on psi_m in the momentum law as
        // a fraction of ln(z/z0); 0 and 1 (the defaults) turn them off
        pp.queryAdd("if_stability_wind_floor",if_stability_wind_floor);
        pp.queryAdd("if_psi_cap_factor",if_psi_cap_factor);
        if (!(if_stability_wind_floor >= zero)) {
            amrex::Abort("erf.if_stability_wind_floor must be >= 0 [m/s] (0 leaves the wind as it is)");
        }
        if (!(if_psi_cap_factor > zero && if_psi_cap_factor <= one)) {
            amrex::Abort("erf.if_psi_cap_factor must be in (0, 1] (1 caps psi at ln(z/z0))");
        }
        // The cell average down over immersed-forcing levels (erf-model/ERF#4918): the plain
        // average (the default), or one weighted by the fluid mass of the fine cells
        pp.queryAdd("if_average_down", if_average_down);
        if (if_average_down != "plain" && if_average_down != "fluid_mass") {
            amrex::Abort("erf.if_average_down must be plain or fluid_mass");
        }
        // (queryAdd puts the default into the table, so the test is on the value, not on contains)
        if (if_average_down == "fluid_mass" &&
            terrain_type != TerrainType::ImmersedForcing && buildings_type != BuildingsType::ImmersedForcing)
        {
            amrex::Abort("erf.if_average_down = fluid_mass needs erf.terrain_type = ImmersedForcing "
                         "or erf.buildings_type = ImmersedForcing");
        }
        // Warn at start-up where a lateral coarse/fine boundary crosses the immersed surface so
        // that the levels disagree: spread of the fine solid fractions under a coarse edge cell
        // at least this (0 turns the check off)
        pp.queryAdd("if_cf_mismatch_warning", if_cf_mismatch_warning);
        if (!(if_cf_mismatch_warning >= zero && if_cf_mismatch_warning <= one)) {
            amrex::Abort("erf.if_cf_mismatch_warning must be in [0, 1] (a spread of solid fractions; 0 turns it off)");
        }
        // Log-law share of the wall coarse face's flux among its fine sub-faces at the lateral
        // coarse/fine faces (erf-model/ERF#4920): the immersed surface or a surface-layer bottom
        // (checked in check_params); constant dz, not EB terrain. The roughness length defaults
        // to erf.if_z0 with immersed forcing and to the inputs' erf.most.z0 or erf.zlo.most.z0
        // (0.1 m when neither is given) otherwise.
        pp.queryAdd("cf_loglaw_fill", cf_loglaw_fill);
        const bool z0_given = pp.query("cf_loglaw_z0", cf_loglaw_z0);
        if (z0_given && !cf_loglaw_fill) {
            amrex::Abort("erf.cf_loglaw_z0 belongs to erf.cf_loglaw_fill = true");
        }
        if (cf_loglaw_fill) {
            if (mesh_type != MeshType::ConstantDz) {
                amrex::Abort("erf.cf_loglaw_fill needs a constant-dz mesh");
            }
            if (terrain_type == TerrainType::EB) {
                amrex::Abort("erf.cf_loglaw_fill is not implemented for erf.terrain_type = EB "
                             "(its wall is not the bottom of the domain)");
            }
            if (!z0_given) {
                if (terrain_type == TerrainType::ImmersedForcing || buildings_type == BuildingsType::ImmersedForcing) {
                    cf_loglaw_z0 = if_z0;
                } else {
                    // the roughness given in the inputs for the bottom surface layer
                    // (erf.most.z0, or erf.zlo.most.z0 when the faces are set one by one); a
                    // roughness a land model supplies later is not seen here
                    cf_loglaw_z0 = amrex::Real(0.1);
                    pp.query("most.z0", cf_loglaw_z0);
                    pp.query("zlo.most.z0", cf_loglaw_z0);
                }
            }
            if (!(cf_loglaw_z0 > zero)) {
                amrex::Abort("erf.cf_loglaw_z0 must be > 0 [m] (the roughness length of the log-law coarse/fine fill)");
            }
        }
        // Say so when a bound is set and the wall law's stability estimate is in use
        if ((terrain_type == TerrainType::ImmersedForcing || buildings_type == BuildingsType::ImmersedForcing) &&
            (if_use_most || if_surf_temp_flux != amrex::Real(1e-8) || if_Olen_in != amrex::Real(1e-8)) &&
            (if_psi_cap_factor < one || if_stability_wind_floor > zero))
        {
            amrex::Print() << "Immersed wall law: erf.if_psi_cap_factor = " << if_psi_cap_factor
                           << ", erf.if_stability_wind_floor = " << if_stability_wind_floor
                           << " m/s (WRF-style bounds; the defaults 1 and 0 turn them off)\n";
        }
 
        if ((if_init_surf_temp > zero        && if_surf_temp_flux != amrex::Real(1e-8)) ||
            (if_init_surf_temp > zero        && if_Olen_in        != amrex::Real(1e-8)) ||
            (if_Olen_in != amrex::Real(1e-8) && if_surf_temp_flux != amrex::Real(1e-8)))
        {
            amrex::Error("Can only specify one of init_surf_temp, surf_temp_flux, or Olen");
        }
        if (if_Olen_in == zero &&
            (terrain_type == TerrainType::ImmersedForcing || buildings_type == BuildingsType::ImmersedForcing)) {
            amrex::Abort("erf.if_Olen = 0 is not an Obukhov length: give a nonzero length [m] (a large |L| is near neutral)");
        }
 
        // The immersed forcing's "almost fluid" threshold (erf-model/ERF#4226). It is recorded
        // under eb2.small_volfrac here, in ReadParameters, before ERF::initializeEB runs, so the
        // EB build (which covers cells below it), the blanking and the forcing kernels all use
        // the same value, and a deck without the key gets 0.005 rather than AMReX's 1e-14.
        if (terrain_type == TerrainType::ImmersedForcing || buildings_type == BuildingsType::ImmersedForcing) {
            amrex::ParmParse pp_eb2("eb2");
            pp_eb2.queryAdd("small_volfrac", if_small_volfrac);
            if (!(if_small_volfrac >= zero && if_small_volfrac < one)) {
                amrex::Abort("eb2.small_volfrac must be in [0, 1) with immersed forcing");
            }
        }
 
        // Implicit drag of the immersed forcing inside the anelastic projection (erf-model/ERF#4919
        // for terrain, #4887 for buildings): each face momentum becomes (m* - grad phi) / (1 + dt rate),
        // with rate the drag rate of the source term on that face, and the solid is not frozen
        // afterwards, so the projected field stays divergence free in the solid too
        pp.queryAdd("if_implicit_projection", if_implicit_projection);
        if (if_implicit_projection) {
            const bool if_terrain   = (terrain_type   == TerrainType::ImmersedForcing);
            const bool if_buildings = (buildings_type == BuildingsType::ImmersedForcing);
            if (!if_terrain && !if_buildings) {
                amrex::Abort("erf.if_implicit_projection needs erf.terrain_type = ImmersedForcing or "
                             "erf.buildings_type = ImmersedForcing");
            }
            if (if_terrain && if_buildings) {
                amrex::Abort("erf.if_implicit_projection is not implemented for erf.terrain_type = ImmersedForcing "
                             "together with erf.buildings_type = ImmersedForcing");
            }
            if (mesh_type != MeshType::ConstantDz) {
                amrex::Abort("erf.if_implicit_projection needs a constant-dz mesh");
            }
            for (int lev = 0; lev <= max_level; ++lev) {
                if (!anelastic[lev]) {
                    amrex::Abort("erf.if_implicit_projection needs erf.anelastic = 1 on every level");
                }
            }
#ifndef ERF_USE_FFT
            amrex::Abort("erf.if_implicit_projection needs a build with FFT (its GMRES solve is FFT-preconditioned)");
#endif
        }
 
        // Flag to do MOST rotations with terrain
        pp.queryAdd("use_rotate_surface_flux",use_rotate_surface_flux);
        if (use_rotate_surface_flux) {
            AMREX_ASSERT_WITH_MESSAGE(terrain_type != TerrainType::None,"MOST stress rotations are only valid with terrain!");
        }
 
        // Which external forcings?
        abl_driver_type = ABLDriverType::None; // Default: no ABL driver for simulating classical fluid dynamics problems
        pp.query_enum_case_insensitive("abl_driver_type",abl_driver_type);
        pp.queryAdd("const_massflux_u", const_massflux_u);
        pp.queryAdd("const_massflux_v", const_massflux_v);
        pp.queryAdd("const_massflux_tau", const_massflux_tau);
        // The control divides by tau; tau <= 0 gives no source or one that drives the flow
        //    away from the target
        if (!(const_massflux_tau > amrex::Real(0))) {
            amrex::Abort("erf.const_massflux_tau must be > 0 (got " +
                         std::to_string(const_massflux_tau) + ")");
        }
        pp.queryAdd("const_massflux_layer_lo", const_massflux_layer_lo);
        pp.queryAdd("const_massflux_layer_hi", const_massflux_layer_hi);
 
        // Which type of inflow turbulent generation
        // Real-data initialization uses parent controls unless explicitly overridden.
        perturbation_auto = (init_type == InitType::WRFInput || init_type == InitType::Metgrid);
        pp.queryAdd("perturbation_auto", perturbation_auto);
        pp.queryAdd("perturbation_auto_period", perturbation_auto_period);
        pp.queryAdd("perturbation_auto_pbl_fraction", perturbation_auto_pbl_fraction);
        if (perturbation_auto &&
            (!(perturbation_auto_period > zero) || !std::isfinite(perturbation_auto_period) ||
             !(perturbation_auto_pbl_fraction > zero && perturbation_auto_pbl_fraction <= one))) {
            amrex::Abort("Automatic perturbations require perturbation_auto_period > 0 and "
                         "0 < perturbation_auto_pbl_fraction <= 1");
        }
        pert_type.resize(max_level+1);
        for (int lev = 0; lev <= max_level; ++lev) {
            pert_type[lev] = PerturbationType::None; // Default
            query_one_or_per_level_enum_case_insensitive(
                pp, "perturbation_type", pert_type[lev], lev, max_level);
        }
 
        amrex::Vector<amrex::Real> abl_pressure_grad_in = {zero, zero, zero};
        pp.queryarr("abl_pressure_grad",abl_pressure_grad_in);
        for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_pressure_grad[i] = abl_pressure_grad_in[i];
 
        amrex::Vector<amrex::Real> abl_geo_forcing_in = {zero, zero, zero};
        if(pp.queryarr("abl_geo_forcing",abl_geo_forcing_in)) {
            amrex::Print() << "Specified abl_geo_forcing: (";
            for (int i = 0; i < AMREX_SPACEDIM; ++i) {
                abl_geo_forcing[i] = abl_geo_forcing_in[i];
                amrex::Print() << abl_geo_forcing[i] << " ";
            }
            amrex::Print() << ")" << std::endl;
        }
 
        if (use_coriolis)
        {
            build_coriolis_forcings_const_lat(pp_prefix);
        }
 
        pp.queryAdd("add_custom_rhotheta_forcing", custom_rhotheta_forcing);
        pp.queryAdd("add_custom_moisture_forcing", custom_moisture_forcing);
        pp.queryAdd("add_custom_w_subsidence", custom_w_subsidence);
        pp.queryAdd("add_do_theta_advection", do_theta_advection);  // If true, apply custom subsidence to (rho*theta) when add_custom_w_subsidence is used
        pp.queryAdd("add_do_mom_advection", do_mom_advection);      // If true, apply custom subsidence to momentum when add_custom_w_subsidence is used
        pp.queryAdd("add_custom_geostrophic_profile", custom_geostrophic_profile);
        pp.queryAdd("custom_forcing_uses_primitive_vars", custom_forcing_prim_vars);
        pp.queryAdd("spatial_rhotheta_forcing", spatial_rhotheta_forcing);
        pp.queryAdd("spatial_moisture_forcing", spatial_moisture_forcing);
 
        pp.queryAdd("nudging_from_input_sounding", nudging_from_input_sounding);
 
        // Nudging towards observations at stations; the stations and the
        // parameters are read by the ObsNudging class, under erf.obs_nudging
        pp.queryAdd("nudging_from_observations", nudging_from_observations);
        if (nudging_from_observations && terrain_type == TerrainType::EB) {
            amrex::Abort("erf.nudging_from_observations does not support erf.terrain_type = EB; "
                         "use a terrain-fitted mesh or immersed forcing");
        }
 
        int nudging_u_z1_set = pp.queryAdd("nudging_u_z1", nudging_u_z1);
        int nudging_u_z2_set = pp.queryAdd("nudging_u_z2", nudging_u_z2);
        pp.queryAdd("nudging_t_z1", nudging_t_z1);
        pp.queryAdd("nudging_t_z2", nudging_t_z2);
        pp.queryAdd("nudging_q_z1", nudging_q_z1);
        pp.queryAdd("nudging_q_z2", nudging_q_z2);
 
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nudging_u_z1 <= nudging_u_z2,
            "nudging_u_z1 must be less than or equal to nudging_u_z2");
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nudging_t_z1 <= nudging_t_z2,
            "nudging_t_z1 must be less than or equal to nudging_t_z2");
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nudging_q_z1 <= nudging_q_z2,
            "nudging_q_z1 must be less than or equal to nudging_q_z2");
 
        pp.queryAdd("nudging_u", nudging_u);
        pp.queryAdd("nudging_t", nudging_t);
        pp.queryAdd("nudging_q", nudging_q);
 
        pp.queryAdd("large_scale_forcing", large_scale_forcing);
 
        if (nudging_from_input_sounding && large_scale_forcing) {
            if (!nudging_u) {
                amrex::Print() << "WARNING: erf.nudging_u = false has no effect when erf.large_scale_forcing "
                                  "= true; momentum nudging is then controlled solely by large_scale_forcing.\n";
            }
            if (nudging_u_z1_set || nudging_u_z2_set) {
                amrex::Print() << "WARNING: erf.nudging_u_z1/erf.nudging_u_z2 are ignored when "
                                  "erf.large_scale_forcing = true; momentum nudging is then applied "
                                  "at all heights.\n";
            }
        }
 
        have_geo_wind_profile = (!abl_geo_wind_table.empty() || custom_geostrophic_profile);
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(!abl_geo_wind_table.empty() && custom_geostrophic_profile),
            "Should not have both abl_geo_wind_table and custom_geostrophic_profile set.");
 
        pp.queryAdd("Ave_Plane", ave_plane);
 
        // Use numerical diffusion?
        pp.queryAdd("num_diff_coeff",num_diff_coeff);
        AMREX_ASSERT_WITH_MESSAGE(( (num_diff_coeff >= zero) && (num_diff_coeff <= one) ),
                                  "Numerical diffusion coefficient must be between 0 & one");
        use_num_diff = (num_diff_coeff > 0);
        if (use_num_diff) {
            amrex::Print() << "6th-order numerical diffusion turned on with coefficient = "
                << num_diff_coeff << std::endl;
            num_diff_coeff *= amrex::Math::powi<-6>(two);
        }
 
            advChoice.init_params(pp_prefix);
           diffChoice.init_params(pp_prefix);
        dampingChoice.init_params(pp_prefix);
         spongeChoice.init_params(pp_prefix);
             ebChoice.init_params(pp_prefix);
 
        turbChoice.resize(max_level+1);
        for (int lev = 0; lev <= max_level; lev++) {
            turbChoice[lev].init_params(lev,max_level,pp_prefix);
        }
 
        if (perturbation_auto) {
            // Automatic forcing has a parent and a three-dimensional LES closure.
            // A horizontal Smagorinsky closure paired with a PBL scheme is mesoscale.
            for (int lev = 0; lev <= max_level; ++lev) {
                const auto& tc = turbChoice[lev];
                if (lev == 0 || tc.les_type == LESType::None || tc.smag2d ||
                    tc.pbl_type != PBLType::None) {
                    pert_type[lev] = PerturbationType::None;
                }
            }
        }
 
        // WRF's Smagorinsky-2D slope limiter divides K_h by the terrain slope factor of the
        // mesh, which is 1 everywhere without a terrain-fitted mesh; and the EB closure
        // (ComputeTurbulentViscosityLES_EB) applies neither limit.
        bool smag2d_requested = false;
        bool smag2d_somewhere = false;
        for (int lev = 0; lev <= max_level; lev++) {
            smag2d_requested = smag2d_requested || turbChoice[lev].smag2d_limits_requested;
            smag2d_somewhere = smag2d_somewhere || turbChoice[lev].smag2d;
        }
        if (smag2d_requested && !smag2d_somewhere) {
            amrex::Error("erf.smag2d_slope_limiter and erf.smag2d_kh_cap apply only to "
                         "erf.les_type = Smagorinsky2D, which no level uses");
        }
        for (int lev = 0; lev <= max_level; lev++) {
            const bool fitted = (terrain_type == TerrainType::StaticFittedMesh ||
                                 terrain_type == TerrainType::MovingFittedMesh);
            if (turbChoice[lev].smag2d_slope_limiter && !fitted) {
                amrex::Error("erf.smag2d_slope_limiter requires erf.terrain_type = StaticFittedMesh "
                             "or MovingFittedMesh");
            }
            if (turbChoice[lev].smag2d_kh_cap > zero && terrain_type == TerrainType::EB) {
                amrex::Error("erf.smag2d_kh_cap is not supported with erf.terrain_type = EB");
            }
        }
 
        // erf.radiation_model = TwoStream selects the two-stream model; its
        // erf.radiation.* sub-options are read only then.
        radChoice.init_params(rad_type == RadiationType::TwoStream, pp_prefix);
        if (rad_type == RadiationType::RRTMGP) {
            rrtmgpSolarChoice.read(pp_prefix);
        }
        topographicChoice.read(pp_prefix, rrtmgpSolarChoice);
        if (topographicChoice.active()) {
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(terrain_type == TerrainType::StaticFittedMesh,
                "Topographic radiation requires erf.terrain_type = StaticFittedMesh");
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rad_type == RadiationType::RRTMGP,
                "Topographic radiation currently requires erf.radiation_model = RRTMGP");
        }
 
        // Noah-MP integrates on the downwelling shortwave, downwelling longwave and solar
        // zenith angle a radiation model writes into it. Only RRTMGP and TwoStream register
        // those outputs with the land model (rad_feeds_lsm, ERF::init_stuff), so under any
        // other choice they are never written and Noah-MP runs with no radiative forcing at
        // all -- the inputs it reads fall back to zero. That is a legitimate thing to want
        // in a test, and rarely what anyone running a land model means, so say it at
        // start-up.
        if (lsm_type == LandSurfaceType::NOAHMP && !rad_feeds_lsm()) {
            amrex::Print() << "WARNING: erf.land_surface_model = NOAHMP with erf.radiation_model = "
                           << amrex::getEnumNameString(rad_type)
                           << ": only RRTMGP and TwoStream supply the downwelling shortwave, "
                              "downwelling longwave and solar zenith angle Noah-MP needs. Noah-MP "
                              "will run with no radiative forcing (zero SWDOWN, GLW and COSZEN).\n";
        }
 
        vert_implicit_fac.resize(max_level+1);
        for (int lev = 0; lev <= max_level; lev++) {
            vert_implicit_fac[lev].resize(3);
            vert_implicit_fac[lev][0] = one;
            vert_implicit_fac[lev][1] = one;
            vert_implicit_fac[lev][2] = zero;
        }
 
        // YSU PBL: use consistent coriolis frequency
        for (int lev = 0; lev <= max_level; lev++) {
            if (turbChoice[lev].pbl_ysu_use_consistent_coriolis) {
                if (use_coriolis) {
                    turbChoice[lev].pbl_ysu_coriolis_freq = coriolis_factor * sinphi;
                    if (lev == 0) {
                        amrex::Print() << "YSU PBL using ERF coriolis frequency: " << turbChoice[lev].pbl_ysu_coriolis_freq << std::endl;
                    }
                } else {
                    amrex::Error("YSU cannot use ERF coriolis frequency if not using coriolis");
                }
            }
        }
        // MRF
        for (int lev = 0; lev <= max_level; lev++) {
            if (turbChoice[lev].pbl_ysu_use_consistent_coriolis) {
                if (use_coriolis) {
                    turbChoice[lev].pbl_ysu_coriolis_freq = coriolis_factor * sinphi;
                    if (lev == 0) {
                        amrex::Print() << "MRF PBL using ERF coriolis frequency: " << turbChoice[lev].pbl_ysu_coriolis_freq << std::endl;
                    }
                } else {
                    amrex::Error("MRF cannot use ERF coriolis frequency if not using coriolis");
                }
            }
        }
 
        // What's the strength of the bdy nudging?
        pp.queryAdd("bdy_nudge_factor",bdy_nudge_factor);
 
        // Optionally use the dry-air density reconstructed from wrfbdy.
        pp.queryAdd("use_wrf_bdy_density", use_wrf_bdy_density);
 
        // Optionally ingest cloud water and cloud ice from wrfinput/wrfbdy.
        pp.queryAdd("use_wrf_bdy_qc_qi", use_wrf_bdy_qc_qi);
 
        if (!use_real_bcs || init_type != InitType::WRFInput) {
            use_wrf_bdy_density = false;
        }
 
        pp.queryAdd("bdy_rho_nudge_factor", bdy_rho_nudge_factor);
 
        // Which approach to nudging the moist variables in the boundary region?
        pp.queryAdd("bdy_moist_nudge_type",bdy_moist_nudge_type);
        if (bdy_moist_nudge_type < 0 || bdy_moist_nudge_type > 3) {
            amrex::Error("erf.bdy_moist_nudge_type must be one of 0, 1, 2, or 3");
        }
        if (bdy_moist_nudge_type == 3 && !use_wrf_bdy_qc_qi) {
            amrex::Error("erf.bdy_moist_nudge_type = 3 requires erf.use_wrf_bdy_qc_qi = true");
        }
        if (use_wrf_bdy_qc_qi) {
            if (moisture_type == MoistureType::None) {
                amrex::Error("erf.use_wrf_bdy_qc_qi requires an active moisture model");
            }
            if (!use_real_bcs || init_type != InitType::WRFInput) {
                amrex::Error("erf.use_wrf_bdy_qc_qi requires WRFInput real boundary conditions");
            }
            if (moisture_indices.qv < RhoQ1_comp || moisture_indices.qv > RhoQ11_comp ||
                moisture_indices.qc < RhoQ1_comp || moisture_indices.qc > RhoQ11_comp) {
                amrex::Error("erf.use_wrf_bdy_qc_qi requires active qv and qc component indices");
            }
            const int optional_moisture_comps[] = {moisture_indices.qi,
                                                   moisture_indices.qr,
                                                   moisture_indices.qs,
                                                   moisture_indices.qg};
            for (const int comp : optional_moisture_comps) {
                if (comp != -1 && (comp < RhoQ1_comp || comp > RhoQ11_comp)) {
                    amrex::Error("erf.use_wrf_bdy_qc_qi found an invalid active moisture component index");
                }
            }
        }
 
        // Are me transporting the scalar component?
        pp.queryAdd("transport_scalar",transport_scalar);
 
        for (int lev = 0; lev <= max_level; ++lev) {
            if (turbChoice[lev].uses_eamxx_shoc()) {
                use_eamxx_shoc = true;
            }
            if (turbChoice[lev].uses_native_shoc()) {
                use_native_shoc = true;
                read_shoc_transport_modes(shoc_transport_mode, shoc_momentum_transport);
            }
        }
 
        // Set up the vertical implicit diffusion solve
        int n_impfac = pp.countval("vert_implicit_fac");
        AMREX_ALWAYS_ASSERT(n_impfac == 0 || n_impfac == 1 || n_impfac==3);
        if (n_impfac == 1) {
            amrex::Real fac_in;
            pp.get("vert_implicit_fac", fac_in);
            for (int lev = 0; lev <= max_level; lev++) {
                for (int i=0; i<3; ++i) {
                    vert_implicit_fac[lev][i] = fac_in;
                }
            }
        } else if (n_impfac == 3) {
            for (int lev = 0; lev <= max_level; lev++) {
                pp.getarr("vert_implicit_fac", vert_implicit_fac[lev]);
            }
        }
 
        // Record whether the user explicitly requested these
        const bool user_set_imp_thermal  = pp.contains("implicit_thermal_diffusion");
        const bool user_set_imp_moisture = pp.contains("implicit_moisture_diffusion");
        const bool user_set_imp_ke       = pp.contains("implicit_ke_diffusion");
        const bool user_set_imp_momentum = pp.contains("implicit_momentum_diffusion");
 
        // If true (default), include implicit contributions to vertical thermal diffusion
        pp.queryAdd("implicit_thermal_diffusion", implicit_thermal_diffusion);
 
        // If true (default), include implicit contributions to vertical moisture diffusion
        pp.queryAdd("implicit_moisture_diffusion", implicit_moisture_diffusion);
 
        // If true (default), include implicit contributions to vertical KE diffusion
        pp.queryAdd("implicit_ke_diffusion", implicit_ke_diffusion);
 
        // If true (default), include implicit contributions in tau13, tau23,
        pp.queryAdd("implicit_momentum_diffusion", implicit_momentum_diffusion);
 
        // If true (default false), also solve the compact terrain-metric part of the vertical
        // diffusion, K_h M d2/dz2 (ERF_TerrainImplicitMetric.H), implicitly (terrain-fitted
        // meshes only)
        pp.queryAdd("implicit_terrain_metric", implicit_terrain_metric);
 
        // The default is true (i.e. that it happens before the acoustic substepping).
        pp.queryAdd("implicit_before_substep", implicit_before_substep);
        for (int lev = 0; lev <= max_level; lev++) {
            if ( (substepping_type[lev] == SubsteppingType::None) && !implicit_before_substep) {
                amrex::Print() << "implicit_before_substep cannot be false without substepping; setting to true." << "\n";
                implicit_before_substep = true;
            }
        }
 
        // Handle shoc case -- SHOC owns whatever vertical diffusion the host does not
        if (uses_shoc_family()) {
            const bool host_owns_scalars  = host_owns_vertical_scalar_diffusion();
            const bool host_owns_momentum = host_owns_vertical_momentum_diffusion();
            if (!host_owns_scalars &&
                ( (user_set_imp_thermal  && implicit_thermal_diffusion ) ||
                  (user_set_imp_moisture && implicit_moisture_diffusion) ||
                  (user_set_imp_ke       && implicit_ke_diffusion      ) ) ) {
                amrex::Print() << "Ignoring erf.implicit_{thermal,moisture,ke}_diffusion = true; "
                               << "SHOC owns vertical scalar diffusion." << "\n";
            }
            if (!host_owns_momentum && user_set_imp_momentum && implicit_momentum_diffusion) {
                amrex::Print() << "Ignoring erf.implicit_momentum_diffusion = true; "
                               << "SHOC owns vertical momentum diffusion." << "\n";
            }
 
            implicit_thermal_diffusion  = (implicit_thermal_diffusion  && host_owns_scalars);
            implicit_moisture_diffusion = (implicit_moisture_diffusion && host_owns_scalars);
            implicit_ke_diffusion       = (implicit_ke_diffusion       && host_owns_scalars);
            implicit_momentum_diffusion = (implicit_momentum_diffusion && host_owns_momentum);
 
            // Make the resulting configuration visible in the run log
            if (!host_owns_scalars && !host_owns_momentum) {
                amrex::Print() << "Turning off vertical implicit solve for all components since SHOC owns diffusion." << "\n";
            } else {
                amrex::Print() << "SHOC hands vertical diffusion of"
                               << (host_owns_scalars  ? " scalars"  : "")
                               << (host_owns_momentum ? " momentum" : "")
                               << " to the host; keeping the vertical implicit solve on (theta = "
                               << implicit_thermal_diffusion
                               << ", moisture = " << implicit_moisture_diffusion
                               << ", tke = "      << implicit_ke_diffusion
                               << ", momenta = "  << implicit_momentum_diffusion << ")." << "\n";
            }
        } // uses_shoc_family
 
        // Do not allow implicit vertical diff if user turns it off
        bool do_vert_implicit = true;
        pp.queryAdd("vert_implicit", do_vert_implicit);
        if (!do_vert_implicit) {
            for (int lev = 0; lev <= max_level; lev++) {
                amrex::Print() << "Turning off native vertical implicit solve from vert_implicit flag at level " << lev << std::endl;
                vert_implicit_fac[lev][0] = zero;
                vert_implicit_fac[lev][1] = zero;
                vert_implicit_fac[lev][2] = zero;
            }
        }
 
        // The implicit vertical diffusion is only second order in an anelastic run when
        // the two stages are the midpoint method and the solve is done in the first stage
        // alone: S^* = S^n + (dt/2) F(S^*) and S^{n+1} = S^n + dt F(S^*) give
        // S^* = (S^n + S^{n+1})/2, i.e. the implicit midpoint rule.  Solving again in the
        // second stage would reduce the pair to backward Euler, and with the RK2 (Heun)
        // stages there is no stage that plays the role of the midpoint at all.
        for (int lev = 0; lev <= max_level; lev++) {
            if (anelastic[lev]) {
                if (anelastic_type[lev] == AnelasticType::RK2) {
                    amrex::Print() << "Turning off native vertical implicit solve flag because anelastic at level " << lev
                                   << " with RK2 scheme." << std::endl;
                    vert_implicit_fac[lev][0] = zero;
                    vert_implicit_fac[lev][1] = zero;
                    vert_implicit_fac[lev][2] = zero;
                } else if (anelastic_type[lev] == AnelasticType::MidPoint) {
                    if (do_vert_implicit && vert_implicit_fac[lev][0] != one) {
                        amrex::Print() << "Forcing implicit factor 1 in first stage of MidPoint to ensure second-order accuracy." << std::endl;
                        vert_implicit_fac[lev][0] = one;
                    }
                    if (vert_implicit_fac[lev][1] > zero) {
                        amrex::Print() << "Turning off the second-stage vertical implicit solve at level " << lev
                                       << " because anelastic with MidPoint solves in the first stage only."
                                       << std::endl;
                    }
                    // The anelastic integrator only takes two stages, so [2] is never used
                    vert_implicit_fac[lev][1] = zero;
                    vert_implicit_fac[lev][2] = zero;
                }
            }
        }
 
        // Do not allow implicit vertical diff with EB
        if (terrain_type == TerrainType::EB) {
            for (int lev = 0; lev <= max_level; lev++) {
                if (do_vert_implicit) {
                    amrex::Print() << "Implicit diffusion is not supported with EB; turning off at level " << lev << "\n";
                    vert_implicit_fac[lev][0] = zero;
                    vert_implicit_fac[lev][1] = zero;
                    vert_implicit_fac[lev][2] = zero;
                }
            }
        }
 
        // The EB diffusion path adds no buoyancy, shear or dissipation source for
        // a prognostic TKE closure
        if (terrain_type == TerrainType::EB) {
            for (int lev = 0; lev <= max_level; lev++) {
                if (turbChoice[lev].use_keqn) {
                    amrex::Error("erf.terrain_type = EB cannot be used with erf.les_type = Deardorff "
                                 "or erf.rans_type = kEqn: the EB diffusion path has no TKE sources");
                }
            }
        }
 
        // Do not allow implicit vertical diff when all the components are false
        if (!implicit_thermal_diffusion  &&
            !implicit_moisture_diffusion &&
            !implicit_ke_diffusion       &&
            !implicit_momentum_diffusion) {
            amrex::Print() << "Thermal, moisture, KE, and momentum diffusion are all turned off; turning off vertical implicit solve." << std::endl;
            for (int lev = 0; lev <= max_level; lev++) {
                vert_implicit_fac[lev][0] = zero;
                vert_implicit_fac[lev][1] = zero;
                vert_implicit_fac[lev][2] = zero;
            }
        }
 
        // Which type of multilevel coupling
        coupling_type = CouplingType::TwoWay; // Default
        pp.query_enum_case_insensitive("coupling_type",coupling_type);
 
        // Test for hybrid (compressible + anelastic) -- in this case we must use one-way coupling
        if (any_anelastic && any_compress) {
            coupling_type = CouplingType::OneWay;
        }
 
        // Which type of windfarm model
        windfarm_type = WindFarmType::None; // Default
        pp.query_enum_case_insensitive("windfarm_type",windfarm_type);
 
        windfarm_loc_type = WindFarmLocType::None;
        pp.query_enum_case_insensitive("windfarm_loc_type",windfarm_loc_type);
 
        pp.queryAdd("windfarm_loc_table",  windfarm_loc_table);
        pp.queryAdd("windfarm_spec_table", windfarm_spec_table);
        pp.queryAdd("windfarm_blade_table", windfarm_blade_table);
        pp.queryAdd("windfarm_airfoil_tables", windfarm_airfoil_tables);
        pp.queryAdd("windfarm_spec_table_extra", windfarm_spec_table_extra);
 
        // Sampling distance upstream of the turbine to find the
        // incoming free stream velocity as a factor of the diameter of the
        // turbine. ie. the sampling distance will be this number multiplied
        // by the diameter of the turbine
        pp.queryAdd("sampling_distance_by_D", sampling_distance_by_D);
        pp.queryAdd("turb_disk_angle_from_x", turb_disk_angle);
 
        pp.queryAdd("windfarm_x_shift",windfarm_x_shift);
        pp.queryAdd("windfarm_y_shift",windfarm_y_shift);
        // Test if time averaged data is to be output
        pp.queryAdd("time_avg_vel",time_avg_vel);
 
        // Interval means and resolved turbulent moments for plotfile output
        pp.queryAdd("compute_mean_vars", compute_mean_vars);
        pp.queryAdd("mean_vars_reset_mode", mean_vars_reset_mode);
        pp.queryAdd("mean_vars_reset_time", mean_vars_reset_time);
 
        if (compute_mean_vars &&
            mean_vars_reset_mode != "plotfile" && mean_vars_reset_mode != "time") {
            amrex::Error("erf.mean_vars_reset_mode must be either 'plotfile' or 'time'");
        }
        if (compute_mean_vars && mean_vars_reset_mode == "time" &&
            mean_vars_reset_time < zero) {
            amrex::Error("erf.mean_vars_reset_time must be non-negative when "
                         "erf.mean_vars_reset_mode = 'time'");
        }
 
        pp.queryAdd("hindcast_lateral_forcing", hindcast_lateral_forcing);
 
        if (hindcast_lateral_forcing) {
            pp.queryAdd("hindcast_boundary_data_dir", hindcast_boundary_data_dir);
 
            if(hindcast_boundary_data_dir.empty()) {
                amrex::Error("ERROR: Missing input parameter 'erf.hindcast_boundary_data_dir' for boundary data for lateral forcing");
            }
            pp.queryAdd("hindcast_data_interval_in_hrs", hindcast_data_interval_in_hrs);
            if(hindcast_data_interval_in_hrs < zero) {
                amrex::Error("ERROR: Input parameter 'erf.hindcast_data_interval_in_hrs' which is the time interval between the "
                             "data files is either missing or set to less than zero");
            }
            pp.queryAdd("hindcast_lateral_sponge_strength", hindcast_lateral_sponge_strength);
            pp.queryAdd("hindcast_lateral_sponge_length", hindcast_lateral_sponge_length);
 
            pp.queryAdd("hindcast_zhi_sponge_length", hindcast_zhi_sponge_length);
            pp.queryAdd("hindcast_zhi_sponge_strength", hindcast_zhi_sponge_strength);
 
            pp.queryAdd("hindcast_zhi_sponge_damping", hindcast_zhi_sponge_damping);
 
            if (hindcast_lateral_forcing and hindcast_lateral_sponge_strength < zero) {
                amrex::Error("ERROR: Missing input parameter 'erf.hindcast_lateral_sponge_strength' or it is specified to be less than zero");
            }
 
            if (hindcast_lateral_forcing and hindcast_lateral_sponge_length < zero) {
                amrex::Error("ERROR: Missing input parameter 'erf.hindcast_lateral_sponge_length' or it is specified to be less than zero");
            }
 
            if (hindcast_zhi_sponge_damping and hindcast_zhi_sponge_length < zero) {
                amrex::Error("ERROR: Missing input parameter 'erf.hindcast_zhi_sponge_length' or it is specified to be less than zero");
            }
 
            if (hindcast_zhi_sponge_damping and hindcast_zhi_sponge_strength < zero) {
                amrex::Error("ERROR: Missing input parameter 'erf.hindcast_zhi_sponge_strength' or it is specified to be less than zero");
            }
        }
 
        pp.queryAdd("hindcast_surface_bcs", hindcast_surface_bcs);
        if(hindcast_surface_bcs) {
            pp.queryAdd("hindcast_surface_data_dir", hindcast_surface_data_dir);
        }
 
        pp.queryAdd("io_hurricane_eye_tracker", io_hurricane_eye_tracker);
        if (io_hurricane_eye_tracker) {
            pp.queryAdd("hurricane_eye_latitude", hurricane_eye_latitude);
            pp.queryAdd("hurricane_eye_longitude", hurricane_eye_longitude);
            if(hurricane_eye_latitude == amrex::Real(-1e10) or hurricane_eye_longitude == amrex::Real(-1e10)) {
                amrex::Error("ERROR: You are using 'erf.io_hurricane_eye_tracker' to write out the files that track the eye of the hurricane"
                              " but have not provided the initial location of the eye of the hurricane to be tracked. There has to be two"
                              " options in the inputs - erf.hurricane_eye_latitude and erf.hurricane_eye_longitude that gives an approximate"
                              " location of the eye in the initial condition");
            }
        }
 
        pp.queryAdd("is_init_for_ensemble", is_init_for_ensemble);
        if(is_init_for_ensemble) {
            amrex::ParmParse pp_ens("ensemble");
            pp_ens.queryAdd("n_members", n_ensemble);
            if(n_ensemble < 2) {
                amrex::Abort("You are using an ensemble run. There needs to be at least 2 ensemble members. "
                             "erf.n_ensemble must be >=2.");
            }
            pp_ens.queryAdd("coarse_bckgnd_data_file", coarse_bckgnd_data_file);
            if (coarse_bckgnd_data_file.empty()) {
                amrex::Abort("coarse_bckgnd_data_file is empty! For ensmeble simulations, there needs to be a coarse background file which "
                             "contains the data which will be interpolated onto the fine mesh. There has to a entry ensemble.coarse_bckgnd_data_file "
                             "which contains the filename in the inputs.");
            }
            pp_ens.queryAdd("ens_pert_amplitude", ens_pert_amplitude);
            if(ens_pert_amplitude <= 0.0) {
                amrex::Error("You are using initialization for ensemble simulations using the inputs option "
                             "ensemble.is_init_for_ensemble=true. In this case, there has to be an option "
                             "ensemble.ens_pert_amplitude which is the value of the amplitude of the perturbation  "
                             "to be added to the background state and has to be greater than 0.0");
            }
            pp_ens.queryAdd("ens_pert_correlated_radius", ens_pert_correlated_radius);
            if(ens_pert_correlated_radius <= 0.0) {
                amrex::Error("You are using initialization with spatially correlated perturbations using the inputs option "
                             "ensemble.is_init_for_ensemble=true. In this case, there has to be an option "
                             "ensemble.ens_pert_correlated_radius which is the value of the the spatial correlation radius, "
                             "and has to be greater than 0.0");
            }
        }
    }
 
    /**
     * @brief Validate solver options after geometry and boundary conditions are known.
     * @param max_level Maximum AMR level configured for the run.
     * @param geom_vect Geometry for each AMR level.
     * @param phys_bc_type Physical boundary-condition types.
     */
    void check_params (int max_level, const amrex::Vector<amrex::Geometry>& geom_vect, amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2> phys_bc_type)
    {
        // The weighted average down acts only where covered coarse cells are averaged down
        if (if_average_down == "fluid_mass") {
            if (max_level > 0 && coupling_type == CouplingType::TwoWay) {
                amrex::Print() << "Immersed forcing: erf.if_average_down = fluid_mass (the cell average down "
                               << "over refined levels is weighted by the fluid mass of the fine cells)\n";
            } else {
                amrex::Warning("erf.if_average_down = fluid_mass has no effect without refined levels "
                               "(amr.max_level > 0) and erf.coupling_type = TwoWay");
            }
        }
        if (if_fraction_stress) {
            // The face-above spacing factor is bounded only while z0 < dz/4 (then d2 > d1 for
            // every solid fraction and the fluid centroid of a full wall cell sits above 2 z0)
            for (int lev = 0; lev <= max_level; ++lev) {
                const amrex::Real dz = geom_vect[lev].CellSize(2);
                if (!(if_z0 < amrex::Real(0.25) * dz)) {
                    amrex::Abort("erf.if_wall_form = fraction_stress needs erf.if_z0 < dz/4 on every level (level " +
                                 std::to_string(lev) + ": z0 = " + std::to_string(if_z0) + " m, dz = " +
                                 std::to_string(dz) + " m)");
                }
            }
        }
        if (if_fraction_stress && (if_wall_tke || if_flat_wall_dissipation)) {
            bool any_keqn = false;
            for (int lev = 0; lev <= max_level; ++lev) {
                any_keqn = any_keqn || (turbChoice[lev].rans_type == RANSType::kEqn);
            }
            if (!any_keqn) {
                amrex::Abort("erf.if_wall_tke and erf.if_flat_wall_dissipation act on the k-eqn RANS closure: "
                             "set erf.rans_type = kEqn on a level, or turn them off (erf.if_flat_wall_dissipation = false)");
            }
        }
        if (if_fraction_stress) {
            amrex::Print() << "Immersed forcing: erf.if_wall_form = fraction_stress (wall stress in the wall cell, "
                           << "the wall face closed to diffusion; face spacing "
                           << ((if_wall_face_spacing == 2) ? "log" : (if_wall_face_spacing == 1) ? "physical" : "grid")
                           << ", wall TKE " << (if_wall_tke ? "on" : "off")
                           << ", flat-wall dissipation " << (if_flat_wall_dissipation ? "on" : "off") << ")\n";
        }
        // The log-law coarse/fine fill acts only on refined levels
        if (cf_loglaw_fill) {
            const bool immersed = (terrain_type == TerrainType::ImmersedForcing ||
                                   buildings_type == BuildingsType::ImmersedForcing);
            if (!immersed && phys_bc_type[amrex::Orientation(amrex::Direction::z, amrex::Orientation::low)] != ERF_BC::surface_layer) {
                amrex::Abort("erf.cf_loglaw_fill needs an immersed surface (erf.terrain_type or "
                             "erf.buildings_type = ImmersedForcing) or zlo.type = surface_layer");
            }
            if (max_level > 0) {
                amrex::Print() << "Coarse/fine faces: erf.cf_loglaw_fill = true (the wall coarse face's flux is shared "
                               << "among its fine sub-faces by the integral of ln(1 + z/z0), z0 = " << cf_loglaw_z0 << " m)\n";
            } else {
                amrex::Warning("erf.cf_loglaw_fill has no effect without refined levels (amr.max_level > 0)");
            }
        }
        if (perturbation_auto && any_perturbation() &&
            (terrain_type == TerrainType::EB || terrain_type == TerrainType::ImmersedForcing)) {
            amrex::Abort("Automatic perturbations currently support flat, stretched, and fitted-terrain meshes; "
                         "immersed/embedded terrain requires a ground-relative parent PBL-height conversion");
        }
#if 0
        // Warn for PBL models and moisture - these may not yet be compatible
        for (int lev = 0; lev <= max_level; lev++) {
            if ((moisture_type != MoistureType::None) && (turbChoice[lev].pbl_type != PBLType::None)) {
                amrex::Warning("\n*** WARNING: Moisture may not yet be compatible with PBL models, \n    proceed with caution ***");
            }
        }
#endif
        //
        // Buoyancy type check
        //
        for (int lev = 0; lev <= max_level; lev++) {
            if (buoyancy_type[lev] != 1 && buoyancy_type[lev] != 2 && buoyancy_type[lev] != 3 && buoyancy_type[lev] != 4) {
                amrex::Error("buoyancy_type must be 1, 2, 3 or 4");
            }
        }
 
        if (!use_lagged_delta_rt && !(terrain_type == TerrainType::MovingFittedMesh)) {
            amrex::Error("Can't turn off lagged_delta_rt when terrain not moving");
        }
 
        //
        // Wind farm checks
        //
        if (windfarm_type==WindFarmType::SimpleAD and sampling_distance_by_D < zero) {
             amrex::Error("To use simplified actuator disks, you need to provide a variable"
                           " erf.sampling_distance_by_D in the inputs which specifies the upstream"
                           " distance as a factor of the turbine diameter at which the incoming free stream"
                           " velocity will be computed at.");
        }
        if ( (windfarm_type==WindFarmType::SimpleAD ||
              windfarm_type==WindFarmType::GeneralAD ) && turb_disk_angle < zero) {
            amrex::Error("To use simplified actuator disks, you need to provide a variable"
                          " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the"
                          " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction"
                          " will have turb_disk_angle value of 90 deg.");
        }
        if (windfarm_loc_type == WindFarmLocType::lat_lon and (windfarm_x_shift < zero or windfarm_y_shift < zero)) {
            amrex::Error("You are using windfarms with latitude-logitude option to position the turbines."
                         " For this you should provide the inputs erf.windfarm_x_shift and"
                         " erf.windfarm_y_shift which are the values by which the bounding box of the"
                         " windfarm is shifted from the x and the y axes.");
        }
 
 
        if ( (const_massflux_u != 0) && !(geom_vect[0].isPeriodic(0)) ) {
            amrex::Error("Constant mass flux (in x) should be used with periodic boundaries");
        }
 
        if ( (const_massflux_v != 0) && !(geom_vect[0].isPeriodic(1)) ) {
            amrex::Error("Constant mass flux (in y) should be used with periodic boundaries");
        }
 
        // Mesoscale diffusion -- test if LES appropriate with this dx,dy
        for (int lev = 0; lev <= max_level; lev++) {
            if ((geom_vect[lev].CellSize(0) > amrex::Real(2000.)) || (geom_vect[lev].CellSize(1) > amrex::Real(2000.)))
            {
                if ( (turbChoice[lev].les_type == LESType::Smagorinsky) && !turbChoice[lev].smag2d ) {
                    amrex::Warning("Should use 2-D Smagorinsky for mesoscale resolution");
                } else if (turbChoice[lev].les_type == LESType::Deardorff) {
                    amrex::Warning("Should not use Deardorff LES for mesoscale resolution");
                }
            }
        }
 
        // Turn off implicit solve if we have no diffusion
        bool l_use_kturb = turbChoice[0].use_kturb;
        for (int lev = 1; lev <= max_level; lev++) {
            l_use_kturb = (l_use_kturb || turbChoice[lev].use_kturb);
        }
        bool l_use_diff    = ( (diffChoice.molec_diff_type != MolecDiffType::None) || l_use_kturb );
 
        for (int lev = 0; lev <= max_level; lev++) {
            bool l_implicit_diff = (vert_implicit_fac[lev][0] > zero ||
                                    vert_implicit_fac[lev][1] > zero ||
                                    vert_implicit_fac[lev][2] > zero);
            if (l_implicit_diff && !l_use_diff) {
                amrex:: Print() << "No molecular or turbulent diffusion, turning off implicit solve at level " << lev << std::endl;
                vert_implicit_fac[lev][0] = zero;
                vert_implicit_fac[lev][1] = zero;
                vert_implicit_fac[lev][2] = zero;
            }
        }
 
        // Two-way coupling refluxes the diffusive fluxes of the slow scalars, but not on an EB mesh
        if (coupling_type == CouplingType::TwoWay && terrain_type == TerrainType::EB &&
            max_level > 0 && l_use_diff) {
            amrex::Warning("erf.coupling_type = TwoWay with erf.terrain_type = EB does not reflux the "
                           "diffusive fluxes: diffusion does not conserve the scalars across "
                           "coarse/fine boundaries");
        }
 
        // The implicit terrain metric is a change to the vertical implicit solves on a
        // terrain-fitted mesh; refuse it where it could not do anything.  This comes after
        // every block that can turn the implicit solve off.
        if (implicit_terrain_metric) {
            if (mesh_type != MeshType::VariableDz) {
                amrex::Abort("erf.implicit_terrain_metric = true requires a terrain-fitted mesh with "
                             "non-flat terrain (erf.terrain_type = StaticFittedMesh or MovingFittedMesh); "
                             "on a flat or a purely stretched mesh there are no terrain-metric terms");
            }
            // In a stage without the implicit solve the whole averaged metric term is explicit
            // for that stage.  On the steep ridge test the option gained little with 1 1 0
            // and ran worse than without it with anelastic MidPoint (measured before the exact
            // split: 15 s instead of 30 s), while 1 1 1 recovered the full time step.  The rule
            // is measured, not derived: every stage must solve implicitly.
            for (int lev = 0; lev <= max_level; lev++) {
                for (int i = 0; i < 3; ++i) {
                    if (vert_implicit_fac[lev][i] <= zero) {
                        amrex::Abort("erf.implicit_terrain_metric = true needs the vertical implicit solve "
                                     "in every Runge-Kutta stage at every level, but erf.vert_implicit_fac "
                                     "is not positive in stage " + std::to_string(i+1) + " at level " +
                                     std::to_string(lev) + ".  Set erf.vert_implicit_fac = 1 1 1; the "
                                     "option cannot be used where the implicit solve is turned off "
                                     "(erf.vert_implicit = false, every erf.implicit_*_diffusion = false, "
                                     "no molecular or turbulent diffusion, or anelastic, which solves in "
                                     "its first stage only or not at all)");
                    }
                }
            }
            // The momentum solves carry most of the stiffness; with them explicit the option
            // treats only the scalar metric terms implicitly
            if (!implicit_momentum_diffusion) {
                amrex::Print() << "WARNING: erf.implicit_terrain_metric = true with "
                               << "erf.implicit_momentum_diffusion = false: the momentum terrain-metric "
                               << "terms stay explicit; only the scalars are treated implicitly" << std::endl;
            }
        } else if (mesh_type == MeshType::VariableDz) {
            // On a terrain-fitted mesh the projected horizontal stresses carry K_h, and their
            // terrain-metric part K_h h^2 d2/dz2 is explicit in any stage without the implicit
            // solve.  On steep slopes that limits the time step (steep test ridge: 45 s before
            // #4214 was fixed, 15 s after, at the default 1 1 0).  Say so once at start-up.
            bool explicit_stage_with_kh = false;
            for (int lev = 0; lev <= max_level; lev++) {
                if (turbChoice[lev].les_type == LESType::None) { continue; }
                for (int i = 0; i < 3; ++i) {
                    if (vert_implicit_fac[lev][i] <= zero) { explicit_stage_with_kh = true; }
                }
            }
            if (explicit_stage_with_kh) {
                amrex::Print() << "NOTE: terrain-fitted mesh with an LES closure (horizontal eddy "
                               << "viscosity K_h) and an explicit vertical-diffusion stage "
                               << "(erf.vert_implicit_fac has a zero): the terrain-metric term "
                               << "K_h h^2 d2/dz2 stays explicit there and limits the time step on "
                               << "steep slopes.  To relax it, use erf.vert_implicit_fac = 1 1 1 with "
                               << "erf.implicit_terrain_metric = true, or erf.smag2d_slope_limiter = true "
                               << "for Smagorinsky2D; erf.diffusive_dt_check reports the limit."
                               << std::endl;
            }
        }
 
        // Countergradient validation: ensure required implicit diffusion is enabled
        for (int lev = 0; lev <= max_level; lev++) {
            if (turbChoice[lev].enable_mrf_countergradient) {
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                    vert_implicit_fac[lev][0] > zero ||
                    vert_implicit_fac[lev][1] > zero ||
                    vert_implicit_fac[lev][2] > zero,
                    "enable_mrf_countergradient requires vert_implicit_fac > 0");
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                    implicit_thermal_diffusion,
                    "enable_mrf_countergradient requires implicit_thermal_diffusion = true");
                if (moisture_type != MoistureType::None) {
                    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                        implicit_moisture_diffusion,
                        "enable_mrf_countergradient with moisture requires implicit_moisture_diffusion = true");
                }
            }
        }
 
        for (int lev = 0; lev <= max_level; lev++) {
            turbChoice[lev].check_params(phys_bc_type);
        }
    }
 
    /**
     * @brief Print the configured solver, physics, and transport options.
     * @param max_level Maximum AMR level configured for the run.
     * @param pp_prefix ParmParse prefix used to display nested choices.
     */
    void display (int max_level, std::string pp_prefix)
    {
        amrex::Print() << "SOLVER CHOICE: " << std::endl;
        if (perturbation_auto) {
            amrex::Print() << "Automatic LES-nest perturbations: parent update period = "
                           << perturbation_auto_period << " s, PBL height fraction = "
                           << perturbation_auto_pbl_fraction << std::endl;
        }
        for (int lev = 0; lev <= max_level; lev++) {
            amrex::Print() << "At level " << lev << " : " << std::endl;
            if (anelastic[lev]) {
                amrex::Print() << "   anelastic with no substepping, using the "
                               << ((anelastic_type[lev] == AnelasticType::MidPoint) ? "MidPoint" : "RK2")
                               << " integrator" << std::endl;
            } else {
               if (substepping_type[lev] == SubsteppingType::None) {
                   amrex::Print() << "   compressible with no substepping" << std::endl;
               } else if (substepping_type[lev] == SubsteppingType::Implicit) {
                   amrex::Print() << "   compressible with implicit substepping" << std::endl;
               }
           }
           if (fixed_density[lev]) {
               amrex::Print() << "   and fixed density" << std::endl;
           }
        }
 
        for (int lev = 0; lev <= max_level; lev++) {
            amrex::Print() << "vert_implicit_fac at level " << lev << " : "
                           << vert_implicit_fac[lev][0] << " "
                           << vert_implicit_fac[lev][1] << " "
                           << vert_implicit_fac[lev][2];
            if (vert_implicit_fac[lev][0] > zero ||
                vert_implicit_fac[lev][1] > zero ||
                vert_implicit_fac[lev][2] > zero)
            {
                amrex::Print() << " (theta = " << implicit_thermal_diffusion
                               << ", moisture = " << implicit_moisture_diffusion
                               << ", tke = " << implicit_ke_diffusion
                               << ", momenta = " << implicit_momentum_diffusion;
                if (implicit_terrain_metric) {
                    amrex::Print() << "; terrain metric implicit";
                }
                amrex::Print() << ")";
            }
        }
        amrex::Print() << std::endl;
        amrex::Print() << "use_coriolis                : " << use_coriolis << std::endl;
        amrex::Print() << "use_gravity                 : " << use_gravity << std::endl;
 
        if (moisture_type == MoistureType::SAM) {
            amrex::Print() << "Moisture Model: SAM" << std::endl;
        } else if (moisture_type == MoistureType::SAM_NoIce) {
            amrex::Print() << "Moisture Model: SAM No Ice" << std::endl;
        } else if (moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
            amrex::Print() << "Moisture Model: SAM No Precip No Ice" << std::endl;
        } else if (moisture_type == MoistureType::Morrison) {
            amrex::Print() << "Moisture Model: Morrison" << std::endl;
        } else if (moisture_type == MoistureType::Morrison_NoIce) {
            amrex::Print() << "Moisture Model: Morrison_NoIce" << std::endl;
        } else if (moisture_type == MoistureType::WSM6) {
            amrex::Print() << "Moisture Model: WSM6" << std::endl;
        } else if (moisture_type == MoistureType::WDM6) {
            amrex::Print() << "Moisture Model: WDM6" << std::endl;
        } else if (moisture_type == MoistureType::Kessler) {
            amrex::Print() << "Moisture Model: Kessler" << std::endl;
        } else if (moisture_type == MoistureType::Kessler_NoRain) {
            amrex::Print() << "Moisture Model: Kessler No Rain" << std::endl;
        } else if (moisture_type == MoistureType::SatAdj) {
            amrex::Print() << "Moisture Model: Saturation Adjustment" << std::endl;
        } else if (moisture_type == MoistureType::SuperDroplets) {
            amrex::Print() << "Moisture Model: SuperDroplets" << std::endl;
        } else if (moisture_type == MoistureType::MoistNoCondensation) {
            amrex::Print() << "Moisture Model: MoistNoCondensation" << std::endl;
        } else if (moisture_type == MoistureType::SBM) {
            amrex::Print() << "Moisture Model: SBM mapped spectral advection (M3)" << std::endl;
        } else {
            amrex::Print() << "Moisture Model: None" << std::endl;
        }
 
        if (terrain_type == TerrainType::StaticFittedMesh) {
            amrex::Print() << "Terrain Type: StaticFittedMesh" << std::endl;
        } else if (terrain_type == TerrainType::MovingFittedMesh) {
            amrex::Print() << "Terrain Type: MovingFittedMesh" << std::endl;
        } else if (terrain_type == TerrainType::EB) {
            amrex::Print() << "Terrain Type: EB" << std::endl;
            ebChoice.display();
        } else if (terrain_type == TerrainType::ImmersedForcing) {
            amrex::Print() << "Terrain Type: ImmersedForcing" << std::endl;
        } else {
            amrex::Print() << "Terrain Type: None" << std::endl;
        }
 
        amrex::Print() << "Flat fitted terrain: "
                       << (flat_terrain ? "true" : "false") << std::endl;
 
        if (buildings_type == BuildingsType::ImmersedForcing) {
            amrex::Print() << "Buildings Type: ImmersedForcing" << std::endl;
        } else {
            amrex::Print() << "Buildings Type: None" << std::endl;
        }
 
        if (mesh_type == MeshType::ConstantDz) {
            amrex::Print() << "   Mesh Type: ConstantDz" << std::endl;
        } else if (mesh_type == MeshType::StretchedDz) {
            amrex::Print() << "   Mesh Type: StretchedDz" << std::endl;
        } else if (mesh_type == MeshType::VariableDz) {
            amrex::Print() << "   Mesh Type: VariableDz" << std::endl;
        } else {
            amrex::Error("No mesh_type set!");
        }
 
        amrex::Print() << "ABL Driver Type: " << std::endl;
        if (abl_driver_type == ABLDriverType::None) {
            amrex::Print() << "    None" << std::endl;
        } else if (abl_driver_type == ABLDriverType::PressureGradient) {
            amrex::Print() << "    Pressure Gradient "
                           << amrex::RealVect(abl_pressure_grad[0],abl_pressure_grad[1],abl_pressure_grad[2])
                           << std::endl;
        } else if (abl_driver_type == ABLDriverType::GeostrophicWind) {
            amrex::Print() << "    Geostrophic Wind "
                           << amrex::RealVect(abl_geo_forcing[0],abl_geo_forcing[1],abl_geo_forcing[2])
                           << std::endl;
        }
 
        if (max_level > 0) {
            amrex::Print() << "Coupling Type: " << std::endl;
            if (coupling_type == CouplingType::TwoWay) {
                amrex::Print() << "    Two-way" << std::endl;
            } else if (coupling_type == CouplingType::OneWay) {
                amrex::Print() << "    One-way" << std::endl;
            }
        }
 
        if (rad_type == RadiationType::RRTMGP) {
            amrex::Print() << "Radiation Model: RRTMGP" << std::endl;
        } else if (rad_type == RadiationType::Simple) {
            amrex::Print() << "Radiation Model: Simple" << std::endl;
        } else if (rad_type == RadiationType::TwoStream) {
            amrex::Print() << "Radiation Model: TwoStream" << std::endl;
        } else {
            amrex::Print() << "Radiation Model: None" << std::endl;
        }
 
        amrex::Print() << "Gradp_type                  : " << gradp_type << std::endl;
 
        for (int lev = 0; lev <= max_level; lev++) {
            amrex::Print() << "Buoyancy_type at level " << lev << " : " << buoyancy_type[lev] << std::endl;
        }
 
            advChoice.display(pp_prefix);
           diffChoice.display();
         spongeChoice.display();
        dampingChoice.display();
 
        for (int lev = 0; lev <= max_level; lev++) {
            turbChoice[lev].display(lev);
        }
    }
 
    /**
     * @brief Build constant-latitude Coriolis forcing parameters.
     * @param pp_prefix ParmParse prefix for the ERF input namespace.
     */
    void build_coriolis_forcings_const_lat (std::string pp_prefix)
    {
        amrex::ParmParse pp(pp_prefix);
 
        // Read the rotational time period (in seconds)
        double rot_time_period = 86400.0;
        pp.queryAdd("rotational_time_period", rot_time_period);
 
        coriolis_factor = static_cast<amrex::Real>(4.0 * PI / rot_time_period);
 
        pp.queryAdd("coriolis_3d", coriolis_3d);
 
        // Convert to radians
        amrex::Real latitude_for_coriolis = amrex::Real(90.0);
        pp.queryAdd("latitude", latitude_for_coriolis);
        latitude_for_coriolis *= (PI/amrex::Real(180.));
        sinphi = std::sin(latitude_for_coriolis);
 
        // Note: cosphi defaults to 0
        if (coriolis_3d) {
            cosphi = std::cos(latitude_for_coriolis);
        }
 
        amrex::Print() << "Coriolis frequency, f = " << coriolis_factor * sinphi << " 1/s" << std::endl;
 
        if (abl_driver_type == ABLDriverType::GeostrophicWind) {
            // Read in the geostrophic wind -- we only use this to construct
            //     the forcing term so no need to keep it
            amrex::Vector<amrex::Real> abl_geo_wind(3);
            pp.queryarr("abl_geo_wind",abl_geo_wind);
 
            if(!pp.queryAdd("abl_geo_wind_table",abl_geo_wind_table)) {
                abl_geo_forcing = {
                    -coriolis_factor * (abl_geo_wind[1]*sinphi - abl_geo_wind[2]*cosphi),
                     coriolis_factor *  abl_geo_wind[0]*sinphi,
                    -coriolis_factor *  abl_geo_wind[0]*cosphi
                };
            } else {
                amrex::Print() << "NOTE: abl_geo_wind_table provided, ignoring input abl_geo_wind" << std::endl;
            }
        }
    }
 
    /**
     * @brief Read one integer value or one value per AMR level.
     * @param max_level Maximum AMR level configured for the run.
     * @param string_to_read ParmParse key to read.
     * @param vec_to_fill Destination vector with one entry per level.
     * @param default_int Default value used when the key is absent.
     * @return True if the key was present in the inputs.
     */
    bool read_int_string (int max_level, const char* string_to_read,
                          amrex::Vector<int>& vec_to_fill, int default_int)
    {
        amrex::ParmParse pp("erf");
        int nvals = pp.countval(string_to_read);
        AMREX_ALWAYS_ASSERT(nvals == 0 || nvals == 1 || nvals >= max_level+1);
        amrex::Vector<int> temp; temp.resize(nvals);
        pp.queryarr(string_to_read,temp);
 
        if (vec_to_fill.size() < max_level+1) {
            vec_to_fill.resize(max_level+1);
        }
 
        if (nvals == 0) {
            for (int i = 0; i <= max_level; ++i) vec_to_fill[i] = default_int;
        } else if (nvals == 1) {
            for (int i = 0; i <= max_level; ++i) vec_to_fill[i] = temp[0];
        } else {
            for (int i = 0; i <= max_level; ++i) vec_to_fill[i] = temp[i];
        }
        return (nvals > 0);
    }
 
    inline static
    InitType init_type = InitType::None;
 
    inline static
    SoundingType sounding_type = SoundingType::Ideal;
 
    inline static
    TerrainType terrain_type = TerrainType::None;
 
    inline static
    bool flat_terrain = false;
 
    inline static
    BuildingsType buildings_type = BuildingsType::None;
 
    inline static
    bool buildings_fc_factories = false;
 
    inline static
    bool use_real_bcs = false;
 
    inline static
    MeshType mesh_type = MeshType::ConstantDz;
 
    /**
     * @brief Override the globally selected mesh type.
     * @param new_mesh_type Mesh type to store.
     */
    static
    void set_mesh_type (MeshType new_mesh_type)
    {
        mesh_type = new_mesh_type;
    }
 
    AdvChoice     advChoice;
    DiffChoice    diffChoice;
    DampingChoice dampingChoice;
    SpongeChoice  spongeChoice;
    amrex::Vector<TurbChoice>  turbChoice;
    EBChoice      ebChoice;
    RadChoice     radChoice;
    RRTMGPSolarChoice rrtmgpSolarChoice;
    TopographicShadingChoice topographicChoice;
 
    int         force_stage1_single_substep = 1;
 
    amrex::Vector<SubsteppingType> substepping_type;
    amrex::Vector<int> anelastic;
    amrex::Vector<AnelasticType> anelastic_type;
    amrex::Vector<int> fixed_density;
    amrex::Vector<int> project_initial_velocity;
    amrex::Vector<int> buoyancy_type;
    amrex::Vector<int> is_land;
 
    // do some extra CFL diagnostics for compressible with substepping
    bool substepping_diag = false;
 
    // time off-centering coefficient, > 0 for forward weighting (i.e., bias
    // towards the future time step)
    amrex::Real beta_s = amrex::Real(0.1);
 
    // if true, beta_s varies on the w-faces as max(beta_s_min, min(beta_s_max, h_xi^2 + h_eta^2))
    // and beta_s itself is not used (implicit substepping only)
    bool use_variable_beta_s = false;
    amrex::Real beta_s_min = amrex::Real(0.1);
    amrex::Real beta_s_max = amrex::Real(0.8);
 
    // divergence damping coefficient: how far (rho theta) is projected forward
    // in the horizontal momentum equations of the acoustic substepping
    amrex::Real beta_d = amrex::Real(0.1);
 
    // This controls the time-centering of the *vertical* differences in the diffusive term for
    // theta, u and v
    // 0: fully explicit
    // 1: fully implicit
    amrex::Vector<amrex::Vector<amrex::Real>> vert_implicit_fac; // one value per RK stage
 
    // if any vert_implicit_fac > 0, then the following apply:
    bool implicit_thermal_diffusion  = true;
    bool implicit_moisture_diffusion = true;
    bool implicit_ke_diffusion       = true;
    bool implicit_momentum_diffusion = true;
    bool implicit_before_substep     = true;
    // On a terrain-fitted mesh, also treat the compact terrain-metric part of the vertical
    // diffusion (K_h M d2/dz2, ERF_TerrainImplicitMetric.H) implicitly, in every implicit solve
    bool implicit_terrain_metric     = false;
 
    int         ncorr               = 1;
    std::string wall_dist_type      = "poisson";
    TerrainPoissonSolver  terrain_poisson_solver  = TerrainPoissonSolver::gmres_fft;
#ifdef AMREX_USE_FLOAT
    amrex::Real poisson_abstol      = amrex::Real(1e-6);
    amrex::Real poisson_reltol      = amrex::Real(1e-6);
#else
    amrex::Real poisson_abstol      = amrex::Real(1e-8);
    amrex::Real poisson_reltol      = amrex::Real(1e-8);
#endif
 
    bool        test_mapfactor         = false;
    amrex::Real test_mapfactor_amp     = amrex::Real(0);
 
    int         gradp_type             = 0;
 
    bool        use_pert_pres_gradient = true;
 
    // Specify what additional physics/forcing modules we use
    bool        use_gravity            = false;
    bool        use_coriolis           = false;
    bool        coriolis_3d            = true;
 
    bool        four_stream_radiation  = false;
 
    // Specify whether to apply other various source terms on substep only
    bool        immersed_forcing_substep  = false;
    bool        forest_substep            = false;
 
    // immersed forcing parameters
    amrex::Real if_Cd_momentum        = amrex::Real(50.0);
    amrex::Real if_Cd_scalar          = amrex::Real(10.0);
    // Use a point-implicit (linearly-implicit) formulation of the immersed-forcing
    // drag. Defaults to false (explicit forward-Euler source).
    bool        if_implicit_drag      = false;
    bool        if_implicit_projection = false;  // erf.if_implicit_projection: solid drag (terrain or buildings) inside the anelastic projection
    // immersed forcing MOST parameters.
    amrex::Real if_z0                    = amrex::Real(0.1);   // [m]
    amrex::Real if_surf_temp_flux        = amrex::Real(1e-8);  // [K m/s]
    amrex::Real if_init_surf_temp        = zero;               // [K]
    amrex::Real if_surf_heating_rate     = zero;               // [K/hr]
    amrex::Real if_Olen_in               = amrex::Real(1e-8);  // [m]
    bool        if_use_most              = false;
    bool        if_snap_partial_cells    = false;
    bool        if_stability_correction  = false;
    bool        if_fraction_stress       = false;  // erf.if_wall_form = fraction_stress
    int         if_wall_face_spacing     = 2;      // erf.if_wall_face_spacing: 0 grid, 1 physical, 2 log (default)
    bool        if_wall_tke              = false;  // erf.if_wall_tke: wall value of the k-eqn TKE in the wall cell
    amrex::Real if_wall_tke_time_factor  = amrex::Real(1.0); // erf.if_wall_tke_time_factor: 1 pins, > 1 relaxes over that many steps
    bool        if_flat_wall_dissipation = false;  // erf.if_flat_wall_dissipation (default: = if_fraction_stress)
    amrex::Real if_ws_floor              = 0.001;  // [m/s]
    amrex::Real if_damp_alpha            = 0.5;    // [-]
    amrex::Real if_stability_wind_floor  = amrex::Real(0.0);  // [m/s]
    amrex::Real if_psi_cap_factor        = amrex::Real(1.0);  // [-]
    amrex::Real if_small_volfrac         = amrex::Real(0.005);  // [-] eb2.small_volfrac
    std::string if_average_down          = "plain";  // or fluid_mass
    amrex::Real if_cf_mismatch_warning   = amrex::Real(0.75); // erf.if_cf_mismatch_warning: c/f surface mismatch threshold
    bool        cf_loglaw_fill           = false;    // erf.cf_loglaw_fill: log-law share of the wall coarse face's flux
    amrex::Real cf_loglaw_z0             = amrex::Real(-1.0); // erf.cf_loglaw_z0 [m]; < 0: resolved from if_z0 / most.z0
 
    // This defaults to true but can be set to false for moving terrain cases only
    bool        use_lagged_delta_rt   = true;
 
    // Flag to rebalance states from wrfinput and wrfbdy
    bool        rebalance_wrf_input   = true;
 
    // Constants
    amrex::Real gravity;
    amrex::Real c_p = Cp_d; // specific heat at constant pressure for dry air [J/(kg-K)]
    amrex::Real rdOcp;
 
    // Staggered z levels for vertical grid stretching
    amrex::Real grid_stretching_ratio   = 0;
    amrex::Real dz0                     = zero;
 
    // Coriolis forcing
    amrex::Real coriolis_factor = zero;
    amrex::Real cosphi          = zero;
    amrex::Real sinphi          = zero;
 
    // User-specified forcings in problem definition
    bool custom_rhotheta_forcing  = false;
    bool custom_moisture_forcing  = false;
    bool custom_w_subsidence      = false;
    bool do_theta_advection       = true;   // Apply custom vertical subsidence to (rho*theta) equation when using custom w subsidence
    bool do_mom_advection         = true;  // Apply custom vertical subsidence to momentum equations when using custom w subsidence
    bool custom_geostrophic_profile = false;
    bool custom_forcing_prim_vars = false;
    bool spatial_rhotheta_forcing = false;
    bool spatial_moisture_forcing = false;
 
    // Do we use source terms to nudge the solution towards
    //    the time-varying data provided in input sounding files?
    bool nudging_from_input_sounding = false;
 
    // Do we nudge the solution towards observations at stations (met masts,
    //    lidars)?  See ERF_ObsNudging.H.
    bool nudging_from_observations = false;
 
    amrex::Real nudging_u_z1 = amrex::Real(-1.0e36);
    amrex::Real nudging_u_z2 = amrex::Real(1.0e36);
    amrex::Real nudging_t_z1 = zero;
    amrex::Real nudging_t_z2 = amrex::Real(10000.0);
    amrex::Real nudging_q_z1 = zero;
    amrex::Real nudging_q_z2 = amrex::Real(10000.0); // TODO: reasonable defaults here?
 
    bool nudging_u = true;      // Nudge u,v when nudging_from_input_sounding is true
    bool nudging_t = true;      // Nudge potential temperature when nudging_from_input_sounding is true
    bool nudging_q = true;      // Nudge water vapor when nudging_from_input_sounding is true
 
    bool large_scale_forcing = false;
    std::string lsf_file = "";
    amrex::Real lsf_timescale = zero;
 
    // MOST stress rotations
    bool use_rotate_surface_flux = false;
 
    // Should we use SHOC?
    bool use_eamxx_shoc = false;
    bool use_native_shoc = false;
 
    // Native shoc transport mode
    ShocTransportMode     shoc_transport_mode     = ShocTransportMode::StateUpdate;
    ShocMomentumTransport shoc_momentum_transport = ShocMomentumTransport::StateUpdate;
 
    // Transport the passive scalar?
    bool transport_scalar = true;
 
    // User wishes to output time averaged velocity fields
    bool time_avg_vel = false;
 
    // User wishes to output interval means and resolved turbulent moments
    bool compute_mean_vars = false;
 
    // Reset interval means after each plotfile or once at a requested time
    std::string mean_vars_reset_mode = "plotfile";
    amrex::Real mean_vars_reset_time = -one;
 
    // Building forcing and its landmask use the same horizontal resolution cutoff.
    [[nodiscard]] bool use_building_immersed_forcing (const amrex::Geometry& geom) const noexcept
    {
        return buildings_type == BuildingsType::ImmersedForcing &&
               std::sqrt(geom.CellSize(0) * geom.CellSize(1)) <=
                   amrex::Real(50.0) * (one + amrex::Real(1.e-6));
    }
 
    /**
     * @brief Query whether any turbulent perturbation mode is enabled on a level.
     * @param lev AMR level index.
     * @return True if source, direct, CPM, or CPM_W perturbations are enabled.
     */
    bool use_perturbation (int lev) const
    {
        return pert_type[lev] == PerturbationType::Source ||
               pert_type[lev] == PerturbationType::Direct ||
               pert_type[lev] == PerturbationType::CPM ||
               pert_type[lev] == PerturbationType::CPM_W;
    }
 
    /**
     * @brief Query whether direct turbulent perturbations are enabled on a level.
     * @param lev AMR level index.
     * @return True if direct, CPM, or CPM_W perturbations are enabled.
     */
    bool use_direct_perturbation (int lev) const
    {
        return pert_type[lev] == PerturbationType::Direct ||
               pert_type[lev] == PerturbationType::CPM ||
               pert_type[lev] == PerturbationType::CPM_W;
    }
 
    /**
     * @brief Query whether source-term turbulent perturbations are enabled on a level.
     * @param lev AMR level index.
     * @return True if source perturbations are enabled.
     */
    bool use_source_perturbation (int lev) const
    {
        return pert_type[lev] == PerturbationType::Source;
    }
 
    /**
     * @brief Query whether vertical-velocity CPM perturbations are enabled on a level.
     * @param lev AMR level index.
     * @return True if CPM_W perturbations are enabled.
     */
    bool use_wvel_perturbation (int lev) const
    {
        return pert_type[lev] == PerturbationType::CPM_W;
    }
 
    /**
     * @brief Query whether any SHOC-family PBL scheme is active.
     * @return True if native or EAMxx SHOC is active.
     */
    bool uses_shoc_family () const noexcept
    {
        return (use_eamxx_shoc || use_native_shoc);
    }
 
    /**
     * @brief Test whether the host owns momentum diffusion.
     * @return True if ERF dycore does vertical diffusion of momentum.
     */
    bool host_owns_vertical_momentum_diffusion () const noexcept
    {
        if (use_eamxx_shoc)   { return false; }
        if (!use_native_shoc) { return true;  }
        return shoc_disables_momentum_transport(shoc_momentum_transport);
    }
 
    /**
     * @brief Test whether the host owns scalar diffusion.
     * @return True if ERF dycore does vertical diffusion of scalars.
     */
    bool host_owns_vertical_scalar_diffusion () const noexcept
    {
        if (use_eamxx_shoc)   { return false; }
        if (!use_native_shoc) { return true;  }
        return false;
    }
 
    /**
     * @brief Query whether the active PBL scheme suppresses microphysics condensation.
     * @return True when the active PBL scheme owns condensation handling.
     */
    bool pbl_suppresses_microphysics_condensation () const noexcept
    {
        return uses_shoc_family();
    }
 
    /**
     * @brief Query whether any level uses a turbulent perturbation mode.
     * @return True if at least one configured level has source, direct, CPM, or CPM_W perturbations.
     */
    bool any_perturbation () const
    {
        for (const auto& type : pert_type) {
            if (type == PerturbationType::Source ||
                type == PerturbationType::Direct ||
                type == PerturbationType::CPM ||
                type == PerturbationType::CPM_W) {
                return true;
            }
        }
        return false;
    }
 
    // Prefer a PBL scheme's height even when the immediate parent is LES.
    // An all-LES hierarchy instead uses the nearest LES surface-layer diagnostic.
    int perturbation_pblh_level (int lev) const
    {
        for (int source = lev-1; source >= 0; --source) {
            if (turbChoice[source].pbl_type != PBLType::None) { return source; }
        }
        for (int source = lev-1; source >= 0; --source) {
            const auto& tc = turbChoice[source];
            if (tc.les_type != LESType::None && !tc.smag2d) { return source; }
        }
        return -1;
    }
 
    // Type of perturbation
    amrex::Vector<PerturbationType> pert_type;
    bool perturbation_auto{false};
    amrex::Real perturbation_auto_period{amrex::Real(600.0)}; // seconds between parent-based updates
    amrex::Real perturbation_auto_pbl_fraction{amrex::Real(0.9)};
 
    // Numerical diffusion
    bool use_num_diff{false};
    amrex::Real num_diff_coeff{zero};
 
    CouplingType coupling_type;
    MoistureType moisture_type;
    int sbm_nbins{4};
    int sbm_moment_mode{1};
    int sbm_cloud_rain_split{-1};
    bool sbm_test_native_qc_write_fault{false};
    bool sbm_test_donor_cfl_fault{false};
    amrex::Vector<amrex::Real> sbm_edges;
    amrex::Vector<amrex::Real> sbm_pivots;
    amrex::Vector<amrex::Real> sbm_fixture_initial_state;
    WindFarmType windfarm_type;
    WindFarmLocType windfarm_loc_type;
    LandSurfaceType lsm_type;
    UrbanType urban_type;
    amrex::Vector<int> urban_enabled_lev;
    RadiationType rad_type;
 
    /** @brief True for the models that run through the IRadiation interface
     *  (RRTMGP, Simple); the two-stream model has its own class. */
    [[nodiscard]] bool rad_uses_interface () const {
        return (rad_type == RadiationType::RRTMGP) || (rad_type == RadiationType::Simple);
    }
 
    /** @brief True when the radiation model's outputs feed the land model
     *  (downwelling shortwave, downwelling longwave, solar zenith angle). RRTMGP
     *  writes them for any land model through lsm_output_ptrs; the two-stream
     *  model writes the broadband three for Noah-MP only
     *  (TwoStreamRadiation::write_land_forcing), since SLM wants the visible /
     *  near-infrared direct / diffuse split it does not compute. Both register
     *  the outputs in ERF::init_stuff and write them in ERF::advance_radiation;
     *  every site that encodes that rule must use this predicate so wiring
     *  another model in (e.g. Simple, which already runs through the interface)
     *  is one edit. */
    [[nodiscard]] bool rad_feeds_lsm () const {
        return (rad_type == RadiationType::RRTMGP) ||
               (rad_type == RadiationType::TwoStream && lsm_type == LandSurfaceType::NOAHMP);
    }
 
    // Does an external ocean coupler supply SST? See the query in init_params.
    bool use_coupled_sst = false;
 
    ABLDriverType abl_driver_type;
    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> abl_pressure_grad;
    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> abl_geo_forcing;
    std::string abl_geo_wind_table;
    bool have_geo_wind_profile {false};
 
    bool variable_coriolis{false};
 
    int ave_plane {2};
 
    // Microphysics params
    MoistureComponentIndices moisture_indices;
 
    bool moisture_tight_coupling {false};
 
    std::string windfarm_loc_table, windfarm_spec_table, windfarm_spec_table_extra;
    std::string windfarm_blade_table, windfarm_airfoil_tables;
    amrex::Real sampling_distance_by_D = -one;
    amrex::Real turb_disk_angle  = -one;
    amrex::Real windfarm_x_shift = -one;
    amrex::Real windfarm_y_shift = -one;
 
    // Nudging factor for bdy sponging : F1 = 1/(nudge_factor * dt)
    amrex::Real bdy_nudge_factor = amrex::Real(10.0);
 
    // Use density reconstructed from standard WRF lateral boundary data.
    bool use_wrf_bdy_density = true;
 
    // Ingest WRF cloud water and active cloud ice at initialization and lateral boundaries.
    bool use_wrf_bdy_qc_qi = false;
 
    // Density-specific Davies factor; non-positive means use bdy_nudge_factor.
    amrex::Real bdy_rho_nudge_factor = amrex::Real(-1.0);
 
    // Which approach to nudging the moist variables in the boundary region?
    int bdy_moist_nudge_type = 1;
 
    // Average the input z-face heights onto the nodes rather than reconstructing them?
    bool avg_grid_faces_to_nodes {false};
 
    // When the ERF mesh is reconstructed rather than taken from the file
    // (avg_grid_faces_to_nodes = false), build zlevels_stag from the domain-mean
    // layer-thickness profile of the wrfinput file.  Set false to recover the older
    // construction: a single geometric stretch whose first layer is the thickest first
    // layer in the file and whose ratio is solved for so the column reaches z_top.
    bool wrfinput_zlevels_from_file {true};
 
    // Interpolate atmospheric state from coarse when creating finer level at later time?
    bool interp_atmos_from_coarse {false};
 
    // Use forest canopy model?
    bool do_forest_drag {false};
 
    // Fixed-leaf-temperature canopy heat exchange (disabled by default)
    bool        forest_biophysics       {false};
    bool        forest_biophysics_heat  {false};
    amrex::Real forest_leaf_theta_fixed {-one};
    // Canopy terms of the k-eqn RANS model (ERF_CanopyTurbStruct.H)
    CanopyTurbulence canopyTurb;
 
    // Enforce constant mass flux?
    amrex::Real const_massflux_u = zero;
    amrex::Real const_massflux_v = zero;
    amrex::Real const_massflux_tau = one;
    amrex::Real const_massflux_layer_lo = amrex::Real(-bogus_large_value);
    amrex::Real const_massflux_layer_hi = amrex::Real( bogus_large_value);
    int massflux_klo {0}; // these are updated in ERF.cpp
    int massflux_khi {0};
 
    std::string hindcast_boundary_data_dir, hindcast_surface_data_dir;
    amrex::Real hindcast_data_interval_in_hrs = -one;
    bool hindcast_lateral_forcing = false;
    bool hindcast_surface_bcs = false;
    amrex::Real hindcast_lateral_sponge_strength = -one, hindcast_lateral_sponge_length = -one;
    amrex::Real hindcast_zhi_sponge_strength = -one, hindcast_zhi_sponge_length = -one;
    bool hindcast_zhi_sponge_damping = false;
 
    bool io_hurricane_eye_tracker = false;
    amrex::Real hurricane_eye_latitude = -1e10, hurricane_eye_longitude = -1e10;
 
    bool is_init_for_ensemble = false;
    int n_ensemble = -1;
    amrex::Real ens_pert_correlated_radius = -1.0;
    amrex::Real ens_pert_amplitude = -1.0;
    std::string coarse_bckgnd_data_file;
};
 
/**
 * @var SolverChoice::init_type
 * @brief Initial-condition source selected for the run.
 * @var SolverChoice::sounding_type
 * @brief Interpretation used for input sounding profiles.
 * @var SolverChoice::terrain_type
 * @brief Terrain or immersed-boundary representation.
 * @var SolverChoice::flat_terrain
 * @brief Permit fixed-index planar averages after validating a flat fitted mesh.
 * @var SolverChoice::buildings_type
 * @brief Building representation.
 * @var SolverChoice::buildings_fc_factories
 * @brief Use native face-centered EB factories for immersed buildings on a static fitted mesh.
 * @var SolverChoice::use_real_bcs
 * @brief Whether real-data lateral boundary conditions are used.
 * @var SolverChoice::mesh_type
 * @brief Vertical mesh representation.
 * @var SolverChoice::advChoice
 * @brief Advection-related options.
 * @var SolverChoice::diffChoice
 * @brief Diffusion-related options.
 * @var SolverChoice::dampingChoice
 * @brief Damping-related options.
 * @var SolverChoice::spongeChoice
 * @brief Sponge-layer options.
 * @var SolverChoice::turbChoice
 * @brief Turbulence options for each AMR level.
 * @var SolverChoice::ebChoice
 * @brief Embedded-boundary options.
 * @var SolverChoice::force_stage1_single_substep
 * @brief Flag forcing a single substep in the first RK stage.
 * @var SolverChoice::substepping_type
 * @brief Acoustic substepping type for each AMR level.
 * @var SolverChoice::anelastic
 * @brief Per-level flag selecting anelastic dynamics.
 * @var SolverChoice::anelastic_type
 * @brief Per-level two-stage scheme used by the anelastic integrator.
 * @var SolverChoice::fixed_density
 * @brief Per-level flag selecting fixed-density treatment; always 1 at an anelastic level.
 * @var SolverChoice::project_initial_velocity
 * @brief Per-level flag for projecting the initial velocity.
 * @var SolverChoice::buoyancy_type
 * @brief Per-level buoyancy formulation selector.
 * @var SolverChoice::is_land
 * @brief Per-level land/water selector.
 * @var SolverChoice::substepping_diag
 * @brief Whether to emit extra CFL diagnostics for compressible substepping.
 * @var SolverChoice::beta_s
 * @brief Time off-centering coefficient for forward weighting.
 * @var SolverChoice::use_variable_beta_s
 * @brief Whether beta_s varies on the w-faces with the terrain slope.
 * @var SolverChoice::beta_s_min
 * @brief Lower bound of the spatially varying beta_s.
 * @var SolverChoice::beta_s_max
 * @brief Upper bound of the spatially varying beta_s.
 * @var SolverChoice::beta_d
 * @brief Divergence damping coefficient used in the acoustic substepping.
 * @var SolverChoice::vert_implicit_fac
 * @brief Per-level, per-stage implicit vertical diffusion factors.
 * @var SolverChoice::implicit_thermal_diffusion
 * @brief Whether implicit vertical thermal diffusion is included.
 * @var SolverChoice::implicit_moisture_diffusion
 * @brief Whether implicit vertical moisture diffusion is included.
 * @var SolverChoice::implicit_ke_diffusion
 * @brief Whether implicit vertical TKE diffusion is included.
 * @var SolverChoice::implicit_momentum_diffusion
 * @brief Whether implicit vertical momentum diffusion is included.
 * @var SolverChoice::implicit_before_substep
 * @brief Whether implicit diffusion is applied before acoustic substepping.
 * @var SolverChoice::implicit_terrain_metric
 * @brief Whether the compact terrain-metric part of the vertical diffusion is solved implicitly.
 * @var SolverChoice::ncorr
 * @brief Number of projection correction iterations.
 * @var SolverChoice::terrain_poisson_solver
 * @brief Projection solver on a terrain-fitted mesh: gmres_fft (default), mlmg or gmres_mlmg.
 * @var SolverChoice::poisson_abstol
 * @brief Absolute tolerance for the Poisson solve.
 * @var SolverChoice::poisson_reltol
 * @brief Relative tolerance for the Poisson solve.
 * @var SolverChoice::test_mapfactor
 * @brief Whether to use test map scale factors.
 * @var SolverChoice::test_mapfactor_amp
 * @brief Amplitude of the spatially varying test map scale factors on level 0 (0: not used).
 * @var SolverChoice::gradp_type
 * @brief Terrain-fitted horizontal pressure-gradient formulation.
 * @var SolverChoice::use_pert_pres_gradient
 * @brief Whether momentum equations use perturbational pressure gradients.
 * @var SolverChoice::use_gravity
 * @brief Whether gravitational forcing is enabled.
 * @var SolverChoice::use_coriolis
 * @brief Whether Coriolis forcing is enabled.
 * @var SolverChoice::coriolis_3d
 * @brief Whether the three-dimensional Coriolis terms are included.
 * @var SolverChoice::four_stream_radiation
 * @brief Whether the four-stream radiation approximation is enabled.
 * @var SolverChoice::immersed_forcing_substep
 * @brief Whether immersed-forcing source terms are applied only during substeps.
 * @var SolverChoice::forest_substep
 * @brief Whether canopy source terms are applied only during substeps.
 * @var SolverChoice::if_Cd_momentum
 * @brief Immersed-forcing drag coefficient for momentum.
 * @var SolverChoice::if_Cd_scalar
 * @brief Immersed-forcing drag coefficient for scalars.
 * @var SolverChoice::if_z0
 * @brief Immersed-forcing roughness length [m].
 * @var SolverChoice::if_surf_temp_flux
 * @brief Immersed-forcing surface temperature flux [K m/s].
 * @var SolverChoice::if_init_surf_temp
 * @brief Initial immersed-forcing surface temperature [K].
 * @var SolverChoice::if_fraction_stress
 * @brief erf.if_wall_form = fraction_stress: the wall stress of ERF_ImmersedWallCell.H in the wall
 *        cell of each column, with the wall face closed to diffusion; false is the legacy MOST
 *        relaxation at fixed 0.5/1.5 dz (the default, unchanged).
 * @var SolverChoice::if_wall_face_spacing
 * @brief Spacing the face above the fraction-stress wall cell differences over: 0 the cell
 *        height, 1 the centroid spacing, 2 (default) the log spacing (ib_wall::face_above_factor).
 * @var SolverChoice::if_wall_tke
 * @brief Hold the k-eqn TKE of the fraction-stress wall cell at its flat-ground wall value
 *        (ib_wall::wall_tke), as erf.dirichlet_k does over flat ground.
 * @var SolverChoice::if_wall_tke_time_factor
 * @brief 1 (the default) pins the wall cell after every stage and in the implicit TKE solve;
 *        a larger value relaxes toward the wall value as a source over that many steps.
 * @var SolverChoice::if_flat_wall_dissipation
 * @brief Partial fraction-stress wall cells dissipate TKE with the length scale of flat ground's
 *        first cell; on by default with the fraction-stress law.
 * @var SolverChoice::if_surf_heating_rate
 * @brief Immersed-forcing surface heating rate [K/hr].
 * @var SolverChoice::if_Olen_in
 * @brief Input Obukhov length for immersed-forcing MOST [m].
 * @var SolverChoice::if_use_most
 * @brief Whether immersed-forcing MOST is enabled.
 * @var SolverChoice::if_average_down
 * @brief Cell average down over immersed-forcing levels: "plain" (the default, the plain cell
 *        average) or "fluid_mass" (the fine cells weighted by their fluid mass, coarse cells solid
 *        on their own level kept; see if_average_down in ERF_Utils.H).
 * @var SolverChoice::if_cf_mismatch_warning
 * @brief Spread of the fine solid fractions under a coarse cell along a lateral coarse/fine face
 *        at or above which the start-up warning counts the cell (the levels describe the surface
 *        differently there); 0 turns the check off.
 * @var SolverChoice::cf_loglaw_fill
 * @brief Share the flux of the first coarse face above a wall among its fine sub-faces at the
 *        lateral coarse/fine faces by ln(1 + z/z0) at their heights above the wall, instead of the
 *        linear interpolation; the immersed surface or a surface-layer bottom, constant dz.
 * @var SolverChoice::cf_loglaw_z0
 * @brief Roughness length [m] of that share; by default erf.if_z0 with immersed forcing, else
 *        the inputs' erf.most.z0 or erf.zlo.most.z0 (0.1 m when neither is given).
 * @var SolverChoice::if_small_volfrac
 * @brief Immersed-forcing solid fraction below which a cell or face counts as fluid; read from
 *        eb2.small_volfrac (default 0.005), the same value the EB build uses.
 * @var SolverChoice::if_snap_partial_cells
 * When true, the immersed forcing reads the blanking snapped to solid (1) or
 * fluid (0) at half, so a building from a height map becomes the same
 * staircase of whole cells an exact box is: the wall law (and the thermal
 * conditions) sit on the boundary solid cells, roofs included, the drag on
 * the interior cells, and nothing on the sliver cells of the reader's
 * one-cell ramp (1-20 percent solid at the corners), on which the raw
 * fractions make the wall law and the drag disagree and grow a vertical
 * checkerboard over hours. The staircase puts full-strength forcing on every
 * wall face, so the switch also selects the point-implicit drag form
 * (if_implicit_drag), without which the explicit rate at rims and corners
 * drives the density negative within seconds. The momentum kernels take the
 * blanking of a face from the two cells it joins (solid if either is), so
 * the staircase is the same with and without the face-centred blanking of
 * the build; a face between a solid and a fluid cell is wall-normal and gets
 * the interior drag toward zero (no penetration) whichever row it lies in,
 * and a face between two solid cells the roof or wall law of its row (with
 * the full log-law target; the partial-cell weight of the target is not
 * applied) or the interior drag. Buildings only (erf.buildings_type = ImmersedForcing): the
 * terrain kernels keep the raw fractions, whose wall law is weighted by the
 * fluid fraction of the face. Default false, the raw fractions, for backward
 * compatibility (on an exact box a boundary cell then carries the wall law
 * and the interior drag together, and the roof only the drag); set it for
 * height-map buildings.
 * @var SolverChoice::if_stability_correction
 * @brief Whether immersed-forcing stability corrections are enabled.
 * @var SolverChoice::if_ws_floor
 * @brief Wind-speed floor [m/s] in the per-step change limit of the building wall-law forcing.
 * @var SolverChoice::if_damp_alpha
 * @brief Immersed-forcing damping coefficient.
 * @var SolverChoice::if_stability_wind_floor
 * @brief Wind-speed floor [m/s] behind the friction velocity of the wall law's stability estimate (0: none).
 * @var SolverChoice::if_psi_cap_factor
 * @brief Cap on psi_m in the immersed momentum wall law as a fraction of ln(z/z0) (1: none).
 * @var SolverChoice::use_lagged_delta_rt
 * @brief Whether to use lagged delta terms for moving terrain.
 * @var SolverChoice::rebalance_wrf_input
 * @brief Whether to rebalance states from wrfinput and wrfbdy data.
 * @var SolverChoice::gravity
 * @brief Effective gravitational acceleration.
 * @var SolverChoice::c_p
 * @brief Specific heat at constant pressure for dry air [J/(kg-K)].
 * @var SolverChoice::rdOcp
 * @brief Ratio of dry-air gas constant to c_p.
 * @var SolverChoice::grid_stretching_ratio
 * @brief Vertical grid stretching ratio.
 * @var SolverChoice::dz0
 * @brief Initial vertical grid spacing for stretched-grid setup.
 * @var SolverChoice::coriolis_factor
 * @brief Twice the planetary rotation rate used for Coriolis forcing.
 * @var SolverChoice::cosphi
 * @brief Cosine of the latitude used for Coriolis forcing.
 * @var SolverChoice::sinphi
 * @brief Sine of the latitude used for Coriolis forcing.
 * @var SolverChoice::custom_rhotheta_forcing
 * @brief Whether custom rho-theta forcing is enabled.
 * @var SolverChoice::custom_moisture_forcing
 * @brief Whether custom moisture forcing is enabled.
 * @var SolverChoice::custom_w_subsidence
 * @brief Whether custom vertical subsidence is enabled.
 * @var SolverChoice::do_theta_advection
 * @brief Whether custom vertical subsidence is applied to rho-theta.
 * @var SolverChoice::do_mom_advection
 * @brief Whether custom vertical subsidence is applied to momentum.
 * @var SolverChoice::custom_geostrophic_profile
 * @brief Whether a custom geostrophic wind profile is used.
 * @var SolverChoice::custom_forcing_prim_vars
 * @brief Whether custom forcing operates on primitive variables.
 * @var SolverChoice::spatial_rhotheta_forcing
 * @brief Whether spatially varying rho-theta forcing is enabled.
 * @var SolverChoice::spatial_moisture_forcing
 * @brief Whether spatially varying moisture forcing is enabled.
 * @var SolverChoice::nudging_from_input_sounding
 * @brief Whether solution fields are nudged toward input sounding data.
 * @var SolverChoice::nudging_from_observations
 * @brief Whether solution fields are nudged toward observations at stations.
 * @var SolverChoice::use_rotate_surface_flux
 * @brief Whether MOST surface fluxes are rotated with terrain.
 * @var SolverChoice::use_eamxx_shoc
 * @brief Whether any level uses the EAMxx SHOC PBL scheme.
 * @var SolverChoice::use_native_shoc
 * @brief Whether any level uses the native SHOC PBL scheme.
 * @var SolverChoice::transport_scalar
 * @brief Whether the passive scalar component is transported.
 * @var SolverChoice::time_avg_vel
 * @brief Whether time-averaged velocity fields are output.
 * @var SolverChoice::pert_type
 * @brief Turbulent perturbation type for each AMR level.
 * @var SolverChoice::perturbation_auto
 * @brief Select LES nests, using parent winds and a PBL ancestor's height,
 * or the nearest LES ancestor's diagnostic when no PBL scheme is present.
 * Defaults to true for WRFInput and Metgrid initialization, false otherwise.
 * @var SolverChoice::perturbation_auto_period
 * @brief Seconds between automatic parent-wind and PBL-height updates.
 * @var SolverChoice::perturbation_auto_pbl_fraction
 * @brief Fraction of the selected ancestor's PBL height occupied by perturbations.
 * @var SolverChoice::use_num_diff
 * @brief Whether sixth-order numerical diffusion is enabled.
 * @var SolverChoice::num_diff_coeff
 * @brief Numerical diffusion coefficient after input scaling.
 * @var SolverChoice::coupling_type
 * @brief Multilevel coupling strategy.
 * @var SolverChoice::moisture_type
 * @brief Moisture or microphysics model.
 * @var SolverChoice::windfarm_type
 * @brief Wind-farm model.
 * @var SolverChoice::windfarm_loc_type
 * @brief Coordinate system for wind-farm turbine locations.
 * @var SolverChoice::lsm_type
 * @brief Land-surface model.
 * @var SolverChoice::rad_type
 * @brief Radiation model.
 * @var SolverChoice::abl_driver_type
 * @brief Atmospheric boundary-layer driver type.
 * @var SolverChoice::abl_pressure_grad
 * @brief Applied ABL pressure-gradient forcing vector.
 * @var SolverChoice::abl_geo_forcing
 * @brief Applied geostrophic-wind forcing vector.
 * @var SolverChoice::abl_geo_wind_table
 * @brief Path to a time-varying geostrophic wind table.
 * @var SolverChoice::have_geo_wind_profile
 * @brief Whether a geostrophic wind profile has been configured.
 * @var SolverChoice::variable_coriolis
 * @brief Whether spatially varying Coriolis forcing is enabled.
 * @var SolverChoice::ave_plane
 * @brief Averaging plane index used by diagnostics.
 * @var SolverChoice::moisture_indices
 * @brief Index map of the moisture data carried by the active scheme: conserved-state
 *        components for the species and qmoist slots for the moist diagnostics.
 * @var SolverChoice::moisture_tight_coupling
 * @brief Whether moisture updates use tight coupling.
 * @var SolverChoice::windfarm_loc_table
 * @brief Wind-farm location table path.
 * @var SolverChoice::windfarm_spec_table
 * @brief Wind-farm turbine specification table path.
 * @var SolverChoice::windfarm_spec_table_extra
 * @brief Additional wind-farm turbine specification table path.
 * @var SolverChoice::windfarm_blade_table
 * @brief Wind-farm blade table path.
 * @var SolverChoice::windfarm_airfoil_tables
 * @brief Wind-farm airfoil table paths.
 * @var SolverChoice::sampling_distance_by_D
 * @brief Turbine inflow sampling distance normalized by rotor diameter.
 * @var SolverChoice::turb_disk_angle
 * @brief Turbine disk angle from the x axis [deg].
 * @var SolverChoice::windfarm_x_shift
 * @brief Wind-farm x-coordinate shift.
 * @var SolverChoice::windfarm_y_shift
 * @brief Wind-farm y-coordinate shift.
 * @var SolverChoice::bdy_nudge_factor
 * @brief Boundary sponge nudging factor.
 * @var SolverChoice::use_wrf_bdy_qc_qi
 * @brief Whether WRF cloud water and active cloud ice are ingested at real boundaries.
 * @var SolverChoice::bdy_moist_nudge_type
 * @brief Moist-variable nudging strategy in boundary regions.
 * @var SolverChoice::do_forest_drag
 * @brief Whether forest canopy drag is enabled.
 * @var SolverChoice::const_massflux_u
 * @brief Target constant mass flux in the x direction.
 * @var SolverChoice::const_massflux_v
 * @brief Target constant mass flux in the y direction.
 * @var SolverChoice::const_massflux_tau
 * @brief Relaxation time scale for constant-mass-flux forcing.
 * @var SolverChoice::const_massflux_layer_lo
 * @brief Lower height of the constant-mass-flux forcing layer.
 * @var SolverChoice::const_massflux_layer_hi
 * @brief Upper height of the constant-mass-flux forcing layer.
 * @var SolverChoice::massflux_klo
 * @brief Lower vertical index for constant-mass-flux forcing.
 * @var SolverChoice::massflux_khi
 * @brief Upper vertical index for constant-mass-flux forcing.
 * @var SolverChoice::hindcast_boundary_data_dir
 * @brief Directory containing hindcast lateral boundary data.
 * @var SolverChoice::hindcast_surface_data_dir
 * @brief Directory containing hindcast surface boundary data.
 * @var SolverChoice::hindcast_data_interval_in_hrs
 * @brief Time interval between hindcast data files [hr].
 * @var SolverChoice::hindcast_lateral_forcing
 * @brief Whether hindcast lateral forcing is enabled.
 * @var SolverChoice::hindcast_surface_bcs
 * @brief Whether hindcast surface boundary conditions are enabled.
 * @var SolverChoice::hindcast_lateral_sponge_strength
 * @brief Strength of the hindcast lateral sponge layer.
 * @var SolverChoice::hindcast_lateral_sponge_length
 * @brief Length of the hindcast lateral sponge layer.
 * @var SolverChoice::hindcast_zhi_sponge_strength
 * @brief Strength of the upper hindcast sponge layer.
 * @var SolverChoice::hindcast_zhi_sponge_length
 * @brief Length of the upper hindcast sponge layer.
 * @var SolverChoice::hindcast_zhi_sponge_damping
 * @brief Whether upper-boundary hindcast sponge damping is enabled.
 * @var SolverChoice::io_hurricane_eye_tracker
 * @brief Whether hurricane-eye tracking output is enabled.
 * @var SolverChoice::hurricane_eye_latitude
 * @brief Initial hurricane-eye latitude.
 * @var SolverChoice::hurricane_eye_longitude
 * @brief Initial hurricane-eye longitude.
 * @var SolverChoice::is_init_for_ensemble
 * @brief Whether initialization is for an ensemble simulation.
 * @var SolverChoice::n_ensemble
 * @brief Number of ensemble members.
 * @var SolverChoice::ens_pert_correlated_radius
 * @brief Spatial correlation radius for ensemble perturbations.
 * @var SolverChoice::ens_pert_amplitude
 * @brief Amplitude of ensemble initialization perturbations.
 * @var SolverChoice::coarse_bckgnd_data_file
 * @brief Coarse background data file used for ensemble initialization.
 */
#endif