ERF_SurfaceLayer.H

Go to the documentation of this file.

#ifndef ERF_SURFACELAYER_H
#define ERF_SURFACELAYER_H
 
#include <fstream>
#include <sstream>
#include <string>
 
#include "AMReX_Geometry.H"
#include "AMReX_ParmParse.H"
#include "AMReX_FArrayBox.H"
#include "AMReX_MultiFab.H"
#include "AMReX_iMultiFab.H"
#include "AMReX_MFInterpolater.H"
 
#include "ERF_IndexDefines.H"
#include "ERF_Constants.H"
#include "ERF_MOSTAverage.H"
#include "ERF_PlanarBoundary.H"
#include "ERF_MOSTStress.H"
#include "ERF_MOSTUtils.H"
#include "ERF_EBMOSTStress.H"
#include "ERF_SurfaceDiagnosticSource.H"
#include "ERF_TerrainMetrics.H"
#include "ERF_PBLHeight.H"
#include "ERF_SurfaceModel.H"
#include "ERF_MicrophysicsUtils.H"
#include "ERF_EB.H"
 
namespace erf_surface_layer {
 
inline bool
planar_sources_supported_for_terrain (TerrainType terrain_type,
                                       bool use_sst,
                                       bool use_tsk,
                                       bool use_coupled_sst,
                                       bool has_lsm_tsurf,
                                       bool has_lsm_fluxes,
                                       bool has_custom_roughness)
{
    return terrain_type != TerrainType::EB ||
           !(use_sst || use_tsk || use_coupled_sst || has_lsm_tsurf ||
             has_lsm_fluxes || has_custom_roughness);
}
 
} // namespace erf_surface_layer
 
 
/** Abstraction layer for different surface layer schemes (e.g. MOST, Cd)
 *
 * van der Laan, P., Kelly, M. C., & Sørensen, N. N. (2017). A new k-epsilon
 * model consistent with Monin-Obukhov similarity theory. Wind Energy,
 * 20(3), 479–amrex::Real(489.) https://doi.org/amrex::Real(10.1002)/we.2017
 *
 * Consistent with Dyer (1974) formulation from page 57, Chapter 2, Modeling
 * the vertical ABL structure in Modelling of Atmospheric Flow Fields,
 * Demetri P Lalas and Corrado F Ratto, January 1996,
 * https://doi.org/amrex::Real(10.1142)/amrex::Real(2975.)
 */
