ERF.H

Go to the documentation of this file.

#ifndef ERF_H_
#define ERF_H_
 
#include <string>
#include <limits>
#include <memory>
#include <set>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include <AMReX_AmrCore.H>
#include <AMReX_BCRec.H>
#include <AMReX_Print.H>
 
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_ParmParse.H>
#include <AMReX_MultiFabUtil.H>
#include <AMReX_FillPatchUtil.H>
#include <AMReX_VisMF.H>
#include <AMReX_PhysBCFunct.H>
#include <AMReX_YAFluxRegister.H>
#include <AMReX_ErrorList.H>
#include <AMReX_EBFabFactory.H>
#include <AMReX_EBMultiFabUtil.H>
#include <AMReX_TableData.H>
 
#include <ERF_EB.H>
#include "Prob/ERF_CloudChamberBudget.H"
#include "Prob/ERF_CloudChamber.H"
 
#ifdef ERF_USE_FFT
#include <AMReX_FFT_Poisson.H>
#endif
 
#ifdef AMREX_MEM_PROFILING
#include <AMReX_MemProfiler.H>
#endif
 
#include "ERF_ProbCommon.H"
 
#include <ERF_NumericalConstants.H>
#include <ERF_IndexDefines.H>
#include <ERF_DataStruct.H>
#include <ERF_BaseStateParams.H>
#include <ERF_TurbPertStruct.H>
#include <ERF_InputSoundingData.H>
#include <ERF_LargeScaleForcingData.H>
#include <ERF_InputSpongeData.H>
#include <ERF_SurfaceLayer.H>
#include <ERF_Derive.H>
#include <ERF_ReadBndryPlanes.H>
#include <ERF_WriteBndryPlanes.H>
#include <ERF_MRI.H>
#include <ERF_PhysBCFunct.H>
#include <ERF_WallScalarBC.H>
#include <ERF_FillPatcher.H>
#include <ERF_SampleData.H>
#include <ERF_StationSampler.H>
#include <ERF_ObsNudging.H>
#include <ERF_ForestDrag.H>
#include <ERF_IBFaceSet.H>
#include <ERF_Provenance.H>
 
#ifdef ERF_USE_PARTICLES
#include "ERF_ParticleData.H"
#endif
 
#include "ERF_EulerianMicrophysics.H"
#include "ERF_LagrangianMicrophysics.H"
#include "ERF_LandSurface.H"
#include "ERF_Urban.H"
#include "ERF_TwoStreamRadiation.H"
 
#ifdef ERF_USE_WINDFARM
#include "ERF_WindFarm.H"
#endif
#include <ERF_SurfaceModel.H>
 
#include <ERF_RadiationModels.H>
 
#ifdef ERF_USE_EAMXX_SHOC
#include "ERF_ShocInterface.H"
#endif
 
#include "ERF_ShocDriver.H"
#include "AuxiliaryState/ERF_AuxiliaryInertTracer.H"
 
#include <iostream>
 
#ifdef AMREX_LAZY
#include <AMReX_Lazy.H>
#endif
 
#ifndef AMREX_USE_MPI
using amrex::MPI_COMM_WORLD;
using amrex::MPI_Comm;
#endif
 
#ifdef ERF_USE_MULTIBLOCK
class MultiBlockContainer;
#endif
 
// Describes one component of a 2D plotfile; defined in
// IO/ERF_Plotfile2DSampledLevel.H, which ERF.H does not need to include.
namespace plotfile2d { struct Plotfile2DOutputDescriptor; }
namespace erf_sbm { class SBMStateManager; class SBMTransport; }
 
/**
 * Enum of possible interpolation types between coarse/fine
*/
AMREX_ENUM(StateInterpType,
    FullState, Perturbational
);
 
/**
 * Main class in ERF code, instantiated from main.cpp
*/
 
