REMORA_NCPlotFile.cpp

Go to the documentation of this file.

#include <iomanip>
#include <iostream>
#include <map>
#include <string>
#include <vector>
#include <ctime>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include <AMReX_Utility.H>
#include <AMReX_buildInfo.H>
#include <AMReX_ParmParse.H>
 
#include "REMORA.H"
#include "REMORA_DateClock.H"
#include "REMORA_NCInterface.H"
#include "REMORA_NCPlotFile.H"
#include "REMORA_IndexDefines.H"
 
using namespace amrex;
 
namespace {
/**
 * \brief Was this 2D field named in remora.plot_vars_2d?
 *
 * Mirrors the name test the AMReX plotfile writer does when it walks plot_var_names_2d,
 * so both writers answer to the same input key.
 */
bool
plot_2d_var_requested (const amrex::Vector<std::string>& names, const std::string& nm)
{
    for (int i = 0; i < names.size(); ++i) {
        if (names[i] == nm) { return true; }
    }
    return false;
}
 
/**
 * \brief Accumulates each rank's hyperslabs and writes them with collective I/O.
 *
 * PnetCDF's collective APIs have to be called by every rank that opened the file,
 * the same number of times and in the same order. REMORA's write loops are MFIter
 * loops, so the number of hyperslabs is per-rank -- a plain put_all() inside those
 * loops would mismatch and hang. ncmpi_put_varn_*_all() takes all of a rank's
 * subarrays for one variable in a single call and tolerates counts that differ
 * between ranks, including zero, which is exactly the shape we need.
 *
 * So the loops hand their hyperslabs here instead of writing, and flush() issues
 * one collective call per variable at the end. Data is copied into m_entries on
 * add(), so callers may reuse or destroy their staging buffer immediately -- unlike
 * the nonblocking iput() path, where PnetCDF aliases buffers larger than 4 KB and
 * only reads them at wait_all().
 */
class VarnCollector
{
public:
    //! Stage one hyperslab of `varid`. `start`/`count` are in NetCDF dimension
    //! order and `dptr` holds product(count) values.
    void add (int varid,
              const std::vector<MPI_Offset>& start,
              const std::vector<MPI_Offset>& count,
              const amrex::Real* dptr)
    {
        AMREX_ASSERT(start.size() == count.size());
        MPI_Offset nelems = 1;
        for (auto c : count) { nelems *= c; }
 
        Entry& e = m_entries[varid];
        e.starts.push_back(start);
        e.counts.push_back(count);
        e.data.insert(e.data.end(), dptr, dptr + nelems);
    }
 
    /**
     * \brief Write everything staged so far. Collective: all ranks must call it.
     *
     * Ranks own different boxes, so they stage different variables -- a rank with
     * no boxes stages nothing at all. The set of variables to write therefore has
     * to be agreed on before any collective call, which is what the reduction over
     * the touched-variable mask below does. Iterating that agreed set in ascending
     * varid order gives every rank the same call sequence.
     */
    void flush (const ncutils::NCFile& ncf)
    {
        const int nvars = ncf.num_variables();
        if (nvars <= 0) { return; }
 
        std::vector<int> touched(nvars, 0);
        for (const auto& kv : m_entries) {
            if (kv.first >= 0 && kv.first < nvars) { touched[kv.first] = 1; }
        }
#ifdef AMREX_USE_MPI
        MPI_Allreduce(MPI_IN_PLACE, touched.data(), nvars, MPI_INT, MPI_MAX,
                      amrex::ParallelContext::CommunicatorSub());
#endif
 
        for (int varid = 0; varid < nvars; ++varid) {
            if (!touched[varid]) { continue; }
 
            auto it = m_entries.find(varid);
            if (it == m_entries.end()) {
                // This rank has nothing for this variable, but the call is
                // collective, so it still has to take part with zero subarrays.
                ncutils::NCVar { ncf.ncid, varid }.put_varn_all(
                    0, nullptr, nullptr, static_cast<const amrex::Real*>(nullptr));
                continue;
            }
 
            Entry& e = it->second;
            const int num = static_cast<int>(e.starts.size());
 
            // ncmpi_put_varn_* wants arrays of pointers into the start/count rows.
            std::vector<MPI_Offset*> start_ptrs(num);
            std::vector<MPI_Offset*> count_ptrs(num);
            for (int i = 0; i < num; ++i) {
                start_ptrs[i] = e.starts[i].data();
                count_ptrs[i] = e.counts[i].data();
            }
 
            ncutils::NCVar { ncf.ncid, varid }.put_varn_all(
                num, start_ptrs.data(), count_ptrs.data(), e.data.data());
        }
 
        m_entries.clear();
    }
 
    /**
     * \brief Handle with the same put() shape as ncutils::NCVar, staging instead
     *        of writing, so the write loops read the same as they always did.
     */
    class StagedVar
    {
    public:
        StagedVar (VarnCollector& coll, int varid) : m_coll(coll), m_varid(varid) {}
 
        void put (const amrex::Real* dptr,
                  const std::vector<MPI_Offset>& start,
                  const std::vector<MPI_Offset>& count) const
        {
            m_coll.add(m_varid, start, count, dptr);
        }
 
    private:
        VarnCollector& m_coll;
        int m_varid;
    };
 
    //! Handle for staging writes to the named variable.
    StagedVar var (const ncutils::NCFile& ncf, const std::string& name)
    {
        return StagedVar(*this, ncf.var(name).varid);
    }
 
private:
    struct Entry {
        std::vector<std::vector<MPI_Offset>> starts;
        std::vector<std::vector<MPI_Offset>> counts;
        std::vector<amrex::Real> data;
    };
 
    //! Keyed by varid and held in a std::map so iteration is in varid order.
    std::map<int, Entry> m_entries;
};
} // namespace
 
/**
 * @param which_step   current step for output
 * @param plotMF       cell-centered plot variables at this level
 * @param lev          level to write
 * @param mf_omega     omega at w points, or nullptr when it is not plotted
 * @param mf_wvel      w at w points, or nullptr when it is not plotted
 */
void REMORA::WriteNCPlotFile(int which_step, MultiFab const* plotMF, int lev,
                             MultiFab const* mf_omega, MultiFab const* mf_wvel) {
    // A refined level is written to its own file, _d02 and up, the way ROMS names a nested
    // grid. Its dimensions come from boxes_at_level, so a level covering nx by ny cells gets
    // xi_rho = nx+2, xi_u = nx+1 and so on, exactly as ROMS shapes a child grid.
    int which_subdomain = 0;
    int which_step_in_chunk = -1;
 
    // Create filename
    std::string plt_string;
    std::string plotfilename;
    if (REMORA::write_history_file) {
        plotfilename = plot_file_name + "_his";
    } else {
        plotfilename = Concatenate(plot_file_name, which_step, file_min_digits);
    }
    // If chunking, concatenate with which file we're in
    if (REMORA::write_history_file and REMORA::chunk_history_file) {
        int which_chunk = history_count / REMORA::steps_per_history_file;
        plotfilename = Concatenate(plotfilename, which_chunk, file_min_digits);
        which_step_in_chunk = history_count - which_chunk * REMORA::steps_per_history_file;
    }
 
    // Set the full IO path for NetCDF output
    std::string FullPath = plotfilename;
    if (lev == 0) {
        const std::string &extension = amrex::Concatenate("_d", lev + 1, 2);
        FullPath += extension + ".nc";
    } else {
        const std::string &extension = amrex::Concatenate("_d", lev + 1 + which_subdomain, 2);
        FullPath += extension + ".nc";
    }
 
    //
    // Check if this file/directory already exists and if so,
    //       have the IOProcessor move the existing
    //       file/directory to filename.old
    //
    if ((!REMORA::write_history_file) || (history_count == 0) ||
        (which_step_in_chunk == 0)) {
        if (amrex::ParallelDescriptor::IOProcessor()) {
            if (amrex::FileExists(FullPath)) {
                std::string newoldname(FullPath + ".old." + amrex::UniqueString());
                amrex::Print() << "WriteNCPlotFile:  " << FullPath << " exists.  Renaming to:  " << newoldname << std::endl;
                if (std::rename(FullPath.c_str(), newoldname.c_str())) {
                    amrex::Abort("WriteNCPlotfile:: std::rename failed");
                }
            }
        }
        ParallelDescriptor::Barrier();
    }
 
    bool is_history;
 
    if (REMORA::write_history_file) {
        is_history = true;
 
        if (history_file_in_use.size() <= lev) {
            history_file_in_use.resize(lev+1);
            history_record_shift.resize(lev+1, 0);
        }
        const long long record = chunk_history_file
                               ? history_count % steps_per_history_file
                               : history_count;
 
        // On the first write to a file, fix where its records start. A restarted
        // run continues an existing file at the restart time, else starts at 0.
        if (FullPath != history_file_in_use[lev]) {
            history_file_in_use[lev] = FullPath;
            long long start = record;
            if (!restart_chkfile.empty()) {
                start = amrex::FileExists(FullPath)
                      ? history_restart_record(lev, FullPath) : 0;
            }
            history_record_shift[lev] = record - start;
        }
        bool write_header = !(amrex::FileExists(FullPath));
 
        auto ncf =
                write_header ?
                        ncutils::NCFile::create(FullPath, NC_CLOBBER|NC_64BIT_DATA, amrex::ParallelContext::CommunicatorSub(), MPI_INFO_NULL) :
                        ncutils::NCFile::open(FullPath, NC_WRITE, amrex::ParallelContext::CommunicatorSub(), MPI_INFO_NULL);
 
        amrex::Print() << "Writing into level " << lev << " NetCDF history file " << FullPath << std::endl;
 
        WriteNCPlotFile_which(lev, which_subdomain, plotMF, mf_omega, mf_wvel, write_header, ncf,
                              is_history, record - history_record_shift[lev]);
 
    } else {
 
        is_history = false;
        bool write_header = true;
 
        // Open new netcdf file to write data
        auto ncf = ncutils::NCFile::create(FullPath, NC_CLOBBER|NC_64BIT_DATA, amrex::ParallelContext::CommunicatorSub(), MPI_INFO_NULL);
        amrex::Print() << "Writing level " << lev << " NetCDF plot file " << FullPath << std::endl;
 
        WriteNCPlotFile_which(lev, which_subdomain, plotMF, mf_omega, mf_wvel, write_header, ncf, is_history);
    }
}
 
