Add regional Thwaites mesh to global_ocean test group

compass/ocean/mesh/cull.py

@@ -305,6 +305,17 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
                                 geojson_filename='land_coverage.geojson',
                                 mask_filename='land_mask.nc')
 
+    # Hook point: allow modification of land mask before culling
+    # This enables regional culling and other mask modifications
+    if os.path.exists('modify_land_mask_hook.py'):
+        logger.info('Running land mask modification hook...')
+        import importlib.util
+        spec = importlib.util.spec_from_file_location("hook", "modify_land_mask_hook.py")
+        hook = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(hook)
+        if hasattr(hook, 'modify_land_mask'):
+            hook.modify_land_mask(logger=logger)
+
     dsBaseMesh = xr.open_dataset('base_mesh.nc')
     dsLandMask = xr.open_dataset('land_mask.nc')
 

compass/ocean/tests/global_ocean/mesh/__init__.py

@@ -12,6 +12,8 @@
 from compass.ocean.tests.global_ocean.mesh.fris02to60 import FRIS02to60BaseMesh
 from compass.ocean.tests.global_ocean.mesh.fris04to60 import FRIS04to60BaseMesh
 from compass.ocean.tests.global_ocean.mesh.fris08to60 import FRIS08to60BaseMesh
+from compass.ocean.tests.global_ocean.mesh.thwaites01to60 import Thwaites01to60BaseMesh
+from compass.ocean.tests.global_ocean.mesh.thwaites01to60.cull_mesh import ThwaitesCullMeshStep
 from compass.ocean.tests.global_ocean.mesh.kuroshio import KuroshioBaseMesh
 from compass.ocean.tests.global_ocean.mesh.qu import (
     IcosMeshFromConfigStep,
@@ -138,6 +140,8 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
             base_mesh_step = FRIS04to60BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['FRIS08to60', 'FRISwISC08to60']:
             base_mesh_step = FRIS08to60BaseMesh(self, name=name, subdir=subdir)
+        elif mesh_name in ['Thwaites01to60', 'ThwaitesWISC01to60']:
+            base_mesh_step = Thwaites01to60BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name.startswith('Kuroshio'):
             base_mesh_step = KuroshioBaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['WC14', 'WCwISC14']:
@@ -191,10 +195,17 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
 
             self.add_step(smoothed_topo)
 
-        self.add_step(CullMeshStep(
-            test_case=self, base_mesh_step=base_mesh_step,
-            with_ice_shelf_cavities=self.with_ice_shelf_cavities,
-            unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
+        # Use custom cull step for Thwaites meshes (supports regional culling)
+        if mesh_name.startswith('Thwaites'):
+            self.add_step(ThwaitesCullMeshStep(
+                test_case=self, base_mesh_step=base_mesh_step,
+                with_ice_shelf_cavities=self.with_ice_shelf_cavities,
+                unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
+        else:
+            self.add_step(CullMeshStep(
+                test_case=self, base_mesh_step=base_mesh_step,
+                with_ice_shelf_cavities=self.with_ice_shelf_cavities,
+                unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
 
     def configure(self, config=None):
         """

compass/ocean/tests/global_ocean/mesh/thwaites01to60/README.md

@@ -0,0 +1,141 @@
+# Thwaites01to60 Mesh
+
+Variable-resolution global ocean mesh with 200m refinement at Thwaites grounding zone.
+
+## Current Status
+
+### ✅ Implemented:
+1. **Base mesh generation** (`__init__.py`)
+   - `Thwaites01to60BaseMesh` class
+   - **Grounding line extraction from BedMachine** (automatic, no geojson needed)
+   - Calculates flotation criterion: thickness > -bed × (ρ_ocean/ρ_ice)
+   - Identifies grounding line as boundary between grounded/floating ice
+   - Creates signed distance field for mesh refinement
+   - Falls back to Gaussian approximation if BedMachine unavailable
+   - Configurable resolutions: 200m GZ, 1km cavity, 3km shelf, 8km far-field
+
+2. **Configuration** (`thwaites01to60.cfg`)
+   - Resolution parameters
+   - Regional domain settings (enabled by default for Amundsen sector)
+   - Thin film parameters (disabled by default)
+
+3. **Regional culling** ✅ **FULLY INTEGRATED**
+   - `ThwaitesCullMeshStep` custom cull step (`cull_mesh.py`)
+   - Hook-based integration into standard culling workflow
+   - Modifies land mask to mark cells outside domain as land
+   - Supports geojson polygon or lat-lon bounding box
+   - **Working and ready to test**
+
+4. **Registration**
+   - Mesh registered in `compass/ocean/tests/global_ocean/mesh/__init__.py`
+   - Custom CullMeshStep wired in for Thwaites meshes
+   - Ready for `compass list` and `compass setup`
+
+### 🚧 TODO:
+
+#### High Priority:
+1. ~~**Integrate regional culling into CullMeshStep**~~ ✅ **DONE**
+   - ✅ Added hook point in base `cull.py`
+   - ✅ Created `ThwaitesCullMeshStep` with hook file generation
+   - ✅ Wired into Mesh test case for Thwaites meshes
+   - ✅ Ready for testing
+
+2. ~~**Replace Gaussian refinement with grounding-line-based refinement**~~ ✅ **DONE**
+   - ✅ Extracts grounding line directly from BedMachine
+   - ✅ Calculates flotation criterion automatically
+   - ✅ Uses `signed_distance_from_geojson` for accurate GZ band
+   - ✅ No manual geojson file required
+   - ⚠️ Requires BedMachine file in bathymetry database
+
+#### Medium Priority:
+3. **Add thin film support**
+   - Modify land mask to keep cells beneath grounded ice
+   - Based on height above flotation from BedMachine
+   - Requires integration with `remap_topography` step
+
+4. **Testing**
+   - Set up and run basic mesh generation
+   - Validate cell widths and resolution
+   - Test regional culling with lat-lon bounds
+   - Test with/without thin film
+
+#### Low Priority:
+5. **Clean up copied FRIS files**
+   - Remove unused geojson files (atlantic.geojson, fris_*.geojson, etc.)
+   - Remove or adapt FRIS-specific namelists if needed
+
+## Usage
+
+### Basic mesh (no regional culling):
+```bash
+compass list | grep Thwaites
+compass setup -t global_ocean/mesh/Thwaites01to60 -w $WORK
+compass run $WORK
+```
+
+### With regional domain (enabled by default):
+Regional culling is **enabled by default** for Amundsen sector (76-73°S, 116-98°W).
+
+To change the domain, edit `thwaites01to60.cfg`:
+```ini
+[thwaites01to60]
+# Option A: lat-lon bounds (currently enabled)
+lat_min = -76.0
+lat_max = -73.0
+lon_min = -116.0
+lon_max = -98.0
+
+# Option B: geojson polygon
+# regional_domain_geojson = amundsen_domain.geojson
+```
+
+To disable regional culling (full global mesh), comment out all domain options.
+
+## Implementation Notes
+
+### Mesh Refinement:
+The implementation now extracts the grounding line directly from BedMachine:
+
+1. **Load BedMachine data** (ice thickness, bed elevation)
+2. **Calculate flotation criterion**:
+   ```python
+   # Ice is grounded when: thickness > -bed × (ρ_ocean/ρ_ice)
+   thickness_flotation = -bed * (1028.0 / 918.0)  # where bed < 0
+   grounded = thickness > thickness_flotation
+   ```
+3. **Identify grounding line**: Boundary between grounded and floating ice
+4. **Create feature collection**: Convert GL cells to geojson polygon
+5. **Calculate signed distance**: Use `signed_distance_from_geojson`
+6. **Apply refinement**:
+   - GZ band: ±15 km from GL → 200m resolution
+   - Cavity: 100 km into cavity → 1 km resolution  
+   - Shelf: 200 km into cavity → 3 km resolution
+   - Far-field: everywhere else → 8 km resolution
+
+**Fallback**: If BedMachine is unavailable, uses Gaussian approximation around Thwaites location (75°S, 106°W).
+
+### Regional Culling:
+**Fully implemented and integrated!**
+
+The `ThwaitesCullMeshStep` (in `cull_mesh.py`) creates a hook file that:
+1. Is executed by the modified `cull.py` after land mask creation
+2. Loads base mesh and land mask
+3. Determines cells inside/outside regional domain
+4. Marks outside cells as "land"
+5. Standard culling then removes these cells
+
+The hook runs between lines 307 and 308 of `compass/ocean/mesh/cull.py`,
+after `land_mask.nc` is created but before culling begins.
+
+### Thin Film:
+Will require:
+1. BedMachine topography (from `remap_topography` step)
+2. Compute height above flotation for grounded ice
+3. Cells with HAF < threshold (e.g., 30m) kept as ocean
+4. Add `thinFilmMask` field for diagnostics
+
+## References
+
+- Design document: `/Users/trhille/Documents/Antarctica/thwaites_grounding_zone_intrusion/Grounding_zone_intrusion_simulation_design.pdf`
+- Implementation plan: `/Users/trhille/Documents/Antarctica/thwaites_grounding_zone_intrusion/QUICK_START.md`
+- FRIS mesh (template): `compass/ocean/tests/global_ocean/mesh/fris01to60/`