class ERF
    : public amrex::AmrCore
{
public:
 
    ////////////////
    // public member functions
 
    // constructor - reads in parameters from inputs file
    //             - sizes multilevel arrays and data structures
    ERF ();
    ~ERF () override;
 
    void ERF_shared ();
 
    // Declare a default move constructor so we ensure the destructor is
    // not called when we return an object of this class by value
    ERF (ERF&&)  noexcept = delete;
 
    // Declare a default move assignment operator
    ERF& operator= (ERF&& other)  noexcept = delete;
 
    // Delete the copy constructor
    ERF (const ERF& other) = delete;
    //
    // Delete the copy assignment operator
    ERF& operator= (const ERF& other) = delete;
 
    // Advance solution to final time
    void Evolve ();
    double EvolveOneStep (double time, double dt_request);
 
    void ConfigureDriverAtmosToOceanCoupling (bool use_coupling_driver,
                                              bool use_two_way_coupling,
                                              bool use_state_contract);
    void SetDriverAtmosToOceanStateContract (bool use_state_contract);
    // Driver currently consumes BoxArray/DistributionMapping only. If we later
    // need exact copy semantics, promote this to expose IndexType/ngrow
    // explicitly or return the underlying MultiFab pointer/reference instead.
    void GetOceanToAtmosSurfaceLayout (amrex::BoxArray& ba,
                                       amrex::DistributionMapping& dm);
    // Returns pointers to level-0 lon_m/lat_m coordinate MultiFabs (populated
    // by the WRF-grid init paths, e.g. ERF_InitFromWRFInput.cpp/ERF_InitFromMetgrid.cpp)
    // used to build ocean->atmos remap weights. Sets both to nullptr if
    // unavailable (e.g. idealized/non-WRF configurations).
    //
    // Staggering contract: both MultiFabs carry CELL index type and hold
    // mass-point values whichever init path filled them -- init_from_wrfinput
    // reads WRF's XLONG / XLAT and init_from_metgrid reads XLONG_M / XLAT_M.
    // Ghost cells at the domain boundary are clamped copies of the edge cell,
    // not extrapolations. These are not corners: a caller building a corner
    // mesh must derive the corners itself. The name says cell centers so that a
    // caller written against the old edge-staggered contract fails to compile
    // instead of silently building its mesh half a cell off.
    void GetOceanToAtmosCellCenterLonLat (const amrex::MultiFab*& lon_cc,
                                          const amrex::MultiFab*& lat_cc) const;
 
    // Returns the level-0, current-time land mask (1 = land, 0 = water) used
    // by the driver to exclude land cells from atmos->ocean remap stencils.
    // Sets lmask to nullptr if unavailable.
    void GetLandMask (const amrex::iMultiFab*& lmask) const;
 
    // True when a radiation model is *configured* at level 0. A coupling driver
    // needs this to know whether the SWrad/LWrad lanes it would hand an ocean
    // model carry real values or just zeros: passing zeros is worse than passing
    // nothing, because the ocean model then suppresses its own radiative forcing
    // in favour of them.
    //
    // Note this is a configuration predicate, not a data-validity one. rad_fluxes
    // is allocated and zero-filled as soon as solverChoice.rad_type != None, so
    // this returns true from level-0 allocation onward - including before the
    // first AdvanceRadiation, when the fluxes are still exactly zero. A caller
    // that needs "the lanes hold real radiation" must additionally know that a
    // radiation step has run.
    [[nodiscard]] bool HasRadiation () const {
        return (!rad_fluxes.empty() && rad_fluxes[0] != nullptr);
    }
 
    // Whether ERF can pack cloud fraction on the atmos->ocean Cloud lane. The
    // pack fills that lane only when moisture is active AND the moisture model
    // carries cloud liquid or ice, so this mirrors that condition exactly; a
    // driver that handed the lane over regardless would be passing its own
    // prefill as though it were an ERF value.
    //
    // Like HasRadiation, this is a configuration predicate rather than a
    // data-validity one: it says the lane will be filled, not that any moisture
    // step has run yet.
    [[nodiscard]] bool HasCloudWater () const {
        return (solverChoice.moisture_type != MoistureType::None) &&
               (solverChoice.moisture_indices.qc != -1 ||
                solverChoice.moisture_indices.qi != -1);
    }
 
    // Which atmos->ocean lanes PackAtmosphericStates actually writes here, as
    // 0/1 in that routine's slot order: the 9-lane state view when
    // state_contract, else the 8-lane flux view.
    //
    // A lane reported 0 is one the pack leaves untouched, so handing it over
    // would pass the driver's own prefill off as an atmospheric value.
    // HasRadiation()/HasCloudWater() answer this for two lanes; this covers all
    // of them, and lives beside the pack so the gates cannot drift apart.
    void GetFillableAtmosToOceanLanes (bool state_contract,
                                       amrex::Vector<int>& can_fill) const;
 
    // weight_mf/index_mf are keyed by WeightFamily slot order
    // {SCALAR_CENTERED, TAUX_FACE, TAUY_FACE, WIND_FACE_TO_CENTER} as built by
    // the driver's ERFRemoraMultiBlockContainer::BuildRemapWeights(); each
    // lane below picks the pair matching its own destination stagger.
    // dst_mskr/dst_msku/dst_mskv are REMORA wet/dry masks (1 = wet, 0 = land)
    // on the rho/u-face/v-face destination layouts. dst_msku/dst_mskv are applied
    // to the TauX/TauY lanes; dst_mskr is applied to the rho-centred *scalar* lanes
    // (Patm, Tair, cloud, and the scalar flux lanes) but deliberately NOT to
    // Uwind/Vwind, which are also rho-centred - so the wind lanes carry live
    // atmospheric values over REMORA land while Patm/Tair carry fallbacks there.
    // REMORA re-masks on its own side, so this asymmetry is currently harmless.
    // Pass nullptr to skip masking.
    void PackAtmosphericStates (amrex::Vector<amrex::MultiFab*>& states, double time,
                                const amrex::Vector<const amrex::MultiFab*>& weight_mf,
                                const amrex::Vector<const amrex::iMultiFab*>& index_mf,
                                int max_stencil_size,
                                const amrex::MultiFab* dst_mskr = nullptr,
                                const amrex::MultiFab* dst_msku = nullptr,
                                const amrex::MultiFab* dst_mskv = nullptr);
    // Apply ocean-model sea-surface temperature to the ERF lower boundary.
    //
    // state[0] carries SST in Kelvin on the layout reported by
    // GetOceanToAtmosSurfaceLayout. erf_coverage is the per-cell coverage flag on
    // that same layout: nonzero where the driver's remap found at least one
    // ocean donor for that ERF cell.
    //
    // Coverage is a contract, not a hint. Where erf_coverage is zero the ocean
    // supplied nothing and state[0] holds the remap's zero fill, which is
    // indistinguishable from a computed zero; ERF keeps its own wrflowinp SST
    // there instead. Passing erf_coverage = nullptr therefore means "no cell is
    // covered", not "every cell is covered" -- a driver that has coverage
    // information must pass it.
    void ApplyOceanSurfaceState (const amrex::Vector<amrex::MultiFab*>& state,
                                 double time,
                                 const amrex::iMultiFab* erf_coverage = nullptr);
 
    // Tag cells for refinement
    void ErrorEst (int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow) override;
 
    // Tag cells in a given box
    void read_box_for_refinement   (std::string& ref_prefix, int& lev_for_box, amrex::RealBox& real_box);
    void update_box_for_refinement (std::string& ref_prefix, int& lev_for_box, amrex::RealBox& real_box, const double time);
 
    // Refinement criteria for hurricanes
    void HurricaneTracker (int lev,
                           double time,
                           const amrex::MultiFab& cc_vel,
                           const amrex::Real velmag_threshold,
                           amrex::TagBoxArray* tags = nullptr);
 
    bool FindInitialEye (int lev,
                         const amrex::MultiFab& cc_vel,
                         const amrex::Real velmag_threshold,
                         amrex::Real& eye_x, amrex::Real& eye_y);
 
    // Hurricane track line output as VTKPolyline
    amrex::Vector<std::array<amrex::Real, 2>> hurricane_track_xy;
    amrex::Vector<std::array<amrex::Real, 2>> hurricane_eye_track_xy,
                                              hurricane_eye_track_latlon,
                                              hurricane_maxvel_vs_time,
                                              hurricane_minpressure_vs_time,
                                              hurricane_tracker_circle;
    amrex::Vector<amrex::MultiFab> weather_forecast_data_1, weather_forecast_data_2;
    amrex::Vector<amrex::Vector<amrex::MultiFab>> forecast_state_1,
                                                  forecast_state_2,
                                                  forecast_state_interp;
    amrex::Vector<amrex::MultiFab> surface_state_1,
                                   surface_state_2,
                                   surface_state_interp;
 
    amrex::Vector<amrex::Vector<amrex::MultiFab>> bckgnd_state;
 
    std::string MakeVTKFilename (int nstep);
    std::string MakeVTKFilename_TrackerCircle (int nstep);
    std::string MakeVTKFilename_EyeTracker_xy (int nstep);
    std::string MakeFilename_EyeTracker_latlon (int nstep);
    std::string MakeFilename_EyeTracker_maxvel (int nstep);
    std::string MakeFilename_EyeTracker_minpressure (int nstep);
    void WriteVTKPolyline (const std::string& filename,
                           amrex::Vector<std::array<amrex::Real, 2>>& points_xy);
 
    void WriteLinePlot (const std::string& filename,
                        amrex::Vector<std::array<amrex::Real, 2>>& points_xy);
 
    void HurricaneEyeTracker (const SolverChoice& solverChoice,
                              const amrex::MultiFab& mf_cc_vel);
 
    void HurricaneEyeTrackerInitial (const SolverChoice& sc,
                                     const amrex::Geometry& geom,
                                     const amrex::Vector<amrex::MultiFab>& S_data,
                                     const amrex::Real& hurricane_eye_latitude,
                                     const amrex::Real& hurricane_eye_longitude);
 
    void HurricaneEyeTrackerNotInitial (const SolverChoice& sc,
                                        const amrex::Geometry& geom,
                                        const amrex::MultiFab& mf_cc_vel,
                                        const amrex::Vector<amrex::MultiFab>& S_data);
 
    void ComputeGlobalMinLocation (const SolverChoice& sc,
                                   const amrex::Geometry& geom,
                                   const amrex::Vector<amrex::MultiFab>& S_data,
                                   amrex::Real* d_val_min_ptr,
                                   int* d_i_min_ptr,
                                   int* d_j_min_ptr,
                                   amrex::Real& global_val_min,
                                   int& global_i_min,
                                   int& global_j_min);
 
    void HurricaneMaxVelTracker (const amrex::Geometry& geom,
                                 const amrex::MultiFab& mf_cc_vel,
                                 const double& time);
 
    void HurricaneMinPressureTracker (MoistureType moisture_type,
                                      const amrex::Geometry& geom,
                                      const amrex::MultiFab& mf_cons_var,
                                      const double& time);
 
    void HurricaneTrackerCircle ();
 
    void ReadStormTrackerRestart ();
 
    // Initialize multilevel data
    void InitData ();
 
    // Initialize multilevel data before MultiBlock
    void InitData_pre ();
 
    // Initialize multilevel data after MultiBlock
    void InitData_post ();
 
    // Interpolate 2D arrays lat_m, lon_m, sinPhi, cosPhi
    void Interp2DArrays (int lev, const amrex::BoxArray& my_ba2d, const amrex::DistributionMapping& my_dm);
 
    // Fill the map factors at lev by interpolation from lev-1
    void interp_mapfac_from_coarse (int lev);
 
    // Fill the surface pressure at lev by interpolation from lev-1
    void interp_psfc_from_coarse (int lev);
 
    // Fill the land mask and the land / soil types at lev by injection from lev-1
    void interp_land_masks_from_coarse (int lev);
 
    // True if this level read its own initialization file, so its map factors
    //     came from that file rather than from its parent
    bool has_own_init_file (int lev) const;
 
    void WriteMyEBSurface ();
 
    // Compute the divergence -- whether EB, no-terrain, flat terrain or general terrain
    void compute_divergence (int lev, amrex::MultiFab& rhs,
                             amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> rho0_u_const,
                             amrex::MultiFab const& mf_mx, amrex::MultiFab const& mf_my,
                             amrex::MultiFab const& mf_vx, amrex::MultiFab const& mf_uy,
                             amrex::MultiFab const& ax_in, amrex::MultiFab const& ay_in,
                             amrex::MultiFab const& dJ_in,
                             amrex::Geometry const& geom_at_lev);
 
    // Project the velocities to be divergence-free -- this is only relevant if anelastic == 1
    void project_initial_velocity (int lev, double time, double dt);
    void project_momenta (int lev, double l_time, double l_dt, amrex::Vector<amrex::MultiFab >& vars);
 
    // Project the velocities to be divergence-free with a thin body
    void project_velocity_tb (int lev, double dt, amrex::Vector<amrex::MultiFab >& vars);
 
    // Calculate wall distance by solving a Poisson equation
    void poisson_wall_dist (int lev);
 
    void make_subdomains (const amrex::BoxList& ba, amrex::Vector<amrex::BoxArray>& bins);
 
#ifdef ERF_USE_FFT
    void build_fft_solvers (int lev);
    void solve_with_fft (int lev, int isub, const amrex::Box& subdomain,
                         amrex::MultiFab& rhs, amrex::MultiFab& p,
                         amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>& fluxes);
#endif
    void solve_with_gmres (int lev, const amrex::Box& subdomain,
                           amrex::MultiFab& rhs, amrex::MultiFab& p,
                           amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>& fluxes,
                           amrex::MultiFab& ax_sub, amrex::MultiFab& ay_sub,
                           amrex::MultiFab& az_sub,
                           amrex::MultiFab& /*dJ_sub*/, amrex::MultiFab& znd_sub);
 
    void ImposeBCsOnPhi (int lev, amrex::MultiFab& phi, const amrex::Box& subdomain);
 
    // Init (NOT restart or regrid)
    void init_only (int lev, double time);
 
    // Restart
    void restart ();
 
    // Check for NaNs and other bad stuff
    void check_state_for_nans (amrex::MultiFab const& S);
    void check_vels_for_nans (amrex::MultiFab const& xvel, amrex::MultiFab const& yvel, amrex::MultiFab const& zvel);
    void check_for_negative_theta (amrex::MultiFab& S);
    void check_for_low_temp       (amrex::MultiFab& S);
 
    void compute_max_pressure_gradient_diagnostic (int lev);
 
    // Max of (buoyancy - vertical gradient of perturbational pressure);
    //    this is identically zero if the state is in HSE and buoyancy_type = 1
    void compute_max_buoyancy_gradp_diagnostic (int lev);
 
    void check_mesh_type (int lev);
 
    // Is it time to write a plotfile or checkpoint?
    // Note that "cur_time" is always in double precision
    bool writeNow (double cur_time, const int nstep, const int plot_int,
                   const double plot_per, const double dt_0, double& last_file_time);
 
    //! True when a level cannot sweep radiation itself -- no box on it holds a whole
    //! column -- and so takes its radiation fields by interpolation from its parent.
    [[nodiscard]] bool rad_level_needs_interpolation (int lev) const;
 
    //! Fill a newly built level's two-stream SEB state from its parent
    void fill_seb_from_coarse (int lev);
 
    //! With erf.radiation.seb_surface_layer_uses_skin, hand the two-stream surface energy
    //! balance's skin temperature of level lev to the zlo surface layer before it
    //! computes its fluxes
    void set_surface_layer_skin (int lev);
 
    //! Start-up checks of the two-stream balance's coupling with the zlo surface layer:
    //! warnings for a scalar default the surface layer's flux replaces and for rotated
    //! surface fluxes, aborts for erf.radiation.seb_surface_layer_uses_skin conflicts
    void check_seb_surface_layer ();
 
    //! check_seb_surface_layer has run
    bool m_seb_surface_layer_checked = false;
 
    //! Per level: the "takes its radiation from its parent, no skin coupling" warning
    //! of set_surface_layer_skin has been printed since the level last coupled
    amrex::Vector<int> m_skin_uncoupled_warned;
 
    //! The zlo surface layer's applied heat and moisture fluxes for the two-stream
    //! surface energy balance's H and LE, or nullptr where it applies none
    void seb_surface_layer_fluxes (int lev,
                                   const amrex::MultiFab*& sens_flux,
                                   const amrex::MultiFab*& laten_flux) const;
 
    //! The conservative surface heat and moisture fluxes (z faces) applied to level
    //! lev, or nullptr where there is none. The sensible_heat_flux/latent_heat_flux
    //! 2D outputs and the two-stream surface energy balance both read them here.
    void surface_flux_sources (int lev,
                               const amrex::MultiFab*& sens_flux,
                               const amrex::MultiFab*& laten_flux) const;
 
    // Called after every level 0 timestep
    void post_timestep (int nstep, double time, double dt_lev);
 
    // Do whatever writing needs to be done at intermediate times
    void WriteAtIntermediateTime (int nstep, double time);
 
    // Do whatever writing needs to be done after time-stepping is complete
    void WriteAtFinalTime ();
 
    // Diagnostics
    void sum_integrated_quantities (double time);
    void sum_derived_quantities (double time);
    void sum_energy_quantities (double time);
 
    void write_1D_profiles (double time);
    void write_1D_profiles_stag (double time);
 
    amrex::Real cloud_fraction (double time);
 
    // Fill the physical boundary conditions for cell-centered velocity (diagnostic only)
    void FillBdyCCVels (amrex::MultiFab& mf_cc_vel, amrex::Geometry& lev_geom);
 
    void sample_points (int lev, double time, amrex::IntVect cell, amrex::MultiFab& mf);
    void sample_lines (int lev, double time, amrex::IntVect cell, amrex::MultiFab& mf);
 
    void derive_diag_profiles (double time,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_u   , amrex::Gpu::HostVector<amrex::Real>& h_avg_v,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_w   , amrex::Gpu::HostVector<amrex::Real>& h_avg_rho,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_th  , amrex::Gpu::HostVector<amrex::Real>& h_avg_ksgs,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_Kmv , amrex::Gpu::HostVector<amrex::Real>& h_avg_Khv,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_qv  , amrex::Gpu::HostVector<amrex::Real>& h_avg_qc,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_qr  ,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_wqv , amrex::Gpu::HostVector<amrex::Real>& h_avg_wqc,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_wqr ,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_qi  , amrex::Gpu::HostVector<amrex::Real>& h_avg_qs,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_qg  ,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_uu  , amrex::Gpu::HostVector<amrex::Real>& h_avg_uv,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_uw,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_vv  , amrex::Gpu::HostVector<amrex::Real>& h_avg_vw,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_ww,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_uth , amrex::Gpu::HostVector<amrex::Real>& h_avg_vth,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_wth, amrex::Gpu::HostVector<amrex::Real>& h_avg_thth,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_ku, amrex::Gpu::HostVector<amrex::Real>& h_avg_kv,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_kw,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_p,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_pu, amrex::Gpu::HostVector<amrex::Real>& h_avg_pv,
                               amrex::Gpu::HostVector<amrex::Real>& h_avg_pw, amrex::Gpu::HostVector<amrex::Real>& h_avg_wthv);
    void derive_diag_profiles_stag (double time,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_u   , amrex::Gpu::HostVector<amrex::Real>& h_avg_v,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_w   , amrex::Gpu::HostVector<amrex::Real>& h_avg_rho,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_th  , amrex::Gpu::HostVector<amrex::Real>& h_avg_ksgs,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_Kmv , amrex::Gpu::HostVector<amrex::Real>& h_avg_Khv,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_qv  , amrex::Gpu::HostVector<amrex::Real>& h_avg_qc,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_qr  ,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_wqv , amrex::Gpu::HostVector<amrex::Real>& h_avg_wqc,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_wqr ,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_qi  , amrex::Gpu::HostVector<amrex::Real>& h_avg_qs,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_qg  ,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_uu  , amrex::Gpu::HostVector<amrex::Real>& h_avg_uv,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_uw,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_vv  , amrex::Gpu::HostVector<amrex::Real>& h_avg_vw,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_ww,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_uth , amrex::Gpu::HostVector<amrex::Real>& h_avg_vth,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_wth, amrex::Gpu::HostVector<amrex::Real>& h_avg_thth,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_ku, amrex::Gpu::HostVector<amrex::Real>& h_avg_kv,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_kw,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_p,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_pu, amrex::Gpu::HostVector<amrex::Real>& h_avg_pv,
                                    amrex::Gpu::HostVector<amrex::Real>& h_avg_pw, amrex::Gpu::HostVector<amrex::Real>& h_avg_wthv);
 
    void derive_stress_profiles (amrex::Gpu::HostVector<amrex::Real>& h_avg_tau11, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau12,
                                 amrex::Gpu::HostVector<amrex::Real>& h_avg_tau13, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau22,
                                 amrex::Gpu::HostVector<amrex::Real>& h_avg_tau23, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau33,
                                 amrex::Gpu::HostVector<amrex::Real>& h_avg_hfx3,  amrex::Gpu::HostVector<amrex::Real>& h_avg_q1fx3,
                                 amrex::Gpu::HostVector<amrex::Real>& h_avg_q2fx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_diss);
    void derive_stress_profiles_stag (amrex::Gpu::HostVector<amrex::Real>& h_avg_tau11, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau12,
                                      amrex::Gpu::HostVector<amrex::Real>& h_avg_tau13, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau22,
                                      amrex::Gpu::HostVector<amrex::Real>& h_avg_tau23, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau33,
                                      amrex::Gpu::HostVector<amrex::Real>& h_avg_hfx3,  amrex::Gpu::HostVector<amrex::Real>& h_avg_q1fx3,
                                      amrex::Gpu::HostVector<amrex::Real>& h_avg_q2fx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_diss);
 
    void derive_forcing_profiles_stag (amrex::Gpu::HostVector<amrex::Real>& h_avg_ttend,  amrex::Gpu::HostVector<amrex::Real>& h_avg_qtend,
                                       amrex::Gpu::HostVector<amrex::Real>& h_avg_wsub,   amrex::Gpu::HostVector<amrex::Real>& h_avg_thtend,
                                       amrex::Gpu::HostVector<amrex::Real>& h_avg_qhtend, amrex::Gpu::HostVector<amrex::Real>& h_avg_tvtend,
                                       amrex::Gpu::HostVector<amrex::Real>& h_avg_qvtend, amrex::Gpu::HostVector<amrex::Real>& h_avg_qcvtend,
                                       amrex::Gpu::HostVector<amrex::Real>& h_avg_tnudge, amrex::Gpu::HostVector<amrex::Real>& h_avg_qnudge,
                                       amrex::Gpu::HostVector<amrex::Real>& h_avg_unudge, amrex::Gpu::HostVector<amrex::Real>& h_avg_vnudge);
 
    // Perform the volume-weighted sum
    amrex::Real
    volWgtSumMF (int lev, const amrex::MultiFab& mf, int comp,
                 const amrex::MultiFab& dJ, const amrex::MultiFab& mfx, const amrex::MultiFab& mfy,
                 bool finemask, bool local = true);
 
    void
    volWgtColumnSum (int lev, const amrex::MultiFab& mf, int comp,
                     amrex::MultiFab& mf_2d, const amrex::MultiFab& dJ);
 
    // Decide if it is time to take an action
    static bool is_it_time_for_action (int nstep, double time, double dt,
                                       int action_interval, double action_per);
 
    // Make a new level using provided BoxArray and DistributionMapping and
    // fill with interpolated coarse level data.
    // overrides the pure virtual function in AmrCore
    void MakeNewLevelFromCoarse (int lev, amrex::Real time, const amrex::BoxArray& ba,
                                 const amrex::DistributionMapping& dm) override;
 
    // Remake an existing level using provided BoxArray and DistributionMapping and
    // fill with existing fine and coarse data.
    // overrides the pure virtual function in AmrCore
    void RemakeLevel (int lev, amrex::Real time, const amrex::BoxArray& ba,
                      const amrex::DistributionMapping& dm) override;
 
    // Delete level data
    // overrides the pure virtual function in AmrCore
    void ClearLevel (int lev) override;
 
    // Make a new level from scratch using provided BoxArray and DistributionMapping.
    // Only used during initialization.
    // overrides the pure virtual function in AmrCore
    void MakeNewLevelFromScratch (int lev, amrex::Real time, const amrex::BoxArray& ba,
                                  const amrex::DistributionMapping& dm) override;
 
    // compute dt from CFL considerations
    double estTimeStep (int lev, long& dt_fast_ratio) const;
 
    // largest explicit eddy-diffusion rates (1/s) on level lev, momentum and scalars, for the
    // diffusive time-step check (erf.diffusive_dt_check / erf.diffusive_dt_limit)
    void ComputeDiffusiveRates (int lev, amrex::Real& rate_mom, amrex::Real& rate_scal) const;
 
    // print the largest terrain slope factor alpha = h dx/dz on level lev: always when force is
    // true (new grids), otherwise only if it changed by more than 1 % since the last report
    void ReportTerrainSlopeFactor (int lev, bool force);
 
#ifdef ERF_USE_WW3_COUPLING
    //amrex::Print() <<  " About to call send_to_ww3 from ERF.H" << std::endl;
    void send_to_ww3 (int lev);
    //amrex::Print() <<  " About to call read_waves from ERF.H" << std::endl;
    void read_waves (int lev);
    //void send_to_ww3 (int lev);
    //void read_waves (int lev);
 
    //void send_to_ww3 (int lev);
#endif
 
    // Interface for advancing the data at one level by one "slow" timestep
    void advance_dycore (int level,
                         amrex::Vector<amrex::MultiFab>& state_old,
                         amrex::Vector<amrex::MultiFab>& state_new,
                         amrex::MultiFab& xvel_old,  amrex::MultiFab& yvel_old,  amrex::MultiFab& zvel_old,
                         amrex::MultiFab& xvel_new,  amrex::MultiFab& yvel_new,  amrex::MultiFab& zvel_new,
                         amrex::MultiFab& source,    amrex::MultiFab& xmom_src,
                         amrex::MultiFab& ymom_src,  amrex::MultiFab& zmom_src,
                         amrex::MultiFab& buoyancy,  amrex::Geometry fine_geom,
                         double dt, double time);
 
    void advance_microphysics (int lev,
                               amrex::MultiFab& cons_in,
                               const double& dt_advance,
                               const int& iteration,
                               const double& time);
 
    void advance_lsm (int lev,
                      amrex::MultiFab& cons_in,
                      amrex::MultiFab& xvel_in,
                      amrex::MultiFab& yvel_in,
                      const double& time,
                      const double& dt_advance);
 
    void advance_urban (int lev,
                        amrex::MultiFab& cons_in,
                        amrex::MultiFab& u_in,
                        amrex::MultiFab& v_in,
                        const amrex::Real& dt_advance,
                        const amrex::Geometry& geom,
                        const amrex::MultiFab* z_phys_nd,
                        amrex::MultiFab& eddyDiffs_in);
 
    void advance_radiation (int lev,
                            amrex::MultiFab& cons_in,
                            const double& dt_advance);
 
#ifdef ERF_USE_EAMXX_SHOC
    void compute_shoc_tendencies (int lev,
                                  amrex::MultiFab* cons,
                                  amrex::MultiFab* xvel,
                                  amrex::MultiFab* yvel,
                                  amrex::MultiFab* zvel,
                                  amrex::Real* w_subsid,
                                  amrex::MultiFab* tau13,
                                  amrex::MultiFab* tau23,
                                  amrex::MultiFab* hfx3,
                                  amrex::MultiFab* qfx3,
                                  amrex::MultiFab* eddyDiffs,
                                  amrex::MultiFab* z_phys_nd,
                                  const double& dt_advance);
#endif
 
    void compute_native_shoc_tendencies (int lev,
                                         amrex::MultiFab* cons,
                                         amrex::MultiFab* xvel,
                                         amrex::MultiFab* yvel,
                                         amrex::MultiFab* zvel,
                                         amrex::Real* w_subsid,
                                         amrex::MultiFab* tau13,
                                         amrex::MultiFab* tau23,
                                         amrex::MultiFab* hfx3,
                                         amrex::MultiFab* qfx3,
                                         amrex::MultiFab* eddyDiffs,
                                         amrex::MultiFab* z_phys_nd,
                                         const double& dt_advance);
 
    void build_fine_mask (int lev, amrex::MultiFab& fine_mask);
 
    void MakeHorizontalAverages ();
    void MakeDiagnosticAverage (amrex::Vector<amrex::Real>& h_havg, amrex::MultiFab& S, int n);
    void derive_upwp (amrex::Vector<amrex::Real>& h_havg);
 
    // Write plotfile to disk
    bool Write3DPlotFile  (int which, PlotFileType plotfile_type, amrex::Vector<std::string> plot_var_names);
    void Write2DPlotFile  (int which, PlotFileType plotfile_type, amrex::Vector<std::string> plot_var_names);
    void ResetIntervalMeansAfter3DPlotfileBatch (int plotfiles_written,
                                                 bool interval_diagnostic_consumed);
    void ResetIntervalMeans ();
 
    // Write subvolume of data into a different "plotfile"
    void WriteSubvolume (int isub, amrex::Vector<std::string> subvol_var_names);
 
    void WriteMultiLevelPlotfileWithTerrain (const std::string &plotfilename,
                                             int nlevels,
                                             const amrex::Vector<const amrex::MultiFab*> &mf,
                                             const amrex::Vector<const amrex::MultiFab*> &mf_nd,
                                             const amrex::Vector<std::string> &varnames,
                                             const amrex::Vector<amrex::Geometry>& my_geom,
                                             double time,
                                             const amrex::Vector<int> &level_steps,
                                             const amrex::Vector<amrex::IntVect>& my_ref_ratio,
                                             const std::string &versionName = "HyperCLaw-V1.1",
                                             const std::string &levelPrefix = "Level_",
                                             const std::string &mfPrefix = "Cell",
                                             const amrex::Vector<std::string>& extra_dirs = amrex::Vector<std::string>()) const;
 
 
    void WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile,
                                                int nlevels,
                                                const amrex::Vector<amrex::BoxArray> &bArray,
                                                const amrex::Vector<std::string> &varnames,
                                                const amrex::Vector<amrex::Geometry>& my_geom,
                                                double time,
                                                const amrex::Vector<int> &level_steps,
                                                const amrex::Vector<amrex::IntVect>& my_ref_ratio,
                                                const std::string &versionName,
                                                const std::string &levelPrefix,
                                                const std::string &mfPrefix) const;
 
    // Filling the plot variables, and the station sampling that reuses it.  Split
    // out of the plotfile writers so that a station column and the corresponding
    // plotfile component are produced the same way.  Private for the same reason
    // canonicalizePlot3DVariables above is: this is output machinery rather than
    // something ERF offers a caller.
 
    // Per-level scratch needed to evaluate the 3D plot variables.  Held by the
    // caller so that the work is done once for a set of variables and then
    // reused for every level that is filled.
    struct Plot3DScratch
    {
        amrex::Vector<amrex::MultiFab> mf_cc_vel;   // cell-centered velocity, 1 ghost
        amrex::Vector<amrex::MultiFab> mf_cc_tau;   // cell-centered stress tensor
        amrex::Vector<amrex::MultiFab> mf_cc_fx;    // cell-centered SFS fluxes
    };
 
    void FillPlot2DVars (int lev,
                         const amrex::Vector<std::string>& plot_var_names,
                         const amrex::Vector<std::string>& varnames,
                         const amrex::Vector<plotfile2d::Plotfile2DOutputDescriptor>& output_descriptors,
                         amrex::MultiFab& mf_dst,
                         int ncomp_mf);
 
    void FillPlot3DVars (int lev,
                         const amrex::Vector<std::string>& plot_var_names,
                         Plot3DScratch& scratch,
                         amrex::MultiFab& mf_dst,
                         int ncomp_mf,
                         amrex::Real tnew);
 
    // max_lev caps the levels whose scratch is built, -1 meaning every level: a
    // caller that reads only the coarse levels should not pay for the fine ones.
    // sync_solution allows the microphysics average-down, which the plotfile
    // writer wants but a diagnostic that runs at every step must not do, since
    // it changes the coarse solution.
    void BuildPlot3DScratch (const amrex::Vector<std::string>& plot_var_names,
                             Plot3DScratch& scratch,
                             int max_lev = -1,
                             bool sync_solution = true);
 
    // Station time-series output
    void init_stations ();
    void resolve_station_positions ();
    void resolve_station_stencils ();
    void sample_stations (amrex::Real time);
    void flush_stations () const;
    // The nodes a station's terrain elevation is read from, at k = the bottom of
    // the domain: z_phys_nd, or with immersed-forcing terrain the immersed surface
    const amrex::MultiFab& station_ground (int lev);
 
    // Nudging towards observations at stations
    void init_obs_nudging ();
 
    void erf_enforce_hse (int lev,
                          amrex::MultiFab& dens, amrex::MultiFab& pres, amrex::MultiFab& pi,
                          amrex::MultiFab& th, amrex::MultiFab& qv,
                          std::unique_ptr<amrex::MultiFab>& z_cc);
 
