Extending the Model

This guide covers three extension points: nodal attributes (spatially varying parameters), physics contributions (terms in the GWCE and momentum equations), and meteorological file formats.

Adding a Nodal Attribute 

Nodal attributes are per-node fields loaded from the mesh file or set to a constant value from configuration. Examples include Manning’s n and the primitive weighting parameter tau0.

All nodal attributes are accessed via the NodalAttributeId enum. Strings appear only at the I/O boundary (reading from NetCDF); the internal API is entirely enum-based for O(1) access with no string hashing at runtime.

Step 1: Register the Enum 

Add an entry to NodalAttributeId in geometry/NodalAttribute.hpp:

enum class NodalAttributeId : int {
  ManningsN = 0,
  Tau0 = 1,
  // ... the other existing attributes ...
  AdvectionState = 6,
  AppleSauce = 7,  // New attribute
};

/// Number of defined nodal attributes (update when adding new entries)
static constexpr size_t NUM_NODAL_ATTRIBUTES = 8;
static_assert(static_cast<int>(NodalAttributeId::AppleSauce) + 1 ==
                  NUM_NODAL_ATTRIBUTES,
              "Update NUM_NODAL_ATTRIBUTES when adding new attributes");

Add the string mapping in attribute_name() in the same file:

constexpr auto attribute_name(const NodalAttributeId id) -> const char* {
  switch (id) {
    case NodalAttributeId::ManningsN:
      return "mannings_n_at_sea_floor";
    case NodalAttributeId::Tau0:
      return "primitive_weighting_in_continuity_equation";
    // ... the other existing attributes ...
    case NodalAttributeId::AppleSauce:
      return "apple_sauce_coefficient";
  }
  Kokkos::abort("attribute_name: unhandled NodalAttributeId");
  return "";
}

The switch has no default, so -Wswitch flags a missing case at compile time. The string returned by attribute_name() must match the NetCDF variable name in the mesh file.

Step 2: Initialize the Attribute 

Add an initialization method to io/input/NodalAttributeInitializer.cpp. The existing initialize_scalar_attribute() helper handles the common pattern of reading from mesh or using a constant:

void NodalAttributeInitializer::initialize_apple_sauce() {
  constexpr double default_value = 1.0;
  constexpr double valid_min = 0.0;
  constexpr double valid_max = 100.0;

  initialize_scalar_attribute(
      Geometry::NodalAttributeId::AppleSauce,
      m_config.physics.apple_sauce_coeff, default_value,
      "m", "apple sauce coefficient", valid_min, valid_max);
}

Call it from initialize():

void NodalAttributeInitializer::initialize() {
  m_mesh->initialize_attributes(m_fields.num_nodes());
  initialize_manning_n();
  initialize_tau0();
  // ... the other attributes ...
  initialize_advection_state();
  initialize_apple_sauce();  // New
  initialize_initial_water_level();
  initialize_wetdry_state();
}

Step 3: Access in Kernels 

Access the attribute in any kernel via the enum:

const auto apple_sauce =
    mesh->attributes().get(Geometry::NodalAttributeId::AppleSauce);

Kokkos::parallel_for(..., KOKKOS_LAMBDA(const int i) {
  const auto val = apple_sauce(i);  // Works for constant or per-node
});

get() returns a const NodalAttributeAccessor& from a pre-built cache array: a single array index, no lookup. If the attribute was not registered, get() throws.

Step 4: Add Configuration 

Add the configuration parameter to the appropriate config struct in core/config/PhysicsConfig.hpp and parse it in core/config/ModelConfigurationFactory.cpp. A value that is either a constant or read from the mesh is a std::optional<double>: parse_mesh_or_numeric() turns the string "mesh" into std::nullopt (read per-node values from the mesh file) and a number into the constant. Add the key to core/config/cocoa_config.schema.yaml as well (mesh_or_numeric); the schema rejects keys it does not list.

Step 5: Update Tests 

Register the new attribute in test fixtures that exercise the kernels using it:

mesh->attributes().register_constant(
    NodalAttributeId::AppleSauce, 1.0, "m");

Vector (Multi-Component) Attributes 

Some attributes have multiple values per node (e.g., internal tide friction stores 3 tensor components). These use vector storage. Read 2D variables [num_nodes, num_components] with NodalAttributeReader::read_vector_attribute() or read_vector_attribute_distributed():

auto attr = m_mesh->is_parallel()
    ? NodalAttributeReader::read_vector_attribute_distributed(
          reader, variable_name, m_mesh->partition(), 0.0, "1/s")
    : NodalAttributeReader::read_vector_attribute(
          reader, variable_name, num_nodes, 0.0, "1/s");

The two-argument accessor returns one component:

const auto accessor =
    mesh->attributes().get(NodalAttributeId::MyVectorAttr);
Kokkos::parallel_for(..., KOKKOS_LAMBDA(const int i) {
  const auto xx = accessor(i, 0);  // First component
  const auto yy = accessor(i, 1);  // Second component
  const auto xy = accessor(i, 2);  // Third component
});

If an attribute may be either scalar (1D in the mesh file) or vector (2D), detect the shape with reader.get_variable_dimensions() and register accordingly, as initialize_internal_tide_friction() does.

Optional Attributes 

An optional attribute applies only when the user enables it. Skip its registration when it is disabled, and check has() before accessing it:

// In initialization: skip if disabled
if (config_value == 0.0) {
  return;  // Don't register -> has() returns false
}

// In kernel: check before applying
if (!mesh->attributes().has(NodalAttributeId::MyOptionalAttr)) {
  return;
}
const auto accessor = mesh->attributes().get(NodalAttributeId::MyOptionalAttr);
// ... apply

See internal_tide_friction for a complete example of an optional vector attribute.

Adding a Physics Contribution 

Physics contributions add terms to the GWCE (continuity) and momentum equation solvers. The framework uses C++20 concepts to enforce compile-time contracts while keeping physics math out of the assembly kernel lambdas.

Architecture 

GWCE contributions use nested Input structs constructed inline at the call site. Each contribution’s Input holds only the data that contribution needs, and that data dies when compute() returns, which keeps peak register pressure on GPU down.

Momentum contributions share one MomentumElementContextT<Elem> struct (MomentumElementContext is the triangle alias).

Terms are written once for both element types: a term whose only dependence on the element type is its vertex count is templated on the element-type tag (Geometry::Triangle, Geometry::Quadrilateral), as GwceTau0T<Elem> is, and each element kernel instantiates it.

Fig. 55 Physics contribution architecture for GWCE and Momentum

Concepts

Three C++20 concepts in numeric/cg/terms/TermsConcepts.hpp define the contracts. The momentum concepts are checked per element type Elem:

Concept	Required Signature	Use Case
`GwcePhysicsTerm`	`static compute(const T::Input&)` `-> Types::Vec2<double>`	Returns (jx, jy) contribution to both x- and y-equation GWCE RHS
`MomentumElementTerm<T, Elem>`	`static compute(const MomentumElementContextT<Elem>&)` `-> Types::Vec2<double>`	Returns one (x, y) contribution that reaches every vertex of the element
`MomentumNodalTerm<T, Elem>`	`static compute(const MomentumElementContextT<Elem>&)` `-> Kokkos::Array<Types::Vec2<double>, Elem::num_nodes>`	Returns per-vertex (x, y) contributions to both x- and y-momentum

Step-by-Step Guide 

Fig. 56 Workflow for adding a new physics contribution

The following steps use a hypothetical physics term called “AppleSauce” as an example.

Step 1: Data Storage

If the new term needs per-node fields (loaded from input or computed each timestep), add a sub-container to data/PhysicsData.hpp:

/// AppleSauce forcing sub-container
struct AppleSauceFields {
  DoubleView coefficient;
  DoubleView intensity;
};

Add a private member, allocate in the constructor, and expose with const/non-const accessors:

// In the public section:
[[nodiscard]] auto apple_sauce() -> AppleSauceFields& {
  return m_apple_sauce;
}
[[nodiscard]] auto apple_sauce() const -> const AppleSauceFields& {
  return m_apple_sauce;
}

// In the private section:
AppleSauceFields m_apple_sauce;

Follow the pattern of FrictionFields, LateralStressFields, and TidePotentialFields.

Step 2: GWCE Contribution

Create numeric/cg/terms/AppleSauce.hpp. Each GWCE contribution nests an Input struct containing only the data it needs, and provides a static compute(const Input&) method:

#pragma once
#include <Kokkos_Core.hpp>
#include "cocoa_kernel/geometry/ElementGather.hpp"
#include "cocoa_kernel/geometry/ElementType.hpp"
#include "cocoa_kernel/numeric/cg/terms/TermsConcepts.hpp"

namespace Cocoa::Numeric::GwceRhsTerms {

template <typename Elem>
struct GwceAppleSauceT {
  using Pack = Kokkos::Array<double, Elem::num_nodes>;

  struct Input {
    Pack apple_sauce;
    // Add only the fields this contribution needs
  };

  KOKKOS_INLINE_FUNCTION
  static auto compute(const Input& input) -> Types::Vec2<double> {
    const auto as_avg = Geometry::avg(input.apple_sauce);
    // Physics math goes here
    const auto jx = /* ... */;
    const auto jy = /* ... */;
    return {jx, jy};
  }
};

static_assert(GwcePhysicsTerm<GwceAppleSauceT<Geometry::Triangle>>);
static_assert(GwcePhysicsTerm<GwceAppleSauceT<Geometry::Quadrilateral>>);

}  // namespace Cocoa::Numeric::GwceRhsTerms

The static_assert lines check the concept contract for each element type at compile time, so a mismatched signature fails there with a readable error rather than inside a template instantiation at the call site. A term whose quadrilateral form is a different expression from the triangle one (rather than the same expression over four vertices) gets its own struct beside the triangle one (see GwceTau0SpatialGradient and GwceQuadTau0SpatialGradient).

Step 3: Momentum Contribution

Add the momentum term to numeric/cg/terms/AppleSauce.hpp (alongside the GWCE term), in namespace Cocoa::Numeric::MomentumTerms. Choose the concept based on whether the result is the same at every vertex or varies per vertex.

Uniform across vertices: follow MomentumAdvection in GwceAdvection.hpp:

struct MomentumAppleSauce {
  template <typename Elem>
  KOKKOS_INLINE_FUNCTION static auto compute(
      const MomentumElementContextT<Elem>& ctx) -> Types::Vec2<double> {
    // Physics math: same value scattered to every vertex
    return {/* x */, /* y */};
  }
};
static_assert(MomentumElementTerm<MomentumAppleSauce, Geometry::Triangle>);
static_assert(
    MomentumElementTerm<MomentumAppleSauce, Geometry::Quadrilateral>);

Different per vertex: follow MomentumLateralStress in LateralStressTerm.hpp:

struct MomentumAppleSauce {
  template <typename Elem>
  KOKKOS_INLINE_FUNCTION static auto compute(
      const MomentumElementContextT<Elem>& ctx)
      -> Kokkos::Array<Types::Vec2<double>, Elem::num_nodes> {
    // Physics math: different value per vertex
    return {/* one Vec2 per vertex */};
  }
};
static_assert(MomentumNodalTerm<MomentumAppleSauce, Geometry::Triangle>);
static_assert(
    MomentumNodalTerm<MomentumAppleSauce, Geometry::Quadrilateral>);

Step 4: Wire Into GWCE Kernel

In numeric/cg/continuity/GwceVectorAssemblyKernels.hpp, the views the RHS assembly reads are gathered on the host into RhsNodalFields (built by rhs_nodal_fields<Elem>()) and captured by value; the per-element body is assemble_element_rhs<Elem>, which rhs_assembly_kernel runs over each type’s wet list. Make three changes:

a) Include the header and add the view to RhsNodalFields and its builder, alongside the other physics views:

// At file top:
#include "cocoa_kernel/numeric/cg/terms/AppleSauce.hpp"

// In RhsNodalFields:
Types::Scalar1DView apple_sauce;

// In rhs_nodal_fields(), alongside the other physics views:
.apple_sauce = physics.apple_sauce().coefficient,

b) Call the contribution in assemble_element_rhs, constructing the Input inline with aggregate initialization and accumulating into the local result:

using GwceRhsTerms::GwceAppleSauceT;
jxy += GwceAppleSauceT<Elem>::compute(
    {.apple_sauce =
         Geometry::gather_element_nodes(fields.apple_sauce, nodes)});

The result accumulates into the Types::Vec2<double> jxy variable, and the Input temporary (with any data gathered into it) dies after compute() returns.

Step 5: Wire Into Momentum Kernel

In numeric/cg/momentum/MomentumRhsKernels.cpp, the same pattern: the views are gathered into MomentumRhsFields (built by momentum_rhs_fields()), and the per-element body is accumulate_element_momentum_rhs<Elem>, which momentum_rhs_element_kernel launches.

a) Include and add the view to MomentumRhsFields and its builder:

#include "cocoa_kernel/numeric/cg/terms/AppleSauce.hpp"

// In MomentumRhsFields:
Types::Scalar1DView apple_sauce;

b) Load into context. Add a Pack apple_sauce{}; member to MomentumElementContextT (numeric/cg/terms/MomentumElementContext.hpp) and gather it with the others:

ctx.apple_sauce = Geometry::gather_element_nodes(fields.apple_sauce, ctx.nm);

c) Call the contribution. For MomentumElementTerm, add to the base contribution:

base_contrib += MomentumAppleSauce::compute(ctx);

For MomentumNodalTerm, add in the per-vertex scatter loops alongside lateral stress (each element type has its own loop):

const auto apple = MomentumAppleSauce::compute(ctx);
for (size_t i = 0; i < Elem::num_nodes; ++i) {
  const auto contrib =
      base_contrib + lateral[i] + baroclinic_contrib[i] + apple[i];
  Kokkos::atomic_add(&rhs_x(ctx.nm[i]), contrib.x);
  Kokkos::atomic_add(&rhs_y(ctx.nm[i]), contrib.y);
}

Step 6: Populate the Data

How the data gets into AppleSauceFields depends on its nature:

Static nodal attribute (loaded once from mesh file):

Register a NodalAttributeId and use NodalAttributeInitializer as described in Adding a Nodal Attribute above.

Computed each timestep (like lateral stress):

Write a computation kernel that populates AppleSauceFields. Call it from run_preprocessing (numeric/cg/continuity/GwcePreprocessingKernels.hpp), which GwceVectorAssembler::assemble runs before the RHS assembly:

// In run_preprocessing:
COCOA_PROFILE_REGION("Preprocessing");
initialize_rhs_fields(model_fields, rhs_1d, model_fields.num_nodes());
elemental_preprocessing_kernel<Geometry::Triangle>(model_fields, config);
elemental_preprocessing_kernel<Geometry::Quadrilateral>(model_fields, config);
sum_lateral_stress(model_fields, comm_manager);
apple_sauce_kernel(model_fields, config);  // New kernel
nodal_postprocessing_kernel(model_fields);

External forcing (like wind or atmospheric pressure):

Add a forcing source to ForcingManager that updates the fields each timestep using the appropriate ramp factor (m_tidal_ramp for tidal forcing, m_meteorological_ramp for meteorological forcing).

Step 7: Ghost Sum (if needed)

A ghost sum is only needed if the field is assembled from element contributions via atomics. Elements have no ghosts, so each rank holds only its own elements’ share of a node on a partition cut, owned or ghost. Before any rank reads the field, sum every rank’s share into the owner and copy the total back to the ghosts with CommunicationManager::sum_and_distribute_fields. Copying the owner’s value to the ghosts alone (exchange_scalar_fields) is not enough: the owner’s value is itself only a partial sum.