class SurfaceLayer
{
 
public:
  /**
   * Construct the surface-layer interface.
   *
   * @param[in]     face                orientation of face (for wall geometries)
   * @param[in]     geom                geometry for all AMR levels
   * @param[in,out] use_rot_surface_flux whether to use rotated surface fluxes
   * @param[in]     a_pp_prefix         ParmParse prefix used by MOST averages
   * @param[in]     Qv_prim             primitive water-vapor fields by level
   * @param[in]     z_phys_nd           nodal physical-height fields by level
   * @param[in]     zlevels_stag        nominal staggered z levels by level
   * @param[in]     a_mesh_type         mesh type
   * @param[in]     a_terrain_type      terrain representation
   * @param[in]     a_turb_choice       turbulence-model options
   * @param[in]     a_rdOcp             configured Rd/cp exponent for T-theta conversions
   * @param[in]     start_low_time      first available low-boundary-data time
   * @param[in]     final_low_time      final available low-boundary-data time
   * @param[in]     low_time_interval   low-boundary-data time interval
   * @param[in]     eb_vec              optional embedded-boundary geometry data
   * @param[in]     surf_model          optional weighted surface-model interface
   * @param[in]     a_has_immersed_buildings  whether building immersed forcing is configured
   */
  explicit SurfaceLayer (amrex::Orientation face,
                         const amrex::Vector<amrex::Geometry>& geom,
                         bool& use_rot_surface_flux,
                         std::string a_pp_prefix,
                         amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
                         amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd,
                         const amrex::Vector<amrex::Vector<amrex::Real>>& zlevels_stag,
                         const MeshType& a_mesh_type,
                         const TerrainType& a_terrain_type,
                         const TurbChoice& a_turb_choice,
                         amrex::Real a_rdOcp,
                         double start_low_time,
                         double final_low_time,
                         double low_time_interval = 0.0,
                         const amrex::Vector<const eb_*>& eb_vec = {},
                         SurfaceModel* surf_model = nullptr,
                         bool a_has_immersed_buildings = false)
  : m_face(face),
    m_pp_prefix(a_pp_prefix),
    m_geom(geom),
    m_rotate(use_rot_surface_flux),
    m_start_low_time(start_low_time),
    m_final_low_time(final_low_time),
    m_low_time_interval(low_time_interval),
    m_eb_vec(eb_vec),
    m_terrain_type(a_terrain_type),
    m_has_immersed_buildings(a_has_immersed_buildings),
    m_rdOcp(a_rdOcp),
    m_ma(face, geom, (z_phys_nd[0] != nullptr), a_pp_prefix, a_mesh_type, a_terrain_type,
         zlevels_stag, eb_vec),
    m_surf_model(surf_model)
  {
      // We have a moisture model if Qv_prim is a valid pointer
      use_moisture = (Qv_prim[0].get());
 
      // Keep standalone SurfaceLayer instances compatible with the historical
      // single-face behavior. ERF replaces this with the complete active-face
      // set after constructing all SurfaceLayer objects.
      m_surface_layer_faces[static_cast<int>(face)] = 1;
 
      // Get roughness
      amrex::ParmParse pp(a_pp_prefix);
      m_z0_given = pp.contains("most.z0");
      pp.queryAdd("most.z0", z0_const);
 
      // Specify how to compute the flux
      if (use_rot_surface_flux) {
          flux_type = FluxCalcType::ROTATE;
      } else {
          std::string flux_string_in;
          std::string flux_string{"moeng"};
          auto read_flux = pp.queryAdd("surface_layer.flux_type", flux_string_in);
          if (read_flux) {
              flux_string = amrex::toLower(flux_string_in);
          }
          if (flux_string == "moeng") {
              flux_type = FluxCalcType::MOENG;
          } else if (flux_string == "rico") {
              flux_type = FluxCalcType::RICO;
          } else if (flux_string == "bulk_coeff") {
              flux_type = FluxCalcType::BULK_COEFF;
          } else if (flux_string == "custom") {
              flux_type = FluxCalcType::CUSTOM;
          } else {
              amrex::Abort("Undefined MOST flux type!");
          }
      }
 
      if (flux_type == FluxCalcType::BULK_COEFF ||
          flux_type == FluxCalcType::CUSTOM ||
          flux_type == FluxCalcType::RICO) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow(),
              "BULK_COEFF, CUSTOM, and RICO surface-layer fluxes are supported only on the z-low face.");
      }
 
      // Include w* to handle free convection (Beljaars 1995, QJRMS)
      pp.queryAdd("most.include_wstar", m_include_wstar);
 
      std::string pblh_string_in;
      std::string pblh_string{"none"};
      auto read_pblh = pp.queryAdd("most.pblh_calc", pblh_string_in);
      if (read_pblh) {
          pblh_string = amrex::toLower(pblh_string_in);
      }
      if (pblh_string == "none") {
          pblh_type = PBLHeightCalcType::None;
      } else if (pblh_string == "mynn25") {
          pblh_type = PBLHeightCalcType::MYNN25;
      } else if (pblh_string == "mynnedmf") {
          pblh_type = PBLHeightCalcType::MYNN25;
      } else if (pblh_string == "ysu") {
          pblh_type = PBLHeightCalcType::YSU;
      } else if (pblh_string == "mrf") {
          pblh_type = PBLHeightCalcType::MRF;
      } else {
          amrex::Abort("Undefined PBLH calc type!");
      }
 
      if (pblh_type != PBLHeightCalcType::None || m_include_wstar) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow(),
              "MOST PBL-height calculation and wstar correction are supported only on the z-low face.");
      }
 
      // The w* correction is computed from the PBL height, so it needs a scheme that
      // actually diagnoses one.  With pblh_calc = "none" the pblh MultiFab keeps the
      // bogus_large_value it was initialized with and calc_wstar turns that into a
      // convective velocity scale of ~1e50, which destroys the surface fluxes.
      if (m_include_wstar && (pblh_type == PBLHeightCalcType::None)) {
          amrex::Abort("erf.most.include_wstar requires a PBL height: set "
                       "erf.most.pblh_calc (MYNN25 is the only scheme implemented)");
      }
 
      if (m_surf_model) {
          use_surface_model = true;
          amrex::Print() << " Using fluxes from surface model!" << std::endl;
          surf_model_fluxes = m_surf_model->are_fluxes();
      }
 
      // Get surface temperature, as a potential temperature (most.pot_surf_temp) or an
      // absolute temperature (most.abs_surf_temp); most.surf_temp is the deprecated
      // spelling of most.pot_surf_temp.  These are read with query, not queryAdd, so
      // that whether a key was given is the query return value: queryAdd would add a
      // key that is absent, and these keys are also read by ERF_InputSoundingData.H and
      // by the SurfaceLayer of every face, so the first reader would make the key look
      // given to the rest.  A key that is given must be positive; it is never ignored.
      {
          amrex::Real pot_st = amrex::Real(-1.);
          amrex::Real abs_st = amrex::Real(-1.);
          amrex::Real old_st = amrex::Real(-1.);
          const bool has_pot = pp.query("most.pot_surf_temp", pot_st);
          const bool has_abs = pp.query("most.abs_surf_temp", abs_st);
          const bool has_old = pp.query("most.surf_temp",     old_st);
          if (int(has_pot) + int(has_abs) + int(has_old) > 1) {
              amrex::Abort("Specify only one of most.pot_surf_temp, most.abs_surf_temp "
                           "and the deprecated most.surf_temp");
          }
          if ((has_pot && !(pot_st > amrex::Real(0))) ||
              (has_abs && !(abs_st > amrex::Real(0))) ||
              (has_old && !(old_st > amrex::Real(0)))) {
              amrex::Abort("most.pot_surf_temp, most.abs_surf_temp and most.surf_temp "
                           "must be positive temperatures [K]");
          }
          if (has_old) {
              amrex::Print() << "WARNING: most.surf_temp is deprecated and is read as a surface "
                             << "potential temperature; use most.pot_surf_temp (potential) or "
                             << "most.abs_surf_temp (absolute) instead" << std::endl;
              pot_st = old_st;
          }
          if (has_abs) {
              surf_temp = abs_st;
              surf_temp_is_absolute = true;
          } else {
              surf_temp = pot_st;
          }
      }
      const bool erf_st = (surf_temp > amrex::Real(0));
      // With an absolute temperature this is only the placeholder until update_fluxes
      // converts it with the surface pressure.
      if (erf_st) { default_land_surf_temp  = surf_temp; }
 
      if (surf_temp_is_absolute) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow() && m_terrain_type != TerrainType::EB,
              "most.abs_surf_temp is supported only on the z-low face without EB terrain");
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              flux_type != FluxCalcType::CUSTOM && flux_type != FluxCalcType::RICO &&
              flux_type != FluxCalcType::BULK_COEFF,
              "most.abs_surf_temp is not supported with custom, rico or bulk_coeff fluxes; "
              "use most.pot_surf_temp");
      }
 
      // Get surface moisture
      bool erf_sq = false;
      if (use_moisture) {
          pp.queryAdd("most.surf_moist", surf_moist);
          erf_sq = (surf_moist >= amrex::Real(0));
      }
      if (erf_sq) { default_land_surf_moist = surf_moist; }
 
      // Custom type user must specify the fluxes
      if (flux_type == FluxCalcType::CUSTOM ||
          flux_type == FluxCalcType::RICO) {
          theta_type = ThetaCalcType::HEAT_FLUX;
          moist_type = MoistCalcType::MOISTURE_FLUX;
          pp.get("most.ustar", custom_ustar);
          pp.get("most.tstar", custom_tstar);
          pp.get("most.qstar", custom_qstar);
          pp.queryAdd("most.rhosurf", custom_rhosurf);
          if (custom_qstar != 0) {
              AMREX_ASSERT_WITH_MESSAGE(use_moisture,
              "Specified custom MOST qv flux without moisture model!");
          }
          amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
                         << custom_ustar << " " << custom_tstar << " "
                         << custom_qstar << std::endl;
 
      // Bulk transfer coefficient (must specify coeffs and surface values)
      } else if (flux_type == FluxCalcType::BULK_COEFF) {
          pp.get("most.Cd", m_Cd);
          pp.get("most.Ch", m_Ch);
          pp.get("most.Cq", m_Cq);
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(erf_st,
              "bulk_coeff fluxes need most.pot_surf_temp");
          default_land_surf_temp = surf_temp;
          pp.get("most.surf_moist", default_land_surf_moist);
          amrex::Print() << "Using specified Cd, Ch, Cq for MOST = "
                         << m_Cd << " " << m_Ch << " "
                         << m_Cq << std::endl;
 
      // Specify surface temperature/moisture or surface flux
      } else {
          if (erf_st) {
              theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
              pp.queryAdd("most.surf_heating_rate", surf_heating_rate); // [K/h]
 
              // Modify rate to be in units of K / s rather than K / hr
              surf_heating_rate /= amrex::Real(3600.0);        // [K/s]
 
              if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
                  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
              }
          } else {
              pp.queryAdd("most.surf_temp_flux", surf_temp_flux);
 
              if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
                  amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
              }
              if (std::abs(surf_temp_flux) >
                  std::numeric_limits<amrex::Real>::epsilon()) {
                  theta_type = ThetaCalcType::HEAT_FLUX;
              } else {
                  theta_type = ThetaCalcType::ADIABATIC;
              }
          }
 
          if (erf_sq) {
              moist_type = MoistCalcType::SURFACE_MOISTURE;
          } else {
              pp.queryAdd("most.surf_moist_flux", surf_moist_flux);
              if (std::abs(surf_moist_flux) >
                  std::numeric_limits<amrex::Real>::epsilon()) {
                  moist_type = MoistCalcType::MOISTURE_FLUX;
              } else {
                  moist_type = MoistCalcType::ADIABATIC;
              }
          }
      }
 
      if (theta_type == ThetaCalcType::HEAT_FLUX ||
          (theta_type == ThetaCalcType::ADIABATIC &&
           moist_type == MoistCalcType::ADIABATIC)) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow(),
              "HEAT_FLUX and ADIABATIC surface-layer fluxes are supported only on the z-low face.");
      }
 
      if (flux_type == FluxCalcType::RICO)
      {
          pp.queryAdd("most.rico.theta_z0", rico_theta_z0);
          pp.queryAdd("most.rico.qsat_z0", rico_qsat_z0);
      }
 
      // Make sure the inputs file doesn't try to use most.roughness_type
      std::string bogus_input;
      if (pp.queryAdd("most.roughness_type", bogus_input) > 0) {
          amrex::Abort("most.roughness_type is deprecated; use "
                       "most.roughness_type_land and/or most.roughness_type_sea");
      }
 
      // Specify how to compute the surface flux over land (if there is any)
      std::string rough_land_string_in;
      std::string rough_land_string{"constant"};
      auto read_rough_land =
          pp.queryAdd("most.roughness_type_land", rough_land_string_in);
      if (read_rough_land) {
          rough_land_string = amrex::toLower(rough_land_string_in);
      }
      if (rough_land_string == "constant") {
          rough_type_land = RoughCalcType::CONSTANT;
      } else {
          amrex::Abort("Undefined MOST roughness type for land!");
      }
 
      // Allow for smooth-flow limit?
      pp.queryAdd("most.smooth_flow_viscosity", smooth_flow_visc);
      amrex::Print() << "The smooth-flow limit will be included for variable roughness models over sea: " << smooth_flow_visc << "\n";
 
      // Specify how to compute the surface flux over sea (if there is any)
      std::string rough_sea_string_in;
      std::string rough_sea_string{"charnock"};
      auto read_rough_sea = pp.queryAdd("most.roughness_type_sea", rough_sea_string_in);
      if (read_rough_sea) {
          rough_sea_string = amrex::toLower(rough_sea_string_in);
      }
      if (rough_sea_string == "charnock") {
          rough_type_sea = RoughCalcType::CHARNOCK;
          pp.queryAdd("most.charnock_constant", cnk_a);
          if (cnk_a > 0) {
              amrex::Print() << "If there is water, Charnock relation with C_a="
                             << cnk_a << " will be used" << std::endl;
          } else {
              amrex::Print() << "If there is water, Charnock relation with variable "
                  "Charnock parameter (COARE3.0) will be used" << std::endl;
          }
      } else if (rough_sea_string == "coare3.0") {
          rough_type_sea = RoughCalcType::CHARNOCK;
          amrex::Print() << "If there is water, Charnock relation with variable "
              "Charnock parameter (COARE3.0) will be used" << std::endl;
          cnk_a = -1;
      } else if (rough_sea_string == "donelan") {
          rough_type_sea = RoughCalcType::DONELAN;
      } else if (rough_sea_string == "modified_charnock") {
          rough_type_sea = RoughCalcType::MODIFIED_CHARNOCK;
          pp.queryAdd("most.modified_charnock_depth", depth);
          if (depth < amrex::Real(10.0) || depth > amrex::Real(100.0) ) {
              amrex::Print() << "Specified depth of " << depth
                             << " is outside the valid range of [10 100], resetting to bounds now."
                             << std::endl;
          }
          // Limiter based upon the fit range in Jiménez & Dudhia
          depth = amrex::min(amrex::max(depth,amrex::Real(10.0)),amrex::Real(100.0));
      } else if (rough_sea_string == "wave_coupled") {
          rough_type_sea = RoughCalcType::WAVE_COUPLED;
      } else if (rough_sea_string == "constant") {
          rough_type_sea = RoughCalcType::CONSTANT;
      } else {
          amrex::Abort("Undefined MOST roughness type for sea!");
      }
 
      // use skin temperature instead of sea-surface temperature
      // (wrfinput data may have lower resolution SST data)
      pp.queryAdd("most.ignore_sst", m_ignore_sst);
 
      // If we're using the RANS k model, then we need to update the dirichlet
      // wall value of k (written into the first cell above the wall and held
      // there through the step) based on the instantaneous u* and θ*; the turbulence modeling
      // choices can vary per level but for now, assume that if specified then
      // all levels are using the same RANS model.
      m_update_k_rans = (a_turb_choice.rans_type == RANSType::kEqn &&
                         a_turb_choice.dirichlet_k == true);
      if (m_update_k_rans) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow(),
              "RANS surface-layer k updates are supported only on the z-low face.");
      }
      if (m_update_k_rans) {
          inv_Cmu2 = one / (a_turb_choice.Cmu0 * a_turb_choice.Cmu0);
          theta_ref = a_turb_choice.theta_ref;
      }
 
  } // constructor
 
  /**
   * Allocate and initialize surface-layer data for one AMR level.
   *
   * @param[in]     lev           level index
   * @param[in]     nlevs         number of AMR levels
   * @param[in]     mfv           conserved and velocity MultiFabs for this level
   * @param[in]     Theta_prim    primitive potential-temperature field
   * @param[in]     Qv_prim       primitive water-vapor field
   * @param[in]     Qr_prim       primitive rain-water field
   * @param[in]     z_phys_nd     nodal physical-height field
   * @param[in]     Hwave         wave-height field
   * @param[in]     Lwave         wavelength field
   * @param[in]     eddyDiffs     eddy-diffusivity field
   * @param[in]     lsm_data      land-surface-model data fields
   * @param[in]     lsm_data_name names for lsm_data entries
   * @param[in]     lsm_flux      land-surface-model flux fields
   * @param[in]     lsm_flux_name names for lsm_flux entries
   * @param[in]     sst_lev       sea-surface-temperature data by time
   * @param[in]     tsk_lev       skin-temperature data by time
   * @param[in]     lmask_lev     land-mask data by time
   */
  void make_SurfaceLayer_at_level (const int& lev,
                                   int nlevs,
                                   const amrex::Vector<amrex::MultiFab*>& mfv,
                                   std::unique_ptr<amrex::MultiFab>& Theta_prim,
                                   std::unique_ptr<amrex::MultiFab>& Qv_prim,
                                   std::unique_ptr<amrex::MultiFab>& Qr_prim,
                                   std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                                   amrex::MultiFab* Hwave,
                                   amrex::MultiFab* Lwave,
                                   amrex::MultiFab* eddyDiffs,
                                   amrex::Vector<amrex::MultiFab*> lsm_data,
                                   amrex::Vector<std::string>      lsm_data_name,
                                   amrex::Vector<amrex::MultiFab*> lsm_flux,
                                   amrex::Vector<std::string>      lsm_flux_name,
                                   amrex::Vector<std::unique_ptr<amrex::MultiFab>>& sst_lev,
                                   amrex::Vector<std::unique_ptr<amrex::MultiFab>>& tsk_lev,
                                   amrex::Vector<std::unique_ptr<amrex::iMultiFab>>& lmask_lev)
  {
      // Update MOST Average
      m_ma.make_MOSTAverage_at_level(lev, mfv,
                                     Theta_prim, Qv_prim, Qr_prim,
                                     z_phys_nd);
 
      // Get CC vars
      amrex::MultiFab& mf = *(mfv[0]);
 
      amrex::ParmParse pp(m_pp_prefix);
 
      // Do we have a time-varying surface roughness that needs to be saved?
      if (lev == 0) {
          const int nghost = 0; // ghost cells not included
          int lmask_min = lmask_min_reduce(*lmask_lev[0].get(), nghost);
          amrex::ParallelDescriptor::ReduceIntMin(lmask_min);
 
          m_var_z0 = (lmask_min < 1) & (rough_type_sea != RoughCalcType::CONSTANT);
          if (m_var_z0) {
              AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                  m_face.coordDir() == 2 && m_face.isLow(),
                  "Variable sea roughness surface-layer fluxes are supported only on the z-low face.");
 
              std::string rough_sea_string{"charnock"};
              pp.queryAdd("most.roughness_type_sea", rough_sea_string);
              amrex::Print() << "Variable sea roughness (type " << rough_sea_string
                             << ")" << std::endl;
          }
      }
 
      if (m_eddyDiffs_lev.size() < lev+1) {
          m_Hwave_lev.resize(nlevs);
          m_Lwave_lev.resize(nlevs);
          m_eddyDiffs_lev.resize(nlevs);
 
          m_lsm_data_lev.resize(nlevs);
          m_lsm_flux_lev.resize(nlevs);
 
          m_sst_lev.resize(nlevs);
          m_tsk_lev.resize(nlevs);
          m_lmask_lev.resize(nlevs);
 
          m_coupled_sst_lev.resize(nlevs, nullptr);
          m_coupled_sst_valid_lev.resize(nlevs, nullptr);
 
          // Size the MOST params for all levels
          z_0.resize(nlevs);
          u_star.resize(nlevs);
          w_star.resize(nlevs);
          t_star.resize(nlevs);
          q_star.resize(nlevs);
          t_surf.resize(nlevs);
          q_surf.resize(nlevs);
          surface_diagnostic_source.resize(nlevs);
          olen.resize(nlevs);
          pblh.resize(nlevs);
          m_planar_bndry.resize(nlevs);
      }
 
      // Get pointers to SST,TSK and LANDMASK data
      int nt_tot_sst = sst_lev.size();
      m_sst_lev[lev].resize(nt_tot_sst);
      for (int nt(0); nt < nt_tot_sst; ++nt) {
          m_sst_lev[lev][nt] = sst_lev[nt].get();
      }
      int nt_tot_tsk = static_cast<int>(tsk_lev.size());
      m_tsk_lev[lev].resize(nt_tot_tsk);
      for (int nt(0); nt < nt_tot_tsk; ++nt) {
          m_tsk_lev[lev][nt] = tsk_lev[nt].get();
      }
      int nt_tot_lmask = static_cast<int>(lmask_lev.size());
      m_lmask_lev[lev].resize(nt_tot_lmask);
      for (int nt(0); nt < nt_tot_lmask; ++nt) {
          m_lmask_lev[lev][nt] = lmask_lev[nt].get();
      }
 
      // Get pointers to wave data
      m_Hwave_lev[lev] = Hwave;
      m_Lwave_lev[lev] = Lwave;
      m_eddyDiffs_lev[lev] = eddyDiffs;
 
      // Text-file driven surface forcing modes. The file always contains at
      // least time(day) and absolute sst(K); the value is normalized to the
      // canonical SurfaceLayer theta field at the z-low surface.
      pp.queryAdd("most.use_sfc_fluxes", m_use_sfc_fluxes);
      pp.queryAdd("most.use_sfc_sst", m_use_sfc_sst);
      if (m_use_sfc_fluxes && m_use_sfc_sst) {
          amrex::Abort("Only one of most.use_sfc_fluxes and most.use_sfc_sst may be enabled");
      }
      if (m_use_sfc_fluxes || m_use_sfc_sst) {
          if (m_use_sfc_fluxes) {
              AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!use_surface_model, "Cannot use surface flux forcing with LSM/Urban model!");
              // The file's H and LE already carry the density, and the stress
              // takes the local density, so a surface density would count it
              // twice in H and LE (see set_sfc_params_from_file_fluxes).
              // The value is tested, not the key: flux_type = custom has already
              // queryAdd-ed most.rhosurf with its default of 0 (use the first cell).
              amrex::Real rhosurf_in = zero;
              pp.query("most.rhosurf", rhosurf_in);
              AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rhosurf_in <= zero,
                  "most.rhosurf cannot be combined with most.use_sfc_fluxes");
          }
          if (m_terrain_type == TerrainType::EB) {
              amrex::Abort("Text-file surface forcing is not supported with EB terrain");
          }
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow(),
              "most.use_sfc_fluxes and most.use_sfc_sst are supported only on the z-low face.");
 
          // load sensible and latent heat fluxes from sfc to prescribe
          std::string sfc_file = "";
          pp.queryAdd("most.sfc_file", sfc_file);
          if (sfc_file.empty()) {
              amrex::Abort("most.sfc_file must be set when using text-file surface forcing");
          }
 
          // sfc contains: time(day) sst(K) H(W/m2) LE(W/m2) TAU(m2/s2)
          sfc = read_cols(sfc_file, 1);
 
          const int min_cols = m_use_sfc_fluxes ? 5 : 2;
          if (static_cast<int>(sfc.size()) < min_cols) {
              amrex::Abort("Surface forcing file does not contain the required number of columns");
          }
 
          // shift time column in days to be relative to current elapsed time
          const amrex::Real start_day = sfc[0][0];
          for (int i = 0; i < static_cast<int>(sfc[0].size()); ++i) {
              sfc[0][i] = 86400.0 * (sfc[0][i] - start_day);
          }
 
          if (m_use_sfc_sst) {
              if (theta_type == ThetaCalcType::HEAT_FLUX || surf_heating_rate != zero) {
                  amrex::Abort("most.use_sfc_sst cannot be combined with prescribed heat flux or surf_heating_rate");
              }
              if (moist_type == MoistCalcType::MOISTURE_FLUX) {
                  amrex::Abort("most.use_sfc_sst cannot be combined with prescribed moisture flux");
              }
              theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
              amrex::Print() << "Using MOST with prescribed SST from most.sfc_file '" << sfc_file << "' over sea" << std::endl;
          }
 
          if (m_use_sfc_fluxes) {
              flux_type = FluxCalcType::CUSTOM;
              amrex::Print() << "Using MOST with prescribed time-varying surface fluxes from '" << sfc_file << "'" << std::endl;
          }
      }
 
      // The constructor's check of most.abs_surf_temp sees only most.flux_type; the
      // modes known here replace the prescribed temperature too.  The surface-flux file
      // makes the fluxes custom, and the SST/TSK and sfc-file SST fills use the
      // prescribed value as their land temperature, which is a potential temperature there.
      if (surf_temp_is_absolute) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_use_sfc_fluxes,
              "most.abs_surf_temp cannot be combined with most.use_sfc_fluxes");
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              (m_sst_lev[lev].empty() || !m_sst_lev[lev][0]) && !m_use_sfc_sst,
              "most.abs_surf_temp cannot be combined with SST/TSK or sfc-file SST input");
      }
 
      // By default an LSM field named "tsurf", "theta" or "t_surf" is copied into
      // t_surf unchanged. With this option the LSM skin temperature ("tsurf" from
      // SLM, "t_sfc" from Noah-MP) is taken as an absolute temperature and
      // converted to potential temperature at the surface pressure.
      pp.queryAdd("most.use_lsm_skin_temperature", m_use_lsm_skin_temperature);
 
      // Get pointers to LSM data and Fluxes
      int ndata = static_cast<int>(lsm_data.size());
      int nflux = static_cast<int>(lsm_flux.size());
      m_lsm_data_name.resize(ndata);
      m_lsm_data_lev[lev].resize(ndata);
      m_lsm_flux_name.resize(nflux);
      m_lsm_flux_lev[lev].resize(nflux);
      for (int n(0); n < ndata; ++n) {
          m_lsm_data_name[n]     = lsm_data_name[n];
          m_lsm_data_lev[lev][n] = lsm_data[n];
          const std::string lc_name = amrex::toLower(lsm_data_name[n]);
          // With most.use_lsm_skin_temperature, "tsurf" (SLM) and "t_sfc" (Noah-MP)
          // are absolute skin temperatures and are converted to potential
          // temperature in get_lsm_tsurf; "theta" and "t_surf" are taken to be
          // potential temperature already.
          if (lc_name == "tsurf" || (m_use_lsm_skin_temperature && lc_name == "t_sfc") ||
              lc_name == "theta" || lc_name == "t_surf") {
              m_has_lsm_tsurf  = true;
              m_lsm_tsurf_indx = n;
              m_lsm_tsurf_is_theta = (lc_name == "theta" || lc_name == "t_surf");
          }
      }
      if (m_use_sfc_sst && m_use_coupled_sst) {
          amrex::Abort("most.use_sfc_sst cannot be combined with an ocean LSM t_surf input");
      }
      int n_valid_lsm_flux = 0;
      bool has_soil_t_flux = false;
      for (int n(0); n < nflux; ++n) {
          m_lsm_flux_name[n]     = lsm_flux_name[n];
          m_lsm_flux_lev[lev][n] = lsm_flux[n];
          if (m_lsm_flux_lev[lev][n]) { ++n_valid_lsm_flux; }
          if (amrex::toLower(m_lsm_flux_name[n]) == "soil_t_flux") {
              has_soil_t_flux = true;
          }
      }
      AMREX_ALWAYS_ASSERT((n_valid_lsm_flux==0 || n_valid_lsm_flux>=4 ||
                           (n_valid_lsm_flux==1 && has_soil_t_flux)));
      if (n_valid_lsm_flux>=4) { m_has_lsm_fluxes = true; }
 
      const bool use_sst = (!m_sst_lev[lev].empty() && m_sst_lev[lev][0]);
      const bool use_tsk = (!m_tsk_lev[lev].empty() && m_tsk_lev[lev][0]);
 
      // Check if there is a user-specified roughness file to be read
      std::string fname;
      bool read_z0 = false;
      if ( (flux_type == FluxCalcType::MOENG) ||
           (flux_type == FluxCalcType::ROTATE)) {
          int count = pp.countval("most.roughness_file_name");
          if (count > 1) {
              AMREX_ALWAYS_ASSERT(count >= lev+1);
              pp.query("most.roughness_file_name", fname, lev);
              read_z0 = true;
          } else if (count == 1) {
              if (lev == 0) {
                  pp.queryAdd("most.roughness_file_name", fname);
              } else {
                  // we will interpolate from the coarsest level
                  fname = "";
              }
              read_z0 = true;
          }
          // else use z0_const
      }
      if (read_z0) {
          AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
              m_face.coordDir() == 2 && m_face.isLow(),
              "Custom MOST roughness is supported only on the z-low face.");
      }
 
      // LSM flux arrays are planar and compute_sfc_params_from_lsm_fluxes only
      // writes the k=0 slab, while EB MOST consumes surface parameters at
      // arbitrary cut-cell k; no valid planar-to-cut-cell mapping exists yet.
      if (!erf_surface_layer::planar_sources_supported_for_terrain(
              m_terrain_type, use_sst, use_tsk, m_use_coupled_sst,
              m_has_lsm_tsurf, m_has_lsm_fluxes, read_z0)) {
          amrex::Abort(
              "EB SurfaceLayer does not support planar SST/TSK, coupled SST, LSM surface "
              "temperature or fluxes, or custom/file-driven roughness; no mapping exists "
              "from those planar inputs to arbitrary EB cut cells.");
      }
 
      // Attributes for MFs and FABs
      //--------------------------------------------------------
      // Create a 2D ba for planar terrain, 3D for EB terrain
      const int dir = m_face.coordDir();
      int sm_index;
      if (m_face.isLow()) {
          sm_index = m_geom[lev].Domain().smallEnd(dir);
      } else {
          sm_index = m_geom[lev].Domain().bigEnd(dir);
      }
 
      amrex::BoxArray ba = mf.boxArray();
      amrex::BoxArray ba_flux;
      amrex::IntVect ng{1,1,0};
 
      // The lateral-wall implementation uses the two-dimensional land-mask
      // layout to identify grids on the requested face.  Lateral walls are
      // intentionally restricted to grids spanning the complete z domain;
      // z-decomposed and partial-height grids are not supported.
      if (dir != 2) {
          const int dom_lo_z = m_geom[lev].Domain().smallEnd(2);
          const int dom_hi_z = m_geom[lev].Domain().bigEnd(2);
          for (int ibox = 0; ibox < ba.size(); ++ibox) {
              AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                  ba[ibox].smallEnd(2) == dom_lo_z && ba[ibox].bigEnd(2) == dom_hi_z,
                  "Surface layer boundaries on x/y faces require Cartesian grids that "
                  "span the full level z domain; partial-height refined grids and grids "
                  "decomposed in z are not supported. Set erf.max_grid_size_z accordingly.");
          }
      }
 
      if (m_terrain_type == TerrainType::EB) {
          // Use full 3D BoxArray for EB terrain
          ba_flux = ba;
          ng = amrex::IntVect{1,1,1};  // Include z ghost cells
      } else {
          // Collapse to 2D for planar terrain
          amrex::BoxList bl2d = ba.boxList();
          for (auto& b : bl2d) {
            b.setRange(dir,sm_index);
          }
          ba_flux = amrex::BoxArray(std::move(bl2d));
          // Lateral faces need state-width ghosts in the tangential z
          // direction.  A z face only needs the original one-cell x/y halo;
          // kernels using the wider state/mask boxes clip to the FAB bounds.
          if (dir == 2) {
              ng = amrex::IntVect{1,1,0};
          } else {
              ng = mf.nGrowVect();
              ng[dir] = 0;
          }
      }
 
      const amrex::DistributionMapping& dm = mf.DistributionMap();
      const int ncomp = 1;
 
      // Surface copies of the planar boxes (see fill_planar_boundary)
      // PlanarBoundary handles duplicate boxes created by a z-split. The x/y
      // layouts use selective lateral exchange instead.
      if (m_terrain_type != TerrainType::EB && dir == 2) {
          const int ksurface = m_face.isLow()
              ? m_geom[lev].Domain().smallEnd(2)
              : m_geom[lev].Domain().bigEnd(2);
          m_planar_bndry[lev].define(ba, ba_flux, dm, ksurface, m_face.isLow());
      }
 
      // Z0 heights FAB
      //--------------------------------------------------------
      z_0[lev].define(ba_flux, dm, ncomp, ng);
      z_0[lev].setVal(z0_const);
      if (read_z0) {
          read_custom_roughness(lev, fname);
      }
 
      // 2D MFs for U*, T*, T_surf
      //--------------------------------------------------------
      u_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      u_star[lev]->setVal(bogus_large_value);
 
      w_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      w_star[lev]->setVal(bogus_large_value);
 
      t_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      t_star[lev]->setVal(zero); // default to neutral
 
      q_star[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      q_star[lev]->setVal(zero); // default to dry
 
      olen[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      olen[lev]->setVal(bogus_large_value);
 
      pblh[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      pblh[lev]->setVal(bogus_large_value);
 
      t_surf[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      t_surf[lev]->setVal(default_land_surf_temp);
 
      q_surf[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      q_surf[lev]->setVal(default_land_surf_moist);
 
      const amrex::iMultiFab& surface_mask = *m_lmask_lev[lev][0];
      AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
          surface_mask.boxArray().size() == mf.boxArray().size(),
          "Surface-layer mask and state must have the same number of boxes.");
      AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
          surface_mask.DistributionMap() == mf.DistributionMap(),
          "Surface-layer mask and state must have identical ownership.");
      AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
          surface_mask.boxArray().size() == u_star[lev]->boxArray().size(),
          "Surface-layer mask and parameters must have the same number of boxes.");
      AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
          surface_mask.DistributionMap() == u_star[lev]->DistributionMap(),
          "Surface-layer mask and parameters must have identical ownership.");
 
      set_immersed_surface_params(lev);
 
      surface_diagnostic_source[lev] = std::make_unique<amrex::MultiFab>(ba_flux, dm, ncomp, ng);
      surface_diagnostic_source[lev]->setVal(
          surface_diagnostics::to_plot_value(surface_diagnostics::SurfaceDiagnosticSource::Missing));
 
      // TODO: Do we want an enum struct for indexing?
 
      if (use_sst || use_tsk || m_has_lsm_tsurf || m_use_coupled_sst || use_surface_model) {
          // Valid SST, TSK, LSM or coupled-ocean data; t_surf set before computing
          // fluxes (avoids extended lambda capture) Note that land temp will be set
          // from m_tsk_lev while sea temp will be set from m_sst_lev
          theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
 
          // Pathways in fill_tsurf_with_sst_and_tsk
          amrex::Print() << "Using MOST with specified surface temperature ";
          if (m_has_lsm_tsurf && !use_sst && !use_tsk) {
              amrex::Print() << "(LSM: " << m_lsm_data_name[m_lsm_tsurf_indx] << ")";
          } else if (!use_sst && !use_tsk) {
              amrex::Print() << "(land: T0, sea: none)";
          } else {
              // NOTE: SST from the LOW file populates TSK in update_sst_tsk.
              //       So if we have TSK, it contains everything and has been
              //       sanity checked for valid SST values.
              if (use_tsk) { m_ignore_sst = true; }
              if (use_tsk) {
                  amrex::Print() << "(land: TSK, ";
              } else {
                  amrex::Print() << "(land: T0, ";
              }
              if (use_tsk && !use_sst) {
                  amrex::Print() << "sea: TSK)";
              } else {
                  amrex::Print() << "sea: SST)";
                  AMREX_ALWAYS_ASSERT(m_sst_lev[lev][0]);
              }
          }
          // The coupler is layered on top of whatever the above selected: it
          // overwrites only the water cells it actually covers, so the pathway
          // named above remains the value for land and for uncovered water.
          if (m_use_coupled_sst) {
              amrex::Print() << " + coupled ocean SST where covered";
          }
          amrex::Print() << std::endl;
      }
  }
 
  /**
   * Set benign surface scales in immersed columns (land mask 2), including ghosts.
   * Called after allocation, checkpoint restoration, and surface-scale updates.
   * Preserve PBL height, roughness, surface temperature, and surface moisture.
   *
   * @param[in] lev level index
   */
  void set_immersed_surface_params (const int& lev);
 
  /**
   * Update surface fluxes and related surface-layer state.
   *
   * @param[in]     lev                      level index
   * @param[in]     elapsed_time             current elapsed simulation time
   * @param[in]     elapsed_time_since_start_low elapsed time relative to low-data start
   * @param[in,out] cons_in                  conserved state used by the flux update
   * @param[in]     z_phys_nd                nodal physical-height field
   * @param[in]     walldist                 wall-distance field
   * @param[in]     max_iters                maximum MOST iteration count
   */
  void
  update_fluxes (const int& lev,
                 const double& elapsed_time,
                 const double& elapsed_time_since_start_low,
                 amrex::MultiFab& cons_in,
                 const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                 const std::unique_ptr<amrex::MultiFab>& walldist,
                 int max_iters = 100);
 
  /**
   * Compute MOST fluxes with a selected flux-iteration functor.
   *
   * @param[in]     lev       level index
   * @param[in]     max_iters maximum MOST iteration count
   * @param[in,out] cons_in   conserved state used by the flux computation
   * @param[in]     most_flux flux-iteration functor
   * @param[in]     is_land   whether the land-surface branch is active
   */
  template <typename FluxIter>
  void compute_fluxes (const int& lev,
                       const int& max_iters,
                       amrex::MultiFab& cons_in,
                       const FluxIter& most_flux,
                       bool is_land);
 
  /**
   * Initialize TKE from the current surface friction velocity.
   *
   * @param[in]     lev       level index
   * @param[in,out] cons      conserved state whose TKE component is initialized
   * @param[in]     z_phys_nd nodal physical-height field
   * @param[in]     tkefac    scale factor applied to the initialized TKE
   * @param[in]     zscale    vertical decay scale
   */
  void init_tke_from_ustar (const int& lev,
                            amrex::MultiFab& cons,
                            const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                            const amrex::Real tkefac = one,
                            const amrex::Real zscale = amrex::Real(700.0));
 
  /**
   * Fill the ghost cells of a planar surface-layer MultiFab, and the valid region of
   * its uncomputed copies when the 3D BoxArray is split in z (see PlanarBoundary).
   *
   * @param[in]     lev level index
   * @param[in,out] mf  planar MultiFab to fill
   */
  void fill_planar_boundary (const int& lev, amrex::MultiFab& mf);
 
  /**
   * Sum a planar field over the valid cells of the surface, counting each surface cell once.
   *
   * On a level whose grids are split in the surface-normal direction the planar BoxArray holds
   * one duplicate box per stacked 3D box (see PlanarBoundary), so a plain sum over the planar
   * MultiFab counts every surface cell once per stacked box.  Only the computed surface copies
   * are read here, so the duplicates need not have been filled by fill_planar_boundary, and the
   * gather needs no communication.  On EB terrain the fields are 3D without duplicates, and the
   * lowest plane of the domain is summed.  Only for a surface layer on a z face.
   *
   * @param[in] lev  level index
   * @param[in] mf   cell-centered planar MultiFab on this surface's planar BoxArray
   * @param[in] comp component to sum
   */
  amrex::Real surface_sum (const int& lev, const amrex::MultiFab& mf, int comp = 0) const;
 
  /**
   * Impose surface-layer boundary conditions for planar terrain.
   *
   * @param[in]     lev        level index
   * @param[in]     mfs        state and velocity fields used by the BC computation
   * @param[in,out] Tau_lev    stress fields to fill
   * @param[in,out] xheat_flux x-face heat flux field
   * @param[in,out] yheat_flux y-face heat flux field
   * @param[in,out] zheat_flux z-face heat flux field
   * @param[in,out] xqv_flux   x-face moisture flux field
   * @param[in,out] yqv_flux   y-face moisture flux field
   * @param[in,out] zqv_flux   z-face moisture flux field
   * @param[in]     z_phys     physical-height field
   */
  void impose_SurfaceLayer_bcs (const int& lev,
                                amrex::Vector<const amrex::MultiFab*> mfs,
                                amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev,
                                amrex::MultiFab* xheat_flux,
                                amrex::MultiFab* yheat_flux,
                                amrex::MultiFab* zheat_flux,
                                amrex::MultiFab* xqv_flux,
                                amrex::MultiFab* yqv_flux,
                                amrex::MultiFab* zqv_flux,
                                const amrex::MultiFab* z_phys,
                                amrex::MultiFab* mf_surface_fluxes=nullptr);
 
  /**
   * Impose surface-layer boundary conditions for embedded-boundary terrain.
   *
   * @param[in]     lev        level index
   * @param[in]     mfs        state and velocity fields used by the BC computation
   * @param[in,out] Tau_lev    EB stress fields to fill
   * @param[in,out] xheat_flux x-face heat flux field
   * @param[in,out] yheat_flux y-face heat flux field
   * @param[in,out] zheat_flux z-face heat flux field
   * @param[in,out] xqv_flux   x-face moisture flux field
   * @param[in,out] yqv_flux   y-face moisture flux field
   * @param[in,out] zqv_flux   z-face moisture flux field
   */
  void impose_SurfaceLayer_bcs_EB (const int& lev,
                                amrex::Vector<const amrex::MultiFab*> mfs,
                                amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>& Tau_lev,
                                amrex::MultiFab* xheat_flux,
                                amrex::MultiFab* yheat_flux,
                                amrex::MultiFab* zheat_flux,
                                amrex::MultiFab* xqv_flux,
                                amrex::MultiFab* yqv_flux,
                                amrex::MultiFab* zqv_flux);
 
  /**
   * Compute planar-terrain surface-layer flux boundary conditions.
   *
   * @param[in]     lev        level index
   * @param[in]     mfs        state and velocity fields used by the BC computation
   * @param[in,out] Tau_lev    stress fields to fill
   * @param[in,out] xheat_flux x-face heat flux field
   * @param[in,out] yheat_flux y-face heat flux field
   * @param[in,out] zheat_flux z-face heat flux field
   * @param[in,out] xqv_flux   x-face moisture flux field
   * @param[in,out] yqv_flux   y-face moisture flux field
   * @param[in,out] zqv_flux   z-face moisture flux field
   * @param[in]     z_phys     physical-height field
   * @param[in]     flux_comp  flux-computation functor
   */
  template <typename FluxCalc>
  void compute_SurfaceLayer_bcs (const int& lev,
                                 amrex::Vector<const amrex::MultiFab*> mfs,
                                 amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev,
                                 amrex::MultiFab* xheat_flux,
                                 amrex::MultiFab* yheat_flux,
                                 amrex::MultiFab* zheat_flux,
                                 amrex::MultiFab* xqv_flux,
                                 amrex::MultiFab* yqv_flux,
                                 amrex::MultiFab* zqv_flux,
                                 const amrex::MultiFab* z_phys,
                                 const FluxCalc& flux_comp,
                                 amrex::MultiFab* mf_surface_fluxes=nullptr);
 
  /**
   * Compute embedded-boundary surface-layer flux boundary conditions.
   *
   * @param[in]     lev        level index
   * @param[in]     mfs        state and velocity fields used by the BC computation
   * @param[in,out] Tau_lev    EB stress fields to fill
   * @param[in,out] xheat_flux x-face heat flux field
   * @param[in,out] yheat_flux y-face heat flux field
   * @param[in,out] zheat_flux z-face heat flux field
   * @param[in,out] xqv_flux   x-face moisture flux field
   * @param[in,out] yqv_flux   y-face moisture flux field
   * @param[in,out] zqv_flux   z-face moisture flux field
   * @param[in]     flux_comp  flux-computation functor
   */
  template <typename FluxCalc>
  void compute_SurfaceLayer_bcs_EB (const int& lev,
                                 amrex::Vector<const amrex::MultiFab*> mfs,
                                 amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>& Tau_lev,
                                 amrex::MultiFab* xheat_flux,
                                 amrex::MultiFab* yheat_flux,
                                 amrex::MultiFab* zheat_flux,
                                 amrex::MultiFab* xqv_flux,
                                 amrex::MultiFab* yqv_flux,
                                 amrex::MultiFab* zqv_flux,
                                 const FluxCalc& flux_comp);
 
  /**
   * Derive MOST surface parameters from LSM fluxes.
   *
   * @param[in]     lev     level index
   * @param[in,out] cons_in conserved state used by the surface-parameter computation
   */
  void compute_sfc_params_from_lsm_fluxes (const int& lev,
                                           amrex::MultiFab& cons_in);
 
  /**
   * Set kinematic u*, t*, and q* from the fluxes interpolated from most.sfc_file.
   *
   * @param[in] lev       level index
   * @param[in] cons_in   conserved state supplying the first-cell density
   * @param[in] tflux_in  sensible heat flux H [W/m2]
   * @param[in] qflux_in  latent heat flux LE [W/m2]
   * @param[in] sfc_tau   kinematic surface stress TAU [m2/s2]
   */
  void set_sfc_params_from_file_fluxes (const int& lev,
                                        const amrex::MultiFab& cons_in,
                                        const amrex::Real& tflux_in,
                                        const amrex::Real& qflux_in,
                                        const amrex::Real& sfc_tau);
 
  /**
   * Fill surface temperature from available SST and TSK data.
   *
   * @param[in] lev  level index
   * @param[in] time interpolation time
   */
  void fill_tsurf_with_sst_and_tsk (const int& lev,
                                    const double& time);
 
  /**
   * Fill surface temperature interpolated from time varying SST file
   *
   * @param[in] lev  level index
   * @param[in] time interpolation time
   */
  void fill_tsurf_with_sfc_sst (const int& lev,
                                const double& time,
                                const amrex::MultiFab& cons_in,
                                const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  /**
   * Overwrite surface temperature with coupled ocean SST where the coupler
   * covers the cell.
   *
   * Runs after fill_tsurf_with_sst_and_tsk so that the lower-boundary data is
   * the base layer: land cells, and water cells with no ocean donor, keep the
   * value written there.
   *
   * @param[in] lev level index
   */
  void fill_tsurf_with_coupled_sst (const int& lev,
                                    const amrex::MultiFab& cons_in,
                                    const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  /**
   * Set t_surf from the prescribed absolute surface temperature (most.abs_surf_temp, plus
   * most.surf_heating_rate times the time), converted to the potential temperature the
   * surface layer works in with the pressure at the surface diagnosed from the lowest
   * cell, as for coupled SST. Like update_surf_temp, which it replaces in that case, it
   * writes every cell of the surface.
   *
   * @param[in] lev       level index
   * @param[in] time      elapsed simulation time
   * @param[in] cons_in   conserved state used to evaluate surface pressure
   * @param[in] z_phys_nd nodal physical-height field
   */
  void fill_tsurf_with_abs_surf_temp (const int& lev,
                                      const double& time,
                                      const amrex::MultiFab& cons_in,
                                      const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  /**
   * Overwrite the land surface temperature with an external skin temperature.
   *
   * The skin is an absolute temperature [K]; it is converted to the potential
   * temperature the surface layer works in with the pressure at the surface
   * diagnosed from the lowest cell, as for coupled SST. Runs last among the
   * surface-temperature writers, so on land it replaces what they wrote.
   *
   * @param[in] lev       level index
   * @param[in] cons_in   conserved state used to evaluate surface pressure
   * @param[in] z_phys_nd nodal physical-height field
   */
  void fill_tsurf_with_skin_temperature (const int& lev,
                                         const amrex::MultiFab& cons_in,
                                         const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  /**
   * Overwrite the land surface mixing ratio from the external skin temperature and the
   * moisture field of set_skin_moisture, with q_air the reference-height mixing ratio the
   * flux uses. Without resistances: beta * q_sat(T_skin, p_s) + (1 - beta) * q_air
   * (erf_surface_moisture::surface_mixing_ratio) with beta the soil-water factor, so the
   * moisture flux is beta times the potential one. With resistances:
   * erf_surface_moisture::two_source_surface_mixing_ratio of the aerodynamic resistance
   * (from the last u* and Obukhov length; neutral before the first flux), the canopy
   * resistance with its vapour-deficit factor, the soil resistance and the soil pore air's
   * relative humidity. Needs set_skin_temperature and the averages of this step.
   *
   * @param[in] lev       level index
   * @param[in] cons_in   conserved state (surface pressure and q_air from the lowest cell)
   * @param[in] z_phys_nd nodal heights, or nullptr on a flat mesh
   */
  void fill_qsurf_with_skin_moisture (const int& lev,
                                      const amrex::MultiFab& cons_in,
                                      const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  /**
   * Fill surface moisture from saturation specific humidity.
   *
   * @param[in] lev       level index
   * @param[in] cons_in   conserved state used to evaluate surface pressure
   * @param[in] z_phys_nd nodal physical-height field
   */
  void fill_qsurf_with_qsat (const int& lev,
                             const amrex::MultiFab& cons_in,
                             const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  void set_pblh (const int& lev, const amrex::MultiFab& pblh_in);
 
  /**
   * Updates current time index for interpolating data from SFC/SST file.
   *
   * @param[in] time elapsed time
   */
  void update_sfc_time_index (const amrex::Real& time);
 
  /**
   * Interpolates the SFC/SST data at the given column and time
   *
   * @param[in] time elapsed time
   * @param[in] col  column index of file data
   */
  amrex::Real interpolate_sfc_column (const amrex::Real& time,
                                      int col) const;
 
  /**
   * Fill surface temperature from the LSM surface-temperature field, converted to
   * potential temperature at the surface pressure when the field is an absolute
   * temperature.
   *
   * @param[in] lev       level index
   * @param[in] cons_in   conserved state (surface pressure from the lowest cell)
   * @param[in] z_phys_nd nodal heights, or nullptr on a flat mesh
   */
  void get_lsm_tsurf (const int& lev,
                      const amrex::MultiFab& cons_in,
                      const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
  /**
   * Wrapper around compute_pblh.
   *
   * @param[in]     lev              level index
   * @param[in,out] vars             state variables used by the PBL-height calculation
   * @param[in]     z_phys_cc        cell-centered physical-height field
   * @param[in]     moisture_indices indices for moisture components
   */
  void update_pblh (const int& lev,
                    amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
                    amrex::MultiFab* z_phys_cc,
                    const MoistureComponentIndices& moisture_indices);
 
  /**
   * Set the PBL height of a refined level whose grids do not run from the ground to the top
   * of the domain in every column from the next coarser level, column by column (see
   * compute_pblh).  Public because it launches a device lambda.
   *
   * @param[in] lev Current level (> 0)
   */
  void fill_pblh_from_coarser (const int& lev);
 
  /**
   * The PBL height of the next coarser level, injected onto the grids of this level over all
   * of their cells, ghost cells included.  Public because it launches a device lambda.
   *
   * @param[in]  lev  Current level (> 0)
   * @param[out] fine MultiFab on the grids of pblh[lev]
   */
  void inject_pblh_from_coarser (const int& lev, amrex::MultiFab& fine) const;
 
  /**
   * Compute planetary-boundary-layer height with the selected estimator.
   *
   * @param[in]     lev              level index
   * @param[in,out] vars             state variables used by the PBL-height calculation
   * @param[in]     z_phys_cc        cell-centered physical-height field
   * @param[in]     est              PBL-height estimator functor
   * @param[in]     moisture_indice  indices for moisture components
   */
  template <typename PBLHeightEstimator>
  void compute_pblh (const int& lev,
                     amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
                     amrex::MultiFab* z_phys_cc,
                     const PBLHeightEstimator& est,
                     const MoistureComponentIndices& moisture_indice);
 
  /**
   * Read custom roughness data for one level.
   *
   * @param[in] lev   level index
   * @param[in] fname roughness-data file name
   */
  void read_custom_roughness (const int& lev,
                              const std::string& fname);
 
  /**
   * Update prescribed surface temperature from the configured heating rate.
   *
   * @param[in] time elapsed simulation time
   */
  void update_surf_temp (const double& time)
  {
      // NOTE: this is a whole-domain setVal, so it overwrites the SST/TSK fill
      //       done earlier in update_fluxes. Coupled SST is applied after this
      //       call and therefore still wins on the water cells it covers.
      // An absolute surface temperature is written, converted, by fill_tsurf_with_abs_surf_temp
      if (surf_temp_is_absolute) { return; }
      if (surf_heating_rate != 0) {
          // Use the actual size of t_surf, not m_geom.size(), which is always
          // max_level+1 and so runs past the levels that exist.  t_surf is sized for
          // all levels up front but filled one level at a time, so we also have to
          // skip the entries that have not been allocated yet.
          int nlevs = static_cast<int>(t_surf.size());
          for (int lev = 0; lev < nlevs; lev++) {
              if (!t_surf[lev]) { continue; }
              t_surf[lev]->setVal(surf_temp + surf_heating_rate * static_cast<amrex::Real>(time));
              amrex::Print() << "Surface temp at t=" << time << ": "
                             << surf_temp + surf_heating_rate * time << std::endl;
          }
      }
  }
 
  /**
   * Update MOST-average field pointers.
   *
   * @param[in] lev        level index
   * @param[in] vars_old   old-time state variables
   * @param[in] Theta_prim primitive potential-temperature fields by level
   * @param[in] Qv_prim    primitive water-vapor fields by level
   * @param[in] Qr_prim    primitive rain-water fields by level
   */
  void update_mac_ptrs (const int& lev,
                        amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
                        amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
                        amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
                        amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qr_prim)
  {
      m_ma.update_field_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
  }
 
  /**
   * Return the friction-velocity field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_u_star (const int& lev) { return u_star[lev].get(); }
 
  /**
   * Return the convective velocity scale field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_w_star (const int& lev) { return w_star[lev].get(); }
 
  /**
   * Do we actually compute w*?  If not then the field returned by get_w_star holds
   * only the value it was initialized with, and must not be reported as a diagnostic.
   */
  [[nodiscard]] bool computes_w_star () const { return m_include_wstar; }
 
  /**
   * Do we actually compute the PBL height?  If not then the field returned by get_pblh
   * holds only the value it was initialized with, and must not be reported as a diagnostic.
   */
  [[nodiscard]] bool computes_pblh () const { return (pblh_type != PBLHeightCalcType::None); }
 
  /**
   * Return the temperature scale field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_t_star (const int& lev) { return t_star[lev].get(); }
 
  /**
   * Return the moisture scale field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_q_star (const int& lev) { return q_star[lev].get(); }
 
  /**
   * Return the Obukhov length field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_olen (const int& lev) { return olen[lev].get(); }
 
  /**
   * Return the planetary-boundary-layer-height field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_pblh (const int& lev) { return pblh[lev].get(); }
 
  /**
   * Return a MOST-average field.
   *
   * @param[in] lev  level index
   * @param[in] comp component index
   */
  const amrex::MultiFab* get_mac_avg (const int& lev, int comp)
  {
      return m_ma.get_average(lev, comp);
  }
 
  /**
   * Return whether the MOST averages are filtered in time.
   */
  bool mac_avg_is_time_averaged () const { return m_ma.do_time_averaging(); }
 
  /**
   * Return the number of MOST-average components.
   */
  int get_num_mac_avg () const { return m_ma.get_navg(); }
 
  /**
   * Return whether the time filter at this level holds meaningful history.
   *
   * @param[in] lev level index
   */
  bool mac_avg_is_initialized (const int& lev) const { return m_ma.time_avg_is_initialized(lev); }
 
  /**
   * Declare the time filter at this level to hold meaningful history (restart).
   *
   * @param[in] lev level index
   */
  void set_mac_avg_initialized (const int& lev) { m_ma.set_time_avg_initialized(lev); }
 
  /**
   * Return a MOST-average field for modification (restart).
   *
   * @param[in] lev  level index
   * @param[in] comp component index
   */
  amrex::MultiFab* get_mac_avg_ptr (const int& lev, int comp) { return m_ma.get_average(lev, comp); }
 
  /**
   * Return the filtered plane averages, which hold the filter state for the
   * plane and EB averaging policies.
   *
   * @param[in] lev level index
   */
  amrex::Vector<amrex::Real> get_mac_plane_avg (const int& lev) const { return m_ma.get_plane_average(lev); }
 
  /**
   * Restore the filtered plane averages from a checkpoint; returns false if the
   * checkpoint does not hold what this run expects.
   *
   * @param[in] lev  level index
   * @param[in] pavg filtered plane averages, one per average component
   */
  bool set_mac_plane_avg (const int& lev, const amrex::Vector<amrex::Real>& pavg) { return m_ma.set_plane_average(lev, pavg); }
 
  /**
   * Return the surface-temperature field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_t_surf (const int& lev) { return t_surf[lev].get(); }
 
  /**
   * Take the land surface temperature of level lev from an external absolute
   * skin temperature [K] on the surface-layer layout, or stop doing so (nullptr).
   * update_fluxes reads it; the owner keeps it alive and resets the pointer
   * whenever it reallocates the field.
   *
   * @param[in] lev    level index
   * @param[in] t_skin absolute skin temperature, or nullptr
   */
  void set_skin_temperature (const int& lev, const amrex::MultiFab* t_skin)
  {
      if (lev >= static_cast<int>(m_skin_tsurf_lev.size())) {
          m_skin_tsurf_lev.resize(lev + 1, nullptr);
      }
      m_skin_tsurf_lev[lev] = t_skin;
  }
 
  /**
   * Why this surface layer cannot take its surface temperature from an external
   * skin temperature, or an empty string when it can: the skin sets t_surf, so the
   * surface layer must compute its heat flux from t_surf (a surface temperature given,
   * no heating rate) on a non-EB lower surface without rotated surface fluxes, and no
   * land model may own t_surf.
   */
  [[nodiscard]] std::string skin_temperature_conflict () const
  {
      if (m_terrain_type == TerrainType::EB) {
          return "the surface layer is on EB terrain";
      }
      // Rotated, the flux is split over hfx1/hfx2/hfx3 and the balance, which removes
      // the vertical-face flux hfx3 only, would lose cos(slope) of what the air gains.
      if (m_rotate) {
          return "erf.use_rotate_surface_flux is set";
      }
      if (theta_type != ThetaCalcType::SURFACE_TEMPERATURE) {
          return "the surface layer does not compute its heat flux from a surface "
                 "temperature (set erf.most.pot_surf_temp or erf.most.abs_surf_temp, not a "
                 "flux or custom/bulk fluxes)";
      }
      if (surf_heating_rate != amrex::Real(0)) {
          return "erf.most.surf_heating_rate is set";
      }
      if (m_has_lsm_tsurf || use_surface_model) {
          return "a land-surface or surface model supplies the surface temperature";
      }
      return "";
  }
 
  /**
   * Take the land surface mixing ratio of level lev from the skin set by
   * set_skin_temperature (fill_qsurf_with_skin_moisture), or stop doing so (nullptr).
   * The owner keeps the field alive and resets the pointer whenever it reallocates it.
   *
   * @param[in] lev             level index
   * @param[in] moisture        per-column field in the components of
   *                            erf_surface_moisture::Component, on the surface-layer layout
   * @param[in] use_resistances form q_surf from the soil resistance and pore humidity and,
   *                            where the field's VegetationFraction is positive, the
   *                            canopy's (two_source_surface_mixing_ratio), not from the
   *                            soil-water factor alone:
   *                            the two-stream balance passes true whenever it has a soil type
   * @param[in] hs              the canopy's vapour-deficit parameter [1/(kg/kg)]; 0 without
   *                            vegetation
   */
  void set_skin_moisture (const int& lev, const amrex::MultiFab* moisture,
                          bool use_resistances = false, amrex::Real hs = amrex::Real(0))
  {
      if (lev >= static_cast<int>(m_skin_moisture_lev.size())) {
          m_skin_moisture_lev.resize(lev + 1, nullptr);
          m_skin_moisture_resistances_lev.resize(lev + 1, 0);
          m_skin_moisture_hs_lev.resize(lev + 1, amrex::Real(0));
      }
      m_skin_moisture_lev[lev] = moisture;
      m_skin_moisture_resistances_lev[lev] = use_resistances ? 1 : 0;
      m_skin_moisture_hs_lev[lev] = hs;
  }
 
  /**
   * Why this surface layer cannot take its surface mixing ratio from the skin
   * (set_skin_moisture), or an empty string when it can: everything
   * skin_temperature_conflict asks, plus a moisture model and a moisture flux computed
   * from a surface mixing ratio (erf.most.surf_moist given).
   */
  [[nodiscard]] std::string skin_moisture_conflict () const
  {
      const std::string temperature = skin_temperature_conflict();
      if (!temperature.empty()) { return temperature; }
      if (!use_moisture) {
          return "the run has no moisture model";
      }
      if (moist_type != MoistCalcType::SURFACE_MOISTURE) {
          return "the surface layer does not compute its moisture flux from a surface "
                 "mixing ratio (set erf.most.surf_moist, not a flux)";
      }
      return "";
  }
 
  /**
   * The land roughness length [m] to use when erf.most.z0 is not given (the two-stream
   * balance's, from Noah-MP's tables). Call before the levels are made, which set z0 from
   * it.
   *
   * Taking it also replaces the erf.most.z0 entry that the constructor's queryAdd put in
   * the ParmParse table. job_info (ParmParse::dumpTable) then records the roughness the
   * run used.
   *
   * @return whether the roughness was taken: an erf.most.z0 given in the inputs stands.
   */
  bool set_default_roughness (amrex::Real z0)
  {
      if (m_z0_given) { return false; }
      z0_const = z0;
      amrex::ParmParse pp(m_pp_prefix);
      pp.remove("most.z0");
      pp.add("most.z0", z0_const);
      return true;
  }
 
  /** The constant roughness length [m] the levels take. */
  [[nodiscard]] amrex::Real roughness () const { return z0_const; }
 
  /**
   * Return whether prescribed surface fluxes are active.
   */
  bool use_sfc_fluxes () const { return m_use_sfc_fluxes; }
 
  /**
   * Return whether the surface fluxes are rotated with the terrain slope
   * (erf.use_rotate_surface_flux), splitting them over the x, y and z faces.
   */
  bool rotates_surface_fluxes () const { return m_rotate; }
 
  /**
   * Set the surface-temperature field to a constant value.
   *
   * @param[in] lev   level index
   * @param[in] tsurf surface temperature
   */
  void set_t_surf (const int& lev, const amrex::Real tsurf) { t_surf[lev]->setVal(tsurf); }
 
  /**
   * Return the surface-moisture field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_q_surf (const int& lev) { return q_surf[lev].get(); }
 
  /**
   * Set the surface-moisture field to a constant value.
   *
   * @param[in] lev   level index
   * @param[in] qsurf surface moisture
   */
  void set_q_surf (const int& lev, const amrex::Real qsurf) { q_surf[lev]->setVal(qsurf); }
 
  /**
   * Return the surface-diagnostic provenance field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_surface_diagnostic_source (const int& lev) { return surface_diagnostic_source[lev].get(); }
 
  /**
   * Return the 2D MF of zref heights.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_zref (const int& lev) { return m_ma.get_zref(lev); }
 
  /**
   * Return the minimum reference height for one level.
   *
   * @param[in] lev level index
   */
  amrex::Real get_zref_min (const int& lev) { return (m_ma.get_zref(lev))->min(0); }
 
  /**
   * Return the roughness-height field.
   *
   * @param[in] lev level index
   */
  amrex::MultiFab* get_z0 (const int& lev) { return &z_0[lev]; }
 
  /**
   * Return whether variable sea roughness is active.
   */
  bool have_variable_sea_roughness () { return m_var_z0; }
 
  /**
   * Return the land-mask field.
   *
   * @param[in] lev level index
   */
  amrex::iMultiFab* get_lmask (const int& lev) { return m_lmask_lev[lev][0]; }
 
  /**
   * Compute the minimum land-mask value over valid and optional ghost cells.
   *
   * @param[in] lmask  land-mask field
   * @param[in] nghost number of ghost cells included in the reduction
   */
  int lmask_min_reduce (amrex::iMultiFab& lmask,
                        const int& nghost)
  {
      int lmask_min = amrex::ReduceMin(lmask, nghost, [=] AMREX_GPU_HOST_DEVICE(
                      amrex::Box const& bx, amrex::Array4<int const> const& lm_arr) -> int
                      {
                          int locmin = std::numeric_limits<int>::max();
                          const auto lo = lbound(bx);
                          const auto hi = ubound(bx);
                          for (int j = lo.y; j <= hi.y; ++j) {
                              for (int i = lo.x; i <= hi.x; ++i) {
                                  locmin = std::min(locmin, lm_arr(i, j, 0));
                              }
                          }
                          return locmin;
                      });
 
      return lmask_min;
  }
 
  /**
   * Update one stored sea-surface-temperature pointer.
   *
   * @param[in] lev     level index
   * @param[in] itime   time-slice index
   * @param[in] sst_ptr sea-surface-temperature field pointer
   */
  void update_sst_ptr (const int lev, const int itime, amrex::MultiFab* sst_ptr) {
      AMREX_ALWAYS_ASSERT(itime >= 0 && itime < static_cast<int>(m_sst_lev[lev].size()));
      m_sst_lev[lev][itime] = sst_ptr;
  }
 
  /**
   * Update one stored skin-temperature pointer.
   *
   * @param[in] lev     level index
   * @param[in] itime   time-slice index
   * @param[in] tsk_ptr skin-temperature field pointer
   */
  void update_tsk_ptr (const int lev, const int itime, amrex::MultiFab* tsk_ptr) {
      AMREX_ALWAYS_ASSERT(itime >= 0 && itime < static_cast<int>(m_tsk_lev[lev].size()));
      m_tsk_lev[lev][itime] = tsk_ptr;
  }
 
  /**
   * Hand over the coupled sea-surface temperature and its per-cell coverage flag.
   *
   * Both pointers stay owned by the caller and must remain valid, and on the same
   * layout, until replaced. Passing a null sst_ptr retracts the coupled lane, so
   * the lower-boundary SST/TSK value stands everywhere again.
   *
   * @param[in] lev       level index
   * @param[in] sst_ptr   coupled sea-surface temperature [K]
   * @param[in] valid_ptr per-cell coverage flag; nonzero where the coupler
   *                      supplied a value. Null means no cell is covered.
   */
  void update_coupled_sst_ptr (const int lev,
                               amrex::MultiFab* sst_ptr,
                               amrex::iMultiFab* valid_ptr) {
      m_coupled_sst_lev[lev]   = sst_ptr;
      m_coupled_sst_valid_lev[lev] = valid_ptr;
  }
 
  /**
   * Declare that an ocean coupler will supply SST for this run.
   *
   * Must be called before make_SurfaceLayer_at_level, which uses it to select
   * ThetaCalcType::SURFACE_TEMPERATURE.
   *
   * @param[in] active whether coupled SST is configured
   */
  void set_coupled_sst_active (const bool active) { m_use_coupled_sst = active; }
 
  /**
   * Set the domain faces that use the surface-layer boundary condition.
   *
   * The face set is used to suppress transpose stress writes at shared edges;
   * the face-normal component remains owned by its corresponding face path.
   *
   * @param[in] active_faces whether a specific x/y/z lo/hi face is enabled
   */
  void set_surface_layer_faces (const amrex::GpuArray<int, AMREX_SPACEDIM*2>& active_faces)
  {
      m_surface_layer_faces = active_faces;
  }
 
  /**
   * Reads columns of data from a text file, returning each column in a vector.
   *
   * @param[in] fname       path to text file
   * @param[in] skip_nlines number of lines to skip before reading data (e.g, header lines)
   * @return Vector containing each column in the file as a vector
   */
  static amrex::Vector<amrex::Vector<amrex::Real>>
  read_cols (const std::string& fname, const int skip_nlines = 1);
 
  /**
   * Fills surface layer data for lateral faces. A normal FillBoundary on
   * lateral faces would bring in data from other ranks or interior grids
   * that do not own data on the lateral face.
   *
   * @param[in] lev level index
   * @param[in,out] selected_field if non-null, fill only this field;
   *                              otherwise fill all lateral surface fields
   */
  void fill_lateral_surface_parameter_ghosts (const int& lev,
                                               amrex::MultiFab* selected_field = nullptr);
 
  enum struct FluxCalcType {
    MOENG = 0,  ///< Moeng functional form
    CUSTOM,     ///< Custom constant flux functional form
    BULK_COEFF, ///< Bulk transfer coefficient functional form
    ROTATE,     ///< Terrain rotation flux functional form
    RICO
  };
 
  enum struct ThetaCalcType {
    ADIABATIC = 0,
    HEAT_FLUX,          ///< Heat-flux specified
    SURFACE_TEMPERATURE ///< Surface temperature specified
  };
 
  enum struct MoistCalcType {
    ADIABATIC = 0,
    MOISTURE_FLUX,      ///< Qv-flux specified
    SURFACE_MOISTURE    ///< Surface Qv specified
  };
 
  enum struct RoughCalcType {
    CONSTANT = 0, ///< Constant z0
    CHARNOCK,
    MODIFIED_CHARNOCK,
    DONELAN,
    WAVE_COUPLED
  };
 
  enum struct PBLHeightCalcType { None, MYNN25, YSU, MRF, YSUNew };
 
  FluxCalcType flux_type{FluxCalcType::MOENG};
  ThetaCalcType theta_type{ThetaCalcType::ADIABATIC};
  MoistCalcType moist_type{MoistCalcType::ADIABATIC};
  RoughCalcType rough_type_land{RoughCalcType::CONSTANT};
  RoughCalcType rough_type_sea{RoughCalcType::CHARNOCK};
  PBLHeightCalcType pblh_type{PBLHeightCalcType::None};
 
private:
  // Set in constructor
  amrex::Orientation m_face;
  std::string m_pp_prefix;
  amrex::Vector<amrex::Geometry> m_geom;
  bool m_rotate = false;
  amrex::GpuArray<int, AMREX_SPACEDIM*2> m_surface_layer_faces{};
  double m_start_low_time;
  double m_final_low_time;
  double m_low_time_interval;
 
  bool m_include_wstar = false;
  amrex::Real z0_const{amrex::Real(0.1)};
  // erf.most.z0 was given (set_default_roughness then leaves z0_const alone).
  bool m_z0_given{false};
  amrex::Real default_land_surf_temp{amrex::Real(300.)};
  // Negative sentinels: a surface temperature is in kelvin and a surface moisture is a
  // mixing ratio, so neither can legitimately be negative.  The constructor probes
  // most.pot_surf_temp / most.abs_surf_temp / most.surf_moist with these and tests the
  // resulting value to decide whether the user supplied them.  They were previously left
  // uninitialized, which was safe only because a plain query() leaves them untouched on
  // a miss.
  amrex::Real surf_temp{amrex::Real(-1.)};
  // surf_temp (and surf_heating_rate) is an absolute temperature (most.abs_surf_temp),
  // converted to t_surf's potential temperature with the surface pressure each step.
  bool surf_temp_is_absolute{false};
  amrex::Real surf_heating_rate{0};
  amrex::Real surf_temp_flux{0};
  amrex::Real default_land_surf_moist{zero};
  amrex::Real surf_moist{amrex::Real(-1.)};
  amrex::Real surf_moist_flux{0};
  amrex::Real custom_ustar{0};
  amrex::Real custom_tstar{0};
  amrex::Real custom_qstar{0};
  amrex::Real custom_rhosurf{0}; // use specified value instead of rho from first cell
  bool specified_rho_surf{false};
  amrex::Real cnk_a{amrex::Real(0.0185)};
  bool smooth_flow_visc{true};
  amrex::Real depth{amrex::Real(30.0)};
  amrex::Vector<amrex::MultiFab> z_0;
  bool m_var_z0{false};
 
  amrex::Real rico_theta_z0{amrex::Real(298.0)};
  amrex::Real rico_qsat_z0{amrex::Real(0.001)};
 
  bool use_moisture;
  bool m_has_lsm_fluxes = false;
  bool m_has_lsm_tsurf  = false;
  int  m_lsm_tsurf_indx = -1;
  bool m_lsm_tsurf_is_theta = false;
  bool m_use_lsm_skin_temperature = false;
 
  // variables for sfc/sst forcing from file
  bool m_use_sfc_fluxes = false;
  bool m_use_sfc_sst = false;
  int sfc_time_ind = 0;
  amrex::Vector<amrex::Vector<amrex::Real>> sfc;
  amrex::Real sfc_qflux = zero;
  amrex::Real sfc_tflux = zero;
  amrex::Real sfc_ustar = zero;
 
  // Set from solverChoice.use_coupled_sst before the per-level Define runs.
  // Configuration only: true means an ocean coupler will supply SST on some
  // subset of the water cells, not that it has supplied any yet.
  bool m_use_coupled_sst = false;
 
  amrex::Real m_Cd = zero;
  amrex::Real m_Ch = zero;
  amrex::Real m_Cq = zero;
  bool m_ignore_sst = false;
 
  amrex::Vector<const eb_*> m_eb_vec;
  TerrainType m_terrain_type;
  bool m_has_immersed_buildings = false;
  amrex::Real m_rdOcp = RdoCp;
  MOSTAverage m_ma;
  SurfaceModel* m_surf_model = nullptr;
  bool use_surface_model = false;
  bool surf_model_fluxes = false;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> u_star;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> w_star;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_star;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> q_star;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> olen;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> pblh;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_surf;
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> q_surf;
  // The planar BoxArray is the z-collapse of the 3D BoxArray, so a 3D BoxArray split in z
  // gives duplicate planar boxes (see PlanarBoundary). fill_planar_boundary goes through
  // this for u*, w* (only when it is computed), t*, q*, L, t_surf, q_surf, pblh and the
  // diagnostic source, after each is computed on the surface copies.
  amrex::Vector<PlanarBoundary> m_planar_bndry;
  // Columns for the PBL-height estimator.  The estimator scans each box it is given from its
  // lowest cell to its highest, so the box must start at the ground.  When some box of a
  // level does not (grids stacked in z, or a refined region aloft), compute_pblh runs the
  // estimator on these boxes instead: the runs of cells that start at the ground, each joined
  // into one box that ends where the grids above that ground end.  Built by
  // define_pblh_columns when the grids of the level change.
  struct PBLHColumns
  {
      amrex::BoxArray ba;                 // Grids the columns were built for ...
      amrex::DistributionMapping dm;      // ... and their ranks
      bool needed = false;                // Whether some box does not start at the ground
      bool full_height = true;            // Whether every column runs from the ground to the top
      amrex::BoxArray ba_col;             // The columns; empty on a level entirely aloft
      amrex::BoxArray ba_col2d;           // Their ground plane
      amrex::DistributionMapping dm_col;
      amrex::MultiFab hold_col;           // No data; keep the copy metadata alive
      amrex::MultiFab hold_col2d;
  };
  amrex::Vector<PBLHColumns> m_pblh_columns;
  void define_pblh_columns (const int& lev,
                            const amrex::BoxArray& ba,
                            const amrex::DistributionMapping& dm);
  // Diagnostic-only provenance mask. This records the cell-centered source
  // path used by the SurfaceLayer scalar diagnostic path. It must not feed
  // back into flux calculations, boundary conditions, or model state.
  amrex::Vector<std::unique_ptr<amrex::MultiFab>> surface_diagnostic_source;
 
  amrex::Vector<amrex::Vector<amrex::MultiFab*>>  m_sst_lev;
  amrex::Vector<amrex::Vector<amrex::MultiFab*>>  m_tsk_lev;
  amrex::Vector<amrex::Vector<amrex::iMultiFab*>> m_lmask_lev;
 
  // Coupled SST from an external ocean model, plus the per-cell coverage flag
  // that says whether the coupler actually supplied a value for that cell.
  // Both are owned by ERF and borrowed here, exactly like m_sst_lev. A null
  // m_coupled_sst_lev[lev] means "no coupled SST has arrived yet"; a zero entry in
  // m_coupled_sst_valid_lev means "this cell has no ocean donor", and in both cases the
  // lower-boundary SST/TSK value filled earlier in update_fluxes stands.
  amrex::Vector<amrex::MultiFab*>  m_coupled_sst_lev;
  amrex::Vector<amrex::iMultiFab*> m_coupled_sst_valid_lev;
  // External absolute skin temperature per level (set_skin_temperature), owned
  // by the caller; nullptr where none is supplied.
  amrex::Vector<const amrex::MultiFab*> m_skin_tsurf_lev;
  // What that skin hands over for its surface mixing ratio, with how to form it from that
  // field, per level (set_skin_moisture), so setting one level never changes another's.
  amrex::Vector<const amrex::MultiFab*> m_skin_moisture_lev;
  amrex::Vector<int> m_skin_moisture_resistances_lev;
  amrex::Vector<amrex::Real> m_skin_moisture_hs_lev;
  amrex::Vector<amrex::Vector<amrex::MultiFab*>>  m_lsm_data_lev;
  amrex::Vector<amrex::Vector<amrex::MultiFab*>>  m_lsm_flux_lev;
  amrex::Vector<std::string> m_lsm_data_name;
  amrex::Vector<std::string> m_lsm_flux_name;
  amrex::Vector<amrex::MultiFab*> m_Hwave_lev;
  amrex::Vector<amrex::MultiFab*> m_Lwave_lev;
  amrex::Vector<amrex::MultiFab*> m_eddyDiffs_lev;
 
  bool m_update_k_rans = false;
  amrex::Real inv_Cmu2 = zero;
  amrex::Real theta_ref = zero;
};
 
#endif /* SURFACELAYER_H */