#ifdef ERF_USE_NETCDF
    //! Write a timestep to 1D vertical column output for coupling
    void writeToNCColumnFile (int lev,
                              const std::string& colfile_name, amrex::Real xloc, amrex::Real yloc,
                              double time);
#endif //ERF_USE_NETCDF
 
    void init_from_input_sounding (int lev);
 
    /**
     * The input sounding sampled at the cell centres of level lev. Level 0 is
     * input_sounding_data itself; a finer level gets its own copy, resampled at
     * geom[lev] / zlevels_stag[lev], so that its profile index k is its own cell k
     * and no fine cell blends in the surface line of the sounding.
     */
    const InputSoundingData& sounding_at_level (int lev);
 
    /**
     * The large-scale forcing profiles interpolated to the cell centres of level lev.
     * Level 0 is lsf itself; a finer level gets its own copy, interpolated at
     * geom[lev] / zlevels_stag[lev] against sounding_at_level(lev), so that its profile
     * index k is its own cell k.
     */
    const LargeScaleForcingData& lsf_at_level (int lev);
 
    void init_immersed_forcing (int lev);
 
    // Immersed-boundary surface energy balance (Source/ImmersedBoundarySEB)
    void init_ibseb ();
    void ibseb_check_refined_levels () const;
    void ibseb_check_sun_matches_two_stream () const;
    void ibseb_check_two_stream_provider () const;
    [[nodiscard]] bool ibseb_after_radiation () const {
        return ibseb_params.enable && ibseb_params.radiation == "two_stream";
    }
    void ibseb_set_two_stream_sun (int lev, amrex::Real time);
    void ibseb_write_checkpoint (const std::string& checkpointname, int lev) const;
    void ibseb_report (int nstep, amrex::Real time);
    amrex::Real ibseb_bulk_richardson_height (int lev, const amrex::MultiFab& cons,
                                              const amrex::MultiFab& xvel, const amrex::MultiFab& yvel);
    void ibseb_advance (int lev, amrex::Real time, amrex::Real dt_lev, const amrex::MultiFab& cons,
                        const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& zvel);
 
    void input_sponge (int lev);
 
    void init_from_hse (int lev);
 
    void init_thin_body (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
 
    void FillForecastStateMultiFabs (const int lev,
                                     const std::string& filename,
                                     const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                                     amrex::Vector<amrex::Vector<amrex::MultiFab>>& forecast_state);
 
    void FillSurfaceStateMultiFabs (const int lev,
                                    const std::string& filename,
                                    amrex::Vector<amrex::MultiFab>& surface_state);
 
 
    void WeatherDataInterpolation (const int nlevs,
                                   const double time,
                                   amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd,
                                   bool regrid_forces_file_read);
 
    void SurfaceDataInterpolation (const int nlevs,
                                   const double time,
                                   amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd,
                                   bool regrid_forces_file_read);
 
    void create_background_state_for_ensemble (int lev,
                                               amrex::MultiFab& mf_cc_pert,
                                               amrex::MultiFab& cons_pert,
                                               amrex::MultiFab& xvel_pert,
                                               amrex::MultiFab& yvel_pert,
                                               amrex::MultiFab& zvel_pert);
 
    void create_random_perturbations (const int lev,
                                      amrex::MultiFab& mf_cc_pert);
 
    void apply_gaussian_smoothing_to_perturbations (const int lev,
                                                    amrex::MultiFab& mf_cc_pert);
 
    void ComputeAndWriteEnsemblePerturbations (const int da_iter);
 
    void PerformDataAssimilation (int da_iter);
 
    void SetDirsForPlotfilesAndCheckpointsForDA (const int da_iter, const int ens_no);
 
    void GetEnsembleCheckpointName (int da_iter, int ens_no);
 
    void init_custom (int lev);
 
#ifdef ERF_USE_MULTIBLOCK
    // constructor used when ERF is created by a multiblock driver
    // calls AmrCore constructor -> AmrMesh constructor
    ERF (const amrex::RealBox& rb, int max_level_in,
         const amrex::Vector<int>& n_cell_in, int coord,
         const amrex::Vector<amrex::IntVect>& ref_ratio,
         const amrex::Array<int,AMREX_SPACEDIM>& is_per,
         std::string prefix);
 
    // Advance a block specified number of time steps
    void Evolve_MB (int MBstep, int max_block_step);
 
    // Set parmparse prefix for MultiBlock
    void SetParmParsePrefix (std::string name) { pp_prefix = name; }
 
    // Set 'this' from multiblock container
    void SetMultiBlockPointer (MultiBlockContainer *mbc) { m_mbc = mbc; }
 
    // Public data copy for MB
    std::vector<amrex::Box> domain_p;
 
    // Primary solution data data containers
    amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_new;
    amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_old;
 
    // Pressure gradient at each level
    amrex::Vector<amrex::Vector<amrex::MultiFab> > gradp;
 
    // Velocity time averaged field
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> vel_t_avg;
    amrex::Vector<double> t_avg_cnt;
 
    // Interval first and second moments used by mean/fluctuation plot fields
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> interval_means;
    amrex::Vector<double> t_mean_cnt;
    // The configured time reset is one global, one-shot event.  Keep this
    // separate from the per-level accumulation storage so regridding a level
    // cannot re-arm an event that has already happened.
    int mean_vars_time_reset_done{0};
#endif
 
    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);
 