Fields that are directly loaded from nodal attributes or computed per-node (without element scatter) need neither.

Follow sum_lateral_stress (GwcePreprocessingKernels.cpp), which sums the lateral stress components before nodal_postprocessing_kernel normalizes them:

void sum_apple_sauce(Data::ModelFields& model_fields,
                     Distributed::CommunicationManager& comm_manager) {
  std::array fields = {
    Types::StridedScalar1DView(model_fields.physics().apple_sauce().intensity)
  };
  comm_manager.sum_and_distribute_fields(fields);  // no-op in serial
}

Call it after the AppleSauce computation kernel and before anything that normalizes or reads the field.

Note

Pass several fields in one std::array to combine them in a single round-trip pair rather than calling sum_and_distribute_field once per view. max_and_distribute_fields is the exact reduction for flags and bounds any rank may raise.

Step 8: Configuration and Build

Configuration. Add parameters to core/config/PhysicsConfig.hpp:

struct AppleSauceConfig {
  bool active{false};
  double coefficient{1.0};
};

Add AppleSauceConfig apple_sauce; to PhysicsConfig. Parse in core/config/ModelConfigurationFactory.cpp, and add the keys to core/config/cocoa_config.schema.yaml.

CMakeLists.txt. Add new headers to COCOA_KERNEL_HEADERS in src/cocoa_kernel/CMakeLists.txt:

numeric/cg/terms/AppleSauce.hpp

Reference: Existing Contributions 

File / Struct	Concept	Description
`Coriolis.hpp` / `GwceCoriolis`	`GwcePhysicsTerm`	Coriolis force (jx, jy) with spherical correction
`BottomFrictionTerm.hpp` / `GwceBottomFrictionT`	`GwcePhysicsTerm`	TKM bottom friction tensor (jx, jy)
`Tau0.hpp` / `GwceTau0T`	`GwcePhysicsTerm`	Spatially variable tau0 weighting (jx, jy)
`TidePotential.hpp` / `GwceTidePotentialT`	`GwcePhysicsTerm`	Tide potential gradient (jx, jy)
`WindStressTerm.hpp` / `GwceWindStressT`	`GwcePhysicsTerm`	Kinematic wind stress (jx, jy) via `wind_stress<DragLaw>()`
`AtmPressureGradient.hpp` / `GwceAtmPressureGradientT`	`GwcePhysicsTerm`	Atmospheric pressure gradient (jx, jy)
`BaroclinicGradientTerm.hpp` / `GwceBaroclinicGradientT`	`GwcePhysicsTerm`	Baroclinic pressure gradient (jx, jy)
`Tau0SpatialGradient.hpp` / `GwceTau0SpatialGradient`, `GwceQuadTau0SpatialGradient`	none (`compute` returns `double`)	Spatially varying tau0 gradient, a scalar (triangle and quadrilateral forms)
`LateralStressTerm.hpp` / `GwceLateralStressT`	`GwcePhysicsTerm`	Lateral stress divergence (jx, jy)
`GwceAdvection.hpp` / `GwceAdvectionT`	`GwcePhysicsTerm`	Non-conservative advection (jx, jy)
`PressureGradientTerm.hpp` / `MomentumPressureGradient`	`MomentumElementTerm`	Pressure gradient (x, y) (uniform per element)
`GwceAdvection.hpp` / `MomentumAdvection`	`MomentumElementTerm`	Non-conservative advection (x, y) (uniform)
`AtmPressureGradient.hpp` / `MomentumAtmPressureGradient`	`MomentumElementTerm`	Atmospheric pressure gradient (x, y) (uniform)
`BaroclinicGradientTerm.hpp` / `MomentumBaroclinicGradient`	`MomentumNodalTerm`	Baroclinic pressure gradient (x, y) (per-node)
`LateralStressTerm.hpp` / `MomentumLateralStress`	`MomentumNodalTerm`	Inline lateral stress (x, y) (per-node)

