From 0ed7118706f2d36831d744027d1dbaee95cf08b6 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 11:13:24 -0600
Subject: [PATCH 01/29] modify sequence grid to accept 2d shape and spacing

---
 sequence/_grid.py | 100 +++++++++++++++++++++++++++++++++++++---------
 1 file changed, 81 insertions(+), 19 deletions(-)

diff --git a/sequence/_grid.py b/sequence/_grid.py
index a422259..087f1f3 100644
--- a/sequence/_grid.py
+++ b/sequence/_grid.py
@@ -15,26 +15,60 @@
 class SequenceModelGrid(RasterModelGrid):
     """Create a Landlab ModelGrid for use with Sequence."""
 
-    def __init__(self, n_cols: int, spacing: float = 100.0):
+    def __init__(
+        self, shape: int | tuple[int, int], spacing: float | tuple[float, float] = 100.0
+    ):
         """Create a *Landlab* :class:`~landlab.grid.base.ModelGrid` for use with *Sequence*.
 
         Parameters
         ----------
-        n_cols : int
-            The number of columns.
-        spacing : float, optional
-            The spacing between columns.
+        shape : int or tuple of int
+            The number of columns in the cross-shore direction or, if a ``tuple``,
+            ``(n_rows, n_cols)``  where rows are in the along-shore direction and
+            columns in the cross-shore direction.
+        spacing : float or tuple of float, optional
+            The spacing between columns or, if ``len(shape) == 2``,
+            ``(row_spacing, col_spacing)``.
 
         Examples
         --------
         >>> from sequence import SequenceModelGrid
-        >>> grid = SequenceModelGrid(500, spacing=10.0)
-        >>> grid.shape
-        (3, 500)
-        >>> grid.spacing
-        (1.0, 10.0)
+        >>> grid = SequenceModelGrid(5, spacing=10.0)
+        >>> grid.y_of_row
+        array([  0.,  1.,  2.])
+        >>> grid.x_of_column
+        array([ 0.,  10.,  20.,  30.,  40.])
+
+        >>> grid = SequenceModelGrid((3, 5), spacing=(10000.0, 10.0))
+        >>> grid.y_of_row
+        array([  0.,  10000.,  20000.,  30000., 40000.])
+        >>> grid.x_of_column
+        array([ 0.,  10.,  20.,  30.,  40.])
         """
-        super().__init__((3, n_cols), xy_spacing=(spacing, 1.0))
+        shape = np.atleast_1d(np.asarray(shape, dtype=int))
+        spacing = np.atleast_1d(np.asarray(spacing, dtype=float))
+
+        if (
+            (len(shape) == 1 and len(spacing) != 1) or (len(shape) != 1 and len(spacing) != len(shape))
+        ):
+            raise ValueError(
+                f"spacing dimension ({len(spacing)}) does not match shape"
+                f" dimension ({len(shape)})"
+            )
+
+        if len(shape) == 1:
+            n_rows, n_cols = 1, shape[0]
+        elif len(shape) == 2:
+            n_rows, n_cols = shape
+        else:
+            raise ValueError(f"invalid number of dimensions for grid ({len(shape)})")
+
+        if len(shape) == 1:
+            row_spacing, col_spacing = 1.0, spacing[0]
+        elif len(shape) == 2:
+            row_spacing, col_spacing = spacing
+
+        super().__init__((n_rows + 2, n_cols), xy_spacing=(col_spacing, row_spacing))
 
         self.status_at_node[self.nodes_at_top_edge] = self.BC_NODE_IS_CLOSED
         self.status_at_node[self.nodes_at_bottom_edge] = self.BC_NODE_IS_CLOSED
@@ -42,6 +76,16 @@ def __init__(self, n_cols: int, spacing: float = 100.0):
         self.at_node["sediment_deposit__thickness"] = self.zeros(at="node")
         self.at_grid["sea_level__elevation"] = 0.0
 
+        self.new_field_location("row", size=int(n_rows))
+
+    @property
+    def x_of_column(self) -> NDArray:
+        return self.x_of_node[self.nodes_at_top_edge]
+
+    @property
+    def y_of_row(self) -> NDArray:
+        return self.y_of_node[self.nodes_at_left_edge]
+
     def get_profile(self, name: str) -> NDArray:
         """Return the values of a field along the grid's profile.
 
@@ -55,7 +99,8 @@ def get_profile(self, name: str) -> NDArray:
         values : ndarray
             The values of the field located at the middle row of nodes.
         """
-        return self.at_node[name].reshape(self.shape)[1]
+        return self.at_node[name].reshape(self.shape)[1:-1]
+        # return self.at_node[name].reshape(self.shape)[row]
 
     @classmethod
     def from_toml(cls, filepath: os.PathLike[str]) -> "SequenceModelGrid":
@@ -85,19 +130,36 @@ def from_dict(cls, params: dict) -> "SequenceModelGrid":
         params : dict
             A dictionary that contains the parameters needed to
             create the grid.
+
+        Examples
+        --------
+        >>> from sequence import SequenceModelGrid
+        >>> params = {"shape": 5, "spacing": 10.0}
+        >>> grid = SequenceModelGrid.from_dict(params)
+        >>> grid.y_of_row
+        array([  0.,  1.,  2.])
+        >>> grid.x_of_column
+        array([ 0.,  10.,  20.,  30.,  40.])
+
+        >>> params = {"shape": (3, 5), "spacing": (10000.0, 10.0)}
+        >>> grid = SequenceModelGrid.from_dict(params)
+        >>> grid.y_of_row
+        array([  0.,  10000.,  20000.,  30000., 40000.])
+        >>> grid.x_of_column
+        array([ 0.,  10.,  20.,  30.,  40.])
         """
-        if "n_cols" in params:
-            n_cols = params["n_cols"]
-        elif "shape" in params:
-            n_cols = params["shape"][1]
+        if "shape" in params:
+            shape = params["shape"]
+        elif "n_cols" in params:
+            shape = params["n_cols"]
         else:
-            raise KeyError("n_cols")
+            raise KeyError("shape")
 
         if "spacing" in params:
             spacing = params["spacing"]
         elif "xy_spacing" in params:
-            spacing = np.broadcast_to(params["xy_spacing"], 2)[0]
+            spacing = np.atleast_1d(params["xy_spacing"])[::-1]
         else:
             raise KeyError("spacing")
 
-        return cls(n_cols, spacing=spacing)
+        return cls(shape, spacing=spacing)

From 5d02608ae804b699ee0f12d44cc111096f9c3459 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 11:42:07 -0600
Subject: [PATCH 02/29] modify to find shoreline at each grid row

---
 sequence/shoreline.py | 143 ++++++++++++++++++++++++++----------------
 1 file changed, 88 insertions(+), 55 deletions(-)

diff --git a/sequence/shoreline.py b/sequence/shoreline.py
index fa59340..bd8f1d2 100644
--- a/sequence/shoreline.py
+++ b/sequence/shoreline.py
@@ -44,7 +44,7 @@ class ShorelineFinder(Component):
             "intent": "out",
             "optional": False,
             "units": "m",
-            "mapping": "grid",
+            "mapping": "row",
             "doc": "Position of the shore line",
         },
         "x_of_shelf_edge": {
@@ -52,7 +52,7 @@ class ShorelineFinder(Component):
             "intent": "out",
             "optional": False,
             "units": "m",
-            "mapping": "grid",
+            "mapping": "row",
             "doc": "Position of the shelf edge",
         },
     }
@@ -71,24 +71,31 @@ def __init__(self, grid: SequenceModelGrid, alpha: float = 0.0005):
         super().__init__(grid)
 
         self._alpha = alpha
-        self.grid.at_grid["x_of_shore"] = 0.0
-        self.grid.at_grid["x_of_shelf_edge"] = 0.0
+        self.grid.at_row["x_of_shore"] = np.zeros(self.grid.shape[0] - 2)
+        self.grid.at_row["x_of_shelf_edge"] = np.zeros(self.grid.shape[0] - 2)
 
-        x = self.grid.x_of_node[self.grid.node_at_cell]
-        z = self.grid.at_node["topographic__elevation"][self.grid.node_at_cell]
-        dz = self.grid.at_node["sediment_deposit__thickness"][self.grid.node_at_cell]
+        x = self.grid.x_of_column
+        z = self.grid.get_profile("topographic__elevation")
+        dz = self.grid.get_profile("sediment_deposit__thickness")
         sea_level = self.grid.at_grid["sea_level__elevation"]
 
-        x_of_shore = find_shoreline(x, z, sea_level=sea_level)
-        try:
-            x_of_shelf_edge = find_shelf_edge(
-                x, dz, x_of_shore=x_of_shore, alpha=self._alpha
-            )
-        except ShelfEdgeError:
-            x_of_shelf_edge = np.nan
+        x_of_shores = np.empty(z.shape[0], dtype=float)
+        x_of_shelf_edges = np.empty(z.shape[0], dtype=float)
+
+        for row in range(z.shape[0]):
+            x_of_shore = find_shoreline(x, z[row], sea_level=sea_level)
+            try:
+                x_of_shelf_edge = find_shelf_edge(
+                    x, dz[row], x_of_shore=x_of_shore, alpha=self._alpha
+                )
+            except ShelfEdgeError:
+                x_of_shelf_edge = np.nan
+
+            x_of_shores[row] = x_of_shore
+            x_of_shelf_edges[row] = x_of_shelf_edge
 
-        self.grid.at_grid["x_of_shore"] = x_of_shore
-        self.grid.at_grid["x_of_shelf_edge"] = x_of_shelf_edge
+        self.grid.at_row["x_of_shore"][:] = x_of_shores
+        self.grid.at_row["x_of_shelf_edge"][:] = x_of_shelf_edges
 
     @property
     def alpha(self) -> float:
@@ -97,23 +104,30 @@ def alpha(self) -> float:
 
     def update(self) -> None:
         """Update the component one time step to find the new shoreline."""
-        x = self.grid.x_of_node[self.grid.node_at_cell]
-        z = self.grid.at_node["topographic__elevation"][self.grid.node_at_cell]
-        dz = self.grid.at_node["sediment_deposit__thickness"][self.grid.node_at_cell]
+        x = self.grid.x_of_column
+        z = self.grid.get_profile("topographic__elevation")
+        dz = self.grid.get_profile("sediment_deposit__thickness")
         sea_level = self.grid.at_grid["sea_level__elevation"]
 
-        x_of_shore = find_shoreline(x, z, sea_level=sea_level)
-        try:
-            x_of_shelf_edge = find_shelf_edge(
-                x, dz, x_of_shore=x_of_shore, alpha=self._alpha
-            )
-        except ShelfEdgeError:
-            x_of_shelf_edge = np.nan
-        else:
-            assert x_of_shelf_edge > x_of_shore
+        x_of_shores = np.empty(z.shape[0], dtype=float)
+        x_of_shelf_edges = np.empty(z.shape[0], dtype=float)
 
-        self.grid.at_grid["x_of_shore"] = x_of_shore
-        self.grid.at_grid["x_of_shelf_edge"] = x_of_shelf_edge
+        for row in range(z.shape[0]):
+            x_of_shore = find_shoreline(x, z[row], sea_level=sea_level)
+            try:
+                x_of_shelf_edge = find_shelf_edge(
+                    x, dz[row], x_of_shore=x_of_shore, alpha=self._alpha
+                )
+            except ShelfEdgeError:
+                x_of_shelf_edge = np.nan
+            else:
+                assert x_of_shelf_edge > x_of_shore
+
+            x_of_shores[row] = x_of_shore
+            x_of_shelf_edges[row] = x_of_shelf_edge
+
+        self.grid.at_row["x_of_shore"][:] = x_of_shores
+        self.grid.at_row["x_of_shelf_edge"][:] = x_of_shelf_edges
 
     def run_one_step(self, dt: Optional[float] = None) -> None:
         """Update the component on time step.
