====================
Extending the Model
====================

This guide covers how to extend Cocoa with new model capabilities, including
nodal attributes (spatially varying parameters) and physics contributions
(terms in the GWCE and momentum equations).

.. contents:: On This Page
   :local:
   :depth: 2

.. _adding-a-nodal-attribute:

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``:

.. code-block:: cpp

   enum class NodalAttributeId : int {
     ManningsN = 0,
     Tau0 = 1,
     AppleSauce = 2,  // New attribute
   };

   /// Number of defined nodal attributes (update when adding new entries)
   static constexpr size_t NUM_NODAL_ATTRIBUTES = 3;
   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:

.. code-block:: cpp

   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";
       case NodalAttributeId::AppleSauce:
         return "apple_sauce_coefficient";
       default:
         return "";
     }
   }

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:

.. code-block:: cpp

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

     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()``:

.. code-block:: cpp

   void NodalAttributeInitializer::initialize() {
     m_mesh->initialize_attributes(m_fields.num_nodes());
     initialize_manning_n();
     initialize_tau0();
     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:

.. code-block:: cpp

   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
   });

The ``get()`` method returns a ``const NodalAttributeAccessor&`` from a
pre-built cache array — pure array indexing with zero overhead. If the
attribute was not registered, ``get()`` throws with a clear error message.

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``. Use the mesh sentinel value
(``Sentinels::MESH_ATTRIBUTE``) to indicate the value should be read from
the mesh file rather than using a constant.

Step 5: Update Tests
^^^^^^^^^^^^^^^^^^^^

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

.. code-block:: cpp

   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:

**Reading:** Use ``NodalAttributeReader::read_vector_attribute()`` or
``read_vector_attribute_distributed()`` for 2D variables
``[num_nodes, num_components]``:

.. code-block:: cpp

   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");

**Access:** Use the two-argument accessor to get a specific component:

.. code-block:: cpp

   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
   });

**Auto-detecting dimensionality:** If an attribute may be either scalar
(1D in the mesh file) or vector (2D), use
``reader.get_variable_dimensions()`` to detect the shape and register
accordingly. See ``initialize_internal_tide_friction()`` for an example.

Optional Attributes
^^^^^^^^^^^^^^^^^^^

Some attributes are optional — they should only be applied when enabled
by the user. For these, skip registration when disabled and check
``has()`` before accessing:

.. code-block:: cpp

   // 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
^^^^^^^^^^^^

The physics contribution system is built around two patterns:

**GWCE contributions** use nested ``Input`` structs constructed inline at the
call site. Each contribution's ``Input`` holds only the data that contribution
needs. Data loaded for one contribution's ``Input`` dies after ``compute()``
returns, reducing peak register pressure on GPU.

**Momentum contributions** use a shared ``MomentumElementContext`` struct.

.. graphviz::
   :align: center
   :caption: Physics contribution architecture for GWCE and Momentum

   digraph physics_arch {
     rankdir=TB;
     node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10];
     edge [color="#555555"];
     bgcolor="transparent";

     kernel [label="Assembly Kernel Lambda\n(loads data, constructs Input,\ncalls contributions, assembles RHS)", fillcolor="#E0E0E0", color="#616161"];

     subgraph cluster_gwce {
       label="GWCE Path";
       style="filled,rounded";
       fillcolor="#E3F2FD";
       color="#1E88E5";
       fontname="Helvetica";

       ginput [label="Contribution::Input\n(per-contribution struct)", fillcolor="#BBDEFB", color="#1565C0"];
       gcompute [label="static compute(Input)\n-> Types::Vec2<double> {jx, jy}", fillcolor="#BBDEFB", color="#1565C0"];
     }

     subgraph cluster_mom {
       label="Momentum Path";
       style="filled,rounded";
       fillcolor="#FFF3E0";
       color="#FB8C00";
       fontname="Helvetica";

       mctx [label="MomentumElementContext\n(shared per-node data)", fillcolor="#FFE0B2", color="#E65100"];
       mcompute [label="static compute(ctx)\n-> Types::Vec2<double> {x, y}\nor compute_nodal(ctx)\n-> Array<Types::Vec2<double>, 3>", fillcolor="#FFE0B2", color="#E65100"];
     }

     kernel -> ginput;
     kernel -> mctx;
     ginput -> gcompute [label="data dies\nafter call", fontsize=8, color="#43A047", fontcolor="#43A047"];
     mctx -> mcompute;
   }

Concepts
~~~~~~~~

Four C++20 concepts in ``numeric/terms/TermsConcepts.hpp`` define the contracts:

.. list-table::
   :header-rows: 1
   :widths: 30 40 30

   * - 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
   * - ``GwceScalarTerm``
     - ``static compute(const T::Input&)``
       ``-> double``
     - Returns a scalar GWCE RHS contribution (uniform across all nodes)
   * - ``MomentumElementTerm``
     - ``static compute(const MomentumElementContext&)``
       ``-> Types::Vec2<double>``
     - Returns uniform (x, y) contribution to both x- and y-momentum for all 3 nodes
   * - ``MomentumNodalTerm``
     - ``static compute_nodal(const MomentumElementContext&)``
       ``-> Kokkos::Array<Types::Vec2<double>, 3>``
     - Returns per-node (x, y) contributions to both x- and y-momentum

Step-by-Step Guide
^^^^^^^^^^^^^^^^^^

.. graphviz::
   :align: center
   :caption: Workflow for adding a new physics contribution

   digraph workflow {
     rankdir=TB;
     node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10];
     edge [color="#555555"];
     bgcolor="transparent";

     s1 [label="Step 1: Data Storage\nAdd fields to PhysicsData", fillcolor="#FFEBEE", color="#E53935"];
     s2 [label="Step 2: GWCE Contribution\nCreate numeric/terms/ header\nwith Input struct + compute()", fillcolor="#E3F2FD", color="#1E88E5"];
     s3 [label="Step 3: Momentum Contribution\nCreate numeric/terms/ header\nwith compute() or compute_nodal()", fillcolor="#E3F2FD", color="#1E88E5"];
     s4 [label="Step 4: Wire GWCE Kernel\nCapture view + call contribution\nin GwceVectorAssemblyKernels", fillcolor="#FFF3E0", color="#FB8C00"];
     s5 [label="Step 5: Wire Momentum Kernel\nCapture view + load context\nin MomentumRhsKernels", fillcolor="#FFF3E0", color="#FB8C00"];
     s6 [label="Step 6: Populate Data\nNodal attribute, computed kernel,\nor external forcing", fillcolor="#E8F5E9", color="#43A047"];
     s7 [label="Step 7: Ghost Exchange\nAdd to CommunicationManager\n(if element-assembled)", fillcolor="#F3E5F5", color="#8E24AA"];
     s8 [label="Step 8: Configuration + Build\nAdd config, parse YAML,\nupdate CMakeLists.txt", fillcolor="#E0E0E0", color="#616161"];

     s1 -> s2 -> s3 -> s4 -> s5 -> s6 -> s7 -> s8;
     s2 -> s3 [style=invis, weight=10];
   }

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``:

.. code-block:: cpp

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

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

.. code-block:: cpp

   // 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/terms/AppleSauce.hpp``. Each GWCE
contribution nests an ``Input`` struct containing only the data it needs,
and provides a static ``compute(const Input&)`` method:

.. code-block:: cpp

   #pragma once
   #include <Kokkos_Core.hpp>
   #include "numeric/terms/TermsConcepts.hpp"

   namespace Cocoa::Numeric::GwceRhsTerms {

   struct GwceAppleSauce {
     struct Input {
       Kokkos::Array<Types::Scalar, 3> apple_sauce;
       // Add only the fields this contribution needs
     };

     KOKKOS_INLINE_FUNCTION
     static Types::Vec2<double> compute(const Input& input) {
       const auto as_avg = (input.apple_sauce[0] + input.apple_sauce[1] +
                            input.apple_sauce[2]) / 3.0_s;
       // Physics math goes here
       const auto jx = /* ... */;
       const auto jy = /* ... */;
       return {jx, jy};
     }
   };

   static_assert(GwcePhysicsTerm<GwceAppleSauce>);

   }  // namespace Cocoa::Numeric::GwceRhsTerms

The ``static_assert`` verifies the concept contract at compile time. If
your method signature doesn't match, you get a clear error instead of a
cryptic template failure at the call site.

The nested ``Input`` pattern ensures that data loaded for this contribution
dies when ``compute()`` returns, freeing registers for subsequent
contributions.

Step 3: Momentum Contribution
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Add the momentum term to ``numeric/terms/AppleSauce.hpp`` (alongside the GWCE
term). Choose the concept based on whether the result is uniform across all 3
nodes or varies per node.

**Uniform across nodes** -- follow ``MomentumAdvection.hpp``:

.. code-block:: cpp

   struct MomentumAppleSauce {
     KOKKOS_INLINE_FUNCTION
     static Types::Vec2<double> compute(const MomentumElementContext& ctx) {
       // Physics math: same value scattered to all 3 nodes
       return {/* x */, /* y */};
     }
   };
   static_assert(MomentumElementTerm<MomentumAppleSauce>);

**Different per node** -- follow ``MomentumLateralStress.hpp``:

.. code-block:: cpp

   struct MomentumAppleSauce {
     KOKKOS_INLINE_FUNCTION
     static Kokkos::Array<Types::Vec2<double>, 3> compute_nodal(
         const MomentumElementContext& ctx) {
       // Physics math: different value per node
       return {/* node 0 */, /* node 1 */, /* node 2 */};
     }
   };
   static_assert(MomentumNodalTerm<MomentumAppleSauce>);

Step 4: Wire Into GWCE Kernel
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In ``numeric/continuity/GwceVectorAssemblyKernels.hpp``, make two changes
to ``main_rhs_assembly_kernel``:

**a) Include the header and capture the view** at the top of the function,
alongside the other physics view captures:

.. code-block:: cpp

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

   // Inside the function, alongside other view captures:
   const auto apple_sauce = physics.apple_sauce().coefficient;

**b) Call the contribution** inside the lambda, constructing the ``Input``
inline with aggregate initialization and accumulating into the local result:

.. code-block:: cpp

   jxy += GwceAppleSauce::compute(
       {load_nodes(apple_sauce)});

Add the ``using`` declaration at the top of the lambda:

.. code-block:: cpp

   using GwceRhsTerms::GwceAppleSauce;

That's it--one line to capture the view, one line to call. The ``Input``
temporary (and any data loaded into it) dies after ``compute()`` returns.
The result is accumulated into a ``Types::Vec2<double> jxy`` variable.

Step 5: Wire Into Momentum Kernel
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In ``numeric/momentum/MomentumRhsKernels.cpp``, same three-part pattern
in ``momentum_rhs_element_kernel``:

**a) Include and capture:**

.. code-block:: cpp

   #include "numeric/terms/AppleSauce.hpp"

   // Capture the view
   const auto apple_sauce = physics.apple_sauce().coefficient;

**b) Load into context:**

.. code-block:: cpp

   ctx.apple_sauce = {apple_sauce(ctx.nm[0]), apple_sauce(ctx.nm[1]),
                      apple_sauce(ctx.nm[2])};

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

.. code-block:: cpp

   const auto apple = MomentumAppleSauce::compute(ctx);
   const auto base_contrib = pressure + advect + apple;

For ``MomentumNodalTerm``, add in the per-node loop alongside lateral stress:

.. code-block:: cpp

   const auto apple = MomentumAppleSauce::compute_nodal(ctx);
   for (int i = 0; i < nodes_per_triangle; ++i) {
     const auto contrib = base_contrib + lateral[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 :ref:`adding-a-nodal-attribute` above.

**Computed each timestep** (like lateral stress):

Write a computation kernel that populates ``AppleSauceFields``. Call it from
``GwceVectorAssembler::assemble`` during the preprocessing phase:

.. code-block:: cpp

   // In GwceVectorAssembler::assemble_rhs, preprocessing block:
   {
     const auto preprocess_region =
         Core::KokkosProfileRegion("GWCE::Preprocessing");
     nodal_preprocessing_kernel(model_fields);
     elemental_preprocessing_kernel(model_fields, config);
     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 Exchange (if needed)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Ghost exchange is only needed if the field is **assembled from element
contributions via atomics**. In that case, ghost nodes only see contributions
from local elements and miss contributions from elements on their owning
rank.

Fields that are directly loaded from nodal attributes or computed per-node
(without element scatter) do **not** need ghost exchange.

Add a method to ``CommunicationManager``:

.. code-block:: cpp

   // In CommunicationManager.hpp (public):
   void post_apple_sauce_exchange(Data::PhysicsData& physics);

   // In CommunicationManager.cpp:
   void CommunicationManager::post_apple_sauce_exchange(
       Data::PhysicsData& physics) {
     if (!m_is_parallel) return;
     std::array fields = {
       Types::StridedScalar1DView(physics.apple_sauce().intensity)
     };
     exchange_scalar_fields(fields);
   }

Call from ``TimeStepper`` after the AppleSauce computation kernel, before
the GWCE/momentum assembly that reads the data.

.. note::

   Use the batched ``exchange_scalar_fields`` with ``std::array`` to exchange
   multiple fields in a single MPI round-trip rather than calling
   ``exchange_scalar_field`` once per view.

Step 8: Configuration and Build
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

**Configuration** -- Add parameters to ``core/config/PhysicsConfig.hpp``:

.. code-block:: cpp

   struct AppleSauceConfig {
     bool active{false};
     Types::Scalar coefficient{1.0};
   };

Add ``AppleSauceConfig apple_sauce;`` to ``PhysicsConfig``. Parse in
``core/config/ModelConfigurationFactory.cpp``.

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

.. code-block:: text

   numeric/terms/AppleSauce.hpp

Reference: Existing Contributions
----------------------------------

.. list-table::
   :header-rows: 1
   :widths: 35 20 45

   * - File / Struct
     - Concept
     - Description
   * - ``Coriolis.hpp`` -- ``GwceCoriolis``
     - ``GwcePhysicsTerm``
     - Coriolis force (jx, jy) with spherical correction
   * - ``BottomFrictionTerm.hpp`` -- ``GwceBottomFriction``
     - ``GwcePhysicsTerm``
     - TKM bottom friction tensor (jx, jy)
   * - ``Tau0.hpp`` -- ``GwceTau0``
     - ``GwcePhysicsTerm``
     - Spatially variable tau0 weighting (jx, jy)
   * - ``TidePotential.hpp`` -- ``GwceTidePotential``
     - ``GwcePhysicsTerm``
     - Tide potential gradient (jx, jy)
   * - ``WindStressTerm.hpp`` -- ``GwceWindStress``
     - ``GwcePhysicsTerm``
     - Kinematic wind stress (jx, jy) via ``wind_stress<DragLaw>()``
   * - ``AtmPressureGradient.hpp`` -- ``GwceAtmPressureGradient``
     - ``GwcePhysicsTerm``
     - Atmospheric pressure gradient (jx, jy)
   * - ``LateralStressTerm.hpp`` -- ``GwceLateralStress``
     - ``GwcePhysicsTerm``
     - Lateral stress divergence (jx, jy)
   * - ``GwceAdvection.hpp`` -- ``GwceAdvection``
     - ``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)
   * - ``LateralStressTerm.hpp`` -- ``MomentumLateralStress``
     - ``MomentumNodalTerm``
     - Inline lateral stress (x, y) (per-node)

.. _adding-meteo-format:

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.

Architecture Overview
^^^^^^^^^^^^^^^^^^^^^

.. graphviz::
   :align: center
   :caption: cocoa_meteo reader architecture

   digraph meteo_arch {
     rankdir=TB;
     node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10];
     edge [color="#555555"];
     bgcolor="transparent";

     factory [label="MeteoReaderFactory\ncreate_reader(config)", fillcolor="#E0E0E0", color="#616161"];

     base [label="MeteoReaderBase\n(abstract interface)", fillcolor="#E3F2FD", color="#1E88E5"];

     cf [label="CfNetcdfReader", fillcolor="#BBDEFB", color="#1565C0"];
     owi_a [label="OwiAsciiReader", fillcolor="#BBDEFB", color="#1565C0"];
     owi_n [label="OwiNetcdfReader", fillcolor="#BBDEFB", color="#1565C0"];
     yours [label="YourFormatReader", fillcolor="#C8E6C9", color="#2E7D32"];

     factory -> base [label="dispatches on\nMeteoFormat"];
     base -> cf;
     base -> owi_a;
     base -> owi_n;
     base -> yours [style=dashed, label="new"];
   }

Step 1: Add the Format Enum
^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

.. code-block:: cpp

   enum class MeteoFormat { CF_NETCDF, OWI_ASCII, OWI_NETCDF, 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``:

.. code-block:: cpp

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

   namespace Cocoa::Meteo {

   class YourFormatReader : public MeteoReaderBase {
    public:
     explicit YourFormatReader(const MeteoReaderConfig& config);

     void open() override;
     void close() override;
     [[nodiscard]] auto is_open() const -> bool override;
     [[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
   };

   // Verify concept compliance at compile time
   static_assert(MeteoReader<YourFormatReader>);

   }  // namespace Cocoa::Meteo

The pure virtual methods define the contract:

.. list-table::
   :header-rows: 1
   :widths: 30 70

   * - Method
     - Description
   * - ``open()``
     - Open file(s) and index snapshot positions
   * - ``close()``
     - Release file handles
   * - ``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) for one time step and one domain.

Step 3: Register in the Factory
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

.. code-block:: cpp

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

Step 4: Add Configuration Parsing
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In ``cocoa_core/config/ModelConfigurationFactory.cpp``, add the string-to-enum
mapping for the new format:

.. code-block:: cpp

   if (format_str == "your_format") {
     reader_config.format = Meteo::MeteoFormat::YOUR_FORMAT;
   }

Add any format-specific fields to ``MeteoReaderConfig`` if needed (e.g.,
custom variable names or file organization parameters).

Step 5: Update CMakeLists.txt
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

.. code-block:: text

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

Key Data Structures
^^^^^^^^^^^^^^^^^^^

**MeteoGrid** --- Grid geometry returned by ``grid_info()``:

.. code-block:: cpp

   // 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 these in your reader's ``read_snapshot()`` or let the caller
apply them --- the convention is to return raw values and let the
interpolation pipeline handle scaling.

Reference: Existing Readers
^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. list-table::
   :header-rows: 1
   :widths: 30 15 55

   * - 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).