Adding a Meteorological File Format 

The cocoa_meteo library provides a pluggable reader framework for meteorological forcing data. To add support for a new file format, implement a reader class that derives from MeteoReaderBase and register it in the factory. Multi-source composition (mixing formats, or several instances of the same format) is handled by CompositeMeteoReader, the wrapper MeteoReaderFactory builds over two or more configured sources. A new format reader therefore composes with every existing one without any of the steps below touching the composite.

Architecture Overview 

Fig. 57 cocoa_meteo reader architecture

Step 1: Add the Format Enum 

Add an entry to MeteoFormat in cocoa_meteo/core/MeteoFormat.hpp:

enum class MeteoFormat { CF_NETCDF, OWI_ASCII, OWI_NETCDF, GRIB, YOUR_FORMAT };

Step 2: Create the Reader Class 

Create a new directory under cocoa_meteo/readers/your_format/ and implement a reader that inherits from MeteoReaderBase:

#pragma once
#include "cocoa_meteo/core/MeteoReaderBase.hpp"
#include "cocoa_meteo/core/MeteoReaderConfig.hpp"

namespace Cocoa::Meteo {

class YourFormatReader : public MeteoReaderBase {
 public:
  // Construction opens the file(s), indexes snapshot positions, and fails
  // loud on any format error. There is no separate open() step, so a
  // constructed reader is always usable.
  explicit YourFormatReader(const MeteoReaderConfig& config);

  [[nodiscard]] auto num_domains() const -> size_t override;
  [[nodiscard]] auto num_snapshots() const -> size_t override;
  [[nodiscard]] auto grid_info(size_t domain_idx) const
      -> MeteoGrid override;
  [[nodiscard]] auto time_values() const
      -> const std::vector<double>& override;
  [[nodiscard]] auto read_snapshot(size_t snapshot_idx,
                                   size_t domain_idx)
      -> MeteoFieldSnapshot override;

 private:
  // Format-specific internal state
};

}  // namespace Cocoa::Meteo

The pure virtual methods define the contract:

Method	Description
`num_domains()`	Number of grid domains (1 for single-domain formats)
`num_snapshots()`	Total number of time snapshots in the file
`grid_info(domain_idx)`	Return a `MeteoGrid` (`RegularGrid` or `IrregularGrid`)
`time_values()`	Return snapshot times as seconds since Unix epoch
`read_snapshot(idx, domain)`	Read and return a `MeteoFieldSnapshot` for the given time/domain

The MeteoFieldSnapshot bundles the grid definition with the three forcing fields (pressure, wind_u, wind_v), and the sea-ice concentration where the source carries it, for one time step and one domain. read_snapshot indexes the domain’s own time axis (domain_time_values(domain)) and must throw on an out-of-range index.

MeteoReaderBase also declares six optional virtuals with defaults that are correct for most formats; override only what your format needs:

Method	Default	Override when…
`is_moving_domain(domain_idx)`	`false`	the domain’s grid coordinates change between snapshots (e.g. a storm-following nest); moving domains are re-sent with every snapshot broadcast and have their interpolation weights recomputed.
`domain_name(domain_idx)`	`""` (empty)	the format has a natural identity for the domain (a filename, file pair, or group name) to show in start/stop log messages.
`domain_time_values(domain_idx)`	the shared `time_values()`	domains within one reader instance do not share a single time axis (e.g. a nested grid active only during part of the run).
`domain_source(domain_idx)`	`0`	never override this in a leaf reader. It exists so `CompositeMeteoReader` can report which configured source (child reader) owns a flattened domain, routing the per-source `wind_reduction` flag to the right nodes; a single-format reader always has exactly one source, so the default is always correct.
`domain_has_ice(domain_idx)`	`false`	the format can carry sea-ice concentration; the domain must then carry it at every one of its snapshots.
`domain_has_met_fields(domain_idx)`	`true`	the format can legally omit pressure and wind (an ice-only source, as `CfNetcdfReader` allows).