@@ -163,11 +177,10 @@ def find_shelf_edge_by_curvature(
     return x[np.argmin(curvature)]
 
 
-# def find_shelf_edge(grid, x, wd, x_of_shore, sea_level=0.0, alpha = 0.0005):
 def find_shelf_edge(
-    x: list[float],
-    dz: list[float],
-    x_of_shore: float = 0.0,
+    x: NDArray[np.floating],
+    dz: NDArray[np.floating],
+    x_of_shore: NDArray[np.floating] | float = 0.0,
     alpha: float = 0.0005,
 ) -> float:
     """Find the shelf edge based on deposit thickness.
@@ -218,7 +231,7 @@ def find_shoreline(
     z: NDArray[np.floating],
     sea_level: float = 0.0,
     kind: str = "cubic",
-) -> float:
+) -> NDArray[np.floating] | float:
     """Find the shoreline of a profile.
 
     Parameters
@@ -269,18 +282,27 @@ def find_shoreline(
 
     x, z = np.asarray(x), np.asarray(z)
 
-    try:
-        index_at_shore = find_shoreline_index(x, z, sea_level=sea_level)
-    except ShorelineError:
-        if z[0] < sea_level:
-            x_of_shoreline = x[0]
+    z = np.atleast_2d(z)
+    n_rows = z.shape[0]
+
+    x_of_shorelines = np.empty(n_rows, dtype=float)
+    for row in range(n_rows):
+        try:
+            index_at_shore = find_shoreline_index(x, z[row], sea_level=sea_level)
+        except ShorelineError:
+            if z[0] < sea_level:
+                x_of_shoreline = x[0]
+            else:
+                x_of_shoreline = x[-1]
         else:
-            x_of_shoreline = x[-1]
-    else:
-        func = interpolate.interp1d(x, z - sea_level, kind=kind)
-        x_of_shoreline = bisect(func, x[index_at_shore - 1], x[index_at_shore])
+            func = interpolate.interp1d(x, z[row] - sea_level, kind=kind)
+            x_of_shoreline = bisect(func, x[index_at_shore - 1], x[index_at_shore])
+        x_of_shorelines[row] = x_of_shoreline
 
-    return x_of_shoreline
+    if n_rows == 1:
+        return x_of_shorelines.item()
+    else:
+        return x_of_shorelines
 
 
 def _find_shoreline_polyfit(
@@ -330,7 +352,7 @@ def _find_shoreline_polyfit(
 
 def find_shoreline_index(
     x: NDArray[np.floating], z: NDArray[np.floating], sea_level: float = 0.0
-) -> int:
+) -> NDArray[np.integer] | int:
     """Find the landward-index of the shoreline.
 
     Parameters
@@ -381,14 +403,25 @@ def find_shoreline_index(
     Traceback (most recent call last):
     sequence.errors.ShorelineError: No shoreline found. The profile is all above sea level
     """
-    (below_water,) = np.where(z < sea_level)
-
-    if len(below_water) == 0 or len(below_water) == len(x):
-        raise ShorelineError(
-            "No shoreline found. The profile is all {} sea level".format(
-                "above" if len(below_water) == 0 else "below"
+    z = np.atleast_2d(z)
+    n_rows = z.shape[0]
+
+    index_of_shoreline = np.full(n_rows, -1)
+    for row in range(n_rows):
+        (below_water,) = np.where(z[row] < sea_level)
+
+        if len(below_water) == 0 or len(below_water) == len(x):
+            raise ShorelineError(
+                "No shoreline found. The profile is all {} sea level".format(
+                    "above" if len(below_water) == 0 else "below"
+                )
             )
-        )
-        return None
+            # return None
+        else:
+            index_of_shoreline[row] = below_water[0]
+            # return below_water[0]
+
+    if n_rows == 1:
+        return index_of_shoreline.item()
     else:
-        return below_water[0]
+        return index_of_shoreline

From c365aacac2cd042ba8ff44767123b5be75824713 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 11:44:05 -0600
Subject: [PATCH 03/29] modify flexure to calculate for each row

---
 sequence/sediment_flexure.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/sequence/sediment_flexure.py b/sequence/sediment_flexure.py
index 3dacc29..853fc58 100644
--- a/sequence/sediment_flexure.py
+++ b/sequence/sediment_flexure.py
@@ -30,9 +30,11 @@ def __init__(self, grid, isostasytime: Optional[float] = 7000.0, **kwds: dict):
         else:
             self._isostasy_time = self.validate_isostasy_time(isostasytime)
 
-        super().__init__(grid, rows=1, method="flexure", **kwds)
+        active_rows = np.arange(1, grid.shape[0] - 1)
 
-        self._subsidence_pool = np.zeros(grid.shape[1], dtype=float)
+        super().__init__(grid, rows=active_rows, method="flexure", **kwds)
+
+        self._subsidence_pool = np.zeros((len(active_rows), grid.shape[1]), dtype=float)
 
     @staticmethod
     def validate_isostasy_time(time: float) -> float:
@@ -103,6 +105,8 @@ def calc_dynamic_deflection(
         """
         isostasy_fraction = self.calc_isostasy_fraction(self.isostasy_time, dt)
 
+        isostasy_fraction /= self.grid.shape[0] - 2
+
         self._subsidence_pool[:] += isostatic_deflection
         deflection = self._subsidence_pool[:] * isostasy_fraction
         self._subsidence_pool[:] -= deflection

From d0ca8d2aebe652b254292ca43d57ca8e8ceb47d1 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 11:45:47 -0600
Subject: [PATCH 04/29] add some errors for invalid output files

---
 sequence/errors.py | 24 +++++++++++++++++++++---
 1 file changed, 21 insertions(+), 3 deletions(-)

diff --git a/sequence/errors.py b/sequence/errors.py
index 92ee018..c6b6a45 100644
--- a/sequence/errors.py
+++ b/sequence/errors.py
@@ -29,15 +29,33 @@ def __str__(self) -> str:
         return self._msg
 
 
-class MissingRequiredVariable(SequenceError):
+class OutputValidationError(SequenceError):
+    """Raise this exception if there is something wrong with an output file."""
+
+    pass
+
+
+class InvalidRowError(OutputValidationError):
+    """Raise this exception if a required variable is missing from an output file."""
+
+    def __init__(self, row: int, n_rows: int):
+        self._row = row
+        self._n_rows = n_rows
+
+    def __str__(self) -> str:
+        """Return error message that includes the requested row."""
+        return f"row {self._row!r} is out of bounds for output file with {self._n_rows} rows"
+
+
+class MissingRequiredVariable(OutputValidationError):
     """Raise this exception if a required variable is missing from an output file."""
 
     def __init__(self, name: str):
         self._name = name
 
-    def __str_(self) -> str:
+    def __str__(self) -> str:
         """Return error message that includes the name of the missing variable."""
-        return f"{self._name!r}"
+        return f"requireed variable ({self._name!r}) is missing from output file"
 
 
 class ParameterMismatchError(SequenceError):

From 92af39b40567d5c98a6b34983965f6a0fb2200c9 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 11:48:27 -0600
Subject: [PATCH 05/29] add file name to sequence-plot cli

---
 sequence/cli.py | 26 +++++++++++++++-----------
 1 file changed, 15 insertions(+), 11 deletions(-)

diff --git a/sequence/cli.py b/sequence/cli.py
index bc91e20..87b7d78 100644
--- a/sequence/cli.py
+++ b/sequence/cli.py
@@ -17,7 +17,7 @@
 from numpy.typing import ArrayLike
 from tqdm import tqdm
 
-from .errors import MissingRequiredVariable
+from .errors import OutputValidationError
 from .input_reader import TimeVaryingConfig
 from .logging import LoggingHandler
 from .plot import plot_file, plot_layers
@@ -381,12 +381,20 @@ def setup(set: str) -> None:
 
 
 @sequence.command()
-@click.option("--set", multiple=True, help="Set model parameters")
+@click.option("--set", "-S", multiple=True, help="Set model parameters")
+@click.argument(
+    "netcdf_file",
+    type=click.Path(
+        exists=True, file_okay=True, dir_okay=False, path_type=pathlib.Path
+    ),
+    nargs=1,
+)
 @click.pass_context
-def plot(ctx: Any, set: str) -> None:
+def plot(ctx: Any, set: str, netcdf_file) -> None:
     """Plot a Sequence output file."""
     verbose = ctx.parent.params["verbose"]
     folder = pathlib.Path.cwd()
+    path_to_file = folder / netcdf_file
 
     config = PLOT_KEYWORDS.copy()
 
@@ -398,7 +406,6 @@ def plot(ctx: Any, set: str) -> None:
         )
 
     config.update(**_load_params_from_strings(set))
-
     if verbose and len(config) > 0:
         logger.info(
             os.linesep.join(
@@ -409,14 +416,11 @@ def plot(ctx: Any, set: str) -> None:
             )
         )
 
-    logger.info(f"Plotting {folder / 'sequence.nc'}")
+    logger.info(f"Plotting {path_to_file!s}")
     try:
-        plot_file(folder / "sequence.nc", **config)
-    except MissingRequiredVariable as error:
-        logger.error(
-            f"{folder / 'sequence.nc'}: output file is missing a required variable "
-            f"({error})"
-        )
+        plot_file(path_to_file, **config)
+    except OutputValidationError as error:
+        logger.error(f"{path_to_file!s}: output file is invalid ({error!s})")
         raise click.Abort() from error
 
 

From 17347d0936814ae15cd51b81f8afd84caa76b54a Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 12:13:19 -0600
Subject: [PATCH 06/29] modify to write multiple rows to output file

---
 sequence/output_writer.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/sequence/output_writer.py b/sequence/output_writer.py
index 108dd98..56fa17f 100644
--- a/sequence/output_writer.py
+++ b/sequence/output_writer.py
@@ -4,6 +4,7 @@
 from os import PathLike
 from typing import Iterable, Optional, Union
 
+import numpy as np
 from landlab import Component
 
 from ._grid import SequenceModelGrid
@@ -51,9 +52,9 @@ def __init__(
         self.filepath = filepath
 
         if rows is not None:
-            self._nodes = grid.nodes[(rows,)].flatten()
+            self._rows = np.asarray(rows) - 1
         else:
-            self._nodes = None
+            self._rows = np.arange(grid.shape[0] - 2) + 1
 
         self._time = 0.0
         self._step_count = 0
@@ -74,7 +75,10 @@ def run_one_step(self, dt: Optional[float] = None) -> None:
                 mode="a",
                 time=self._time,
                 names={"node": self.fields},
-                ids={"node": self._nodes},
+                ids={
+                    "row": self._rows,
+                    "column": np.arange(self.grid.shape[1] - 2),
+                },
             )
         self._time += dt
         self._step_count += 1

From ac893b74c8208feaa6eb4caa366d97ad210ace75 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 12:28:13 -0600
Subject: [PATCH 07/29] fix some typing; add some docstrings

---
 sequence/_grid.py | 34 +++++++++++++++++++---------------
 1 file changed, 19 insertions(+), 15 deletions(-)

diff --git a/sequence/_grid.py b/sequence/_grid.py
index 087f1f3..e6a941f 100644
--- a/sequence/_grid.py
+++ b/sequence/_grid.py
@@ -45,28 +45,30 @@ def __init__(
         >>> grid.x_of_column
         array([ 0.,  10.,  20.,  30.,  40.])
         """
-        shape = np.atleast_1d(np.asarray(shape, dtype=int))
-        spacing = np.atleast_1d(np.asarray(spacing, dtype=float))
+        row_col_shape: list[int] = np.atleast_1d(np.asarray(shape, dtype=int))
+        row_col_spacing: list[float] = np.atleast_1d(np.asarray(spacing, dtype=float))
 
-        if (
-            (len(shape) == 1 and len(spacing) != 1) or (len(shape) != 1 and len(spacing) != len(shape))
+        if (len(row_col_shape) == 1 and len(row_col_spacing) != 1) or (
+            len(row_col_shape) != 1 and len(row_col_spacing) != len(row_col_shape)
         ):
             raise ValueError(
-                f"spacing dimension ({len(spacing)}) does not match shape"
-                f" dimension ({len(shape)})"
+                f"spacing dimension ({len(row_col_spacing)}) does not match shape"
+                f" dimension ({len(row_col_shape)})"
             )
 
-        if len(shape) == 1:
-            n_rows, n_cols = 1, shape[0]
-        elif len(shape) == 2:
-            n_rows, n_cols = shape
+        if len(row_col_shape) == 1:
+            n_rows, n_cols = 1, row_col_shape[0]
+        elif len(row_col_shape) == 2:
+            n_rows, n_cols = row_col_shape
         else:
-            raise ValueError(f"invalid number of dimensions for grid ({len(shape)})")
+            raise ValueError(
+                f"invalid number of dimensions for grid ({len(row_col_shape)})"
+            )
 
-        if len(shape) == 1:
-            row_spacing, col_spacing = 1.0, spacing[0]
-        elif len(shape) == 2:
-            row_spacing, col_spacing = spacing
+        if len(row_col_shape) == 1:
+            row_spacing, col_spacing = 1.0, row_col_spacing[0]
+        elif len(row_col_shape) == 2:
+            row_spacing, col_spacing = row_col_spacing
 
         super().__init__((n_rows + 2, n_cols), xy_spacing=(col_spacing, row_spacing))
 
@@ -80,10 +82,12 @@ def __init__(
 
     @property
     def x_of_column(self) -> NDArray:
+        """X-coordinate for each column of the grid."""
         return self.x_of_node[self.nodes_at_top_edge]
 
     @property
     def y_of_row(self) -> NDArray:
+        """Y-coordinate for each row of the grid."""
         return self.y_of_node[self.nodes_at_left_edge]
 
     def get_profile(self, name: str) -> NDArray:

From 8b8ed57172c39dce1f65e33c53c65d5aca1af2e4 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 12:57:59 -0600
Subject: [PATCH 08/29] modify fluvial to only operate on the middle row

---
 sequence/fluvial.py | 38 +++++++++++++++++++++++++++-----------
 1 file changed, 27 insertions(+), 11 deletions(-)

diff --git a/sequence/fluvial.py b/sequence/fluvial.py
index 6469a00..4dd7aed 100644
--- a/sequence/fluvial.py
+++ b/sequence/fluvial.py
@@ -122,6 +122,9 @@ def run_one_step(self, dt: float) -> None:
         dt : float
             The time step to update the component by.
         """
+
+        # return
+
         # Upstream boundary conditions  */
         mud_vol = self.sediment_load * (1.0 - self.sand_frac) / self.sand_frac
         sand_vol = self.sediment_load
@@ -142,18 +145,22 @@ def run_one_step(self, dt: float) -> None:
         # channel_width = sand_vol * self.basin_width / qs / 31536000.
         channel_width = sand_vol / qs
 
-        x = self.grid.x_of_node.reshape(self.grid.shape)[1]
+        # x = self.grid.x_of_node.reshape(self.grid.shape)[1]
+        x = self.grid.x_of_column
         elevation = self.grid.get_profile("topographic__elevation")
+        middle_row = elevation.shape[0] // 2
 
         shore = find_shoreline(
-            x, elevation, sea_level=self.grid.at_grid["sea_level__elevation"]
+            x,
+            elevation[middle_row],
+            sea_level=self.grid.at_grid["sea_level__elevation"],
         )
 
         land = x < shore
         water = x >= shore
         # land = self.grid.x_of_node[self.grid.node_at_cell] < shore
         # slope = np.gradient(z[1, land]) / self.grid.dx
-        slope = np.gradient(elevation) / self.grid.dx
+        slope = np.gradient(elevation[middle_row]) / self.grid.dx
         # slp = np.zeros(1)
         # slp[0] = self.plain_slope
         # slope = np.concatenate((slp,slop))
@@ -163,7 +170,7 @@ def run_one_step(self, dt: float) -> None:
         # )
 
         # channel_depth = np.zeros(self.grid.number_of_cells)
-        channel_depth = np.zeros_like(elevation)  # use upsteam slope
+        channel_depth = np.zeros_like(elevation[middle_row])  # use upsteam slope
         # channel_depth[0] = (
         #    (self.sand_density - 1000.) / 1000. * self.sand_grain / self.plain_slope
         channel_depth[land] = (
@@ -178,13 +185,19 @@ def run_one_step(self, dt: float) -> None:
         # width_cb = channel_width/epsilon
 
         # Original: r_cb = (model.new_height[i]-model.height[i]+model.d_sl);
-        r_cb = self.grid.get_profile("sediment_deposit__thickness") * epsilon
-        dz = self.grid.get_profile("bedrock_surface__increment_of_elevation").copy()
+        r_cb = (
+            self.grid.get_profile("sediment_deposit__thickness")[middle_row] * epsilon
+        )
+        dz = self.grid.get_profile("bedrock_surface__increment_of_elevation")[
+            middle_row
+        ].copy()
         # original: r_b = model.thickness[i];
-        r_b = self.grid.get_profile("sediment_deposit__thickness").copy()
+        r_b = self.grid.get_profile("sediment_deposit__thickness")[middle_row].copy()
         # r_fp = np.zeros(self.grid.shape[1])
-        r_fp = np.zeros_like(elevation)
-        percent_sand = self.grid.get_profile("delta_sediment_sand__volume_fraction")
+        r_fp = np.zeros_like(elevation[middle_row])
+        percent_sand = self.grid.get_profile("delta_sediment_sand__volume_fraction")[
+            middle_row
+        ]
         percent_sand.fill(0.0)
 
         for i in np.where(land)[0]:
@@ -265,7 +278,8 @@ def run_one_step(self, dt: float) -> None:
 
         water_depth = (
             self.grid.at_grid["sea_level__elevation"]
-            - self.grid.at_node["topographic__elevation"]
+            # - self.grid.at_node["topographic__elevation"]
+            - self.grid.get_profile("topographic__elevation")[middle_row]
         )
         add_mud = np.zeros(self.grid.shape[1])
         taper = 1
@@ -287,7 +301,9 @@ def run_one_step(self, dt: float) -> None:
                 if add_mud[i] < 0.0:
                     add_mud[i] = 0.0
 
-        sediment_thickness = self.grid.get_profile("sediment_deposit__thickness")
+        sediment_thickness = self.grid.get_profile("sediment_deposit__thickness")[
+            middle_row
+        ]
         sediment_thickness[water] += add_mud[water]
 
         np.divide(

From 83a5b3fe65a1edb5bbf463a2930af7c8164d0052 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 12:58:54 -0600
Subject: [PATCH 09/29] write layers on row and columns

---
 sequence/netcdf.py | 85 ++++++++++++++++++++++++++++++++++++----------
 1 file changed, 68 insertions(+), 17 deletions(-)

diff --git a/sequence/netcdf.py b/sequence/netcdf.py
index 8c8a9a9..e230275 100644
--- a/sequence/netcdf.py
+++ b/sequence/netcdf.py
@@ -1,4 +1,5 @@
 """Write a `SequenceModelGrid` to a NetCDF file."""
+import contextlib
 import os
 import warnings
 from collections import defaultdict
@@ -6,6 +7,7 @@
 from typing import Any, Iterable, Optional, Sequence, Union
 
 import netCDF4 as nc
+import numpy as np
 from numpy.typing import NDArray
 
 from ._grid import SequenceModelGrid
@@ -43,15 +45,28 @@ def _create_grid_dimension(
 ) -> Any:
     """Create grid dimensions for a netcdf file."""
     if ids is None:
-        # ids = Ellipsis
         ids = slice(None)
 
-    if at not in ("node", "link", "patch", "corner", "face", "cell", "grid"):
+    if at not in (
+        "node",
+        "link",
+        "patch",
+        "corner",
+        "face",
+        "cell",
+        "grid",
+        "row",
+        "column",
+    ):
         raise ValueError(f"unknown grid location {at}")
 
     if at not in root.dimensions:
         if at == "grid":
             number_of_elements = 1
+        elif at == "row":
+            number_of_elements = len(getattr(grid, f"y_of_{at}")[ids])
+        elif at == "column":
+            number_of_elements = len(getattr(grid, f"x_of_{at}")[ids])
         else:
             number_of_elements = len(getattr(grid, f"xy_of_{at}")[ids])
         root.createDimension(at, number_of_elements)
@@ -95,9 +110,16 @@ def _set_grid_coordinates(
     if at == "grid":
         return
 
-    coords = getattr(grid, f"xy_of_{at}")
-    root.variables[f"x_of_{at}"][:] = coords[ids, 0]
-    root.variables[f"y_of_{at}"][:] = coords[ids, 1]
+    if at == "row":
+        coords = getattr(grid, f"y_of_{at}")
+        root.variables[f"y_of_{at}"][:] = coords[ids]
+    elif at == "column":
+        coords = getattr(grid, f"x_of_{at}")
+        root.variables[f"x_of_{at}"][:] = coords[ids]
+    else:
+        coords = getattr(grid, f"xy_of_{at}")
+        root.variables[f"x_of_{at}"][:] = coords[ids, 0]
+        root.variables[f"y_of_{at}"][:] = coords[ids, 1]
 
 
 def _create_field(
@@ -139,6 +161,12 @@ def _set_field(
         )
         names = names & set(grid[at])
 
+    if at == "grid":
+        with contextlib.suppress(KeyError):
+            names.remove("x_of_shelf_edge")
+        with contextlib.suppress(KeyError):
+            names.remove("x_of_shore")
+
     _create_field(root, grid, at=at, names=names)
 
     if "time" in root.dimensions:
@@ -149,7 +177,6 @@ def _set_field(
             else:
                 values = grid[at][name]
             root.variables[_netcdf_var_name(name, at)][n_times - 1, :] = values
-
     else:
         for name in names:
             root.variables[_netcdf_var_name(name, at)][:] = grid[at][name][ids]
@@ -165,16 +192,21 @@ def _create_layers(
         netcdf_name = _netcdf_var_name(name, "layer")
         if netcdf_name not in root.variables:
             root.createVariable(
-                netcdf_name, _netcdf_type(grid.event_layers[name]), ("layer", "cell")
+                netcdf_name,
+                _netcdf_type(grid.event_layers[name]),
+                ("layer", "row", "column"),
             )
 
     netcdf_name = _netcdf_var_name("thickness", "layer")
     if netcdf_name not in root.variables:
-        root.createVariable(netcdf_name, "f8", ("layer", "cell"))
+        root.createVariable(netcdf_name, "f8", ("layer", "row", "column"))
 
 
 def _set_layers(
-    root: Any, grid: SequenceModelGrid, names: Optional[Iterable[str]] = None
+    root: Any,
+    grid: SequenceModelGrid,
+    names: Optional[Iterable[str]] = None,
+    ids: Optional[Union[slice, Iterable[int]]] = None,
 ) -> None:
     """Set values for variables at a grid layers."""
     if isinstance(names, str):
@@ -188,12 +220,15 @@ def _set_layers(
 
     n_layers = grid.event_layers.number_of_layers
     for name in names:
-        root.variables[_netcdf_var_name(name, "layer")][
-            :n_layers, :
-        ] = grid.event_layers[name][:]
-    root.variables[_netcdf_var_name("thickness", "layer")][
-        :n_layers, :
-    ] = grid.event_layers.dz[:]
+        layers = grid.event_layers[name][:, ids].reshape(
+            (-1, grid.shape[0] - 2, grid.shape[1] - 2)
+        )
+        root.variables[_netcdf_var_name(name, "layer")][:n_layers, :, :] = layers
+
+    dz = grid.event_layers.dz[:, ids].reshape(
+        (-1, grid.shape[0] - 2, grid.shape[1] - 2)
+    )
+    root.variables[_netcdf_var_name("thickness", "layer")][:n_layers, :, :] = dz
 
 
 def to_netcdf(
@@ -203,7 +238,7 @@ def to_netcdf(
     format: str = "NETCDF4",
     time: float = 0.0,
     at: Optional[Union[str, Sequence[str]]] = None,
-    ids: Optional[Union[dict[str, Iterable[int]], int, Iterable[int], slice]] = None,
+    ids: Optional[dict[str, NDArray[np.integer]] | int | Iterable[int] | slice] = None,
     names: Optional[
         Union[dict[str, Optional[Iterable[str]]], str, Iterable[str]]
     ] = None,
@@ -267,6 +302,12 @@ def to_netcdf(
             else:
                 names_dict[loc] = list(names_)
 
+    with contextlib.suppress(ValueError):
+        names_dict["grid"].remove("x_of_shore")
+    with contextlib.suppress(ValueError):
+        names_dict["grid"].remove("x_of_shelf_edge")
+    names_dict["row"] = ["x_of_shore", "x_of_shelf_edge"]
+
     ids_dict: dict[str, Union[slice, Iterable[int]]] = defaultdict(lambda: slice(None))
     if not isinstance(ids, dict):
         for loc in at:
@@ -286,6 +327,8 @@ def to_netcdf(
         root.createDimension("time", None)
         root.createVariable("time", "f8", ("time",))
 
+        _set_grid_coordinates(root, grid, at="row", ids=ids_dict["row"])
+        _set_grid_coordinates(root, grid, at="column", ids=ids_dict["column"])
         for loc in at:
             _set_grid_coordinates(root, grid, at=loc, ids=ids_dict[loc])
 
@@ -294,7 +337,15 @@ def to_netcdf(
 
     for loc in at:
         _set_field(root, grid, at=loc, ids=ids_dict[loc], names=names_dict[loc])
+    _set_field(
+        root,
+        grid,
+        at="row",
+        ids=ids_dict["row"],
+        names=["x_of_shore", "x_of_shelf_edge"],
+    )
+
     if with_layers:
-        _set_layers(root, grid, names=None)
+        _set_layers(root, grid, names=None, ids=ids_dict.get("cell", None))
 
     root.close()

From 8732f67faa5ecb2626ebd3daf47ea390560819fd Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 13:00:25 -0600
Subject: [PATCH 10/29] modify plotter to plot a given row

---
 sequence/plot.py | 49 +++++++++++++++++++++++++++++++++++-------------
 1 file changed, 36 insertions(+), 13 deletions(-)

diff --git a/sequence/plot.py b/sequence/plot.py
index dcd2f59..37c7b3f 100644
--- a/sequence/plot.py
+++ b/sequence/plot.py
@@ -11,7 +11,7 @@
 from scipy.interpolate import interp1d
 
 from ._grid import SequenceModelGrid
-from .errors import MissingRequiredVariable
+from .errors import InvalidRowError, MissingRequiredVariable
 
 
 def plot_layers(
@@ -185,9 +185,10 @@ def plot_grid(grid: SequenceModelGrid, **kwds: Any) -> None:
     :func:`plot_file` : Plot a `SequenceModelGrid`'s layers from a file.
     """
     elevation_at_layer = (
-        grid.at_node["bedrock_surface__elevation"][grid.node_at_cell] + grid.at_layer.z
+        grid.get_profile("bedrock_surface__elevation") + grid.at_layer.z
     )
-    x_of_stack = grid.x_of_node[grid.node_at_cell]
+
+    x_of_stack = grid.x_of_column
 
     x_of_shore = grid.at_layer_grid["x_of_shore"].flatten()
     x_of_shelf_edge = grid.at_layer_grid["x_of_shelf_edge"].flatten()
@@ -209,7 +210,9 @@ def plot_grid(grid: SequenceModelGrid, **kwds: Any) -> None:
     )
 
 
-def plot_file(filename: Union[str, PathLike], **kwds: Any) -> None:
+def plot_file(
+    filename: Union[str, PathLike], row: Optional[int] = None, **kwds: Any
+) -> None:
     """Plot a `SequenceModelGrid` from a *Sequence* output file.
 
     Parameters
@@ -224,25 +227,45 @@ def plot_file(filename: Union[str, PathLike], **kwds: Any) -> None:
     """
     kwds.setdefault("title", f"{filename}")
     with xr.open_dataset(filename) as ds:
+        if "row" not in ds.dims:
+            raise MissingRequiredVariable("row")
+
+        if row is None:
+            row = ds.dims["row"] // 2
+        elif (row >= ds.dims["row"]) or (row < -ds.dims["row"]):
+            raise InvalidRowError(row, ds.dims["row"])
+
         try:
-            thickness_at_layer = ds["at_layer:thickness"]
-            x_of_shore = ds["at_grid:x_of_shore"].data.squeeze()
-            x_of_shelf_edge = ds["at_grid:x_of_shelf_edge"].data.squeeze()
-            bedrock = ds["at_node:bedrock_surface__elevation"].data.squeeze()
+            thickness_at_layer = ds["at_layer:thickness"][:, row, :]
+            x_of_shore = ds["at_row:x_of_shore"].data.squeeze()
+            x_of_shelf_edge = ds["at_row:x_of_shelf_edge"].data.squeeze()
+            bedrock = (
+                ds["at_node:bedrock_surface__elevation"]
+                .data.reshape((-1, ds.dims["row"] + 2, ds.dims["column"] + 2))
+                .squeeze()
+            )
             time = ds["time"]
-            time_at_layer = ds["at_layer:age"]
+            time_at_layer = ds["at_layer:age"].data.reshape(
+                (-1, ds.dims["row"], ds.dims["column"])
+            )
         except KeyError as err:
             raise MissingRequiredVariable(str(err)) from err
 
         try:
-            x_of_stack = ds["x_of_cell"].data.squeeze()
+            x_of_stack = (
+                ds["x_of_cell"]
+                .data.reshape((ds.dims["row"], ds.dims["column"]))[row, :]
+                .squeeze()
+            )
         except KeyError:
             x_of_stack = np.arange(ds.dims["cell"])
 
-        elevation_at_layer = bedrock[-1, 1:-1] + np.cumsum(thickness_at_layer, axis=0)
+        elevation_at_layer = bedrock[-1, row, 1:-1] + np.cumsum(
+            thickness_at_layer, axis=0
+        )
 
-    x_of_shore = interp1d(time, x_of_shore)(time_at_layer[:, 0])
-    x_of_shelf_edge = interp1d(time, x_of_shelf_edge)(time_at_layer[:, 0])
+    x_of_shore = interp1d(time, x_of_shore[:, row])(time_at_layer[:, row, 0])
+    x_of_shelf_edge = interp1d(time, x_of_shelf_edge[:, row])(time_at_layer[:, row, 0])
 
     plot_layers(
         elevation_at_layer,

From c5cd1d9ef7030c4b30a23f6ea9a358c19349fc14 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 13:01:40 -0600
Subject: [PATCH 11/29] track shoreline and shelf edge at each row

---
 sequence/sequence_model.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sequence/sequence_model.py b/sequence/sequence_model.py
index 68d389e..2dd3079 100644
--- a/sequence/sequence_model.py
+++ b/sequence/sequence_model.py
@@ -116,8 +116,8 @@ def __init__(
         if output is not None:
             self._components["output"] = OutputWriter(self._grid, **output)
 
-        self.grid.at_grid["x_of_shore"] = np.nan
-        self.grid.at_grid["x_of_shelf_edge"] = np.nan
+        self.grid.at_row["x_of_shore"][:] = np.nan
+        self.grid.at_row["x_of_shelf_edge"][:] = np.nan
         self.grid.at_grid["sea_level__elevation"] = 0.0
         self._n_archived_layers = 0
 

From 2382fc23ce84891c6baa88c60161e0aa6d4a4e85 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 13:05:45 -0600
Subject: [PATCH 12/29] modify to operate with multiple rows

---
 sequence/submarine.py | 108 ++++++++++++++++++++++++++++--------------
 1 file changed, 72 insertions(+), 36 deletions(-)

diff --git a/sequence/submarine.py b/sequence/submarine.py
index 4284b38..eca04a5 100644
--- a/sequence/submarine.py
+++ b/sequence/submarine.py
@@ -168,10 +168,6 @@ def sediment_load(self, value: float) -> None:
         self._ksh = self._load / self._plain_slope
         self.grid.at_grid["sediment_load"] = self._load
 
-    # @property
-    # def k0(self):
-    #     return self._k0
-
     @property
     def load0(self) -> float:
         """Return the sediment load entering the profile."""
@@ -196,7 +192,9 @@ def sea_level(self) -> float:
     def sea_level(self, sea_level: float) -> None:
         self.grid.at_grid["sea_level__elevation"] = sea_level
 
-    def calc_diffusion_coef(self, x_of_shore: float) -> NDArray[np.floating]:
+    def calc_diffusion_coef(
+        self, x_of_shore: NDArray[np.floating] | float
+    ) -> NDArray[np.floating]:
         """Calculate and store diffusion coefficient values.
 
         Parameters
@@ -232,34 +230,70 @@ def calc_diffusion_coef(self, x_of_shore: float) -> NDArray[np.floating]:
         >>> diffusion_coef is grid.at_node["kd"]
         True
         """
-        sea_level = self.grid.at_grid["sea_level__elevation"]
-        water_depth = sea_level - self._grid.at_node["topographic__elevation"]
-
-        under_water = water_depth > 0.0
-        deep_water = water_depth > self._wave_base
-        land = ~under_water
+        x_of_shore = np.atleast_1d(x_of_shore)
 
-        k = self.grid.at_node["kd"]
-
-        x = self.grid.x_of_node
-        b = (self._shoreface_height * self._alpha + self._shelf_slope) * self.grid.dx
-
-        k[under_water] = (
-            self._load
-            * ((x[under_water] - x_of_shore) + self.grid.dx)
-            / (water_depth[under_water] + b)
-        )
+        sea_level = self.grid.at_grid["sea_level__elevation"]
 
-        k[deep_water] *= np.exp(
-            -(water_depth[deep_water] - self._wave_base) / self._wave_base
-        )
-
-        self._load = self._load0 * (1 + sea_level * self._load_sl)
-        self._ksh = self._load / self._plain_slope
-        if self._basin_width > 0.0:
-            k[land] = self._ksh * (self._basin_width + x[land]) / self._basin_width
-        else:
-            k[land] = self._ksh
+        water_depth = sea_level - self._grid.get_profile("topographic__elevation")
+        k = self.grid.get_profile("kd")
+        x = self.grid.x_of_column
+
+        n_rows = k.shape[0]
+
+        assert len(x_of_shore) == n_rows
+
+        for row in range(n_rows):
+            under_water = water_depth[row] > 0.0
+            deep_water = water_depth[row] > self._wave_base
+            land = ~under_water
+
+            b = (
+                self._shoreface_height * self._alpha + self._shelf_slope
+            ) * self.grid.dx
+
+            k[row, under_water] = (
+                self._load
+                * ((x[under_water] - x_of_shore[row]) + self.grid.dx)
+                / (water_depth[row, under_water] + b)
+            )
+
+            k[row, deep_water] *= np.exp(
+                -(water_depth[row, deep_water] - self._wave_base) / self._wave_base
+            )
+
+            self._load = self._load0 * (1 + sea_level * self._load_sl)
+            self._ksh = self._load / self._plain_slope
+
+            if self._basin_width > 0.0:
+                k[row, land] = (
+                    self._ksh * (self._basin_width + x[land]) / self._basin_width
+                )
+            else:
+                k[row, land] = self._ksh
+
+            # TODO: modify diffusion outside of the channel row.
+            if row != n_rows // 2:
+                k[row, land] *= 0.5
+
+            # if row == n_rows // 2:
+            #     if self._basin_width > 0.0:
+            #         k[row, land] = self._ksh * (self._basin_width + x[land]) / self._basin_width
+            #     else:
+            #         k[row, land] = self._ksh
+            # else:  # outside of the channel with low diffusivity
+            #     if self._basin_width > 0.0:
+            #         k[row, land] = self._ksh * (self.grid.dx + x[land]) / self._basin_width
+            #     else:
+            #         k[row, land] = self._ksh * (self.grid.dx) / self._basin_width
+
+            # if self._basin_width > 0.0:
+            #     k[land] = self._ksh * (self._basin_width + x[land]) / self._basin_width
+            # else:
+            #     k[land] = self._ksh
+
+        k = self.grid.at_node["kd"].reshape(self.grid.shape)
+        k[0, :] = k[1, :]
+        k[-1, :] = k[-1, :]
 
         return k
 
@@ -272,20 +306,22 @@ def run_one_step(self, dt: float) -> None:
             Time step to advance component by.
         """
         shore = find_shoreline(
-            self.grid.x_of_node[self.grid.node_at_cell],
-            self.grid.at_node["topographic__elevation"][self.grid.node_at_cell],
+            self.grid.x_of_column,
+            self.grid.get_profile("topographic__elevation"),
             sea_level=self.grid.at_grid["sea_level__elevation"],
         )
 
         self.calc_diffusion_coef(shore)
 
         # set elevation at upstream boundary to ensure proper sediment influx
-        x = self.grid.x_of_node.reshape(self.grid.shape)
+        x = self.grid.x_of_column
         z = self._grid.at_node["topographic__elevation"].reshape(self.grid.shape)
         # k = self._grid.at_node["kd"].reshape(self.grid.shape)
         # z[1, 0] = z[1,1] + self._load / k[1, 0] * (x[1,1]-x[1,0])
-        z[1, 0] = z[1, 1] + self._plain_slope * (x[1, 1] - x[1, 0])
-        # self._load/self._load0)
+        # z[1, 0] = z[1, 1] + self._plain_slope * (x[1, 1] - x[1, 0])
+        z[1, 0] = z[1, 1] + self._plain_slope * (x[1] - x[0])
+
+        z[1:-1, 0] = z[1:-1, 1] + self._plain_slope * (x[1] - x[0])
 
         z_before = self.grid.at_node["topographic__elevation"].copy()
 

From 44de4af03969807ab8a5a4f79084a2bc0d0d3ab4 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 13:09:47 -0600
Subject: [PATCH 13/29] remove some lint

---
 sequence/fluvial.py | 4 ----
 sequence/plot.py    | 2 ++
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/sequence/fluvial.py b/sequence/fluvial.py
index 4dd7aed..d7b4bb6 100644
--- a/sequence/fluvial.py
+++ b/sequence/fluvial.py
@@ -122,10 +122,6 @@ def run_one_step(self, dt: float) -> None:
         dt : float
             The time step to update the component by.
         """
-
-        # return
-
-        # Upstream boundary conditions  */
         mud_vol = self.sediment_load * (1.0 - self.sand_frac) / self.sand_frac
         sand_vol = self.sediment_load
         qs = (
diff --git a/sequence/plot.py b/sequence/plot.py
index 37c7b3f..0862736 100644
--- a/sequence/plot.py
+++ b/sequence/plot.py
@@ -219,6 +219,8 @@ def plot_file(
     ----------
     filename : path-like
         Path to the file to plot.
+    row : int, optional
+        Row to plot. If not provided, plot the middle row.
 
     See Also
     --------

From 59030d829aca78356a09d47fabc41591eb1309d7 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 20:57:08 -0600
Subject: [PATCH 14/29] add properties for number of rows/columns

---
 sequence/_grid.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/sequence/_grid.py b/sequence/_grid.py
index e6a941f..afd846f 100644
--- a/sequence/_grid.py
+++ b/sequence/_grid.py
@@ -85,6 +85,14 @@ def x_of_column(self) -> NDArray:
         """X-coordinate for each column of the grid."""
         return self.x_of_node[self.nodes_at_top_edge]
 
+    @property
+    def number_of_rows(self) -> int:
+        return self.shape[0] - 2
+
+    @property
+    def number_of_columns(self) -> int:
+        return self.shape[1]
+
     @property
     def y_of_row(self) -> NDArray:
         """Y-coordinate for each row of the grid."""

From 677e17405355c60760a90b8ca4f1020149fed22f Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 20:58:29 -0600
Subject: [PATCH 15/29] plot rows of a grid

---
 sequence/plot.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/sequence/plot.py b/sequence/plot.py
index 0862736..b8b3655 100644
--- a/sequence/plot.py
+++ b/sequence/plot.py
@@ -171,7 +171,7 @@ def plot_layers(
     plt.show()
 
 
-def plot_grid(grid: SequenceModelGrid, **kwds: Any) -> None:
+def plot_grid(grid: SequenceModelGrid, row: Optional[int] = None, **kwds: Any) -> None:
     """Plot a :class:`~SequenceModelGrid`.
 
     Parameters
@@ -184,20 +184,23 @@ def plot_grid(grid: SequenceModelGrid, **kwds: Any) -> None:
     :func:`plot_layers` : Plot layers from a 2D array of elevations.
     :func:`plot_file` : Plot a `SequenceModelGrid`'s layers from a file.
     """
+    if row is None:
+        row = grid.number_of_rows // 2
+
     elevation_at_layer = (
-        grid.get_profile("bedrock_surface__elevation") + grid.at_layer.z
+        grid.get_profile("bedrock_surface__elevation")[row, 1:-1] + grid.at_layer.z
     )
 
-    x_of_stack = grid.x_of_column
+    x_of_stack = grid.x_of_column[1:-1]
 
-    x_of_shore = grid.at_layer_grid["x_of_shore"].flatten()
-    x_of_shelf_edge = grid.at_layer_grid["x_of_shelf_edge"].flatten()
+    x_of_shore = grid.at_layer_row["x_of_shore"]
+    x_of_shelf_edge = grid.at_layer_row["x_of_shelf_edge"]
 
     time_at_layer_grid = grid.at_layer_grid["age"].flatten()
     time_at_layer = grid.at_layer["age"]
 
-    x_of_shore = interp1d(time_at_layer_grid, x_of_shore)(time_at_layer[:, 0])
-    x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge)(time_at_layer[:, 0])
+    x_of_shore = interp1d(time_at_layer_grid, x_of_shore[:, row])(time_at_layer[:, 0])
+    x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge[:, row])(time_at_layer[:, 0])
 
     kwds.setdefault("title", f"time = {time_at_layer[-1, 0]} years")
 

From c1f632ab54faaa0cd4208c986eaa8e2b5b4c246f Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 20:59:15 -0600
Subject: [PATCH 16/29] keep track of shoreline, shelf_edge at rows over time

---
 sequence/sequence.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/sequence/sequence.py b/sequence/sequence.py
index 38825b6..5a9bd4b 100644
--- a/sequence/sequence.py
+++ b/sequence/sequence.py
@@ -78,6 +78,7 @@ def __init__(
         self._n_archived_layers = 0
 
         self.grid.at_layer_grid = EventLayers(1)
+        self.grid.at_layer_row = EventLayers(grid.number_of_rows)
 
         if "bedrock_surface__elevation" not in self.grid.at_node:
             self.grid.add_field(
@@ -225,6 +226,14 @@ def add_layer(self, dz_at_cell: NDArray[np.floating]) -> None:
             x_of_shelf_edge = self.grid.at_grid["x_of_shelf_edge"]
         except KeyError:
             x_of_shelf_edge = np.nan
+        try:
+            x_of_shores = self.grid.at_row["x_of_shore"]
+        except KeyError:
+            x_of_shores = np.full(self.grid.number_of_rows, np.nan)
+        try:
+            x_of_shelf_edges = self.grid.at_row["x_of_shelf_edge"]
+        except KeyError:
+            x_of_shelf_edges = np.full(self.grid.number_of_rows, np.nan)
 
         self.grid.at_node["topographic__elevation"][
             self.grid.node_at_cell
@@ -238,6 +247,13 @@ def add_layer(self, dz_at_cell: NDArray[np.floating]) -> None:
             x_of_shore=x_of_shore,
             x_of_shelf_edge=x_of_shelf_edge,
         )
+        self.grid.at_layer_row.add(
+            1.0,
+            age=self.time,
+            sea_level=self.grid.at_grid["sea_level__elevation"],
+            x_of_shore=x_of_shores,
+            x_of_shelf_edge=x_of_shelf_edges,
+        )
 
         if (
             self.grid.event_layers.number_of_layers - self._n_archived_layers

From f3de8ddf2a727fc0aa57664134f9742e21b090b9 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Sat, 21 Oct 2023 21:00:19 -0600
Subject: [PATCH 17/29] wrap long line

---
 sequence/submarine.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/sequence/submarine.py b/sequence/submarine.py
index eca04a5..38c9f19 100644
--- a/sequence/submarine.py
+++ b/sequence/submarine.py
@@ -277,7 +277,9 @@ def calc_diffusion_coef(
 
             # if row == n_rows // 2:
             #     if self._basin_width > 0.0:
-            #         k[row, land] = self._ksh * (self._basin_width + x[land]) / self._basin_width
+            #         k[row, land] = (
+            #             self._ksh * (self._basin_width + x[land]) / self._basin_width
+            #         )
             #     else:
             #         k[row, land] = self._ksh
             # else:  # outside of the channel with low diffusivity

From 07693c5fe8df7e38fcb73e20c28f26615b9a1067 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 23 Oct 2023 13:42:51 -0600
Subject: [PATCH 18/29] add some docstrings

---
 sequence/_grid.py | 6 +++++-
 sequence/plot.py  | 6 +++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/sequence/_grid.py b/sequence/_grid.py
index afd846f..deb4008 100644
--- a/sequence/_grid.py
+++ b/sequence/_grid.py
@@ -78,6 +78,7 @@ def __init__(
         self.at_node["sediment_deposit__thickness"] = self.zeros(at="node")
         self.at_grid["sea_level__elevation"] = 0.0
 
+        # self.new_field_location("row", size=int(n_rows + 2))
         self.new_field_location("row", size=int(n_rows))
 
     @property
@@ -87,10 +88,13 @@ def x_of_column(self) -> NDArray:
 
     @property
     def number_of_rows(self) -> int:
-        return self.shape[0] - 2
+        """Number of rows of nodes."""
+        return self.shape[0]
+        # return self.shape[0] - 2
 
     @property
     def number_of_columns(self) -> int:
+        """Number of columns of nodes."""
         return self.shape[1]
 
     @property
diff --git a/sequence/plot.py b/sequence/plot.py
index b8b3655..46a9f3c 100644
--- a/sequence/plot.py
+++ b/sequence/plot.py
@@ -178,6 +178,8 @@ def plot_grid(grid: SequenceModelGrid, row: Optional[int] = None, **kwds: Any) -
     ----------
     grid : SequenceModelGrid
         The grid to plot.
+    row : int, optional
+        The row of the grid to plot. If not provided, plot the middle row.
 
     See Also
     --------
@@ -200,7 +202,9 @@ def plot_grid(grid: SequenceModelGrid, row: Optional[int] = None, **kwds: Any) -
     time_at_layer = grid.at_layer["age"]
 
     x_of_shore = interp1d(time_at_layer_grid, x_of_shore[:, row])(time_at_layer[:, 0])
-    x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge[:, row])(time_at_layer[:, 0])
+    x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge[:, row])(
+        time_at_layer[:, 0]
+    )
 
     kwds.setdefault("title", f"time = {time_at_layer[-1, 0]} years")
 

From 9991ecc839e0c705ef5a6cd85cf7283f084394a7 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 23 Oct 2023 13:44:45 -0600
Subject: [PATCH 19/29] modify for multiple rows

---
 sequence/shoreline.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sequence/shoreline.py b/sequence/shoreline.py
index bd8f1d2..6bf4725 100644
--- a/sequence/shoreline.py
+++ b/sequence/shoreline.py
@@ -290,7 +290,7 @@ def find_shoreline(
         try:
             index_at_shore = find_shoreline_index(x, z[row], sea_level=sea_level)
         except ShorelineError:
-            if z[0] < sea_level:
+            if z[row, 0] < sea_level:
                 x_of_shoreline = x[0]
             else:
                 x_of_shoreline = x[-1]

From ae430f4ed9339d38595db01bfa4f9caa024b873b Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 23 Oct 2023 13:47:52 -0600
Subject: [PATCH 20/29] update tests for multiple rows

---
 sequence/submarine.py                  |  6 +++---
 tests/test_grid.py                     |  6 ++----
 tests/test_netcdf.py                   |  2 +-
 tests/test_output_writer.py            | 27 +++++++++++++-------------
 tests/test_sequence.py                 | 12 ++++++------
 tests/test_sequence_model/marmara.yaml |  3 +--
 tests/test_subsidence.py               |  4 ++--
 7 files changed, 29 insertions(+), 31 deletions(-)

diff --git a/sequence/submarine.py b/sequence/submarine.py
index 38c9f19..fff093f 100644
--- a/sequence/submarine.py
+++ b/sequence/submarine.py
@@ -208,10 +208,10 @@ def calc_diffusion_coef(
         sea level and three below, two of which are in deep water (below the
         default 60 m wave base).
 
-        >>> from landlab import RasterModelGrid
+        >>> from sequence import SequenceModelGrid
         >>> import numpy as np
 
-        >>> grid = RasterModelGrid((3, 6), xy_spacing=200.0)
+        >>> grid = SequenceModelGrid((1, 6), spacing=(1.0, 200.0))
         >>> z = grid.add_zeros("topographic__elevation", at="node")
         >>> z[6:12] = np.array([3., 3., 1., -1., -85., -85.])
         >>> z.reshape((3, 6))
@@ -297,7 +297,7 @@ def calc_diffusion_coef(
         k[0, :] = k[1, :]
         k[-1, :] = k[-1, :]
 
-        return k
+        return self.grid.at_node["kd"]
 
     def run_one_step(self, dt: float) -> None:
         """Advance component one time step.
diff --git a/tests/test_grid.py b/tests/test_grid.py
index 44edb85..a84b309 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -22,13 +22,11 @@ def test_from_dict():
 def test_from_dict_old_style():
     expected = SequenceModelGrid(500, spacing=5.0)
 
-    actual = SequenceModelGrid.from_dict(
-        {"shape": (13, 500), "xy_spacing": (5.0, 25.0)}
-    )
+    actual = SequenceModelGrid.from_dict({"shape": (1, 500), "xy_spacing": (5.0, 1.0)})
     assert actual.shape == expected.shape
     assert actual.spacing == expected.spacing
 
-    actual = SequenceModelGrid.from_dict({"shape": (13, 500), "xy_spacing": 5.0})
+    actual = SequenceModelGrid.from_dict({"shape": 500, "xy_spacing": 5.0})
     assert actual.shape == expected.shape
     assert actual.spacing == expected.spacing
 
diff --git a/tests/test_netcdf.py b/tests/test_netcdf.py
index 7c75f24..12c4e27 100644
--- a/tests/test_netcdf.py
+++ b/tests/test_netcdf.py
@@ -94,7 +94,7 @@ def test_with_layers(tmpdir):
     grid = SequenceModelGrid(4)
     grid.event_layers.add(10.0, age=0.0, water_depth=np.arange(2))
     with tmpdir.as_cwd():
-        to_netcdf(grid, "test.nc")
+        to_netcdf(grid, "test.nc", ids={"row": [1], "column": [1, 2]})
         ds = xr.open_dataset("test.nc")
     assert ds.dims["layer"] == 1
     assert np.all(ds["at_layer:thickness"] == np.array([[10.0, 10.0]]))
diff --git a/tests/test_output_writer.py b/tests/test_output_writer.py
index dac021b..6869eed 100644
--- a/tests/test_output_writer.py
+++ b/tests/test_output_writer.py
@@ -1,13 +1,13 @@
 import numpy as np
 import xarray as xr
-from landlab import RasterModelGrid
 from pytest import approx, raises
 
+from sequence import SequenceModelGrid
 from sequence.output_writer import OutputWriter
 
 
 def test_no_fields(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         writer = OutputWriter(grid, "test.nc")
         writer.run_one_step()
@@ -19,7 +19,7 @@ def test_no_fields(tmpdir):
 
 
 def test_default_output_interval(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         writer = OutputWriter(grid, "test.nc")
         for _ in range(4):
@@ -32,7 +32,7 @@ def test_default_output_interval(tmpdir):
 
 
 def test_output_interval(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         writer = OutputWriter(grid, "test.nc", interval=2)
         for _ in range(3):
@@ -45,7 +45,7 @@ def test_output_interval(tmpdir):
 
 
 def test_default_clobber_is_false(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         OutputWriter(grid, "test.nc").run_one_step()
         with raises(RuntimeError):
@@ -53,7 +53,7 @@ def test_default_clobber_is_false(tmpdir):
 
 
 def test_clobber_keyword(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         writer = OutputWriter(grid, "test.nc")
         for _ in range(10):
@@ -70,7 +70,7 @@ def test_clobber_keyword(tmpdir):
 
 
 def test_default_is_all_rows(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         OutputWriter(grid, "test.nc").run_one_step()
         with xr.open_dataset("test.nc") as ds:
@@ -78,14 +78,15 @@ def test_default_is_all_rows(tmpdir):
 
 
 def test_rows_keyword(tmpdir):
-    grid = RasterModelGrid((3, 4))
+    grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0))
     with tmpdir.as_cwd():
         OutputWriter(grid, "test.nc", rows=(1,)).run_one_step()
         with xr.open_dataset("test.nc") as ds:
-            assert np.all(ds.x_of_node == approx([0, 1, 2, 3]))
-            assert np.all(ds.y_of_node == approx(1.0))
+            assert np.all(ds.x_of_column == approx([0, 1]))
+            assert np.all(ds.y_of_row == approx(0.0))
 
-        OutputWriter(grid, "test.nc", rows=((0, 2),), clobber=True).run_one_step()
+        # OutputWriter(grid, "test.nc", rows=((0, 2),), clobber=True).run_one_step()
+        OutputWriter(grid, "test.nc", rows=(0, 2), clobber=True).run_one_step()
         with xr.open_dataset("test.nc") as ds:
-            assert np.all(ds.x_of_node == approx([0, 1, 2, 3, 0, 1, 2, 3]))
-            assert np.all(ds.y_of_node == approx([0, 0, 0, 0, 2, 2, 2, 2]))
+            assert np.all(ds.x_of_column == approx([0, 1]))
+            assert np.all(ds.y_of_row == approx([2, 1]))
diff --git a/tests/test_sequence.py b/tests/test_sequence.py
index b04b5fc..9ecdbe1 100644
--- a/tests/test_sequence.py
+++ b/tests/test_sequence.py
@@ -70,15 +70,15 @@ def test_sequence_x_of_shore_missing(grid):
 def test_sequence_x_of_shore(grid):
     seq = Sequence(grid, components=[ShorelineFinder(grid)])
 
-    assert "x_of_shore" in grid.at_grid
-    assert "x_of_shelf_edge" in grid.at_grid
+    assert "x_of_shore" in grid.at_row
+    assert "x_of_shelf_edge" in grid.at_row
 
     assert grid.at_layer_grid.number_of_layers == 1
-    assert grid.at_layer_grid["x_of_shore"][-1] == pytest.approx(20000.0)
-    assert grid.at_layer_grid["x_of_shelf_edge"][-1] > 20000.0
+    assert np.all(grid.at_layer_row["x_of_shore"][-1, :] == pytest.approx(20000.0))
+    assert np.all(grid.at_layer_row["x_of_shelf_edge"][-1, :] > 20000.0)
 
     seq.update()
 
     assert grid.at_layer_grid.number_of_layers == 2
-    assert grid.at_layer_grid["x_of_shore"][-1] == pytest.approx(20000.0)
-    assert grid.at_layer_grid["x_of_shelf_edge"][-1] > 20000.0
+    assert np.all(grid.at_layer_row["x_of_shore"][-1, :] == pytest.approx(20000.0))
+    assert np.all(grid.at_layer_row["x_of_shelf_edge"][-1, :] > 20000.0)
diff --git a/tests/test_sequence_model/marmara.yaml b/tests/test_sequence_model/marmara.yaml
index ae28f82..c3593e1 100644
--- a/tests/test_sequence_model/marmara.yaml
+++ b/tests/test_sequence_model/marmara.yaml
@@ -15,8 +15,7 @@ output:
   interval: 10
   filepath: marmara.nc
   clobber: True
-  rows:
-  - 1
+  rows: [0, 1, 2]
   fields:
   - sediment_deposit__thickness
   - kd
diff --git a/tests/test_subsidence.py b/tests/test_subsidence.py
index 7c7a9af..ce89d4d 100644
--- a/tests/test_subsidence.py
+++ b/tests/test_subsidence.py
@@ -73,7 +73,7 @@ def test_linear_subsidence(tmpdir):
         subsidence.run_one_step(dt=1.0)
         assert_array_equal(
             grid.get_profile("bedrock_surface__increment_of_elevation"),
-            [0.0, 10.0, 20.0, 30.0, 40.0],
+            [[0.0, 10.0, 20.0, 30.0, 40.0]],
         )
 
 
@@ -91,7 +91,7 @@ def test_add_to_existing_subsidence(tmpdir):
         subsidence.run_one_step(dt=1.0)
         assert_array_equal(
             grid.get_profile("bedrock_surface__increment_of_elevation"),
-            [1.0, 11.0, 21.0, 31.0, 41.0],
+            [[1.0, 11.0, 21.0, 31.0, 41.0]],
         )
 
 

From 4a524c56da566783c5c376b87f6c3284f9c882c7 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 23 Oct 2023 13:55:28 -0600
Subject: [PATCH 21/29] drop testing with python 3.9

---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index c968a56..7657d14 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -20,7 +20,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ["3.9", "3.10", "3.11"]
+        python-version: ["3.10", "3.11"]
 
     steps:
       - uses: actions/checkout@v2

From 62ebd07297401b71eb3cb544e4ffbdf7085a32d3 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 12:09:24 -0700
Subject: [PATCH 22/29] fix the default rows to write to the output file

---
 sequence/output_writer.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sequence/output_writer.py b/sequence/output_writer.py
index 56fa17f..9d816f9 100644
--- a/sequence/output_writer.py
+++ b/sequence/output_writer.py
@@ -54,7 +54,7 @@ def __init__(
         if rows is not None:
             self._rows = np.asarray(rows) - 1
         else:
-            self._rows = np.arange(grid.shape[0] - 2) + 1
+            self._rows = np.arange(grid.shape[0] - 2)
 
         self._time = 0.0
         self._step_count = 0

From b1b6e2bfa18395f838a66f06186c74e6bfc361e1 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 12:10:01 -0700
Subject: [PATCH 23/29] update default values for components

---
 sequence/sequence_model.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/sequence/sequence_model.py b/sequence/sequence_model.py
index 2dd3079..fc23b25 100644
--- a/sequence/sequence_model.py
+++ b/sequence/sequence_model.py
@@ -39,13 +39,12 @@ class SequenceModel:
         "flexure",
     )
     DEFAULT_PARAMS = {
-        "grid": {"n_cols": 100, "spacing": 1000.0},
+        "grid": {"shape": (1, 100), "spacing": (1.0, 1000.0)},
         "clock": {"start": 0.0, "stop": 600000.0, "step": 100.0},
         "output": {
             "interval": 10,
             "filepath": "sequence.nc",
             "clobber": True,
-            "rows": [1],
             "fields": ["sediment_deposit__thickness", "bedrock_surface__elevation"],
         },
         "submarine_diffusion": {

From c315df880dd27edcbb495871a17f297b385262d9 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 12:11:30 -0700
Subject: [PATCH 24/29] keep x_of_shore and x_of_shelf_edge arrays 2d

---
 sequence/plot.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sequence/plot.py b/sequence/plot.py
index 46a9f3c..276859d 100644
--- a/sequence/plot.py
+++ b/sequence/plot.py
@@ -246,8 +246,8 @@ def plot_file(
 
         try:
             thickness_at_layer = ds["at_layer:thickness"][:, row, :]
-            x_of_shore = ds["at_row:x_of_shore"].data.squeeze()
-            x_of_shelf_edge = ds["at_row:x_of_shelf_edge"].data.squeeze()
+            x_of_shore = ds["at_row:x_of_shore"].data
+            x_of_shelf_edge = ds["at_row:x_of_shelf_edge"].data
             bedrock = (
                 ds["at_node:bedrock_surface__elevation"]
                 .data.reshape((-1, ds.dims["row"] + 2, ds.dims["column"] + 2))

From 9681f9354dc4200f48d924a907a0b273c9851244 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 12:23:47 -0700
Subject: [PATCH 25/29] fix toml literal block

---
 docs/tutorials/cli/sealevel.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/tutorials/cli/sealevel.md b/docs/tutorials/cli/sealevel.md
index 8922851..e17ee09 100644
--- a/docs/tutorials/cli/sealevel.md
+++ b/docs/tutorials/cli/sealevel.md
@@ -27,7 +27,7 @@ we replace this section with the following,
 
 ```toml
 [sequence.sea_level]
-filepath = sealevel.csv
+filepath = "sealevel.csv"
 ```
 
 With this configuration, *Sequence* will read sea level values from this

From a2f77bed0150c59138e22b9cb70fee6b41341c87 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 12:37:57 -0700
Subject: [PATCH 26/29] fix job name

---
 .github/workflows/notebooks.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/notebooks.yml b/.github/workflows/notebooks.yml
index 1e5fdcf..e30e8e1 100644
--- a/.github/workflows/notebooks.yml
+++ b/.github/workflows/notebooks.yml
@@ -33,5 +33,5 @@ jobs:
           - name: Install dependencies
             run: pip install nox
 
-          - name: Build documentation
+          - name: Test notebooks
             run: nox -s test-notebooks

From ac4dbc185e0a2a232da6671785c051228fc01d42 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 12:39:22 -0700
Subject: [PATCH 27/29] fix default row to plot for plot_grid

---
 sequence/plot.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/sequence/plot.py b/sequence/plot.py
index 276859d..4f78fd0 100644
--- a/sequence/plot.py
+++ b/sequence/plot.py
@@ -187,7 +187,7 @@ def plot_grid(grid: SequenceModelGrid, row: Optional[int] = None, **kwds: Any) -
     :func:`plot_file` : Plot a `SequenceModelGrid`'s layers from a file.
     """
     if row is None:
-        row = grid.number_of_rows // 2
+        row = (grid.number_of_rows - 2) // 2
 
     elevation_at_layer = (
         grid.get_profile("bedrock_surface__elevation")[row, 1:-1] + grid.at_layer.z
@@ -201,8 +201,10 @@ def plot_grid(grid: SequenceModelGrid, row: Optional[int] = None, **kwds: Any) -
     time_at_layer_grid = grid.at_layer_grid["age"].flatten()
     time_at_layer = grid.at_layer["age"]
 
-    x_of_shore = interp1d(time_at_layer_grid, x_of_shore[:, row])(time_at_layer[:, 0])
-    x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge[:, row])(
+    x_of_shore = interp1d(time_at_layer_grid, x_of_shore[:, row].squeeze())(
+        time_at_layer[:, 0]
+    )
+    x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge[:, row].squeeze())(
         time_at_layer[:, 0]
     )
 

From 4e990c9744a0cc5e57ad958dcf108ce915e0ad90 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 13:00:46 -0700
Subject: [PATCH 28/29] add news fragment

---
 news/78.feature | 3 +++
 1 file changed, 3 insertions(+)
 create mode 100644 news/78.feature

diff --git a/news/78.feature b/news/78.feature
new file mode 100644
index 0000000..dd4171d
--- /dev/null
+++ b/news/78.feature
@@ -0,0 +1,3 @@
+
+Added a quasi-2d mode for *Sequence* that keeps track of multiple cross-shore
+rows.

From eef79c7dc4319d74f49780106a4ccdca763d7120 Mon Sep 17 00:00:00 2001
From: mcflugen <mcflugen@gmail.com>
Date: Mon, 20 Nov 2023 15:49:47 -0700
Subject: [PATCH 29/29] update input file docs for sequence.grid section

---
 docs/reference/cli/index.md       |  2 +-
 docs/reference/cli/input_files.md | 27 +++++++++++++++++++++++++--
 2 files changed, 26 insertions(+), 3 deletions(-)

diff --git a/docs/reference/cli/index.md b/docs/reference/cli/index.md
index 2e868de..d8e492c 100644
--- a/docs/reference/cli/index.md
+++ b/docs/reference/cli/index.md
@@ -1,4 +1,4 @@
-# ``sequence``
+# Command-line interface
 
 The *Sequence* model can be run from the command line using the
 *sequence* program.
diff --git a/docs/reference/cli/input_files.md b/docs/reference/cli/input_files.md
index 71d77b9..e0b73a1 100644
--- a/docs/reference/cli/input_files.md
+++ b/docs/reference/cli/input_files.md
@@ -58,11 +58,34 @@ spacing = 1000.0
 ```
 
 In this case we have a grid that represents a 1D profile that consists of
-500 columns (i.e. vertical stacks) of sediment (the *n_cols* parameter).
+100 columns (i.e. vertical stacks) of sediment (the *n_cols* parameter).
 
 The *spacing* parameter is the width of each of your sediment stacks in meters.
 Thus, the length of you domain is the product of the number of columns with
-the spacing (that is, for this example, 500 * 100 m or 50 km).
+the spacing (that is, for this example, 100 * 1000 m or 100 km).
+
+You can also run *Sequence* in a quasi-2d mode in which *Sequence* models parallel
+cross-shore profiles. Each profile will have its own sediment supply that is
+transported into the ocean and then can be diffused between profiles. To set up
+such a case, the grid section can, for example, be modified in the following way,
+
+```toml
+[sequence.grid]
+shape = [3, 100]
+spacing = [10000.0, 1000.0]
+```
+In this case *Sequence* will create three parallel profiles (with an along-shore
+width of 10000.0 m).
+
+:::{note}
+The following is equivalent to the 1D example above,
+
+```toml
+[sequence.grid]
+shape = [1, 100]
+spacing = [1.0, 1000.0]
+```
+:::
 
 (the-output-section)=