#ifdef ERF_USE_NETCDF
    // Version in which we always set the velocity values but don't set theta
    void fill_from_realbdy (const amrex::Vector<amrex::MultiFab*>& mfs,
                            double time,
                            bool cons_only,
                            int icomp_cons,
                            int ncomp_cons,
                            amrex::IntVect ngvect_cons,
                            amrex::IntVect ngvect_vels);
#endif
 
#ifdef ERF_USE_NETCDF
    // Read a wrfinput file at this level.
    //
    // With read_atmos_state = true (the default) this reads the full state: the
    // atmospheric fields (U, V, W, THM, the density triple and the moisture species)
    // as well as the surface fields, and it remaps the state onto the ERF vertical
    // grid, builds the base state and reads the boundary data.
    //
    // With read_atmos_state = false only the surface fields are read -- terrain,
    // surface pressure, surface and skin temperatures, the land/soil masks, the map
    // factors, the lat/lon fields and the LSM variables -- and the terrain and metric
    // terms are built from them.  The atmospheric state is left untouched for the
    // caller to fill by interpolation from the coarse level (see the
    // erf.interp_atmos_from_coarse path in MakeNewLevelFromCoarse), which is also
    // responsible for calling rebuild_base_state_from_wrfinput afterwards.
    void init_from_wrfinput (int lev, amrex::MultiFab& mf_PSFC, bool read_atmos_state = true);
    void init_from_metgrid  (int lev);
    void init_from_ncfile   (int lev);
 
    // Rebuild the WRF base state on the grids of base_state_mf at this level, using the
    // level-0 reference parameters (wrf_bsp) and this level's z_phys_cc.  This is the same
    // construction init_from_wrfinput performs at initialization, factored out so that a
    // level remade by a regrid gets the same base state rather than one interpolated from
    // its parent.  Requires z_phys_cc[lev] and wrf_bsp to be valid.
    void rebuild_base_state_from_wrfinput (int lev, amrex::MultiFab& base_state_mf);
 
    // Re-read the six reference parameters that define the WRF base state from the level-0
    // wrfinput file.  These are read at initialization but are not part of the checkpoint,
    // so on restart wrf_bsp would otherwise be left unset -- and the first regrid that has
    // to rebuild a base state would abort in init_base_state_from_wrfinput.
    void restore_base_state_params_on_restart ();
#endif // ERF_USE_NETCDF
 
#ifdef ERF_USE_WINDFARM
    void init_windfarm (int lev);
    void advance_windfarm (const amrex::Geometry& a_geom,
                           const double& dt_advance,
                           amrex::MultiFab& cons_in,
                           amrex::MultiFab& U_old,
                           amrex::MultiFab& V_old,
                           amrex::MultiFab& W_old,
                           amrex::MultiFab& mf_vars_windfarm,
                           const amrex::MultiFab& mf_Nturb,
                           const amrex::MultiFab& mf_SMark,
                           const double& time);
#endif
 
    void MakeEBGeometry ();
    void make_eb_box ();
    void make_eb_regular ();
 
    // more flexible version of AverageDown() that lets you average down across multiple levels
    void AverageDownTo (int crse_lev,
                        int scomp,
                        int ncomp,
                        bool do_perturbational_and_momenta = true); // NOLINT
 
    /*! \brief Conservation-preserving fine-to-coarse average of RhoTheta and
     *  the active moist conserved components using detJ/mfac weighting, reusing
     *  AverageDownTo with perturbational-state and momentum handling disabled. */
    void AverageDownMoistStateTo (int crse_lev);
 
    // Write checkpoint file to disk
    void WriteCheckpointFile () const;
 
    // Read checkpoint file from disk
    void ReadCheckpointFile ();
 
    // Read only velocities from checkpoint file
    void ReadVelsOnlyFromCheckpointFile (int lev_to_fill, std::string& chkfile);
 
    // read checkpoint file from disk -- called after instantiating m_SurfaceLayer
    void ReadCheckpointFileSurfaceLayer ();
 
    void fill_terrain_surface  (int lev, amrex::MultiFab& z_phys, double elapsed_time);
    void init_zphys            (int lev, double elapsed_time);
    void remake_zphys          (int lev, amrex::Real time, std::unique_ptr<amrex::MultiFab>& temp_zphys_nd);
    void update_terrain_arrays (int lev);
    void make_terrain_blanking (int lev);
 
    bool m_driver_has_atm2ocn_coupling = false;
    bool m_driver_uses_two_way_coupling = false;
    bool m_driver_atm2ocn_uses_state_contract = true;
 
    // Coupled sea-surface temperature and its per-cell coverage flag, on the
    // layout reported by GetOceanToAtmosSurfaceLayout. Allocated on the first
    // ApplyOceanSurfaceState call and handed to the SurfaceLayer as borrowed
    // pointers, so ERF must outlive the SurfaceLayer's use of them. Level 0
    // only: the coupler exchanges on the base grid.
    std::unique_ptr<amrex::MultiFab>  m_coupled_sst;
    std::unique_ptr<amrex::iMultiFab> m_coupled_sst_valid;
 
    // One-shot latch so a misconfigured deck says so once rather than every step.
    bool m_warned_coupled_sst_declined = false;
    // The coupled-SST diagnostics run on the first ApplyOceanSurfaceState call,
    // then only at erf.v >= 2.
    bool m_checked_coupled_sst = false;
    // Out-of-range report: the worst extremes and NaN count already reported;
    // at erf.v >= 2 it repeats only when a later call is worse. Start at the bounds.
    amrex::Real m_warned_coupled_sst_min = amrex::Real(260.0);
    amrex::Real m_warned_coupled_sst_max = amrex::Real(320.0);
    amrex::Long m_warned_coupled_sst_nan = 0;
 
