Fix the IBSEB build break and legacy surface-layer input detection

Docs/sphinx_doc/SurfaceLayer.rst

@@ -307,8 +307,29 @@ For example,
     erf.ylo.most.radius = 0
 
 
+The face name is omitted when ``zlo`` is the only surface-layer face, so an input
+file written before this capability existed keeps working unchanged: with
+``zlo.type = "surface_layer"`` and no other surface-layer face, ``erf.most.*`` and
+``erf.surface_layer.*`` are read as before.
+
+As soon as a second surface-layer face is declared, every face is expected to
+qualify its inputs, ``zlo`` included.  Unqualified inputs are still honored in
+that case -- they are applied to ``zlo``, with a warning asking that they be
+migrated to ``erf.zlo.most.*`` and ``erf.zlo.surface_layer.*`` -- so that adding a
+wall to a working deck cannot silently fall back to the built-in defaults.  Giving
+both spellings at once (``erf.most.z0`` together with ``erf.zlo.most.z0``, say) is
+an error and aborts, since there is no way to tell which one was meant.
+
 Note that not all existing options are supported when using MOST on other faces (such as interpolation and time averaging).
-Currently the `MOENG` flux type is supported on all faces.
+Currently the `MOENG` flux type is supported on all faces; ``bulk_coeff``,
+``custom`` and ``rico`` abort anywhere else.  Of the surface thermodynamic
+pathways, only the prescribed surface temperature (``most.surf_temp``) carries the
+orientation through to the flux, so a prescribed surface heat flux
+(``most.surf_temp_flux``), the fully adiabatic case, and any variable sea
+roughness (``most.roughness_type_sea`` other than ``constant`` over water cells)
+abort on a face other than zlo rather than run with the zlo sign convention.  The
+MOST PBL-height diagnostic (``most.pblh_calc``) and the free-convection correction
+it feeds (``most.include_wstar``) are likewise zlo-only.
 Lateral (x/y) surface-layer boundaries currently require every grid on that
 level to span its full vertical extent.  Grids decomposed in z and
 partial-height refined grids are not supported; configure

Source/ERF.H

@@ -4,6 +4,7 @@
 #include <string>
 #include <limits>
 #include <memory>
+#include <set>
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -797,6 +798,14 @@ public:
 
     std::string pp_prefix {"erf"};
 
+    // The inputs table as the user wrote it, captured at the top of
+    // ReadParameters before anything in ERF can add to it.  ParmParse::queryAdd
+    // inserts the default it was handed whenever a key is absent, so by the time
+    // initialization is under way the table no longer distinguishes a key the
+    // user typed from a key some other query defaulted (see
+    // has_surface_layer_inputs).
+    std::set<std::string> user_specified_inputs;
+
     void fill_from_bndryregs (const amrex::Vector<amrex::MultiFab*>& mfs,
                               double time);
 
@@ -901,7 +910,7 @@ private:
     void ReadParameters ();
     void ParameterSanityChecks ();
 
-    static bool has_surface_layer_inputs (const std::string& prefix);
+    [[nodiscard]] bool has_surface_layer_inputs (const std::string& prefix) const;
 
     // set covered coarse cells/faces to be the average of overlying fine cells/faces
     void AverageDown ();

Source/ERF.cpp

@@ -1266,9 +1266,10 @@ ERF::InitData_post ()
     }
 
     // With multiple surface-layer faces, face-qualified prefixes are normally
-    // required.  Preserve the historical unqualified zlo inputs when a user
-    // adds another surface-layer face, but detect the inputs before any
-    // SurfaceLayer constructor can insert queryAdd defaults into the table.
+    // required.  Preserve the historical unqualified zlo inputs when a user adds
+    // another surface-layer face.  What the user actually wrote is read from the
+    // snapshot taken in ReadParameters; see has_surface_layer_inputs for why the
+    // live ParmParse table cannot be asked this question here.
     bool use_legacy_zlo_prefix = false;
     if (n_faces > 1 &&
         phys_bc_type[Orientation::zlo()] == ERF_BC::surface_layer) {
@@ -2560,6 +2561,12 @@ ERF::init_only (int lev, double elapsed_time)
 void
 ERF::ReadParameters ()
 {
+    // Record the inputs table before anything below can modify it.  Every query
+    // here and in the classes built afterwards is a queryAdd, which inserts the
+    // default it was handed when the key is absent, so this is the only point at
+    // which "the key is in the table" still means "the user wrote it".
+    user_specified_inputs = ParmParse::getEntries(pp_prefix);
+
     std::string prob_name = "Undefined";
     ParmParse pp_pn("erf");
     pp_pn.queryAdd("prob_name", prob_name);
@@ -3729,14 +3736,26 @@ ERF::check_mesh_type(int lev)
    }
 }
 
+/**
+ * Whether the user wrote any surface-layer input under the given prefix.
+ *
+ * This asks the snapshot taken in ReadParameters, not the live ParmParse table.
+ * The live table cannot answer the question: ERF_InputSoundingData.H queries
+ * erf.most.surf_temp and erf.most.surf_moist with queryAdd while reading the
+ * sounding, which runs during InitData_pre and leaves both keys in the table
+ * holding their negative sentinels.  Testing the live table would therefore
+ * report unqualified MOST inputs for every input_sounding run, whether or not
+ * the user wrote any.
+ *
+ * @param[in] prefix the prefix to test, e.g. "erf" or "erf.zlo"
+ */
 bool
-ERF::has_surface_layer_inputs (const std::string& prefix)
+ERF::has_surface_layer_inputs (const std::string& prefix) const
 {
-    const auto entries = ParmParse::getEntries(prefix);
     const std::string most_prefix = prefix + ".most.";
     const std::string surface_layer_prefix = prefix + ".surface_layer.";
 
-    for (const auto& key : entries) {
+    for (const auto& key : user_specified_inputs) {
         if (key.compare(0, most_prefix.size(), most_prefix) == 0 ||
             key.compare(0, surface_layer_prefix.size(), surface_layer_prefix) == 0) {
             return true;

Source/ImmersedBoundarySEB/ERF_IBSEB.cpp

@@ -140,10 +140,14 @@ ERF::ibseb_advance (int lev, Real time, Real dt, const MultiFab& cons,
     const MultiFab* olen2d = nullptr;
     const MultiFab* pblh2d = nullptr;
     Real z_i_bulk = 0.0;
-    if (m_SurfaceLayer && ibseb_params.stability_correction) { olen2d = m_SurfaceLayer->get_olen(lev); }
+    // The surface layer now exists per domain face.  What the wall function wants
+    // here is the ground beneath the buildings, so take zlo; a surface layer on a
+    // lateral or upper wall says nothing about the stability of this column.
+    const auto& ground_sl = m_SurfaceLayer[Orientation(Direction::z, Orientation::low)];
+    if (ground_sl && ibseb_params.stability_correction) { olen2d = ground_sl->get_olen(lev); }
     if (ibseb_params.convective_velocity == "deardorff") {
-        if (m_SurfaceLayer && m_SurfaceLayer->computes_pblh() && ibseb_params.z_i_mode == "pblh") {
-            pblh2d = m_SurfaceLayer->get_pblh(lev);
+        if (ground_sl && ground_sl->computes_pblh() && ibseb_params.z_i_mode == "pblh") {
+            pblh2d = ground_sl->get_pblh(lev);
         }
         z_i_bulk = (ibseb_params.z_i_mode == "fixed") ? ibseb_params.z_i
                                                      : ibseb_bulk_richardson_height(lev, cons, xvel, yvel);