/**
 * Find where a restarted run continues an existing history file: its first
 * record at or after the restart time. A file for a different grid is renamed
 * aside, and the run starts a new one at record 0.
 *
 * @param lev               level the file holds
 * @param path              history file to continue
 * @return                  record to write the restart step to
 */
long long REMORA::history_restart_record (int lev, const std::string& path)
{
    const Box subdomain = (lev == 0) ? geom[lev].Domain()
                                     : boxes_at_level[lev][0];
    const double t = model_time(t_new[lev]);
    long long nrec = 0;
    long long first = 0;
    bool fits;
    {
        auto ncf = ncutils::NCFile::open(
            path, NC_NOWRITE, amrex::ParallelContext::CommunicatorSub(),
            MPI_INFO_NULL);
        auto dim_is = [&] (const std::string& name, long long len) {
            return ncf.has_dim(name) && ncf.dim(name).len() == len;
        };
        fits = dim_is("xi_rho", subdomain.length(0) + 2) &&
               dim_is("eta_rho", subdomain.length(1) + 2) &&
               dim_is("s_rho", subdomain.length(2)) &&
               ncf.has_var("ocean_time");
        if (fits) {
            nrec = ncf.dim("ocean_time").len();
            std::vector<double> times(nrec);
            if (nrec > 0) {
                ncf.var("ocean_time").get_all(times.data(), {0}, {nrec});
            }
            // Skip records from before the restart, allowing half a step
            const double slack = 0.5 * static_cast<double>(dt[lev]);
            while (first < nrec && times[first] < t - slack) { ++first; }
        }
        ncf.close();
    }
 
    if (!fits) {
        if (amrex::ParallelDescriptor::IOProcessor()) {
            std::string newoldname(path + ".old." + amrex::UniqueString());
            amrex::Print() << "WriteNCPlotFile:  " << path
                           << " does not match level " << lev
                           << ".  Renaming to:  " << newoldname << std::endl;
            if (std::rename(path.c_str(), newoldname.c_str())) {
                amrex::Abort("WriteNCPlotfile:: std::rename failed");
            }
        }
        ParallelDescriptor::Barrier();
        return 0;
    }
 
    amrex::Print() << "WriteNCPlotFile:  continuing " << path
                   << " at record " << first;
    if (nrec > first) {
        amrex::Print() << " (" << nrec - first
                       << " later record(s) from the earlier run remain"
                       << " until written over)";
    }
    amrex::Print() << std::endl;
    return first;
}
 
/**
 * @param lev               level of data to output
 * @param which_subdomain   index of subdomain if lev != 0
 * @param plotMF            cell-centered plot variables at this level
 * @param mf_omega          omega at w points, or nullptr when it is not plotted
 * @param mf_wvel           w at w points, or nullptr when it is not plotted
 * @param write_header      whether to write a header
 * @param ncf               netcdf file object
 * @param is_history        whether the file being written is a history file
 * @param history_record    record of the history file to write to
 */