private:
 
    ///////////////////////////
    // private member functions
    ///////////////////////////
 
    // read in some parameters from inputs file
    void ReadParameters ();
    void ParameterSanityChecks ();
 
    [[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 ();
 
    void update_diffusive_arrays (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
 
    void make_flux_register (int lev);
 
    void Construct_ERFFillPatchers (int lev);
 
    void Define_ERFFillPatchers (int lev);
 
    void init1DArrays ();
 
    void init_bcs ();
    void init_phys_bcs (bool& read_prim_theta);
 
    void init_stuff (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
                     amrex::Vector<amrex::MultiFab>& lev_new, amrex::Vector<amrex::MultiFab>& lev_old,
                     amrex::MultiFab& tmp_base_state,
                     std::unique_ptr<amrex::MultiFab>& tmp_zphys_nd);
 
    // Initialize the Turbulent perturbation
    void turbPert_update (const int lev, const double dt);
    void turbPert_amplitude (const int lev);
    void turbPert_amplitude_w (const int lev);
 
    // Initialize the integrator object
    void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf);
 
    // Create the physbcs objects
    void make_physbcs (int lev);
 
    // (Re)build the land surface model on this level's current grids
    void make_lsm_at_level (int lev, bool from_regrid = false,
                            int source_lev = -1, bool source_is_raw_input = false,
                            bool initialize_now = false);
 
    // Initialize the microphysics object
    void initializeMicrophysics (const int&);
 
#ifdef ERF_USE_WINDFARM
    // Initialize the windfarm object
    void initializeWindFarm (const int&);
#endif
 
    // Initialize embedded boundary geometries
    void initializeEB ();
 
    // Compute a vector of new MultiFabs by copying from valid region and filling ghost cells
    //
    // NOTE: FillPatch takes in an empty MF, and returns cell-centered + velocities (not momenta)
    //
    // This one works only at level = 0 (base state does not change)
    void FillPatchCrseLevel (int lev, double time,
                             const amrex::Vector<amrex::MultiFab*>& mfs_vel,
                             bool cons_only=false);
 
    // This one works only at level > 0 (base state does change)
    void FillPatchFineLevel (int lev, double time,
                             const amrex::Vector<amrex::MultiFab*>& mfs_vel,
                             const amrex::Vector<amrex::MultiFab*>& mfs_mom,
                             const amrex::MultiFab& old_base_state,
                             const amrex::MultiFab& new_base_state,
                             bool fillset=true, bool cons_only=false);
 
    // Compute new multifabs by copying data from valid region and filling ghost cells.
    // Unlike FillPatch, FillIntermediatePatch will use the supplied multifabs instead of fine level data.
    // This is to support filling boundary cells at an intermediate time between old/new times
    // on the fine level when valid data at a specific time is already available (such as
    // at each RK stage when integrating between initial and final times at a given level).
    //
    // NOTE: FillIntermediatePatch takes in updated momenta, and returns both updated velocity and momenta
    //
    void FillIntermediatePatch (int lev, double time,
                                const amrex::Vector<amrex::MultiFab*>& mfs_vel,
                                const amrex::Vector<amrex::MultiFab*>& mfs_mom,
                                int ng_cons, int ng_vel, bool cons_only, int icomp_cons, int ncomp_cons);
 
    // Fill all multifabs (and all components) in a vector of multifabs corresponding to the
    // grid variables defined in vars_old and vars_new just as FillCoarsePatch.
    void FillCoarsePatch (int lev, double time);
 
    // advance a level by dt
    // includes a recursive call for finer levels
    void timeStep (int lev, double time, int iteration);
 
    // advance a single level for a single time step
    void Advance (int lev, double time, double dt_lev, int iteration, int ncycle);
 
    //! Initialize HSE
    void initHSE ();
    void initHSE (int lev);
 
    // Fill base_state[lev] by conservative interpolation from base_state[lev-1], then apply
    // the base-state physical bc's.  This is how the part of a fine level that lies inside
    // the domain but outside the fine grids -- and the fine ghost cells -- get their values,
    // so it must run before any construction that overwrites only the fine grids.
    void interp_base_state_from_coarse (int lev);
 
    //! Initialize Rayleigh damping profiles at a level
    void initRayleigh_at_level (const int& lev);
 
    //! Initialize sponge profiles
    void initSponge ();
 
    //! Set Rayleigh mean profiles from input sounding
    void setRayleighRefFromSounding (bool restarting);
 
    //! Set sponge mean profiles from input sounding
    void setSpongeRefFromSounding (bool restarting);
 
    // a wrapper for estTimeStep()
    void ComputeDt (int step = -1, double cur_time_d = 0.0);
 
    // get plotfile name
    [[nodiscard]] std::string PlotFileName (int lev) const;
 
    // set plotfile variable names
    static amrex::Vector<std::string> PlotFileVarNames (amrex::Vector<std::string> plot_var_names);
 
    // set which variables and derived quantities go into plotfiles
    void setPlotVariables (const std::string& pp_plot_var_names, amrex::Vector<std::string>& plot_var_names);
    void canonicalizePlot3DVariables (amrex::Vector<std::string>& plot_var_names);
    void appendPlot3DVariables (const amrex::Vector<std::string>& plot_var_names,
                                amrex::Vector<std::string>& a_plot_var_names);
    void setPlotVariables2D (const std::string& pp_plot_var_names, amrex::Vector<std::string>& plot_var_names);
    // append variables to plot
    void appendPlotVariables (const std::string& pp_plot_var_names, amrex::Vector<std::string>& plot_var_names);
 
    void setSubVolVariables (const std::string& pp_subvol_var_names, amrex::Vector<std::string>& subvol_var_names);
 
#ifdef ERF_USE_NETCDF
    //! Create 1D vertical column output for coupling
    void createNCColumnFile (int lev,
                             const std::string& colfile_name, amrex::Real xloc, amrex::Real yloc);
 
    // Copy from the NC*fabs into the MultiFabs holding the boundary data
    void init_from_wrfbdy (amrex::Vector<amrex::FArrayBox*> x_vel_lateral,
                           amrex::Vector<amrex::FArrayBox*> y_vel_lateral,
                           amrex::Vector<amrex::FArrayBox*> z_vel_lateral,
                           amrex::Vector<amrex::FArrayBox*> T_lateral);
 
    static double start_bdy_time;
    static double final_bdy_time;
 
    static double start_low_time;
    static double final_low_time;
 
    static double bdy_time_interval;
    static double low_time_interval;
 
    // *** *** FArrayBox's for holding the SURFACE data
    // amrex::IArrayBox NC_IVGTYP_fab; // Vegetation type (IVGTYP); Discrete numbers;
    // amrex::FArrayBox NC_z0_fab;     // Surface Roughness, z0 = z0 (IVGTYP)
    // amrex::FArrayBox NC_PSFC_fab;   // Surface pressure
 
    // TODO: Clarify the relation between SST and TSK
    // amrex::FArrayBox NC_SST_fab;    // Sea Surface Temperature; Defined even for land area
    // amrex::FArrayBox NC_TSK_fab;    // Surface Skin Temperature; Appears to be same as SST...
 
    // Vectors (over time) of Vector (over variables) of FArrayBoxs for holding the data read from the wrfbdy NetCDF file
    amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_xlo;
    amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_xhi;
    amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_ylo;
    amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_yhi;
 
    amrex::Vector<amrex::Vector<amrex::FArrayBox>> low_data_zlo;
 
    // Maximum value of terrain read in from wrfinput at level 0
    amrex::Real z_top = zero;
 
#endif // ERF_USE_NETCDF
 
    // Latitude and Longitude on the grid, always at cell centers: the wrfinput
    // path reads XLAT / XLONG and the metgrid path reads XLAT_M / XLONG_M, both
    // of which are mass-point values, so a consumer that wants the mass point
    // can use these directly.
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> lat_m, lon_m;
 
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> sinPhi_m, cosPhi_m; // Coriolis factors on the grid
 
    // Struct for working with the sounding data we take as an input. Its profiles are
    // sampled at the level-0 cell centres; use sounding_at_level for any other level.
    InputSoundingData input_sounding_data;
 
    // The sounding sampled at the cell centres of each level > 0, built on first use
    // from input_sounding_data (index 0 is unused)
    amrex::Vector<std::unique_ptr<InputSoundingData>> input_sounding_data_fine;
 
    // Struct for working with the sponge data we take as an input
    InputSpongeData input_sponge_data;
 
    // Vector (6 planes) of DeviceVectors (ncell in plane) for Dirichlet BC data
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> xvel_bc_data;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> yvel_bc_data;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> zvel_bc_data;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> th_bc_data;
 
    // Planar averages for immersed forcing (computed on slow step, used in slow and fast steps)
    amrex::Vector<amrex::TableData<amrex::Real, 1>> r_plane_avg;   // density planar average per level
    amrex::Vector<amrex::TableData<amrex::Real, 1>> t_plane_avg;   // RhoTheta planar average per level
 
    // Function to read and populate above vectors (if input file exists)
    bool init_Dirichlet_bc_data (const std::string input_file);
 
    // Read the file passed to amr.restart and use it as an initial condition for
    // the current simulation. Supports a different number of components and
    // ghost cells.
    void InitializeFromFile ();
 
    // Initialize the new-time data at a level from the initial_data MultiFab
    void InitializeLevelFromData (int lev, const amrex::MultiFab& initial_data);
 
    // utility to skip to next line in Header
    static void GotoNextLine (std::istream& is);
 
    // Single level functions called by advance()
    void post_update (amrex::MultiFab& state_mf, double time, const amrex::Geometry& geom);
    void fill_rhs (amrex::MultiFab& rhs_mf, const amrex::MultiFab& state_mf, double time, const amrex::Geometry& geom);
 
    ////////////////
    // private data members
 
    std::unique_ptr<ProblemBase> prob = nullptr;  // problem-specific functions
 
    amrex::Vector<int> num_boxes_at_level;                   // how many boxes specified at each level by tagging criteria
    amrex::Vector<int> num_files_at_level;                   // how many wrfinput files specified at each level
    amrex::Vector<amrex::Vector<amrex::Box>> boxes_at_level; //      the boxes specified at each level by tagging criteria
 
    amrex::Vector<int> istep;      // which step?
    amrex::Vector<int> nsubsteps;  // how many substeps on each level?
 
    // One entry per level: set when that level has just been built by the
    // erf.interp_atmos_from_coarse path, whose atmospheric state comes from FillCoarsePatch
    // and is not yet thermodynamically consistent (rho, theta and qv are interpolated
    // independently).  advance_radiation interpolates the radiation fields from the parent
    // rather than running RRTMGP on that state, and clears the entry once it has done so, so
    // exactly one step is skipped per level creation.
    //
    // NOTE: this is deliberately a flag owned by the level's lifetime rather than a test on
    //       istep[lev].  istep is sized once in ReadParameters and is never reset when a
    //       level is removed, so it carries over from a level's previous life: a level that
    //       ran for fifty steps, was de-refined, and was then re-created would have
    //       istep[lev] == 50 and would silently miss the skip.
    //
    // (using int instead of bool to avoid the std::vector<bool> specialization)
    amrex::Vector<int> rad_interp_from_coarse_pending;
 
    // keep track of old time, new time, and time step at each level
    amrex::Vector<double> t_new;
    amrex::Vector<double> t_old;
    amrex::Vector<double> dt;
    amrex::Vector<long> dt_mri_ratio;
 
#ifndef ERF_USE_MULTIBLOCK
    // Array of multifabs to store the solution at each level of refinement;
    // after advancing a level we use "swap".
    amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_new;
    amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_old;
 
    // Pressure gradient at each level
    amrex::Vector<amrex::Vector<amrex::MultiFab> > gradp;
 
    // Velocity time averaged field
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> vel_t_avg;
    amrex::Vector<double> t_avg_cnt;
 
    // Interval first and second moments used by mean/fluctuation plot fields
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> interval_means;
    amrex::Vector<double> t_mean_cnt;
    // Global one-shot time-reset state; interval_means remains per level.
    int mean_vars_time_reset_done{0};
#endif
    amrex::Vector<std::unique_ptr<MRISplitIntegrator<amrex::Vector<amrex::MultiFab> > > > mri_integrator_mem;
 
    // Used for anelastic or when we impose an initial projection
    amrex::Vector<amrex::MultiFab> pp_inc;
 
    // Used only for fast substepping
    amrex::Vector<amrex::MultiFab> lagged_delta_rt;
    amrex::Vector<amrex::MultiFab> avg_xmom;
    amrex::Vector<amrex::MultiFab> avg_ymom;
    amrex::Vector<amrex::MultiFab> avg_zmom;
 
    // Vector over levels of routines to impose physical boundary conditions
    amrex::Vector<std::unique_ptr<ERFPhysBCFunct_cons>> physbcs_cons;
    amrex::Vector<std::unique_ptr<ERFPhysBCFunct_u>> physbcs_u;
    amrex::Vector<std::unique_ptr<ERFPhysBCFunct_v>> physbcs_v;
    amrex::Vector<std::unique_ptr<ERFPhysBCFunct_w>> physbcs_w;
    amrex::Vector<std::unique_ptr<ERFPhysBCFunct_base>> physbcs_base;
 
    // Store primitive variable for MOST BC
    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;
 
    // Scratch space for time integrator
    amrex::Vector<amrex::MultiFab> rU_old;
    amrex::Vector<amrex::MultiFab> rU_new;
    amrex::Vector<amrex::MultiFab> rV_old;
    amrex::Vector<amrex::MultiFab> rV_new;
    amrex::Vector<amrex::MultiFab> rW_old;
    amrex::Vector<amrex::MultiFab> rW_new;
 
    // amrex::Vector<amrex::MultiFab> xmom_crse_rhs;
    // amrex::Vector<amrex::MultiFab> ymom_crse_rhs;
    amrex::Vector<amrex::MultiFab> zmom_crse_rhs;
 
    std::unique_ptr<Microphysics> micro;
    std::unique_ptr<erf_sbm::SBMStateManager> sbm_state_manager;
    // Declared after the state manager so its non-owning layout reference is
    // destroyed before the authoritative SBM owner.
    std::unique_ptr<erf_sbm::SBMTransport> sbm_transport;
    std::unique_ptr<erf_auxiliary::AuxiliaryInertTracer> auxiliary_inert_tracer;
    std::unique_ptr<CloudChamberBudget> cloud_chamber_budget;
    erf_cloud_chamber::Config cloud_chamber_config;
    amrex::Vector<amrex::Vector<amrex::MultiFab*>> qmoist; // (lev,ncomp) rain_accum, snow_accum, graup_accum
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> precip; // SLM precipitation-rate scratch
 
    // Variables for wind farm parametrization models
 
#ifdef ERF_USE_WINDFARM
    std::unique_ptr<WindFarm> windfarm;
    amrex::Vector<amrex::MultiFab> Nturb;
    amrex::Vector<amrex::MultiFab> vars_windfarm; // Fitch: Vabs, Vabsdt, dudt, dvdt, dTKEdt
                                                  // EWP:                 dudt, dvdt, dTKEdt
 
    amrex::Vector<amrex::MultiFab> SMark; // A multifab that holds an integer corresponding
                                                   // to the number of the wind turbine to sample
                                                   // velocity upstream of the turbine
#endif
 
    std::unique_ptr<SurfaceModel> m_SurfaceModel = nullptr;
 
    LandSurface lsm;
    Urban urban;
    amrex::Vector<std::string> lsm_data_name;                // (ncomp) -- e.g., sw_dn_dir/dif_vis/nir, lw_dn, cos_zen
    amrex::Vector<amrex::Vector<amrex::MultiFab*>> lsm_data; // (lev,ncomp) Data for LSM coupling
    amrex::Vector<std::string> lsm_flux_name;                // (ncomp) -- e.g., hfx3, q1fx3, tau13, tau23
    amrex::Vector<amrex::Vector<amrex::MultiFab*>> lsm_flux; // (lev,ncomp) Fluxes from LSM needed by DYCORE
    amrex::Vector<amrex::Vector<amrex::MultiFab*>> urban_data; // (lev,ncomp) Data for urban coupling
    amrex::Vector<std::string> urban_data_name;                // (ncomp) Urban coupling field names
    amrex::Vector<amrex::Vector<amrex::MultiFab*>> urban_flux; // (lev,ncomp) Fluxes from urban model
 
    LargeScaleForcingData lsf;
 
    // The large-scale forcing profiles on each level > 0, built on first use from lsf
    // (index 0 is unused)
    amrex::Vector<std::unique_ptr<LargeScaleForcingData>> lsf_fine;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> nudge_data; // holds nudging inputs: tnudge, qnudge, unudge, vnudge
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> lsf_data; //  holds large scale forcing tendencies: ttend, qtend, wsub
 
    amrex::Vector<std::unique_ptr<IRadiation>> rad; // Radiation model at each level
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> qheating_rates;  // radiation heating rate source terms (SW, LW) (DYCORE)
    // Radiation fluxes (SW up, SW down, LW up, LW down) [W/m^2] at the
    // interfaces (OUTPUT ONLY). Index k holds the lower interface of layer k
    // (k = klo is the surface); the top-of-atmosphere interface is stored in
    // the z-ghost cell k = khi + 1, so the array has one ghost cell in z.
    // Both RRTMGP and the two-stream model fill all nz + 1 levels.
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> rad_fluxes;
 
    amrex::Real solar_declin = 0.0;
    amrex::Real calday = 0.0;
 
    // Two-stream radiation model (erf.radiation_model = TwoStream): owns its
    // surface and SEB state; every call from ERF is behind a test of
    // solverChoice.rad_type == RadiationType::TwoStream.
    TwoStreamRadiation two_stream_rad;
 
#ifdef ERF_USE_EAMXX_SHOC
    amrex::Vector<std::unique_ptr<SHOCInterface>> eamxx_shoc_interface; // EAMxx SHOC model at each level
#endif
 
    amrex::Vector<std::unique_ptr<ShocDriver>> native_shoc_driver; // Native SHOC model at each level
 
    bool plot_rad = false;
    int rad_datalog_int = -1;
 
    // Fillpatcher classes for coarse-fine boundaries
    int cf_width{0};
    int cf_set_width{0};
    amrex::Vector<ERFFillPatcher> FPr_c;
    amrex::Vector<ERFFillPatcher> FPr_u;
    amrex::Vector<ERFFillPatcher> FPr_v;
    amrex::Vector<ERFFillPatcher> FPr_w;
 
    // Diffusive stresses and Smag
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> Tau;
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> Tau_corr; // correction to tau for implicit
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> eddyDiffs_lev;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> SmnSmn_lev;
 
    // Sea Surface Temps, Skin Temperature and Land Masks (lev, ntimes)
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>  sst_lev;
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>  tsk_lev;
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::iMultiFab>>> lmask_lev;
 
    // Land grid types and urban fraction
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::iMultiFab>>> land_type_lev;
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::iMultiFab>>> soil_type_lev;
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>  urb_frac_lev;
 
    // Other SFS terms
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_hfx1_lev, SFS_hfx2_lev, SFS_hfx3_lev;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_diss_lev;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q1fx1_lev, SFS_q1fx2_lev, SFS_q1fx3_lev;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q2fx3_lev;
 
    // Diffusive fluxes for EB
    // Tau_EB[level][component][grid_type]
    amrex::Vector<amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>> Tau_EB;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> hfx3_EB;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> qfx3_EB;
 
    // Grid stretching
    amrex::Vector<amrex::Vector<amrex::Real>> zlevels_stag; // nominal height levels
 
    // Data structures for terrain-fitted coordinates
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_nd;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_cc;
 
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   detJ_cc;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   ax;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   ay;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   az;
 
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_nd_src;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_cc_src;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   detJ_cc_src;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   ax_src;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   ay_src;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   az_src;
 
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_nd_new;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>   detJ_cc_new;
 
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_t_rk;
 
    // Immersed-boundary surface energy balance: inputs and one face set per level
    IBSEBParams ibseb_params;
    amrex::Vector<std::unique_ptr<IBFaceSet>> m_ibseb;
    TwoStreamRadiation::OrbitalCache ibseb_orbit;   ///< Orbit of the year for erf.ibseb.sun_mode = two_stream
    amrex::Real ibseb_z_i_level0 = 0.0;             ///< Level 0's bulk-Richardson depth of its current step, for the refined levels
 
    // Data structures for terrain immersed forcing
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> terrain_blanking;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> terrain_blanking_xface;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> terrain_blanking_yface;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> terrain_blanking_zface;
 
    // Wall distance function
    amrex::Vector<std::unique_ptr<amrex::MultiFab>>  walldist;
 
    // Map scale factors -- vector across levels of vector across types
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> mapfac;
 
    // Mask to be used to identify cells covered by finer level (if level > 0)
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> fine_mask;
 
    amrex::Vector<amrex::Vector<amrex::Real>> stretched_dz_h;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> stretched_dz_d;
 
    amrex::Vector<amrex::MultiFab> base_state;
    amrex::Vector<amrex::MultiFab> base_state_new;
 
    // Wave coupling data
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> Hwave;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> Lwave;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> Hwave_onegrid;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> Lwave_onegrid;
    bool finished_wave = false;
 
    // array of flux registers
    amrex::Vector<std::unique_ptr<amrex::YAFluxRegister>> advflux_reg;
 
    // A BCRec is essentially a 2*DIM integer array storing the boundary
    // condition type at each lo/hi walls in each direction. We have one BCRec
    // for each component of the cell-centered variables and each velocity component.
    amrex::Vector           <amrex::BCRec> domain_bcs_type;
    amrex::Gpu::DeviceVector<amrex::BCRec> domain_bcs_type_d;
 
    // We store these so that we can print them out in the job_info file
    amrex::Array<std::string,2*AMREX_SPACEDIM> domain_bc_type;
 
    // These hold the Dirichlet values at walls which need them ...
    amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals;
 
    // These hold the Neumann values at walls which need them ...
    amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals;
 
    // Host-side, face-local metadata for solid-wall scalar inputs. This is
    // deliberately separate from the values copied into device ghost-fill
    // kernels so explicit zero values remain distinguishable from omission.
    amrex::Array<amrex::Array<erf_wall_scalar_bc::WallScalarBCIntent, AMREX_SPACEDIM*2>, NBCVAR_max>
        m_wall_scalar_bc_intent{};
 
    // Per-face flag: non-reflecting inflow (extrapolate RhoTheta instead of prescribing)
    amrex::Array<bool, AMREX_SPACEDIM*2> m_bc_nonreflecting = {{false}};
    // Per face: whether its Dirichlet theta comes from a dirichlet_file, and in which form
    // (0: not from a file, 1: the file holds rho*theta, 2: the file holds theta)
    amrex::Array<int, AMREX_SPACEDIM*2> m_th_file_face = {{0}};
 
    // These are the "physical" boundary condition types (e.g. "inflow")
    amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2> phys_bc_type;
 
    // These are the masks for the thin immersed body
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> xflux_imask;
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> yflux_imask;
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> zflux_imask;
    //amrex::Vector<std::unique_ptr<amrex::iMultiFab>> overset_imask;
 
    // These are the body forces that result from the thin immersed body
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> thin_xforce;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> thin_yforce;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> thin_zforce;
 
    static int last_plot3d_file_step_1;
    static int last_plot3d_file_step_2;
    static int last_plot2d_file_step_1;
    static int last_plot2d_file_step_2;
    static int last_check_file_step;
 
    static double last_plot3d_file_time_1;
    static double last_plot3d_file_time_2;
    static double last_plot2d_file_time_1;
    static double last_plot2d_file_time_2;
    static double last_check_file_time;
 
    amrex::Vector<int> last_subvol_step;
    amrex::Vector<double> last_subvol_time;
 
    static bool plot_file_on_restart;
 
    // Output control
    const int datwidth = 14;
    const int datprecision = 6;
    const int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ amrex::Real(0.01) ==> min prec = 10
                                  // epoch time for 2030-01-01 00:00:00 is 1,893,456,000 s with dt ~ amrex::Real(0.01) ==> min prec = 12
 
    ////////////////
    // runtime parameters
 
    // Maximum number of steps
    int max_step = -1;
 
    // Start and stop times
    static double start_time;
    static double  stop_time;
 
    bool use_datetime = false;
    const std::string datetime_format = "%Y-%m-%d %H:%M:%S"; // ISO 8601 standard
 
    // if >= 0 we restart from a checkpoint
    std::string restart_chkfile = "";
 
    // Host-only metadata shared by all native output paths. It never enters
    // device code or changes the model state read from a checkpoint.
    erf_provenance::ExecutionProvenance execution_provenance;
 
    // Time step controls
    static amrex::Real cfl;
    static amrex::Real sub_cfl;
    static amrex::Real init_shrink;
    static amrex::Real change_max;
    static double dt_max_initial;
    static double dt_max;
 
    // Explicit eddy-diffusion time-step check (opt-in): warn when dt * rate exceeds
    // diffusive_cfl (diffusive_dt_check) and optionally limit an adaptive dt to it
    static bool diffusive_dt_check;
    static bool diffusive_dt_limit;
    static amrex::Real diffusive_cfl;
    // true until the first ComputeDt of this run has finished (restart hold of the limit)
    bool diffusive_first_call = true;
    // largest Fourier number already reported per level (warnings repeat only when it grows)
    amrex::Vector<amrex::Real> diffusive_fourier_warned;
    // grids on which, and value at which, the terrain slope factor was last reported, per level
    amrex::Vector<amrex::BoxArray> slope_report_grids;
    amrex::Vector<amrex::Real> slope_report_alpha;
 
    // Fixed dt at each level (only used if positive)
    amrex::Vector<double> fixed_dt;
    amrex::Vector<double> fixed_fast_dt;
    static int fixed_mri_dt_ratio;
 
    // how often each level regrids the higher levels of refinement
    // (after a level advances that many time steps)
    int regrid_int = -1;
 
    // Should we regrid level 0 after restart if Ngrids < Nprocs or
    //    max_grid_size has changed?
    bool regrid_level_0_on_restart = false;
 
    //
    // The grids each level's checkpoint was written on, and whether level 0's grids actually
    // moved when it was read. ReadCheckpointFile decides the level-0 grids before it builds
    // the levels, so every field is read on the grids in the file and redistributed onto the
    // grids the run will use; these carry the former. Empty away from a restart.
    //
    amrex::Vector<amrex::BoxArray>            ba_chkfile;
    amrex::Vector<amrex::DistributionMapping> dm_chkfile;
    bool regrid_level_0_on_restart_in_effect = false;
 
    //
    // Read one checkpoint MultiFab and deposit it into dst, which may live on different grids.
    //
    // The source is read into an undefined MultiFab, which VisMF::Read then defines on the
    // BoxArray recorded in the file, so no caller needs to know the checkpoint's grids. Only
    // the VALID cells of the source are offered by default: a file's ghost cells are not
    // guaranteed to agree with its valid data, and ParallelCopy intersects the grown source
    // with the grown destination without any ordering guarantee between overlapping sources,
    // so offering them can land stale values in valid cells (issue 4225). Pass
    // src_ghosts = true only for a field whose halo the checkpoint carries deliberately and
    // whose exterior halo cannot be rebuilt -- the terrain heights are the case that needs it.
    //
    void read_checkpoint_mf (amrex::MultiFab& dst,
                             const std::string& mf_name,
                             int lev,
                             int scomp, int dcomp, int ncomp,
                             const amrex::IntVect& dst_ghosts,
                             bool src_ghosts = false);
 
    //
    // As above, returning false instead of aborting when the file is not there, for fields a
    // older checkpoint may not carry.
    //
    bool read_checkpoint_mf_if_present (amrex::MultiFab& dst,
                                        const std::string& mf_name,
                                        int lev,
                                        int scomp, int dcomp, int ncomp,
                                        const amrex::IntVect& dst_ghosts,
                                        bool src_ghosts = false);
 
    // plotfile prefix and frequency
    std::string plot3d_file_1 {"plt_1_"};
    std::string plot3d_file_2 {"plt_2_"};
    std::string plot2d_file_1 {"plt2d_1_"};
    std::string plot2d_file_2 {"plt2d_2_"};
    std::string subvol_file {"subvol"};
    bool m_expand_plotvars_to_unif_rr = false;
    int m_plot3d_int_1 = -1;
    int m_plot3d_int_2 = -1;
    int m_plot2d_int_1 = -1;
    int m_plot2d_int_2 = -1;
 
    amrex::Vector<int> m_subvol_int;
    amrex::Vector<double> m_subvol_per;
 
    double m_plot3d_per_1 = -1.0;
    double m_plot3d_per_2 = -1.0;
    double m_plot2d_per_1 = -1.0;
    double m_plot2d_per_2 = -1.0;
    bool m_plot_face_vels = false;
    bool m_plot_face_terrain_blanking = false;
 
    bool plot_lsm = false;
    bool plot_surfmodel = false;
 
    // other sampling output control
    int profile_int = -1;
    bool destag_profiles = true;
 
    // Checkpoint type, prefix and frequency
    std::string check_file {"chk"};
    int m_check_int = -1;
    double m_check_per = -1.0;
 
    amrex::Vector<std::string> subvol3d_var_names;
 
    amrex::Vector<std::string> plot3d_var_names_1;
    amrex::Vector<std::string> plot3d_var_names_2;
    amrex::Vector<std::string> plot2d_var_names_1;
    amrex::Vector<std::string> plot2d_var_names_2;
    const amrex::Vector<std::string> cons_names     {"density", "rhotheta", "rhoKE", "rhoadv_0",
                                                     "rhoQ1", "rhoQ2", "rhoQ3",
                                                     "rhoQ4", "rhoQ5", "rhoQ6",
                                                     "rhoQ7", "rhoQ8", "rhoQ9",
                                                     "rhoQ10", "rhoQ11"};
 
    // **************************************************************************************
    // NOTE: The order of variable names here **MUST MATCH THE ORDER** in IO/ERF_Plotfile.cpp
    // **************************************************************************************
    const amrex::Vector<std::string> derived_names {"temp", "theta", "KE", "scalar", "soundspeed",
                                                    "reflectivity","max_reflectivity", "precipitable", "mucape",
                                                    "vorticity_x","vorticity_y","vorticity_z",
                                                    "helicity", "local_helicity", "vort_stretching",
                                                    "magvel", "divU",
                                                    "pres_hse", "dens_hse", "theta_hse", "pi_hse", "qv_hse",
                                                    "pressure", "pert_pres", "pert_dens",
                                                    "buoyancy", "eq_pot_temp", "VPD",
#ifdef ERF_USE_WINDFARM
                                                    "num_turb", "SMark0", "SMark1",
#endif
                                                    "dpdx", "dpdy", "dpdz", "pres_hse_x", "pres_hse_y",
                                                    "z_phys", "detJ", "h_xi", "h_eta", "h_zeta",
                                                    "mapfac", "lat_m", "lon_m",
                                                    // Time averaged velocity
                                                    "u_t_avg", "v_t_avg", "w_t_avg", "umag_t_avg",
                                                    // Interval means and resolved turbulence moments
                                                    "u_mean", "v_mean", "w_mean", "theta_mean",
                                                    "uu_mean", "vv_mean", "ww_mean", "uw_mean", "vw_mean", "wtheta_mean",
                                                    "uu_fluct", "vv_fluct", "ww_fluct", "uw_fluct", "vw_fluct", "wtheta_fluct", "tke_resolved",
                                                    // eddy viscosity
                                                    "nut",
                                                    // eddy diffusivity of momentum
                                                    "Kmv","Kmh",
                                                    // eddy diffusivity of heat
                                                    "Khv","Khh",
                                                    // turbulence lengthscale
                                                    "Lturb",
                                                    // k-eqn RANS: smoothed turbulent Richardson number and stability functions
                                                    "Rt", "cmu", "cmu_prime",
                                                    // native SHOC diagnostics
                                                    "pblh", "shoc_cldfrac", "shoc_ql", "shoc_ql2", "shoc_cond",
                                                    "wqls_sec", "wthv_sec", "w_sec",
                                                    "thl_sec", "qw_sec", "qwthl_sec", "wthl_sec", "wqw_sec", "w3",
                                                    "brunt", "isotropy", "shear_prod", "buoy_prod", "diss_tke",
                                                    // wall distance
                                                    "walldist",
                                                    // dissipation
                                                    "diss",
                                                    // moisture vars
                                                    "moist_density", "qv", "qc", "qi", "qrain", "qsnow", "qgraup",
                                                    "qt", "qn", "qp", "qsat", "nc", "ni", "nr", "ns", "ng", "nn",
                                                    "rain_accum", "snow_accum", "graup_accum",
                                                    "rel_humidity", "condensation_rate",
                                                    // Terrain IB mask
                                                    "terrain_IB_mask",
                                                    "ibseb_nfaces", "ibseb_tskin",
                                                    "ibseb_sw_abs", "ibseb_shadow",
                                                    "ibseb_lw_net", "ibseb_f_sky", "ibseb_H", "ibseb_G",
                                                    // EB variables
                                                    "volfrac"
#ifdef ERF_COMPUTE_ERROR
                                                    // error vars
                                                    ,"xvel_err", "yvel_err", "zvel_err", "pp_err"
#endif
                                                    ,"qsrc_sw", "qsrc_lw",
                                                    "Tau11", "Tau12", "Tau13", "Tau21", "Tau22", "Tau23", "Tau31", "Tau32", "Tau33",
                                                    "hfx1", "hfx2", "hfx3", "q1fx1", "q1fx2", "q1fx3", "q2fx3"
                                                   };
 
    // **************************************************************************************
    // NOTE: This is the set of derived variables that WriteSubvolume knows how to fill; it
    //       is an unordered allow-list, since setSubVolVariables walks "derived_names" and
    //       keeps the entries found here.  That means the names are selected in
    //       "derived_names" order, which is the order WriteSubvolume must compute them in.
    // **************************************************************************************
    const amrex::Vector<std::string> derived_subvol_names {"soundspeed", "temp", "theta", "KE", "scalar",
                                                           "precipitable", "mucape"};
 
    // algorithm choices
    static SolverChoice solverChoice;
 
    // Turbulent perturbation structure
    TurbulentPerturbation turbPert;
 
#ifdef ERF_USE_PARTICLES
    // Particle container with all particle species
    ParticleData particleData;
 
    // variables and functions for tracers particles
    bool  m_use_tracer_particles; /*!< tracer particles that advect with flow with optional sedimentation */
 
    /*! Read tracer and hydro particles parameters */
    void readTracersParams ();
 
    /*! Initialize tracer and hydro particles */
    void initializeTracers ( amrex::ParGDBBase*,
                             const amrex::Vector<std::unique_ptr<amrex::MultiFab>>&,
                             const double time);
 
    /*! Restart tracer and hydro particles */
    void restartTracers ( amrex::ParGDBBase*, const std::string& );
 
    /*! Evolve tracers and hydro particles */
    void evolveTracers ( int,
                        double,
                        amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
                        const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& );
 
#endif
 
#ifdef ERF_USE_MULTIBLOCK
    MultiBlockContainer *m_mbc = nullptr;
#endif
 
    static int verbose;
    static int mg_verbose;
    static bool use_fft;
 
    static int check_for_nans;
 
    // Diagnostic output interval
    static int sum_interval;
    static int pert_interval;
    static double sum_per;
 
    // Write in native AMReX or NetCDF format for each plotfile
    static PlotFileType plotfile3d_type_1;
    static PlotFileType plotfile3d_type_2;
    static PlotFileType plotfile2d_type_1;
    static PlotFileType plotfile2d_type_2;
 
    int file_name_digits = 5;
    bool use_real_time_in_pltname = false;
 
    static StateInterpType interpolation_type;
 
    // NetCDF initialization (wrfinput) file
    static amrex::Vector<amrex::Vector<std::string>> nc_init_file;
    static amrex::Vector<amrex::Vector<int>> have_read_nc_init_file;
 
    // NetCDF initialization (wrfbdy/met_em) file
    static std::string nc_bdy_file;
    int real_width{0};
    bool real_extrap_w{true};
 
    // NetCDF initialization of zlo boundary values
    static std::string nc_low_file;
 
    // Options for vertical interpolation of met_em*.nc data.
    bool        metgrid_debug_quiescent{false};
    bool        metgrid_debug_isothermal{false};
    bool        metgrid_debug_dry{false};
    bool        metgrid_debug_psfc{false};
    bool        metgrid_debug_msf{false};
    bool        metgrid_interp_theta{false};
    bool        metgrid_basic_linear{false};
    bool        metgrid_use_below_sfc{true};
    bool        metgrid_use_sfc{true};
    bool        metgrid_retain_sfc{false};
    amrex::Real metgrid_proximity{amrex::Real(500.0)};
    int         metgrid_order{2};
    int         metgrid_force_sfc_k{6};
 
    // Options for ERF boundary files for real simulations.
    bool        write_erfbdy{false};   // User-settable: write erfbdy during initialization
    bool        use_erfbdy{false};     // User-settable for wrfinput: read erfbdy instead of
                                       //    nc_bdy_file.  Set internally for metgrid and on restart.
    std::string erfbdy_file{"erfbdy"}; // User-settable: erfbdy file name
    int         nvars_erfbdy{0};       // Internal: number of variables in erfbdy
 
    amrex::Vector<amrex::BoxArray> ba1d;
    amrex::Vector<amrex::BoxArray> ba2d;
 
    // Vertical extent of the grid column over each (i,j), stored on ba2d (so the
    // compressed z index is 0) with one ghost cell in x and y so that a box can also
    // read the extent of the column on the other side of each of its faces.
    //     component 0 : smallEnd(2) of the box covering this column
    //     component 1 : bigEnd(2)   of the box covering this column
    // Ghost entries that no box covers keep the "unconstraining" sentinels below.
    // This is well defined only because grids may not be split in the vertical --
    // see ERF::grids_are_split_in_z.  Used by the implicit vertical diffusion solves.
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> column_kextent;
 
    static constexpr int column_kextent_lo_sentinel = std::numeric_limits<int>::min()/2;
    static constexpr int column_kextent_hi_sentinel = std::numeric_limits<int>::max()/2;
 
    void define_column_kextent (int lev, const amrex::BoxArray& ba,
                                const amrex::DistributionMapping& dm);
 
    static bool grids_are_split_in_z (const amrex::BoxArray& ba);
 
    // The six parameters defining the WRF analytic reference profile, plus the layer
    // interfaces derived from them.  These are read from the LEVEL 0 wrfinput file only
    // and then shared by every level: the levels of an AMR hierarchy must be built from
    // one reference profile, or their base states cannot agree across the coarse/fine
    // interface.  (A nested wrfinput file that disagrees produces a warning.)
    BaseStateParams wrf_bsp;
 
    std::unique_ptr<amrex::MultiFab> wrf_C1H;
    std::unique_ptr<amrex::MultiFab> wrf_C2H;
    std::unique_ptr<amrex::MultiFab> wrf_MUB;
    std::unique_ptr<amrex::MultiFab> wrf_RDNW;
    std::unique_ptr<amrex::MultiFab> wrf_PHB; // Base geopotential
 
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> mf_PSFC;
 
    // 1D CDF output (for ingestion in AMR-Wind)
    static int         output_1d_column;
    static int         column_interval;
    static double column_per;
    static amrex::Real column_loc_x;
    static amrex::Real column_loc_y;
    static std::string column_file_name;
 
    // 2D BndryRegister output (for ingestion in AMR-Wind)
    static int         output_bndry_planes;
    static int         bndry_output_planes_interval;
    static double bndry_output_planes_per;
    static double bndry_output_planes_start_time;
 
    // 2D BndryRegister input
    static int          input_bndry_planes;
 
    static int ng_dens_hse;
    static int ng_pres_hse;
 
    // Custom source terms
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> rhotheta_src;
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> rhoqt_src;
 
    amrex::Vector<           amrex::Vector<amrex::Real> > h_w_subsid;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > d_w_subsid;
 
    amrex::Vector<           amrex::Vector<amrex::Real> > h_u_geos;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > d_u_geos;
 
    amrex::Vector<           amrex::Vector<amrex::Real> > h_v_geos;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > d_v_geos;
 
    // Function to read and populate above host vectors (if input file exists)
    void init_geo_wind_profile (const std::string input_file,
                                amrex::Vector<amrex::Real>& u_geos,
                                amrex::Gpu::DeviceVector<amrex::Real>& u_geos_d,
                                amrex::Vector<amrex::Real>& v_geos,
                                amrex::Gpu::DeviceVector<amrex::Real>& v_geos_d,
                                const amrex::Geometry& lgeom,
                                const amrex::Vector<amrex::Real>& zlev_stag);
 
    // This is a vector over levels of vectors across quantities of Vectors
    amrex::Vector<amrex::Vector<amrex::Vector<amrex::Real> > > h_rayleigh_ptrs;
    amrex::Vector<amrex::Vector<amrex::Vector<amrex::Real> > > h_sponge_ptrs;
 
    // These are vectors over levels of Vectors
    amrex::Vector<amrex::Vector<amrex::Real> > h_sinesq_ptrs;
    amrex::Vector<amrex::Vector<amrex::Real> > h_sinesq_stag_ptrs;
 
    // This is a vector over levels of vectors across quantities of DeviceVectors
    amrex::Vector<amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > > d_rayleigh_ptrs;
    amrex::Vector<amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > > d_sponge_ptrs;
 
    // These are vectors over levels of DeviceVectors
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > d_sinesq_ptrs;
    amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > d_sinesq_stag_ptrs;
 
    amrex::Vector<amrex::Real> h_havg_density;
    amrex::Vector<amrex::Real> h_havg_temperature;
    amrex::Vector<amrex::Real> h_havg_pressure;
    amrex::Vector<amrex::Real> h_havg_qv;
    amrex::Vector<amrex::Real> h_havg_qc;
 
    amrex::Gpu::DeviceVector<amrex::Real> d_havg_density;
    amrex::Gpu::DeviceVector<amrex::Real> d_havg_temperature;
    amrex::Gpu::DeviceVector<amrex::Real> d_havg_pressure;
    amrex::Gpu::DeviceVector<amrex::Real> d_havg_qv;
    amrex::Gpu::DeviceVector<amrex::Real> d_havg_qc;
 
    void refinement_criteria_setup ();
 
    std::unique_ptr<WriteBndryPlanes> m_w2d  = nullptr;
    std::unique_ptr<ReadBndryPlanes>  m_r2d  = nullptr;
    amrex::Vector<std::unique_ptr<SurfaceLayer>>     m_SurfaceLayer;
    amrex::Vector<std::unique_ptr<ForestDrag>> m_forest_drag;
 
    //
    // Holds info for dynamically generated tagging criteria
    //
    static amrex::Vector<amrex::AMRErrorTag> ref_tags;
 
    //
    // Name of the refinement indicator that generated each entry of ref_tags.
    // Not every indicator generates a ref_tag (e.g. "storm_tracker" does not),
    // so this must be carried alongside ref_tags rather than re-deriving it
    // from the positional index into erf.refinement_indicators.
    //
    static amrex::Vector<std::string> ref_tag_indicator_names;
 
    amrex::Vector<amrex::Vector<amrex::BoxArray>> subdomains;
 
    amrex::Vector<amrex::Real> dz_min;
 
    static AMREX_FORCE_INLINE
    int
    ComputeGhostCells (const SolverChoice& sc)
    {
        // We never have fewer than 2 ghost cells
        int ngrow = 2;
 
        // This sets the minimum to 3 if using num_diff but we might still need more
        if (sc.use_num_diff) {
            ngrow = 3;
        }
 
        // With EB we always have 4 so no need to do further tests
        if (sc.terrain_type == TerrainType::EB) {
            ngrow = 4;
            return ngrow;
        }
 
        if (
             (sc.advChoice.dycore_horiz_adv_type    == AdvType::Centered_6th)
          || (sc.advChoice.dycore_vert_adv_type     == AdvType::Centered_6th)
          || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Centered_6th)
          || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Centered_6th)
          || (sc.advChoice.dryscal_horiz_adv_type   == AdvType::Centered_6th)
          || (sc.advChoice.dryscal_vert_adv_type    == AdvType::Centered_6th) )
        { ngrow = 3; }
        else if (
             (sc.advChoice.dycore_horiz_adv_type    == AdvType::Upwind_5th)
          || (sc.advChoice.dycore_vert_adv_type     == AdvType::Upwind_5th)
          || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th)
          || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th)
          || (sc.advChoice.dryscal_horiz_adv_type   == AdvType::Upwind_5th)
          || (sc.advChoice.dryscal_vert_adv_type    == AdvType::Upwind_5th) )
        { ngrow = 3; }
        else if (
              (sc.advChoice.dycore_horiz_adv_type    == AdvType::Weno_5)
           || (sc.advChoice.dycore_vert_adv_type     == AdvType::Weno_5)
           || (sc.advChoice.dryscal_horiz_adv_type   == AdvType::Weno_5)
           || (sc.advChoice.dryscal_vert_adv_type    == AdvType::Weno_5)
           || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Weno_5)
           || (sc.advChoice.moistscal_vert_adv_type  == AdvType::Weno_5) )
        { ngrow = 3; }
        else if (
              (sc.advChoice.dycore_horiz_adv_type    == AdvType::Weno_5Z)
           || (sc.advChoice.dycore_vert_adv_type     == AdvType::Weno_5Z)
           || (sc.advChoice.dryscal_horiz_adv_type   == AdvType::Weno_5Z)
           || (sc.advChoice.dryscal_vert_adv_type    == AdvType::Weno_5Z)
           || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Weno_5Z)
           || (sc.advChoice.moistscal_vert_adv_type  == AdvType::Weno_5Z) )
        { ngrow = 3; }
        else if (
              (sc.advChoice.dycore_horiz_adv_type    == AdvType::Weno_7)
           || (sc.advChoice.dycore_vert_adv_type     == AdvType::Weno_7)
           || (sc.advChoice.dryscal_horiz_adv_type   == AdvType::Weno_7)
           || (sc.advChoice.dryscal_vert_adv_type    == AdvType::Weno_7)
           || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Weno_7)
           || (sc.advChoice.moistscal_vert_adv_type  == AdvType::Weno_7) )
        { ngrow = 3; }
        else if (
              (sc.advChoice.dycore_horiz_adv_type    == AdvType::Weno_7Z)
           || (sc.advChoice.dycore_vert_adv_type     == AdvType::Weno_7Z)
           || (sc.advChoice.dryscal_horiz_adv_type   == AdvType::Weno_7Z)
           || (sc.advChoice.dryscal_vert_adv_type    == AdvType::Weno_7Z)
           || (sc.advChoice.moistscal_horiz_adv_type == AdvType::Weno_7Z)
           || (sc.advChoice.moistscal_vert_adv_type  == AdvType::Weno_7Z) )
        { ngrow = 4; }
 
        return ngrow;
    }
 
    AMREX_FORCE_INLINE
    amrex::YAFluxRegister* getAdvFluxReg (int lev)
    {
        return advflux_reg[lev].get();
    }
 
    AMREX_FORCE_INLINE
    std::ostream&
    DataLog (int i)
    {
        return *datalog[i];
    }
 
    AMREX_FORCE_INLINE
    std::ostream&
    DerDataLog (int i)
    {
        return *der_datalog[i];
    }
 
    AMREX_FORCE_INLINE
    int
    NumDataLogs () noexcept
    {
        return static_cast<int>(datalog.size());
    }
 
    AMREX_FORCE_INLINE
    int
    NumDerDataLogs () noexcept
    {
        return static_cast<int>(der_datalog.size());
    }
 
 
    AMREX_FORCE_INLINE
    std::ostream&
    SamplePointLog (int i)
    {
        return *sampleptlog[i];
    }
 
    AMREX_FORCE_INLINE
    int
    NumSamplePointLogs () noexcept
    {
        return static_cast<int>(sampleptlog.size());
    }
 
    AMREX_FORCE_INLINE
    std::ostream&
    SampleLineLog (int i)
    {
        return *samplelinelog[i];
    }
 
    AMREX_FORCE_INLINE
    int
    NumSampleLineLogs () noexcept
    {
        return static_cast<int>(samplelinelog.size());
    }
 
    amrex::IntVect&
    SamplePoint (int i)
    {
        return samplepoint[i];
    }
 
    AMREX_FORCE_INLINE
    int
    NumSamplePoints () noexcept
    {
        return static_cast<int>(samplepoint.size());
    }
 
    amrex::IntVect&
    SampleLine (int i)
    {
        return sampleline[i];
    }
 
    AMREX_FORCE_INLINE
    int
    NumSampleLines () noexcept
    {
        return static_cast<int>(sampleline.size());
    }
 
    static double startCPUTime;
    static double previousCPUTimeUsed;
 
    static double
    getCPUTime ()
    {
      int numCores = amrex::ParallelDescriptor::NProcs();
#ifdef _OPENMP
      numCores = numCores * omp_get_max_threads();
#endif
 
      double T =
        numCores * (amrex::ParallelDescriptor::second() - startCPUTime) +
        previousCPUTimeUsed;
 
      return T;
    }
 
    void setRecordDataInfo (int i, const std::string& filename) // NOLINT
    {
        if (amrex::ParallelDescriptor::IOProcessor())
        {
            datalog[i] = std::make_unique<std::fstream>();
            datalog[i]->open(filename.c_str(),std::ios::out|std::ios::app);
            if (!datalog[i]->good()) {
                amrex::FileOpenFailed(filename);
            }
        }
        amrex::ParallelDescriptor::Barrier("ERF::setRecordDataInfo");
    }
 
    void setRecordDerDataInfo (int i, const std::string& filename) // NOLINT
    {
        if (amrex::ParallelDescriptor::IOProcessor())
        {
            der_datalog[i] = std::make_unique<std::fstream>();
            der_datalog[i]->open(filename.c_str(),std::ios::out|std::ios::app);
            if (!der_datalog[i]->good()) {
                amrex::FileOpenFailed(filename);
            }
        }
        amrex::ParallelDescriptor::Barrier("ERF::setRecordDerDataInfo");
    }
 
    void setRecordEnergyDataInfo (int i, const std::string& filename) // NOLINT
    {
        if (amrex::ParallelDescriptor::IOProcessor())
        {
            tot_e_datalog[i] = std::make_unique<std::fstream>();
            tot_e_datalog[i]->open(filename.c_str(),std::ios::out|std::ios::app);
            if (!tot_e_datalog[i]->good()) {
                amrex::FileOpenFailed(filename);
            }
        }
        amrex::ParallelDescriptor::Barrier("ERF::setRecordEnergyDataInfo");
    }
 
    void setRecordSamplePointInfo (int i, int lev, amrex::IntVect& cell, const std::string& filename) // NOLINT
    {
        amrex::MultiFab dummy(grids[lev],dmap[lev],1,0);
        for (amrex::MFIter mfi(dummy); mfi.isValid(); ++mfi)
        {
            const amrex::Box& bx = mfi.validbox();
            if (bx.contains(cell)) {
                sampleptlog[i] = std::make_unique<std::fstream>();
                sampleptlog[i]->open(filename.c_str(),std::ios::out|std::ios::app);
                if (!sampleptlog[i]->good()) {
                    amrex::FileOpenFailed(filename);
                }
            }
        }
        amrex::ParallelDescriptor::Barrier("ERF::setRecordSamplePointInfo");
    }
 
    void setRecordSampleLineInfo (int i, int lev, amrex::IntVect& cell, const std::string& filename) // NOLINT
    {
        amrex::MultiFab dummy(grids[lev],dmap[lev],1,0);
        for (amrex::MFIter mfi(dummy); mfi.isValid(); ++mfi)
        {
            const amrex::Box& bx = mfi.validbox();
            if (bx.contains(cell)) {
                samplelinelog[i] = std::make_unique<std::fstream>();
                samplelinelog[i]->open(filename.c_str(),std::ios::out|std::ios::app);
                if (!samplelinelog[i]->good()) {
                    amrex::FileOpenFailed(filename);
                }
            }
        }
        amrex::ParallelDescriptor::Barrier("ERF::setRecordSampleLineInfo");
    }
 
    // Data sampler for line and plane output
    int line_sampling_interval = -1;
    int plane_sampling_interval = -1;
    amrex::Real line_sampling_per = -one;
    amrex::Real plane_sampling_per = -one;
    std::unique_ptr<LineSampler> line_sampler = nullptr;
    std::unique_ptr<PlaneSampler> plane_sampler = nullptr;
 
    // Station time-series output
    int station_sampling_interval = -1;
    amrex::Real station_sampling_per = -one;
    std::unique_ptr<StationSampler> station_sampler = nullptr;
    // The 3D and 2D plot variables the stations need, in the order the fill
    // routines produce their components
    amrex::Vector<std::string> station_vars_3d;
    amrex::Vector<std::string> station_vars_2d;
    // The immersed-forcing terrain surface under each level's grids, for the stations
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> station_ib_ground;
 
    // Nudging towards observations at stations; built only when
    // solverChoice.nudging_from_observations is true
    std::unique_ptr<ObsNudging> obs_nudging = nullptr;
 
    amrex::Vector<std::unique_ptr<std::fstream> >       datalog;
    amrex::Vector<std::unique_ptr<std::fstream> >   der_datalog;
    amrex::Vector<std::unique_ptr<std::fstream> > tot_e_datalog;
    amrex::Vector<std::string>       datalogname;
    amrex::Vector<std::string>   der_datalogname;
    amrex::Vector<std::string> tot_e_datalogname;
 
    amrex::Vector<std::unique_ptr<std::fstream> > sampleptlog;
    amrex::Vector<std::string> sampleptlogname;
    amrex::Vector<amrex::IntVect> samplepoint;
 
    amrex::Vector<std::unique_ptr<std::fstream> > samplelinelog;
    amrex::Vector<std::string> samplelinelogname;
    amrex::Vector<amrex::IntVect> sampleline;
 
    //! The filename of the ith datalog file.
    [[nodiscard]] std::string DataLogName (int i) const noexcept { return datalogname[i]; }
    [[nodiscard]] std::string DerDataLogName (int i) const noexcept { return der_datalogname[i]; }
 
    //! The filename of the ith sampleptlog file.
    [[nodiscard]] std::string SamplePointLogName (int i) const noexcept { return sampleptlogname[i]; }
 
    //! The filename of the ith samplelinelog file.
    [[nodiscard]] std::string SampleLineLogName (int i) const noexcept { return samplelinelogname[i]; }
 
    // array of EB objects
    amrex::Vector<std::unique_ptr<eb_>> eb;
 
    [[nodiscard]] eb_ const& get_eb (int lev) const noexcept {
        AMREX_ASSERT(lev >= 0 && lev < eb.size() && eb[lev] != nullptr);
        return *eb[lev];
    }
 
    //
    // NOTE: this dereferences the cell-centered EB factory, which is only built when
    //       we are using EB or immersed forcing -- see ERF::MakeNewLevelFromScratch.
    //       Calling this in any other case is undefined behavior, so we assert here
    //       rather than relying on the (build-dependent) libstdc++ hardening check.
    //
    [[nodiscard]] amrex::EBFArrayBoxFactory const&
    EBFactory (int lev) const noexcept {
        AMREX_ASSERT(lev >= 0 && lev < eb.size() && eb[lev] != nullptr);
        AMREX_ASSERT(eb[lev]->get_const_factory() != nullptr);
        return *(eb[lev]->get_const_factory());
    }
 
    [[nodiscard]] static int nghost_eb_basic ()
    { return 5; }
 
    // We need 5 for doing StateRedistribution; otherwise 4 would be enough
    [[nodiscard]] static int nghost_eb_volume ()
    { return 5; }
 
    [[nodiscard]] static int nghost_eb_full ()
    { return 4; }
 
#ifdef ERF_USE_FFT
    // Vector over levels of Vector over subdomains of solvers
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::FFT::Poisson<amrex::MultiFab>>>> m_3D_poisson;
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::FFT::PoissonHybrid<amrex::MultiFab>>>> m_2D_poisson;
#endif
 
public:
    void writeJobInfo (const std::string& dir,
                       erf_provenance::ArtifactType artifact_type,
                       int artifact_step,
                       double artifact_time_seconds) const;
    static void writeBuildInfo (std::ostream& os);
 
    static void print_banner (MPI_Comm /*comm*/, std::ostream& /*out*/);
    static void print_usage (MPI_Comm /*comm*/, std::ostream& /*out*/);
    static void print_error (MPI_Comm /*comm*/, const std::string& msg);
    static void print_summary (std::ostream&);
    static void print_tpls (std::ostream& /*out*/);
};
 
#endif