Step 3: Register in the Factory 

Add a case to create_reader() in cocoa_meteo/MeteoReaderFactory.cpp:

case MeteoFormat::YOUR_FORMAT:
  return std::make_unique<YourFormatReader>(config);

Step 4: Add Configuration Parsing 

In src/cocoa_kernel/core/config/ForcingParsers.cpp, add the format’s string spellings to the METEO_FORMATS lookup table that parse_meteo_format() reads (bump the table’s size template argument):

constexpr Core::EnumTable<Meteo::MeteoFormat, 9> METEO_FORMATS{
    {{{"cf_netcdf", Meteo::MeteoFormat::CF_NETCDF},
      {"cf", Meteo::MeteoFormat::CF_NETCDF},
      {"netcdf", Meteo::MeteoFormat::CF_NETCDF},
      {"owi_ascii", Meteo::MeteoFormat::OWI_ASCII},
      {"owi", Meteo::MeteoFormat::OWI_ASCII},
      {"owi_netcdf", Meteo::MeteoFormat::OWI_NETCDF},
      {"grib", Meteo::MeteoFormat::GRIB},
      {"grib2", Meteo::MeteoFormat::GRIB},
      {"your_format", Meteo::MeteoFormat::YOUR_FORMAT}}}};

This table is consulted once per configured source, so a run with several domains: entries can mix your_format with any other format for free. Add any format-specific fields to MeteoReaderConfig if needed (e.g., custom variable names or file organization parameters), and parse them alongside the format string in the same file.

Step 5: Update CMakeLists.txt 

Add the new source files to src/cocoa_meteo/CMakeLists.txt:

readers/your_format/YourFormatReader.hpp
readers/your_format/YourFormatReader.cpp

Key Data Structures 

MeteoGrid: grid geometry returned by grid_info():

// Regular rectangular grid
struct RegularGrid {
  double xll, yll;    // Lower-left corner (longitude, latitude)
  double dx, dy;      // Grid spacing [degrees]
  size_t nx, ny;      // Grid dimensions
  bool global;        // True if 360-degree longitude wrap
};

// Irregular grid (explicit coordinates per point)
struct IrregularGrid {
  size_t nx, ny;
  MeteoField longitude;  // [ny * nx] coordinate field
  MeteoField latitude;   // [ny * nx] coordinate field
};

using MeteoGrid = std::variant<RegularGrid, IrregularGrid>;

Use RegularGrid for structured data on a latitude-longitude grid. Use IrregularGrid when grid coordinates vary per point (e.g., moving vortex grids).

MeteoField: 2D field storage [ny * nx].

Data is stored in row-major order: data[j * nx + i] where i is the longitude index and j is the latitude index. This matches the NetCDF [lat, lon] dimension convention.

Scale factors:

The MeteoReaderConfig provides pressure_scale and wind_scale fields. Apply them in your reader’s read_snapshot() before returning, as every existing reader does, by calling apply_scale_factors() from core/MeteoFieldSnapshot.hpp. Nothing downstream scales the snapshot, so a reader that returns raw values yields unscaled winds and pressures.

Reference: Existing Readers 

Reader	Domains	Notes
`CfNetcdfReader`	Single	CF-compliant NetCDF with customizable variable names. Auto-detects global grids.
`OwiAsciiReader`	Multiple	Paired pressure/wind ASCII files with fixed-width headers. Indexes file positions on open for random-access reads.
`OwiNetcdfReader`	Multiple	NetCDF group-based format with support for moving grids (per-timestep coordinates).
`GribReader`	Multiple	GRIB2 via ecCodes, gated behind `cocoa_ENABLE_GRIB`. One domain per configured weather product (priority level); supports moving (storm-following) nests.
`CompositeMeteoReader`	Flattened	Not a format reader; wraps two or more configured sources of any format(s) into one flat domain space, list order = coverage priority.