Every regional MOM6 domain has 4 cardinal open boundaries — the outer north/south/east/west edges. But sometimes the boundary you actually want is inside the domain: a strait behind a headland, a channel between two landmasses, a gap in an otherwise-closed coastline. MOM6 supports this fine — an OBC segment is really just “a straight line, somewhere on the grid” — but building one by hand exposes a few conventions that the 4 cardinal boundaries never make you think about, because they always land on the same easy special case.
This notebook builds the smallest possible domain that needs one: a single
small peninsula, with an interior segment feeding the water behind it. It
uses regional_mom6.segment.Segment directly (via Case.configure_forcings’s
boundaries argument, which accepts a mix of cardinal strings and Segment
objects) rather than the legacy cardinal-only path.
Step 1: A domain with something to cut through¶
A flat, 12x8 T-cell domain, with one 3-row-tall, 4-column-wide peninsula hanging down from the north edge — not spanning the full width, so the ocean stays connected around both sides of it.
from CrocoDash.grid import Grid
from CrocoDash.topo import Topo
NX, NY = 12, 8
RESOLUTION = 1.0 # degrees -- this domain is synthetic, so the value itself doesn't matter
grid = Grid(
lenx=NX * RESOLUTION,
leny=NY * RESOLUTION + 1e-6, # avoid a floating-point rounding grid.ny down by 1
resolution=RESOLUTION,
xstart=-70.0,
ystart=20.0,
name="interior_demo",
)
topo = Topo(grid, min_depth=10.0, git=False)
topo.set_flat(1000.0)
# T-cell indices (0 = south/west). The peninsula: rows 5-7 (the top 3 rows),
# columns 4-7 (the middle 4 of 12) -- land everywhere else stays open ocean.
PENINSULA_ROWS = (5, 7)
PENINSULA_COLS = (4, 7)
depth = topo.depth.values.copy()
depth[PENINSULA_ROWS[0] : PENINSULA_ROWS[1] + 1, PENINSULA_COLS[0] : PENINSULA_COLS[1] + 1] = 0.0
topo.depth = depth
# Defensive, same as any hand-edited mask: removes one-cell-wide channels/lakes
# and diagonal-only land touches, which MOM6's own mask consistency check
# (see Step 5) would otherwise reject. A no-op here, but good habit.
topo.fill_inland_lakes_and_channels()
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6, 4))
ax.imshow(topo.depth > 0, origin="lower", cmap="Blues", vmin=-0.3, vmax=1.3)
ax.set_xlabel("T-column (I index)")
ax.set_ylabel("T-row (J index)")
ax.set_title("interior_demo: ocean (blue) vs. the peninsula (white)")
plt.show()
Step 2: Where exactly does the interior segment go?¶
The natural place for an open boundary here is the peninsula’s own southern base: T-row 4, spanning T-columns 4-7 (right under the land). Water flowing in there feeds the water south of the peninsula, same as any other boundary.
Building this by hand means getting three things right that a cardinal edge never has to worry about:
1. Supergrid vs. native indices. hgrid (what Segment.from_hgrid slices)
is at 2x MOM6’s own resolution — even indices are cell corners, odd indices
are T-cell centers. For T-row k, the supergrid index is 2 * k + 1. Our
strait sits at T-row 4, so index=9; it spans T-columns 4-7, so
index_range=slice(2*4+1, 2*7+2) = slice(9, 16) (7 supergrid points — odd,
which MOM6’s BRUSHCUTTER_MODE requires; Segment checks this for you and
raises a clear error if you get it wrong).
2. T-points vs. U-/V-points. A horizontal line (fixed J, our case)
sits on MOM6’s V-point grid, because the velocity component that flows
through a horizontal boundary is the meridional one, V — this is baked
into MOM6 itself (the Fortran routine that parses "J=..." strings is
literally called setup_v_point_obc). A vertical line (fixed I) would sit
on the U-point grid instead, for the same reason with the zonal
component. (axis="nyp" below means “J is fixed” — a horizontal line.)
3. Which side is the interior. MOM6 requires the domain interior to stay
on the same side all the way around a segment. Our open water is south of
the strait, the peninsula is north — so mom6_index_reverse=True (for a
horizontal line: False means interior-is-north, True means
interior-is-south). Get this backwards and MOM6 forces the wrong side to
land at runtime, since it always trusts this direction over the bathymetry
file.
(For the full, precise version of all of this — including the land-masking
mechanics Segment handles for you automatically — see the “Indexing
conventions” section of the Segment class docstring in
regional_mom6/segment.py.)
import xarray as xr
from regional_mom6.segment import Segment
hgrid = grid._supergrid.to_ds(name=grid.name, author="interior_obc_segments demo")
strait = Segment.from_hgrid(
hgrid,
axis="nyp", # J fixed -> a horizontal line -> V-point grid
index=2 * 4 + 1, # supergrid index for T-row 4 (9)
index_range=slice(2 * 4 + 1, 2 * 7 + 2), # T-columns 4-7 (slice(9, 16))
segment_name="strait",
topo=topo,
mom6_index_reverse=True, # interior (ocean) is SOUTH of this line
)
fig, ax = plt.subplots(figsize=(6, 4))
ax.imshow(topo.depth > 0, origin="lower", cmap="Blues", vmin=-0.3, vmax=1.3)
lo, hi = PENINSULA_COLS
ax.plot([lo, hi], [4, 4], color="crimson", lw=4, solid_capstyle="butt", label="strait (interior segment)")
ax.legend(loc="lower center")
ax.set_xlabel("T-column (I index)")
ax.set_ylabel("T-row (J index)")
ax.set_title("The interior segment sits at the peninsula's southern base")
plt.show()
Step 3: The other 3 boundaries are just cardinal edges¶
South, east, and west are full outer edges — nothing special about them.
Case.configure_forcings’s boundaries list accepts a mix of plain cardinal
strings and Segment objects side by side, in any order.
boundaries = ["south", "east", "west", strait]
Step 4: Build the CESM case¶
Same as the Getting Started tutorial
from here on — a vertical grid, then a Case.
from pathlib import Path
from CrocoDash.vgrid import VGrid
from CrocoDash.case import Case
vgrid = VGrid.uniform(nk=10, depth=1000.0, name="interior_demo")
casename = "interior_obc_segments_demo"
cesmroot = "<CESM>"
inputdir = Path.home() / "croc_input" / casename
caseroot = Path.home() / "croc_cases" / casename
case = Case(
cesmroot=cesmroot,
caseroot=caseroot,
inputdir=inputdir,
ocn_grid=grid,
ocn_vgrid=vgrid,
ocn_topo=topo,
compset="1850_DATM%NYF_SLND_SICE_MOM6%REGIONAL_SROF_SGLC_SWAV",
atm_grid_name="T62",
machine="ubuntu-latest",
project="NCGD0011",
override=True,
)
Step 5: Forcings¶
configure_forcings takes the boundaries list from Step 3 directly — the
interior strait segment goes through exactly the same download/regrid
pipeline as the 3 cardinal ones. This is also the point where MOM6’s mask
consistency requirement gets checked: the T-cell immediately north of the
strait (the peninsula’s own base row) must already be land in the bathymetry
file, or the case would fail at runtime with an “ESMF mesh and MOM6 domain
masks are inconsistent” error. It already is here (Step 1), so nothing more
to do — but this is the check to remember if you build your own interior
segment and hit that error.
This demo doesn’t need real ocean data — a real case would download GLORYS
here, but since this domain doesn’t correspond to anywhere on Earth, we
write a tiny synthetic dataset instead in the same on-disk format
process_forcings expects (the download step is skipped automatically once
those files already exist).
case.configure_forcings(
date_range=["2020-01-01 00:00:00", "2020-01-03 00:00:00"],
boundaries=boundaries,
)
import numpy as np
import pandas as pd
raw_dir = inputdir / "extract_forcings" / "raw_data"
raw_dir.mkdir(parents=True, exist_ok=True)
lon = np.arange(-72.0, -56.0, 0.25)
lat = np.arange(18.0, 32.0, 0.25)
depth_levels = np.array([0, 10, 25, 50, 100, 200, 500, 1000.0])
dates = pd.date_range("2020-01-01", "2020-01-03", freq="D")
nt, nz, ny, nxg = len(dates), len(depth_levels), len(lat), len(lon)
synthetic = xr.Dataset(
{
"uo": (("time", "depth", "latitude", "longitude"), np.zeros((nt, nz, ny, nxg), "float32")),
"vo": (("time", "depth", "latitude", "longitude"), np.zeros((nt, nz, ny, nxg), "float32")),
"thetao": (("time", "depth", "latitude", "longitude"), np.full((nt, nz, ny, nxg), 15.0, "float32")),
"so": (("time", "depth", "latitude", "longitude"), np.full((nt, nz, ny, nxg), 35.0, "float32")),
"zos": (("time", "latitude", "longitude"), np.zeros((nt, ny, nxg), "float32")),
},
coords={"time": dates, "depth": depth_levels, "latitude": lat, "longitude": lon},
)
for var, units in [("uo", "m s-1"), ("vo", "m s-1"), ("thetao", "degC"), ("so", "PSU"), ("zos", "m")]:
synthetic[var].attrs["units"] = units
synthetic["depth"].attrs["units"] = "m"
for name in ("south", "east", "west", "strait"):
synthetic.to_netcdf(raw_dir / f"{name}_unprocessed.2020-01-01_2020-01-03.nc")
synthetic.isel(time=[0]).to_netcdf(raw_dir / "ic_unprocessed.nc")
case.process_forcings()
Step 6: Reading the result¶
process_forcings already wrote the OBC_SEGMENT_00N lines into
user_nl_mom — including one for strait, MOM6’s own numeric encoding of
everything decided in Step 2.
user_nl_mom = caseroot / "user_nl_mom"
for line in user_nl_mom.read_text().splitlines():
if line.strip().startswith("OBC_SEGMENT_00") and "DATA" not in line and "TIMESCALE" not in line:
print(line)
The strait line reads J=5,I=7:4, not J=4,I=7:4 — even though we placed
it at T-row 4. That’s not a mistake: it’s Segment automatically
compensating for an asymmetry in how MOM6 masks land next to a segment like
this one (see Step 2’s docstring pointer for the full mechanics). The number
you get back is the row MOM6 will force to land — here, J=5, the
peninsula’s own base row, exactly as intended — not the segment’s own row.
I=7:4 (descending) is the mom6_index_reverse=True from Step 2, telling
MOM6 which direction the interior side is in.
You never need to do this arithmetic by hand — pick index/index_range to
point at the row/column you want the segment’s own data to come from, and
mom6_index_reverse to say which side is interior; Segment takes care of
the rest.
Step 7: Build and run¶
From here it’s identical to any other CrocoDash case — review user_nl_mom
in caseroot, then:
./case.setup
./case.build
./case.submit