void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, MultiFab const* plotMF,
                                   MultiFab const* mf_omega, MultiFab const* mf_wvel,
                                   bool write_header, ncutils::NCFile &ncf,
                                   bool is_history, long long history_record)
{
    // Number of cells in this "domain" at this level
    std::vector<int> n_cells;
 
    // One level per file, so the geometry records below describe just this level.
    int flev = lev + 1;
 
    Box subdomain;
    if (lev == 0) {
        subdomain = geom[lev].Domain();
    } else {
        subdomain = boxes_at_level[lev][which_subdomain];
    }
 
    int nx = subdomain.length(0);
    int ny = subdomain.length(1);
    int nz = subdomain.length(2);
 
    n_cells.push_back(nx);
    n_cells.push_back(ny);
    n_cells.push_back(nz);
 
    const std::string nt_name = "ocean_time";
    const std::string ndim_name = "num_geo_dimensions";
 
    const std::string flev_name = "FINEST_LEVEL";
 
    const std::string nx_name = "NX";
    const std::string ny_name = "NY";
    const std::string nz_name = "NZ";
 
    const std::string nx_r_name = "xi_rho";
    const std::string ny_r_name = "eta_rho";
    const std::string nz_r_name = "s_rho";
 
    const std::string nx_u_name = "xi_u";
    const std::string ny_u_name = "eta_u";
 
    const std::string nx_v_name = "xi_v";
    const std::string ny_v_name = "eta_v";
 
    const std::string nx_p_name = "xi_psi";
    const std::string ny_p_name = "eta_psi";
    const std::string nz_w_name = "s_w";
 
    if (write_header) {
        ncf.enter_def_mode();
        ncf.put_attr("title", "REMORA data ");
        // The time dimension is unlimited so the record count reflects what was
        // actually written rather than an up-front estimate from max_step/plot_int.
        // Note PnetCDF does not prefill record variables, so every element of a
        // time-varying variable has to be written explicitly (see the loops below).
        ncf.def_dim(nt_name, NC_UNLIMITED);
        ncf.def_dim(ndim_name, AMREX_SPACEDIM);
 
        ncf.def_dim(nx_r_name, nx + 2);
        ncf.def_dim(ny_r_name, ny + 2);
        ncf.def_dim(nz_r_name, nz);
 
        ncf.def_dim(nx_u_name, nx + 1);
        ncf.def_dim(ny_u_name, ny + 2);
 
        ncf.def_dim(nx_v_name, nx + 2);
        ncf.def_dim(ny_v_name, ny + 1);
 
        ncf.def_dim(nx_p_name, nx + 1);
        ncf.def_dim(ny_p_name, ny + 1);
 
        ncf.def_dim(nz_w_name, nz + 1);
 
        ncf.def_dim(flev_name, flev);
 
        ncf.def_dim(nx_name, n_cells[0]);
        ncf.def_dim(ny_name, n_cells[1]);
        ncf.def_dim(nz_name, n_cells[2]);
 
        ncf.def_var("probLo", ncutils::NCDType::Real, { ndim_name });
        ncf.var("probLo").put_attr("long_name","Low side of problem domain in internal AMReX grid");
        ncf.var("probLo").put_attr("units","meter");
        ncf.def_var("probHi", ncutils::NCDType::Real, { ndim_name });
        ncf.var("probHi").put_attr("long_name","High side of problem domain in internal AMReX grid");
        ncf.var("probHi").put_attr("units","meter");
 
        ncf.def_var("Geom.smallend", NC_INT, { flev_name, ndim_name });
        ncf.var("Geom.smallend").put_attr("long_name","Low side of problem domain in index space");
        ncf.def_var("Geom.bigend", NC_INT, { flev_name, ndim_name });
        ncf.var("Geom.bigend").put_attr("long_name","High side of problem domain in index space");
        ncf.def_var("CellSize", ncutils::NCDType::Real, { flev_name, ndim_name });
        ncf.var("CellSize").put_attr("long_name","Cell size on internal AMReX grid");
        ncf.var("CellSize").put_attr("units","meter");
 
        ncf.def_var("theta_s",ncutils::NCDType::Real,{});
        ncf.var("theta_s").put_attr("long_name","S-coordinate surface control parameter");
        ncf.def_var("theta_b",ncutils::NCDType::Real,{});
        ncf.var("theta_b").put_attr("long_name","S-coordinate bottom control parameter");
        ncf.def_var("hc",ncutils::NCDType::Real,{});
        ncf.var("hc").put_attr("long_name","S-coordinate parameter, critical depth");
        ncf.var("hc").put_attr("units","meter");
 
        ncf.def_var("grid",NC_INT, {});
        ncf.var("grid").put_attr("cf_role","grid_topology");
        ncf.var("grid").put_attr("topology_dimension",std::vector({2}));
        ncf.var("grid").put_attr("node_dimensions", "xi_psi eta_psi");
        ncf.var("grid").put_attr("face_dimensions", "xi_rho: xi_psi (padding: both) eta_rho: eta_psi (padding: both)");
        ncf.var("grid").put_attr("edge1_dimensions", "xi_u: xi_psi eta_u: eta_psi (padding: both)");
        ncf.var("grid").put_attr("edge2_dimensions", "xi_v: xi_psi (padding: both) eta_v: eta_psi");
        ncf.var("grid").put_attr("node_coordinates", "x_psi y_psi");
        ncf.var("grid").put_attr("face_coordinates", "x_rho y_rho");
        ncf.var("grid").put_attr("edge1_coordinates", "x_u y_u");
        ncf.var("grid").put_attr("edge2_coordinates", "x_v y_v");
        ncf.var("grid").put_attr("vertical_dimensions", "s_rho: s_w (padding: none)");
 
        ncf.def_var("s_rho",ncutils::NCDType::Real, {nz_r_name});
        ncf.var("s_rho").put_attr("long_name","S-coordinate at RHO-points");
        ncf.var("s_rho").put_attr("field","s_rho, scalar");
 
        ncf.def_var("s_w",ncutils::NCDType::Real, {nz_w_name});
        ncf.var("s_w").put_attr("long_name","S-coordinate at W-points");
        ncf.var("s_w").put_attr("field","s_w, scalar");
 
        ncf.def_var("pm",ncutils::NCDType::Real, {ny_r_name, nx_r_name});
        ncf.var("pm").put_attr("long_name","curvilinear coordinate metric in XI");
        ncf.var("pm").put_attr("units","meter-1");
        ncf.var("pm").put_attr("grid","grid");
        ncf.var("pm").put_attr("location","face");
        ncf.var("pm").put_attr("coordinates","x_rho y_rho");
        ncf.var("pm").put_attr("field","pm, scalar");
 
        ncf.def_var("pn",ncutils::NCDType::Real, {ny_r_name, nx_r_name});
        ncf.var("pn").put_attr("long_name","curvilinear coordinate metric in ETA");
        ncf.var("pn").put_attr("units","meter-1");
        ncf.var("pn").put_attr("grid","grid");
        ncf.var("pn").put_attr("location","face");
        ncf.var("pn").put_attr("coordinates","x_rho y_rho");
        ncf.var("pn").put_attr("field","pn, scalar");
 
        ncf.def_var("f",ncutils::NCDType::Real, {ny_r_name, nx_r_name});
        ncf.var("f").put_attr("long_name","Coriolis parameter at RHO-points");
        ncf.var("f").put_attr("units","second-1");
        ncf.var("f").put_attr("grid","grid");
        ncf.var("f").put_attr("location","face");
        ncf.var("f").put_attr("coordinates","x_rho y_rho");
        ncf.var("f").put_attr("field","coriolis, scalar");
 
        ncf.def_var("x_rho",ncutils::NCDType::Real, {ny_r_name, nx_r_name});
        ncf.var("x_rho").put_attr("long_name","x-locations of RHO-points");
        ncf.var("x_rho").put_attr("units","meter");
        ncf.var("x_rho").put_attr("field","x_rho, scalar");
 
        ncf.def_var("y_rho",ncutils::NCDType::Real, {ny_r_name, nx_r_name});
        ncf.var("y_rho").put_attr("long_name","y-locations of RHO-points");
        ncf.var("y_rho").put_attr("units","meter");
        ncf.var("y_rho").put_attr("field","y_rho, scalar");
 
        ncf.def_var("x_u",ncutils::NCDType::Real, {ny_u_name, nx_u_name});
        ncf.var("x_u").put_attr("long_name","x-locations of U-points");
        ncf.var("x_u").put_attr("units","meter");
        ncf.var("x_u").put_attr("field","x_u, scalar");
 
        ncf.def_var("y_u",ncutils::NCDType::Real, {ny_u_name, nx_u_name});
        ncf.var("y_u").put_attr("long_name","y-locations of U-points");
        ncf.var("y_u").put_attr("units","meter");
        ncf.var("y_u").put_attr("field","y_u, scalar");
 
        ncf.def_var("x_v",ncutils::NCDType::Real, {ny_v_name, nx_v_name});
        ncf.var("x_v").put_attr("long_name","x-locations of V-points");
        ncf.var("x_v").put_attr("units","meter");
        ncf.var("x_v").put_attr("field","x_v, scalar");
 
        ncf.def_var("y_v",ncutils::NCDType::Real, {ny_v_name, nx_v_name});
        ncf.var("y_v").put_attr("long_name","y-locations of V-points");
        ncf.var("y_v").put_attr("units","meter");
        ncf.var("y_v").put_attr("field","y_v, scalar");
 
        ncf.def_var("x_psi",ncutils::NCDType::Real, {ny_p_name, nx_p_name});
        ncf.var("x_psi").put_attr("long_name","x-locations of PSI-points");
        ncf.var("x_psi").put_attr("units","meter");
        ncf.var("x_psi").put_attr("field","x_psi, scalar");
 
        ncf.def_var("y_psi",ncutils::NCDType::Real, {ny_p_name, nx_p_name});
        ncf.var("y_psi").put_attr("long_name","y-locations of PSI-points");
        ncf.var("y_psi").put_attr("units","meter");
        ncf.var("y_psi").put_attr("field","y_psi, scalar");
 
        // Both time stamps run on the model clock, so both carry the reference
        // date remora.time_ref names, as ROMS writes Rclock%string.
        const std::string ref_date_string = remora_ref_date_string(solverChoice.time_ref);
        const std::string ref_calendar = remora_ref_calendar(solverChoice.time_ref);
 
        // Double regardless of Real, like dstart: seconds since the reference
        // date outgrow a float's resolution within days.
        ncf.def_var("ocean_time", NC_DOUBLE, { nt_name });
        ncf.var("ocean_time").put_attr("long_name","time since initialization");
        ncf.var("ocean_time").put_attr("units","seconds since " + ref_date_string);
        ncf.var("ocean_time").put_attr("calendar",ref_calendar);
        ncf.var("ocean_time").put_attr("field","time, scalar, series");
 
        // ROMS DSTART: when the run starts, in days on the model clock. Double
        // regardless of Real -- days since year 1 needs more than a float's
        // seven significant digits, which would round this to the nearest hour.
        ncf.def_var("dstart", NC_DOUBLE, {});
        ncf.var("dstart").put_attr("long_name","time stamp assigned to model initialization");
        ncf.var("dstart").put_attr("units","days since " + ref_date_string);
        ncf.var("dstart").put_attr("calendar",ref_calendar);
 
        ncf.def_var("Cs_r", ncutils::NCDType::Real, {nz_r_name});
        ncf.var("Cs_r").put_attr("long_name", "S-coordinate stretching curves at RHO points");
        ncf.var("Cs_r").put_attr("valid_min",std::vector({-one}));
        ncf.var("Cs_r").put_attr("valid_max",std::vector({zero}));
        ncf.var("Cs_r").put_attr("field","Cs_r, scalar");
 
        ncf.def_var("Cs_w", ncutils::NCDType::Real, {nz_w_name});
        ncf.var("Cs_w").put_attr("long_name", "S-coordinate stretching curves at W points");
        ncf.var("Cs_w").put_attr("valid_min",std::vector({-one}));
        ncf.var("Cs_w").put_attr("valid_max",std::vector({zero}));
        ncf.var("Cs_w").put_attr("field","Cs_w, scalar");
 
        ncf.def_var("h", ncutils::NCDType::Real, { ny_r_name, nx_r_name });
        ncf.var("h").put_attr("long_name","bathymetry at RHO-points");
        ncf.var("h").put_attr("units","meter");
        ncf.var("h").put_attr("grid","grid");
        ncf.var("h").put_attr("location","face");
        ncf.var("h").put_attr("coordinates","x_rho y_rho");
        ncf.var("h").put_attr("field","bath, scalar");
 
        ncf.def_var_fill("zeta", ncutils::NCDType::Real, { nt_name, ny_r_name, nx_r_name }, &netcdf_fill_value);
        ncf.var("zeta").put_attr("long_name","free-surface");
        ncf.var("zeta").put_attr("units","meter");
        ncf.var("zeta").put_attr("time","ocean_time");
        ncf.var("zeta").put_attr("grid","grid");
        ncf.var("zeta").put_attr("location","face");
        ncf.var("zeta").put_attr("coordinates","x_rho y_rho ocean_time");
        ncf.var("zeta").put_attr("field","free-surface, scalar, series");
 
        for (int n = 0; n < ncons; ++n) {
            int comp = -1;
            for (int i = 0; i < plot_var_names_3d.size(); i++) {
                if (plot_var_names_3d[i] == cons_names[n]) comp = i;
            }
            if (comp >= 0) {
                const std::string& nm = cons_names[n];
                ncf.def_var_fill(nm, ncutils::NCDType::Real, { nt_name, nz_r_name, ny_r_name, nx_r_name }, &netcdf_fill_value);
                if (n == Temp_comp) {
                    ncf.var(nm).put_attr("long_name", "potential temperature");
                    ncf.var(nm).put_attr("units", "Celsius");
                    ncf.var(nm).put_attr("field", "temperature, scalar, series");
                } else if (n == Salt_comp) {
                    ncf.var(nm).put_attr("long_name", "salinity");
                    ncf.var(nm).put_attr("field", "salinity, scalar, series");
                } else {
                    ncf.var(nm).put_attr("long_name", nm);
                    ncf.var(nm).put_attr("field", nm + ", scalar, series");
                }
                ncf.var(nm).put_attr("time", "ocean_time");
                ncf.var(nm).put_attr("grid", "grid");
                ncf.var(nm).put_attr("location", "face");
                ncf.var(nm).put_attr("coordinates", "x_rho y_rho s_rho ocean_time");
            }
        }
 
        {
            int comp = -1;
            for (int i = 0; i < plot_var_names_3d.size(); i++) {
                if (plot_var_names_3d[i] == "vorticity") comp = i;
            }
            if (comp >= 0) {
               ncf.def_var_fill("vorticity", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_r_name, nx_r_name }, &netcdf_fill_value);
               ncf.var("vorticity").put_attr("long_name","vorticity");
               ncf.var("vorticity").put_attr("time","ocean_time");
               ncf.var("vorticity").put_attr("grid","grid");
               ncf.var("vorticity").put_attr("location","face");
               ncf.var("vorticity").put_attr("coordinates","x_rho y_rho s_rho ocean_time");
               ncf.var("vorticity").put_attr("field","vorticity, scalar, series");
            }
        } // end vorticity
 
        if (plot_omega) {
            ncf.def_var_fill("omega", ncutils::NCDType::Real, { nt_name, nz_w_name, ny_r_name, nx_r_name }, &netcdf_fill_value);
            ncf.var("omega").put_attr("long_name","S-coordinate vertical momentum component");
            ncf.var("omega").put_attr("units","meter second-1");
            ncf.var("omega").put_attr("time","ocean_time");
            ncf.var("omega").put_attr("grid","grid");
            ncf.var("omega").put_attr("location","face");
            ncf.var("omega").put_attr("coordinates","x_rho y_rho s_w ocean_time");
            ncf.var("omega").put_attr("field","omega, scalar, series");
        }
 
        if (plot_wvel) {
            ncf.def_var_fill("w", ncutils::NCDType::Real, { nt_name, nz_w_name, ny_r_name, nx_r_name }, &netcdf_fill_value);
            ncf.var("w").put_attr("long_name","vertical momentum component");
            ncf.var("w").put_attr("units","meter second-1");
            ncf.var("w").put_attr("time","ocean_time");
            ncf.var("w").put_attr("grid","grid");
            ncf.var("w").put_attr("location","face");
            ncf.var("w").put_attr("coordinates","x_rho y_rho s_w ocean_time");
            ncf.var("w").put_attr("field","w-velocity, scalar, series");
        }
 
        // Output 2D horizontal mixing coefficients if using scaled_to_grid option
        if (solverChoice.horiz_mixing_type == HorizMixingType::scaled_to_grid) {
            ncf.def_var_fill("visc2", ncutils::NCDType::Real, { ny_r_name, nx_r_name }, &netcdf_fill_value);
            ncf.var("visc2").put_attr("long_name","horizontal harmonic viscosity coefficient at RHO-points");
            ncf.var("visc2").put_attr("units","meter2 second-1");
            ncf.var("visc2").put_attr("grid","grid");
            ncf.var("visc2").put_attr("location","face");
            ncf.var("visc2").put_attr("coordinates","x_rho y_rho");
            ncf.var("visc2").put_attr("field","visc2, scalar");
 
            for (int n = 0; n < ncons; ++n) {
                const std::string nm = std::string("diff2_") + cons_names[n];
                ncf.def_var_fill(nm, ncutils::NCDType::Real, { ny_r_name, nx_r_name }, &netcdf_fill_value);
                ncf.var(nm).put_attr("long_name", std::string("horizontal harmonic diffusivity coefficient for ") + cons_names[n] + " at RHO-points");
                ncf.var(nm).put_attr("units","meter2 second-1");
                ncf.var(nm).put_attr("grid","grid");
                ncf.var(nm).put_attr("location","face");
                ncf.var(nm).put_attr("coordinates","x_rho y_rho");
                ncf.var(nm).put_attr("field", nm + ", scalar");
            }
        }
 
        ncf.def_var_fill("u", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_u_name, nx_u_name }, &netcdf_fill_value);
        ncf.var("u").put_attr("long_name","u-momentum component");
        ncf.var("u").put_attr("units","meter second-1");
        ncf.var("u").put_attr("time","ocean_time");
        ncf.var("u").put_attr("grid","grid");
        ncf.var("u").put_attr("location","edge1");
        ncf.var("u").put_attr("coordinates","x_u y_u s_rho ocean_time");
        ncf.var("u").put_attr("field","u-velocity, scalar, series");
 
        ncf.def_var_fill("v", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_v_name, nx_v_name }, &netcdf_fill_value);
        ncf.var("v").put_attr("long_name","v-momentum component");
        ncf.var("v").put_attr("units","meter second-1");
        ncf.var("v").put_attr("time","ocean_time");
        ncf.var("v").put_attr("grid","grid");
        ncf.var("v").put_attr("location","edge2");
        ncf.var("v").put_attr("coordinates","x_v y_v s_rho ocean_time");
        ncf.var("v").put_attr("field","v-velocity, scalar, series");
 
        ncf.def_var_fill("ubar", ncutils::NCDType::Real, { nt_name, ny_u_name, nx_u_name }, &netcdf_fill_value);
        ncf.var("ubar").put_attr("long_name","vertically integrated u-momentum component");
        ncf.var("ubar").put_attr("units","meter second-1");
        ncf.var("ubar").put_attr("time","ocean_time");
        ncf.var("ubar").put_attr("grid","grid");
        ncf.var("ubar").put_attr("location","edge1");
        ncf.var("ubar").put_attr("coordinates","x_u y_u ocean_time");
        ncf.var("ubar").put_attr("field","ubar-velocity, scalar, series");
 
        ncf.def_var_fill("vbar", ncutils::NCDType::Real, { nt_name, ny_v_name, nx_v_name }, &netcdf_fill_value);
        ncf.var("vbar").put_attr("long_name","vertically integrated v-momentum component");
        ncf.var("vbar").put_attr("units","meter second-1");
        ncf.var("vbar").put_attr("time","ocean_time");
        ncf.var("vbar").put_attr("grid","grid");
        ncf.var("vbar").put_attr("location","edge2");
        ncf.var("vbar").put_attr("coordinates","x_v y_v ocean_time");
        ncf.var("vbar").put_attr("field","vbar-velocity, scalar, series");
 
        ncf.def_var("sustr", ncutils::NCDType::Real, { nt_name, ny_u_name, nx_u_name });
        ncf.var("sustr").put_attr("long_name","surface u-momentum stress");
        ncf.var("sustr").put_attr("units","newton meter-2");
        ncf.var("sustr").put_attr("time","ocean_time");
        ncf.var("sustr").put_attr("grid","grid");
        ncf.var("sustr").put_attr("location","edge1");
        ncf.var("sustr").put_attr("coordinates","x_u y_u ocean_time");
        ncf.var("sustr").put_attr("field","surface u-momentum stress, scalar, series");
 
        ncf.def_var("svstr", ncutils::NCDType::Real, { nt_name, ny_v_name, nx_v_name });
        ncf.var("svstr").put_attr("long_name","surface v-momentum stress");
        ncf.var("svstr").put_attr("units","newton meter-2");
        ncf.var("svstr").put_attr("time","ocean_time");
        ncf.var("svstr").put_attr("grid","grid");
        ncf.var("svstr").put_attr("location","edge2");
        ncf.var("svstr").put_attr("coordinates","x_v y_v ocean_time");
        ncf.var("svstr").put_attr("field","surface v-momentum stress, scalar, series");
 
        ncf.def_var("mask_rho", ncutils::NCDType::Real, { nt_name, ny_r_name, nx_r_name });
        ncf.var("mask_rho").put_attr("long_name","mask on RHO-points");
        ncf.var("mask_rho").put_attr("time","ocean_time");
        ncf.var("mask_rho").put_attr("flag_values",std::vector({Real(0.0),Real(1.0)}));
        ncf.var("mask_rho").put_attr("flag_meanings","land water");
 
        ncf.def_var("mask_u", ncutils::NCDType::Real, { nt_name, ny_u_name, nx_u_name });
        ncf.var("mask_u").put_attr("long_name","mask on U-points");
        ncf.var("mask_u").put_attr("time","ocean_time");
        ncf.var("mask_u").put_attr("flag_values",std::vector({Real(0.0),Real(1.0)}));
        ncf.var("mask_u").put_attr("flag_meanings","land water");
 
        ncf.def_var("mask_v", ncutils::NCDType::Real, { nt_name, ny_v_name, nx_v_name });
        ncf.var("mask_v").put_attr("long_name","mask on V-points");
        ncf.var("mask_v").put_attr("time","ocean_time");
        ncf.var("mask_v").put_attr("flag_values",std::vector({Real(0.0),Real(1.0)}));
        ncf.var("mask_v").put_attr("flag_meanings","land water");
 
        // Surface tracer fluxes, one per cell-centered tracer, matching what the AMReX
        // plotfile writer already offers and requested through the same remora.plot_vars_2d
        // key. vec_stflux is ncons wide and unconditionally allocated, so any tracer can be
        // named. Components past temp and salt read zero unless something fills them --
        // Fennel's air-sea gas exchange acts on the tracer directly rather than through this
        // array -- so they show what the solver actually applied, which may be nothing.
        for (int n = 0; n < ncons; ++n) {
            const std::string nm = std::string("stflux_") + cons_names[n];
            if (!plot_2d_var_requested(plot_var_names_2d, nm)) { continue; }
            ncf.def_var_fill(nm, ncutils::NCDType::Real, { nt_name, ny_r_name, nx_r_name }, &netcdf_fill_value);
            ncf.var(nm).put_attr("long_name", std::string("surface flux of ") + cons_names[n]);
            // Kinematic, as ssflux is: the tracer's own units times meter second-1.
            ncf.var(nm).put_attr("units", (n == Temp_comp) ? "Celsius meter second-1"
                                                           : "meter second-1");
            ncf.var(nm).put_attr("time","ocean_time");
            ncf.var(nm).put_attr("grid","grid");
            ncf.var(nm).put_attr("location","face");
            ncf.var(nm).put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var(nm).put_attr("field", nm + ", scalar, series");
        }
 
        if (solverChoice.output_forcing) {
            // Surface air temperature (Celsius)
            ncf.def_var("Tair", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("Tair").put_attr("long_name","surface air temperature");
            ncf.var("Tair").put_attr("units","Celsius");
            ncf.var("Tair").put_attr("time","ocean_time");
            ncf.var("Tair").put_attr("grid","grid");
            ncf.var("Tair").put_attr("location","face");
            ncf.var("Tair").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("Tair").put_attr("field","Tair, scalar, series");
 
            // Surface air pressure (Pascal)
            ncf.def_var("Pair", ncutils::NCDType::Real,{ nt_name, ny_r_name, nx_r_name });
            ncf.var("Pair").put_attr("long_name","surface air pressure");
            ncf.var("Pair").put_attr("units","Pascal");
            ncf.var("Pair").put_attr("time","ocean_time");
            ncf.var("Pair").put_attr("grid","grid");
            ncf.var("Pair").put_attr("location","face");
            ncf.var("Pair").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("Pair").put_attr("field","Pair, scalar, series");
 
            // Surface net heat flux (W/m2)
            ncf.def_var("qnet", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("qnet").put_attr("long_name","surface net heat flux");
            ncf.var("qnet").put_attr("units","watt meter-2");
            ncf.var("qnet").put_attr("time","ocean_time");
            ncf.var("qnet").put_attr("grid","grid");
            ncf.var("qnet").put_attr("location","face");
            ncf.var("qnet").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("qnet").put_attr("field","surface heat flux, scalar, series");
 
            // Surface net salt flux (kinematic)
            ncf.def_var("ssflux", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("ssflux").put_attr("long_name","kinematic surface net salt flux, SALT*(E-P)/rhow");
            ncf.var("ssflux").put_attr("units","meter second-1");
            ncf.var("ssflux").put_attr("time","ocean_time");
            ncf.var("ssflux").put_attr("grid","grid");
            ncf.var("ssflux").put_attr("location","face");
            ncf.var("ssflux").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("ssflux").put_attr("field","surface net salt flux, scalar, series");
 
            // Latent heat flux (W/m2)
            ncf.def_var("latent", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("latent").put_attr("long_name","net latent heat flux");
            ncf.var("latent").put_attr("units","watt meter-2");
            ncf.var("latent").put_attr("time","ocean_time");
            ncf.var("latent").put_attr("grid","grid");
            ncf.var("latent").put_attr("location","face");
            ncf.var("latent").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("latent").put_attr("field","latent heat flux, scalar, series");
 
            // Sensible heat flux (W/m2)
            ncf.def_var("sensible", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("sensible").put_attr("long_name","net sensible heat flux");
            ncf.var("sensible").put_attr("units","watt meter-2");
            ncf.var("sensible").put_attr("time","ocean_time");
            ncf.var("sensible").put_attr("grid","grid");
            ncf.var("sensible").put_attr("location","face");
            ncf.var("sensible").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("sensible").put_attr("field","sensible heat flux, scalar, series");
 
            // Longwave radiation (W/m2)
            ncf.def_var("lwrad", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("lwrad").put_attr("long_name","net longwave radiation flux");
            ncf.var("lwrad").put_attr("units","watt meter-2");
            ncf.var("lwrad").put_attr("time","ocean_time");
            ncf.var("lwrad").put_attr("grid","grid");
            ncf.var("lwrad").put_attr("location","face");
            ncf.var("lwrad").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("lwrad").put_attr("field","longwave radiation, scalar, series");
 
            // Shortwave radiation (W/m2)
            ncf.def_var("swrad", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("swrad").put_attr("long_name","solar shortwave radiation flux");
            ncf.var("swrad").put_attr("units","watt meter-2");
            ncf.var("swrad").put_attr("time","ocean_time");
            ncf.var("swrad").put_attr("grid","grid");
            ncf.var("swrad").put_attr("location","face");
            ncf.var("swrad").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("swrad").put_attr("field","shortwave radiation, scalar, series");
 
            // Evaporation rate (kg m-2 s-1)
            ncf.def_var("evaporation", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("evaporation").put_attr("long_name","evaporation rate");
            ncf.var("evaporation").put_attr("units","kilogram meter-2 second-1");
            ncf.var("evaporation").put_attr("time","ocean_time");
            ncf.var("evaporation").put_attr("grid","grid");
            ncf.var("evaporation").put_attr("location","face");
            ncf.var("evaporation").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("evaporation").put_attr("field","evaporation, scalar, series");
 
            // Rain rate (kg m-2 s-1)
            ncf.def_var("rain", ncutils::NCDType::Real, {nt_name, ny_r_name, nx_r_name });
            ncf.var("rain").put_attr("long_name","rain fall rate");
            ncf.var("rain").put_attr("units","kilogram meter-2 second-1");
            ncf.var("rain").put_attr("time","ocean_time");
            ncf.var("rain").put_attr("grid","grid");
            ncf.var("rain").put_attr("location","face");
            ncf.var("rain").put_attr("coordinates","x_rho y_rho ocean_time");
            ncf.var("rain").put_attr("field","rain, scalar, series");
        }
        // Right now this is hard-wired to {temp, salt, tracer, u, v}
        ncf.put_attr("space_dimension", std::vector<int> { AMREX_SPACEDIM });
//        ncf.put_attr("current_time", std::vector<double> { time });
        ncf.put_attr("CurrentLevel", std::vector<int> { flev });
        ncf.put_attr("DefaultGeometry", std::vector<int> { amrex::DefaultGeometry().Coord() });
 
        ncf.exit_def_mode();
 
        // We are doing single-level writes but it doesn't have to be level 0
        //
        // Write out the header information.
        //
 
        Real dx[AMREX_SPACEDIM];
        for (int i = 0; i < AMREX_SPACEDIM; i++) {
            dx[i] = geom[lev].CellSize()[i];
        }
        const auto *base = geom[lev].ProbLo();
        RealBox rb(subdomain, dx, base);
 
        amrex::Vector<Real> probLo;
        amrex::Vector<Real> probHi;
        for (int i = 0; i < AMREX_SPACEDIM; i++) {
            probLo.push_back(rb.lo(i));
            probHi.push_back(rb.hi(i));
        }
 
        //nc_probLo.par_access(NC_COLLECTIVE);
        // small variable data written by just the master proc
        ncmpi_begin_indep_data(ncf.ncid);
        if (amrex::ParallelDescriptor::IOProcessor()) // only master proc
        {
            auto nc_probLo = ncf.var("probLo");
 
            nc_probLo.put(probLo.data(), { 0 }, { AMREX_SPACEDIM });
 
            auto nc_probHi = ncf.var("probHi");
            //nc_probHi.par_access(NC_COLLECTIVE);
            nc_probHi.put(probHi.data(), { 0 }, { AMREX_SPACEDIM });
 
            amrex::Vector<int> smallend;
            amrex::Vector<int> bigend;
            for (int i = lev; i < flev; i++) {
                smallend.clear();
                bigend.clear();
                for (int j = 0; j < AMREX_SPACEDIM; j++) {
                    smallend.push_back(subdomain.smallEnd(j));
                    bigend.push_back(subdomain.bigEnd(j));
                }
                auto nc_Geom_smallend = ncf.var("Geom.smallend");
                //nc_Geom_smallend.par_access(NC_COLLECTIVE);
                nc_Geom_smallend.put(smallend.data(), { static_cast<long long int>(i - lev), 0 }, { 1,
                AMREX_SPACEDIM });
 
                auto nc_Geom_bigend = ncf.var("Geom.bigend");
                //nc_Geom_bigend.par_access(NC_COLLECTIVE);
                nc_Geom_bigend.put(bigend.data(), { static_cast<long long int>(i - lev), 0 }, { 1,
                AMREX_SPACEDIM });
            }
 
            amrex::Vector<Real> CellSize;
            for (int i = lev; i < flev; i++) {
                CellSize.clear();
                for (Real &j : dx) {
                    CellSize.push_back(amrex::Real(j));
                }
                auto nc_CellSize = ncf.var("CellSize");
                //nc_CellSize.par_access(NC_COLLECTIVE);
                nc_CellSize.put(CellSize.data(), { static_cast<long long int>(i - lev), 0 }, { 1,
                AMREX_SPACEDIM });
            }
            Real hc = solverChoice.tcline;
            Real theta_s = solverChoice.theta_s;
            Real theta_b = solverChoice.theta_b;
            ncf.var("hc").put(&hc);
            ncf.var("theta_s").put(&theta_s);
            ncf.var("theta_b").put(&theta_b);
 
            // remora.start_time in seconds is ROMS DSTART in days.
            double dstart = start_time / 86400.0;
            ncf.var("dstart").put(&dstart);
 
        }
        ncmpi_end_indep_data(ncf.ncid);
 
    } // end if write_header
 
    // Past this point every write is collective. The loops below stage their
    // hyperslabs into the collector rather than writing them, because the number
    // of hyperslabs a rank contributes depends on how many boxes it owns; the
    // single flush() at the end turns each variable into one ncmpi_put_varn_*_all.
    VarnCollector collector;
 
    //
    // We compute the offsets based on location of the box within the domain
    //
    long long local_start_nt = is_history ? history_record : 0;
    long long local_nt = 1; // We write data for only one time
 
    // ocean_time is on the model clock, in double. Written here rather than staged,
    // since the collector holds Real; collective, with only the IO rank contributing.
    {
        const double ocean_time = model_time(t_new[lev]);
        const long long nt_count = amrex::ParallelDescriptor::IOProcessor() ? local_nt : 0;
        ncf.var("ocean_time").put_all(&ocean_time, { local_start_nt }, { nt_count });
    }
 
    // Check whether there are any nans or infs in variables that we will write out
    if (vec_Zt_avg1[lev]->contains_nan() || vec_Zt_avg1[lev]->contains_inf()) {
        amrex::Abort("Found while writing output: zeta contains nan or inf");
    }
    // Check every cell-centered tracer that is actually being written: temperature,
    // salinity, the passive scalars and the biology tracers. plotMF is indexed by
    // position in plot_var_names_3d, not by cons component, so look the name up the
    // same way the write loops below do.
    for (int n = 0; n < ncons; ++n) {
        int comp = -1;
        for (int i = 0; i < plot_var_names_3d.size(); i++) {
            if (plot_var_names_3d[i] == cons_names[n]) comp = i;
        }
        if (comp < 0) { continue; }
        if (plotMF->contains_nan(comp,1) || plotMF->contains_inf(comp,1)) {
            amrex::Abort("Found while writing output: " + cons_names[n] +
                         " contains nan or inf");
        }
    }
    if (xvel_new[lev]->contains_nan() || xvel_new[lev]->contains_inf()) {
        amrex::Abort("Found while writing output: velocity u contains nan or inf");
    }
    if (vec_ubar[lev]->contains_nan(0,1) || vec_ubar[lev]->contains_inf(0,1)) {
        amrex::Abort("Found while writing output: velocity ubar contains nan or inf");
    }
    if (yvel_new[lev]->contains_nan() || yvel_new[lev]->contains_inf()) {
        amrex::Abort("Found while writing output: velocity v contains nan or inf");
    }
    if (vec_vbar[lev]->contains_nan(0,1) || vec_vbar[lev]->contains_inf(0,1)) {
        amrex::Abort("Found while writing output: velocity vbar contains nan or inf");
    }
    if (mf_omega && (mf_omega->contains_nan() || mf_omega->contains_inf())) {
        amrex::Abort("Found while writing output: omega contains nan or inf");
    }
    if (mf_wvel && (mf_wvel->contains_nan() || mf_wvel->contains_inf())) {
        amrex::Abort("Found while writing output: velocity w contains nan or inf");
    }
 
    for (MFIter mfi(*plotMF, false); mfi.isValid(); ++mfi) {
        auto bx = mfi.validbox();
        if (subdomain.contains(bx)) {
            //
            // We only include one grow cell at subdomain boundaries, not internal grid boundaries
            //
            Box tmp_bx(bx);
            if (tmp_bx.smallEnd()[0] == subdomain.smallEnd()[0])
                tmp_bx.growLo(0, 1);
            if (tmp_bx.smallEnd()[1] == subdomain.smallEnd()[1])
                tmp_bx.growLo(1, 1);
            if (tmp_bx.bigEnd()[0] == subdomain.bigEnd()[0])
                tmp_bx.growHi(0, 1);
            if (tmp_bx.bigEnd()[1] == subdomain.bigEnd()[1])
                tmp_bx.growHi(1, 1);
            // amrex::Print() << "    BX " << bx << std::endl;
            // amrex::Print() << "TMP_BX " << tmp_bx << std::endl;
 
            Box tmp_bx_2d(tmp_bx);
            tmp_bx_2d.makeSlab(2, 0);
 
            Box tmp_bx_1d(tmp_bx);
            tmp_bx_1d.makeSlab(0, 0);
            tmp_bx_1d.makeSlab(1, 0);
 
            //
            // These are the dimensions of the data we write for only this box
            //
            long long local_nx = tmp_bx.length()[0];
            long long local_ny = tmp_bx.length()[1];
            long long local_nz = tmp_bx.length()[2];
 
            // We do the "+1" because the offset needs to start at 0
            // Offsets are relative to this level's subdomain: a refined level's boxes
            // start at its own index space, not at zero, while the file is sized to
            // that subdomain.
            long long local_start_x = static_cast<long long>(tmp_bx.smallEnd()[0] - subdomain.smallEnd()[0] + 1);
            long long local_start_y = static_cast<long long>(tmp_bx.smallEnd()[1] - subdomain.smallEnd()[1] + 1);
            long long local_start_z = static_cast<long long>(tmp_bx.smallEnd()[2] - subdomain.smallEnd()[2]);
 
            if (write_header) {
                // Only write out s_rho and s_w at x=0,y=0 to avoid NaNs
                if (bx.contains(subdomain.smallEnd()))
                {
                    {
                        amrex::Vector<amrex::Real> tmp_srho(local_nz);
 
#ifdef AMREX_USE_GPU
                        Gpu::dtoh_memcpy(tmp_srho.data(), s_r.data(), sizeof(amrex::Real)*local_nz);
#else
                        std::memcpy(tmp_srho.data(), s_r.data(), sizeof(amrex::Real)*local_nz);
#endif
                        Gpu::streamSynchronize();
 
                        auto nc_plot_var = collector.var(ncf, "s_rho");
                        //nc_plot_var.par_access(NC_INDEPENDENT);
                        nc_plot_var.put(tmp_srho.data(), { local_start_z }, { local_nz });
                    }
                    {
                        amrex::Vector<amrex::Real> tmp_sw(local_nz+1);
 
#ifdef AMREX_USE_GPU
                        Gpu::dtoh_memcpy(tmp_sw.data(), s_w.data(), sizeof(amrex::Real)*(local_nz+1));
#else
                        std::memcpy(tmp_sw.data(), s_w.data(), sizeof(amrex::Real)*(local_nz+1));
#endif
                        Gpu::streamSynchronize();
 
                        auto nc_plot_var = collector.var(ncf, "s_w");
                        //nc_plot_var.par_access(NC_INDEPENDENT);
                        nc_plot_var.put(tmp_sw.data(), { local_start_z }, { local_nz + 1});
                    }
                    {
                        amrex::Vector<amrex::Real> tmp_Csrho(local_nz);
 
#ifdef AMREX_USE_GPU
                        Gpu::dtoh_memcpy(tmp_Csrho.data(), Cs_r.data(), sizeof(amrex::Real)*(local_nz));
#else
                        std::memcpy(tmp_Csrho.data(), Cs_r.data(), sizeof(amrex::Real)*(local_nz));
#endif
                        Gpu::streamSynchronize();
 
                        auto nc_plot_var = collector.var(ncf, "Cs_r");
                        //nc_plot_var.par_access(NC_INDEPENDENT);
                        nc_plot_var.put(tmp_Csrho.data(), { local_start_z }, { local_nz });
                    }
                    {
                        amrex::Vector<amrex::Real> tmp_Csw(local_nz+1);
 
#ifdef AMREX_USE_GPU
                        Gpu::dtoh_memcpy(tmp_Csw.data(), Cs_w.data(), sizeof(amrex::Real)*(local_nz+1));
#else
                        std::memcpy(tmp_Csw.data(), Cs_w.data(), sizeof(amrex::Real)*(local_nz+1));
#endif
 
                        Gpu::streamSynchronize();
 
                        auto nc_plot_var = collector.var(ncf, "Cs_w");
                        //nc_plot_var.par_access(NC_INDEPENDENT);
                        nc_plot_var.put(tmp_Csw.data(), { local_start_z }, { local_nz + 1});
                    }
                }
 
                {
                    FArrayBox tmp_bathy;
                    tmp_bathy.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    tmp_bathy.template copy<RunOn::Device>((*vec_h[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "h");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
                    nc_plot_var.put(tmp_bathy.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
 
                {
                    FArrayBox tmp_pm;
                    tmp_pm.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    tmp_pm.template copy<RunOn::Device>((*vec_pm[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "pm");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
 
                    nc_plot_var.put(tmp_pm.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
 
                {
                    FArrayBox tmp_pn;
                    tmp_pn.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    tmp_pn.template copy<RunOn::Device>((*vec_pn[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "pn");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
                    nc_plot_var.put(tmp_pn.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
 
                {
                    FArrayBox tmp_f;
                    tmp_f.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    tmp_f.template copy<RunOn::Device>((*vec_fcor[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "f");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
                    nc_plot_var.put(tmp_f.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
 
                {
                    FArrayBox tmp_xr;
                    tmp_xr.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    tmp_xr.template copy<RunOn::Device>((*vec_xr[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "x_rho");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
 
                    nc_plot_var.put(tmp_xr.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
 
                {
                    FArrayBox tmp_yr;
                    tmp_yr.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    tmp_yr.template copy<RunOn::Device>((*vec_yr[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "y_rho");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
                    nc_plot_var.put(tmp_yr.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
            }
 
            {
                FArrayBox tmp_zeta;
                tmp_zeta.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp_zeta.template copy<RunOn::Device>((*vec_Zt_avg1[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "zeta");
                nc_plot_var.put(tmp_zeta.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny,
                        local_nx });
            }
 
            {
                FArrayBox tmp_mask_rho;
                tmp_mask_rho.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp_mask_rho.template copy<RunOn::Device>((*vec_mskr[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "mask_rho");
                nc_plot_var.put(tmp_mask_rho.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny,
                        local_nx });
            }
 
            // stflux_*, one per tracer. Defined above under the same condition.
            for (int n = 0; n < ncons; ++n) {
                const std::string nm = std::string("stflux_") + cons_names[n];
                if (!plot_2d_var_requested(plot_var_names_2d, nm)) { continue; }
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_stflux[lev])[mfi.index()], n, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_var = collector.var(ncf, nm);
                nc_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x },
                                          { local_nt,       local_ny,       local_nx });
            }
 
            if (solverChoice.output_forcing)
            {
                const Real Hscale = solverChoice.rho0 * Cp;
                // The copy and the mult below are both async on the same stream, so the
                // sync has to follow the last of them and precede the host-side put().
                // Tair
                {
                    FArrayBox tmp_Tair;
                    tmp_Tair.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp_Tair.template copy<RunOn::Device>((*vec_Tair[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "Tair");
                    nc_plot_var.put(tmp_Tair.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
                }
                // Pair
                {
                    FArrayBox tmp_Pair;
                    tmp_Pair.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp_Pair.template copy<RunOn::Device>((*vec_Pair[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "Pair");
                    nc_plot_var.put(tmp_Pair.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
                }
                // qnet  (stored °C m/s → write W/m²)
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    // Copy stflux Temp component
                    tmp.template copy<RunOn::Device>(
                        (*vec_stflux[lev])[mfi.index()],
                        Temp_comp,  // source component
                        0,          // dest component
                        1           // number of comps
                    );
 
                    // Convert °C·m/s → W/m²
                    tmp.mult<RunOn::Device>(Hscale);
 
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "qnet");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // ssflux = surface net freshwater flux (kg/m²/s converted to m/s)
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
 
                    // Copy stflux Salt component
                    tmp.template copy<RunOn::Device>(
                        (*vec_stflux[lev])[mfi.index()],
                        Salt_comp, // source component
                        0,         // destination component
                        1          // number of components
                    );
 
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "ssflux");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // latent  (stored °C m/s → write W/m²)
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_lhflx[lev])[mfi.index()], 0, 0, 1);
 
                    // Convert °C·m/s → W/m²
                    tmp.mult<RunOn::Device>(Hscale);
 
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "latent");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // sensible  (stored °C m/s → write W/m²)
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_shflx[lev])[mfi.index()], 0, 0, 1);
 
                    // Convert °C·m/s → W/m²
                    tmp.mult<RunOn::Device>(Hscale);
 
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "sensible");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // lwrad  (stored °C m/s → write W/m²)
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_lrflx[lev])[mfi.index()], 0, 0, 1);
 
                    // Convert °C·m/s → W/m²
                    tmp.mult<RunOn::Device>(Hscale);
 
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "lwrad");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // swrad, note this is stored explicitly as W/m², not degC m/s in REMORA.bulk_flux.cpp
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_srflx[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "swrad");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // evaporation
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_evap[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "evaporation");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
                // rain
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_rain[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "rain");
                    nc_var.put(tmp.dataPtr(),
                            { local_start_nt, local_start_y, local_start_x },
                            { local_nt,       local_ny,       local_nx });
                }
            } // end output forcing
 
            // **************************************************************************
            for (int n = 0; n < ncons; ++n) {
                int comp = -1;
                for (int i = 0; i < plot_var_names_3d.size(); i++) {
                    if (plot_var_names_3d[i] == cons_names[n]) comp = i;
                }
                if (comp >= 0) {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*plotMF)[mfi.index()], comp, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, cons_names[n]);
                    nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_z, local_start_y, local_start_x }, { local_nt,
                            local_nz, local_ny, local_nx });
                }
            }
            // **************************************************************************
 
            // **************************************************************************
            { // Vorticity
                int comp = -1;
                for (int i = 0; i < plot_var_names_3d.size(); i++) {
                    if (plot_var_names_3d[i] == "vorticity") comp = i;
                }
                if (comp >= 0) {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*plotMF)[mfi.index()], comp, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, plot_var_names_3d[comp]);
                    nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_z, local_start_y, local_start_x }, { local_nt,
                            local_nz, local_ny, local_nx });
                } // if vorticity exists in plotMF
            } // end vorticity
            // **************************************************************************
 
            // **************************************************************************
            // omega and w, on the same columns as the tracers but at the nz+1 w points
            // **************************************************************************
            {
                Box tmp_bx_w(tmp_bx);
                tmp_bx_w.surroundingNodes(2);
                std::pair<const char*, MultiFab const*> w_vars[] = {{"omega", mf_omega}, {"w", mf_wvel}};
                for (const auto& wv : w_vars) {
                    if (wv.second == nullptr) { continue; }
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_w, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*wv.second)[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, wv.first);
                    nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_z, local_start_y, local_start_x }, { local_nt,
                            local_nz+1, local_ny, local_nx });
                }
            } // end omega and w
            // **************************************************************************
 
            // **************************************************************************
            // Horizontal mixing coefficients (scaled_to_grid):
            // vertically homogeneous and time-invariant -> write static 2D fields once.
            // **************************************************************************
            if ((solverChoice.horiz_mixing_type == HorizMixingType::scaled_to_grid) && write_header) {
                { // visc2
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_visc2_r[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, "visc2");
                    nc_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
 
                // diff2_*
                for (int n = 0; n < ncons; ++n) {
                    const std::string nm = std::string("diff2_") + cons_names[n];
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_diff2[lev])[mfi.index()], n, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_var = collector.var(ncf, nm);
                    nc_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
            }
 
        } // subdomain
    } // mfi
 
    // Flush at each loop boundary rather than once at the end, so the collector
    // only ever holds one loop's worth of staged data. All ranks reach these
    // points, which is what the collective call requires.
    collector.flush(ncf);
 
    // Writing u (we loop over cons to get cell-centered box)
    for (MFIter mfi(*plotMF, false); mfi.isValid(); ++mfi) {
        Box bx = mfi.validbox();
 
        if (subdomain.contains(bx)) {
            //
            // We only include one grow cell at subdomain boundaries, not internal grid boundaries
            //
            Box tmp_bx(bx);
            tmp_bx.surroundingNodes(0);
            if (tmp_bx.smallEnd()[1] == subdomain.smallEnd()[1])
                tmp_bx.growLo(1, 1);
            if (tmp_bx.bigEnd()[1] == subdomain.bigEnd()[1])
                tmp_bx.growHi(1, 1);
            Box tmp_bx_2d(tmp_bx);
            tmp_bx_2d.makeSlab(2, 0);
 
            //
            // These are the dimensions of the data we write for only this box
            //
            long long local_nx = tmp_bx.length()[0];
            long long local_ny = tmp_bx.length()[1];
            long long local_nz = tmp_bx.length()[2];
 
            // We do the "+1" because the offset needs to start at 0
            long long local_start_x = static_cast<long long>(tmp_bx.smallEnd()[0] - subdomain.smallEnd()[0]);
            long long local_start_y = static_cast<long long>(tmp_bx.smallEnd()[1] - subdomain.smallEnd()[1] + 1);
            long long local_start_z = static_cast<long long>(tmp_bx.smallEnd()[2] - subdomain.smallEnd()[2]);
 
            if (write_header) {
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_xu[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "x_u");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
                    nc_plot_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
                {
                    FArrayBox tmp;
                    tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                    tmp.template copy<RunOn::Device>((*vec_yu[lev])[mfi.index()], 0, 0, 1);
                    Gpu::streamSynchronize();
 
                    auto nc_plot_var = collector.var(ncf, "y_u");
                    //nc_plot_var.par_access(NC_INDEPENDENT);
                    nc_plot_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
            }
 
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*xvel_new[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "u");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_z, local_start_y, local_start_x }, { local_nt,
                        local_nz, local_ny, local_nx });
            }
 
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_ubar[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "ubar");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
            }
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_sustr[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "sustr");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
            }
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_msku[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "mask_u");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
            }
        } // in subdomain
    } // mfi
 
    collector.flush(ncf);
 
    // Writing v (we loop over cons to get cell-centered box)
    for (MFIter mfi(*plotMF, false); mfi.isValid(); ++mfi) {
        Box bx = mfi.validbox();
 
        if (subdomain.contains(bx)) {
            //
            // We only include one grow cell at subdomain boundaries, not internal grid boundaries
            //
            Box tmp_bx(bx);
            tmp_bx.surroundingNodes(1);
            if (tmp_bx.smallEnd()[0] == subdomain.smallEnd()[0])
                tmp_bx.growLo(0, 1);
            if (tmp_bx.bigEnd()[0] == subdomain.bigEnd()[0])
                tmp_bx.growHi(0, 1);
            // amrex::Print() << "    BX " << bx << std::endl;
            // amrex::Print() << "TMP_BX " << tmp_bx << std::endl;
 
            Box tmp_bx_2d(tmp_bx);
            tmp_bx_2d.makeSlab(2, 0);
 
            //
            // These are the dimensions of the data we write for only this box
            //
            long long local_nx = tmp_bx.length()[0];
            long long local_ny = tmp_bx.length()[1];
            long long local_nz = tmp_bx.length()[2];
 
            // We do the "+1" because the offset needs to start at 0
            long long local_start_x = static_cast<long long>(tmp_bx.smallEnd()[0] - subdomain.smallEnd()[0] + 1);
            long long local_start_y = static_cast<long long>(tmp_bx.smallEnd()[1] - subdomain.smallEnd()[1]);
            long long local_start_z = static_cast<long long>(tmp_bx.smallEnd()[2] - subdomain.smallEnd()[2]);
 
            if (write_header) {
                {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_xv[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "x_v");
                //nc_plot_var.par_access(NC_INDEPENDENT);
                nc_plot_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
                {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_yv[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "y_v");
                //nc_plot_var.par_access(NC_INDEPENDENT);
                nc_plot_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
                }
            }
 
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*yvel_new[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "v");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_z, local_start_y, local_start_x }, { local_nt,
                        local_nz, local_ny, local_nx });
            }
 
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_vbar[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "vbar");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
            }
 
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_svstr[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "svstr");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
            }
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_mskv[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "mask_v");
                nc_plot_var.put(tmp.dataPtr(), { local_start_nt, local_start_y, local_start_x }, { local_nt, local_ny, local_nx });
            }
 
        } // in subdomain
    } // mfi
 
    collector.flush(ncf);
 
    for (MFIter mfi(*plotMF, false); mfi.isValid(); ++mfi) {
        Box bx = mfi.validbox();
 
        if (subdomain.contains(bx)) {
            //
            // We only include one grow cell at subdomain boundaries, not internal grid boundaries
            //
            Box tmp_bx(bx);
            tmp_bx.surroundingNodes(0);
            tmp_bx.surroundingNodes(1);
 
            Box tmp_bx_2d(tmp_bx);
            tmp_bx_2d.makeSlab(2, 0);
 
            //
            // These are the dimensions of the data we write for only this box
            //
            long long local_nx = tmp_bx.length()[0];
            long long local_ny = tmp_bx.length()[1];
 
            // We do the "+1" because the offset needs to start at 0
            long long local_start_x = static_cast<long long>(tmp_bx.smallEnd()[0] - subdomain.smallEnd()[0]);
            long long local_start_y = static_cast<long long>(tmp_bx.smallEnd()[1] - subdomain.smallEnd()[1]);
 
            if (write_header) {
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_xp[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "x_psi");
                //nc_plot_var.par_access(NC_INDEPENDENT);
                nc_plot_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
            }
            {
                FArrayBox tmp;
                tmp.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena());
                tmp.template copy<RunOn::Device>((*vec_yp[lev])[mfi.index()], 0, 0, 1);
                Gpu::streamSynchronize();
 
                auto nc_plot_var = collector.var(ncf, "y_psi");
                //nc_plot_var.par_access(NC_INDEPENDENT);
                nc_plot_var.put(tmp.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx });
            }
 
            } // header
        } // in subdomain
    } // mfi
 
    // One collective ncmpi_put_varn_*_all per variable. Collective writes keep
    // numrecs synced in the file header as they go, so no separate
    // ncmpi_end_indep_data() sync is needed.
    collector.flush(ncf);
 
    ncf.close();
 
    REMORA::total_nc_plot_file_step += 1;
}