From 21219ab0148415d597be488c97843e2912bb4ad9 Mon Sep 17 00:00:00 2001
From: Trevor Hillebrand <trhille@lanl.gov>
Date: Fri, 7 Aug 2026 12:18:54 -0700
Subject: [PATCH] Seed AIS bed flood fill from open ocean

The Antarctica mesh-gen step seeded the bed-topography flood fill at the
grid center, which is thick interior ice far from any grounding line. The
high_dist_bed cap zeroes that cell, so gridded_flood_fill returned an empty
mask and use_bed produced no refinement on the full ice sheet.

Add a 'flood_fill_start=auto' mode to build_cell_width that seeds the flood
fill from the deepest ice-free, below-sea-level (open ocean) cell, chosen by
data values so it is independent of grid resolution. Use it in the Antarctica
step. Existing callers using [None, None] or an explicit seed are unchanged.
---
 compass/landice/mesh.py                  | 24 ++++++++++++++++++++----
 compass/landice/tests/antarctica/mesh.py |  6 +-----
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py
index 2a2b86fd38..8d4116e9e4 100644
--- a/compass/landice/mesh.py
+++ b/compass/landice/mesh.py
@@ -616,10 +616,13 @@ def build_cell_width(self, section_name, gridded_dataset,
     gridded_dataset : str
         name of NetCDF file used to define cell spacing
 
-    flood_fill_start : list of ints
-        ``i`` and ``j`` indices used to define starting location for flood
-        fill. Most cases will use ``[None, None]``, which will just start the
-        flood fill in the center of the gridded dataset.
+    flood_fill_start : list of ints or str
+        ``i`` and ``j`` indices used to define starting location for the
+        bed-topography flood fill. Most cases will use ``[None, None]``, which
+        skips the flood fill. Pass ``'auto'`` to seed the flood fill from the
+        deepest ice-free, below-sea-level cell (open ocean), which is the
+        appropriate choice for marine ice sheets and is independent of the
+        gridded dataset's resolution.
 
     Returns
     -------
@@ -659,6 +662,19 @@ def build_cell_width(self, section_name, gridded_dataset,
 
     f.close()
 
+    # Seed the bed flood fill from open ocean (deepest ice-free, below-sea-
+    # level cell) when requested, so it lands in the ocean regardless of grid
+    # resolution rather than on ice-covered cells that get masked out.
+    if isinstance(flood_fill_start, str) and flood_fill_start == 'auto':
+        ocean = np.logical_and(thk == 0.0, topg < 0.0)
+        if not ocean.any():
+            raise ValueError(
+                "flood_fill_start='auto' requires at least one ice-free, "
+                'below-sea-level (ocean) cell in the gridded dataset, but '
+                'none were found.')
+        flood_fill_start = list(np.unravel_index(
+            np.argmin(np.where(ocean, topg, np.inf)), topg.shape))
+
     # Get bounds defined by user, or use bounds from the gridded dataset.
     bnds = get_mesh_config_bounding_box(
         section,
diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py
index d18b80ebcc..877b8efa78 100644
--- a/compass/landice/tests/antarctica/mesh.py
+++ b/compass/landice/tests/antarctica/mesh.py
@@ -1,5 +1,4 @@
 import netCDF4
-import xarray as xr
 from mpas_tools.logging import check_call
 
 from compass.landice.mesh import (
@@ -81,15 +80,12 @@ def run(self):
         else:
             bm_updated_gridded_dataset = source_gridded_dataset
 
-        ds = xr.open_dataset(bm_updated_gridded_dataset)
-        nx, ny = ds.sizes["x1"], ds.sizes["y1"]
-        ds.close()
         logger.info('calling build_cell_width')
         cell_width, x1, y1, geom_points, geom_edges, floodFillMask = \
             build_cell_width(
                 self, section_name=section_name,
                 gridded_dataset=bm_updated_gridded_dataset,
-                flood_fill_start=[nx // 2, ny // 2])
+                flood_fill_start='auto')
 
         # Now build the base mesh and perform the standard interpolation
         build_mali_mesh(