One Way Nesting a Regional Ocean Model: Strait of Gibraltar Example

This notebook demonstrates MOM6-in-CESM domain nesting: running progressively finer child domains whose open-boundary conditions (OBCs) come from the next coarser parent, rather than from a global reanalysis like GLORYS.

Feature: the Strait of Gibraltar, the narrow (~14 km wide) sill connecting the Atlantic and the Mediterranean. Warm, salty Mediterranean water sinks over the Camarinal Sill (~300 m) and cascades down the Gulf of Cadiz slope as Mediterranean Overflow Water, spinning off eddies (“Meddies”), a classic hydraulic-control/overflow process that only shows up once resolution is fine enough to resolve the sill itself.

Each level in this hierarchy refines by 4× (the max recommended child:parent ratio), so the tiny child reaches sub-kilometer, “hyper-resolution”, while every domain stays small enough to build and run on ~10 cores.

Domain hierarchy¶

Parent      (gibraltar_parent)      1/12°   (~7.5-9 km)   -8.0→-2.0 lon,  34.5→37.0 lat   (Gulf of Cadiz + Alboran approach)
│  4× refinement
└── Middle child                     1/48°   (~1.9-2.3 km) -6.0→-5.0 lon,  35.75→36.25 lat  (strait mouth)
    │  4× refinement
    └── Tiny child                   1/192°  (~0.47-0.58 km) -5.9→-5.65 lon, 35.83→35.955 lat  (Camarinal Sill)

Rule of thumb: each child should be no more than 4× finer than its parent (4:1 ratio). This hierarchy now sits right at that maximum every step, rather than the more conservative 2× used previously. The parent alone (1/12°) can’t resolve the 14 km-wide strait at all; the middle child (1/48°) resolves the strait mouth with several points across it; the tiny child (1/192°) zooms in tight enough to resolve the Camarinal Sill itself at sub-kilometer resolution.

Workflow¶

The three domains are built up together, step by step, rather than one domain at a time. Each step below does the same thing for the parent and both children, since the calls are basically identical across domains:

Paths: case/input paths for all three domains
Grids: horizontal grid for all three (created fresh. There’s no existing reference case for Gibraltar, unlike a production run)
Topography: bathymetry for all three, from GEBCO
Vertical grids: hyperbolic vertical grid for all three
CESM cases: create all three cases. Each pinned to a small NTASKS_OCN (~10 cores or fewer)
Forcings: configure OBCs for all three: parent from GLORYS, children from CESM_MOM_OUTPUT.get_mom6_output_data() reading their parent’s nesting output

Setup: imports¶

%load_ext autoreload
%autoreload 2
from pathlib import Path
from CrocoDash.grid import Grid
from CrocoDash.topo import Topo
from CrocoDash.vgrid import VGrid
from CrocoDash.case import Case

Step 1: Paths (parent, middle child, tiny child)¶

CESMROOT = Path("~/work/installs/CROCESM_workshop_2025").expanduser()
PROJECT  = "ncgd0011"
GEBCO    = Path("<GEBCO>")

# --- Parent (Gulf of Cadiz + Alboran approach, 1/12°) ---
PARENT_CASE  = Path("~/croc_cases/gibraltar_parent").expanduser()
PARENT_INPUT = Path("<inputdir>/gibraltar_parent")

# --- Middle child (strait mouth, 1/24°) ---
MID_CASE  = Path("~/croc_cases/gibraltar_mid_child").expanduser()
MID_INPUT = Path("<inputdir>/gibraltar_mid_child")

# --- Tiny child (Camarinal Sill / Tarifa Narrows, 1/48°) ---
TINY_CASE  = Path("~/croc_cases/gibraltar_tiny_child").expanduser()
TINY_INPUT = Path("<inputdir>/gibraltar_tiny_child")

# --- Nesting sources (each child's OBC comes from the next coarser domain's own MOM6 output) ---
# See the note in Step 6 for what this directory should actually contain -- it's not just
# "point this at wherever the parent archives its history output."
PARENT_NESTING_OUTPUT = Path("~/scratch/archive/gibraltar_parent/ocn/hist/z_files").expanduser()
MID_NESTING_OUTPUT    = Path("~/scratch/archive/gibraltar_mid_child/ocn/hist/z_files").expanduser()

Step 2: Horizontal grids (parent, middle child, tiny child)¶

# Parent: 1/12°, Gulf of Cadiz + western Alboran Sea approach to the strait
parent_grid = Grid(
    resolution = 1/12,
    xstart     = -8.0,
    lenx       = 6.0,    # -8.0 → -2.0
    ystart     = 34.5,
    leny       = 2.5,    # 34.5 → 37.0
    name       = "gibraltar_parent",
)

# Middle child: 1/48°, 4× refinement over the parent, zoomed on the strait mouth
mid_grid = Grid(
    resolution = 1/48,
    xstart     = -6.0,
    lenx       = 1.0,    # -6.0 → -5.0
    ystart     = 35.75,
    leny       = 0.5,    # 35.75 → 36.25
    name       = "gibraltar_mid_child",
)

# Tiny child: 1/192°, 4× refinement over the middle child, zoomed on the Camarinal Sill
tiny_grid = Grid(
    resolution = 1/192,
    xstart     = -5.9,
    lenx       = 0.25,   # -5.9 → -5.65
    ystart     = 35.83,
    leny       = 0.125,  # 35.83 → 35.955
    name       = "gibraltar_tiny_child",
)

print(f"Parent grid: {parent_grid}")
print(f"Middle child grid: {mid_grid}")
print(f"Tiny child grid: {tiny_grid}")

Step 3: Bathymetry / topography (parent, middle child, tiny child)¶

# Parent: from GEBCO
parent_topo = Topo(grid=parent_grid, min_depth=10.0)
parent_topo.set_from_dataset(
    bathymetry_path           = GEBCO,
    longitude_coordinate_name = "lon",
    latitude_coordinate_name  = "lat",
    vertical_coordinate_name  = "elevation",
)
parent_topo.depth.plot(figsize=(8, 4))
print(f"Parent topo depth range: {float(parent_topo.depth.min()):.1f} – {float(parent_topo.depth.max()):.1f} m")

# Middle child: from GEBCO
mid_topo = Topo(grid=mid_grid, min_depth=10.0)
mid_topo.set_from_dataset(
    bathymetry_path           = GEBCO,
    longitude_coordinate_name = "lon",
    latitude_coordinate_name  = "lat",
    vertical_coordinate_name  = "elevation",
)
mid_topo.depth.plot(figsize=(6, 4))

# Tiny child: from GEBCO; at 1/192° this should clearly resolve the Camarinal Sill (~300 m)
# rising out of the deeper strait channel on either side
tiny_topo = Topo(grid=tiny_grid, min_depth=10.0)
tiny_topo.set_from_dataset(
    bathymetry_path           = GEBCO,
    longitude_coordinate_name = "lon",
    latitude_coordinate_name  = "lat",
    vertical_coordinate_name  = "elevation",
)
tiny_topo.depth.plot(figsize=(6, 4));

Step 4: Vertical grids (parent, middle child, tiny child)¶

# Parent: 20 levels, hyperbolic stretching
parent_vgrid = VGrid.hyperbolic(
    nk    = 20,
    depth = float(parent_topo.depth.max()),
    ratio = 10.0,
)

# Middle child: 30 levels, hyperbolic stretching
mid_vgrid = VGrid.hyperbolic(
    nk    = 30,
    depth = float(mid_topo.depth.max()),
    ratio = 10.0,
)

# Tiny child: 40 levels, hyperbolic stretching (finer near-surface resolution to capture the overflow layer)
tiny_vgrid = VGrid.hyperbolic(
    nk    = 40,
    depth = float(tiny_topo.depth.max()),
    ratio = 10.0,
)

Step 5: Create the CESM cases (parent, middle child, tiny child)¶

All three grids are small enough that NTASKS_OCN is pinned to 10 cores or fewer. This demo is meant to build and run quickly, not to scale up.

COMPSET = "1850_DATM%NYF_SLND_SICE_MOM6%REGIONAL_SROF_SGLC_SWAV"

parent_case = Case(
    cesmroot      = CESMROOT,
    caseroot      = PARENT_CASE,
    inputdir      = PARENT_INPUT,
    compset       = COMPSET,
    ocn_grid      = parent_grid,
    ocn_topo      = parent_topo,
    ocn_vgrid     = parent_vgrid,
    machine       = "derecho",
    project       = PROJECT,
    ntasks_ocn    = 10,
    override      = True,
    atm_grid_name = "T62"
)

mid_case = Case(
    cesmroot      = CESMROOT,
    caseroot      = MID_CASE,
    inputdir      = MID_INPUT,
    compset       = COMPSET,
    ocn_grid      = mid_grid,
    ocn_topo      = mid_topo,
    ocn_vgrid     = mid_vgrid,
    machine       = "derecho",
    project       = PROJECT,
    ntasks_ocn    = 8,
    override      = True,
    atm_grid_name = "T62"
)

tiny_case = Case(
    cesmroot      = CESMROOT,
    caseroot      = TINY_CASE,
    inputdir      = TINY_INPUT,
    compset       = COMPSET,
    ocn_grid      = tiny_grid,
    ocn_topo      = tiny_topo,
    ocn_vgrid     = tiny_vgrid,
    machine       = "derecho",
    project       = PROJECT,
    ntasks_ocn    = 6,
    override      = True,
    atm_grid_name = "T62"
)

Step 6: Configure and process forcings (parent, middle child, tiny child)¶

The parent gets its OBCs from GLORYS (there’s no existing reference case to reuse for Gibraltar, unlike a production run). Both children instead read their parent’s nesting output as OBCs, via CESM_MOM_OUTPUT.get_mom6_output_data.

Note

: CESM_MOM_OUTPUT (CrocoDash/raw_data_access/datasets/cesm_ocean_output.py) is a dedicated product for native MOM6 output, with MOM6-native variable/coordinate metadata (thetao/so/uo/vo/zos on xh/yh/xq/yq/z_l), split out from the CESM-POP-oriented CESM_POP_OUTPUT product so the two don’t share a mismatched metadata set. It has two access methods: get_mom6_output_data (multi-variable-per- file history/diag_table output, used below) and get_mom6_single_variable_data (if a parent run’s output is instead organized in the CESM single-variable-per-file/tseries convention). CESM_MOM_OUTPUT has no universal default dataset_path (it’s meant to point at any user-supplied directory). Each child passes its own parent’s nesting output directory in via function_overrides={"dataset_path": ...} (raw_data_access_flexible_default_args), instead of hand-editing config.json between configure_forcings() and process_forcings().

What actually needs to be in that directory: get_mom6_output_data’s default variables list is ["zos", "thetao", "so", "uo", "vo"]. Sea surface height, temperature, salinity, and the two horizontal velocity components. That’s the complete set MOM6’s OBC segments need (SSH/TEMP/SALT/U/V in the OBC_SEGMENT_*_DATA namelist, plus the same five for the initial condition). A parent run’s full history output, the .mom6.h.z.* and .mom6.h.native.* files used elsewhere in this notebook. Carries dozens of other diagnostics (h, rhopot0/rhopot2, difvho/difvso, Kv_u/Kv_v, transport terms, mixed-layer diagnostics, etc.) that a child domain has no use for and that just make the directory get_mom6_output_data has to read through far larger than it needs to be, the equivalent single-month file for the nesting demo is ~30 GB for exactly this reason. When you set up a parent run whose purpose is to feed a nested child, add a dedicated diag_table stream restricted to just those 5 fields (3D uo/vo/thetao/so plus 2D zos) and point dataset_path at that stream’s output directory, not at the full-history one. This demo’s parent/mid cases predate that setup and use their full z_files output (which happens to already isolate just the sfc- and z-level files, though still with extra variables beyond the 5 needed). Fine at this toy domain’s file sizes, but not a pattern to copy at production scale.</cell id=“cell-13”>

PARENT_DATES = ["2000-01-01 00:00:00", "2000-01-15 00:00:00"]
MID_DATES    = ["2000-01-01 00:00:00", "2000-01-15 00:00:00"]
TINY_DATES   = ["2000-01-01 00:00:00", "2000-01-15 00:00:00"]

parent_case.configure_forcings(date_range=PARENT_DATES, function_name = "get_glorys_data_from_rda")

# dataset_path is CESM_MOM_OUTPUT's one non-required arg with no universal default (it's
# meant to point at any user-supplied directory), so each child passes its own parent's
# nesting-output directory in via function_overrides -- see the note above for what that
# directory should contain. function_overrides is written straight to config.json, so its
# values must be JSON-serializable (plain str), not a Path.
mid_case.configure_forcings(
    date_range         = MID_DATES,
    product_name       = "cesm_mom_output",
    function_name      = "get_mom6_output_data",
    function_overrides = {"dataset_path": PARENT_NESTING_OUTPUT.as_posix()},
)

tiny_case.configure_forcings(
    date_range         = TINY_DATES,
    product_name       = "cesm_mom_output",
    function_name      = "get_mom6_output_data",
    function_overrides = {"dataset_path": MID_NESTING_OUTPUT.as_posix()},
)

parent_case.process_forcings()

Step 7: Visualize the nested current fields¶

All three domains actually ran to completion (see ~/scratch/archive/gibraltar_*), so we can pull real surface-current output (speed) and land-mask (geolon/geolat/wet) from each domain’s mom6.h.sfc and mom6.h.static history files and paint all three directly on top of one another, in the same map, at each domain’s true geographic location, parent on the bottom, middle child on top of it, tiny child on top of that, with a thin white outline marking where each child’s footprint actually starts. One shared color scale across all three, and the displayed window is cropped in around the strait so the nested children aren’t tiny slivers in the full parent extent. We also animate it over the run so you can watch the Gibraltar jet meander day by day.

from pathlib import Path
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.patches as mpatches
from matplotlib.colors import Normalize

ARCHIVE = Path("~/scratch/archive").expanduser()
DOMAINS = ["gibraltar_parent", "gibraltar_mid_child", "gibraltar_tiny_child"]
LABELS = {"gibraltar_parent": "1/12°", "gibraltar_mid_child": "1/48°", "gibraltar_tiny_child": "1/192°"}
TRIM = 2       # crop the outermost OBC nudging-zone cells (they host single-cell corner artifacts)
N_FRAMES = 12  # common to all three domains (tiny child is the shortest run)

# crop the displayed window so the nested children aren't tiny slivers in the full parent extent
LON_CROP = (-7.3, -3.3)
LAT_CROP = (34.9, 37.0)

sfc = {d: xr.open_dataset(ARCHIVE / d / "ocn/hist" / f"{d}.mom6.h.sfc.2000-01.nc", decode_timedelta=False) for d in DOMAINS}
static = {d: xr.open_dataset(ARCHIVE / d / "ocn/hist" / f"{d}.mom6.h.static.nc") for d in DOMAINS}
sl2 = (slice(TRIM, -TRIM), slice(TRIM, -TRIM))
geo = {d: dict(lon=static[d].geolon.values[sl2], lat=static[d].geolat.values[sl2],
               wet=static[d].wet.values[sl2]) for d in DOMAINS}

all_speed = {d: [] for d in DOMAINS}
for d in DOMAINS:
    for t in range(N_FRAMES):
        speed = sfc[d].isel(time=t).speed.values[sl2]
        all_speed[d].append(np.where(geo[d]["wet"] > 0, speed, np.nan))

# one shared color scale across all three nested domains
pooled = np.concatenate([np.stack(all_speed[d]).ravel() for d in DOMAINS])
norm = Normalize(vmin=0, vmax=np.nanpercentile(pooled, 99))
cmap = plt.get_cmap("RdBu_r")
dates = sfc["gibraltar_parent"].time.values[:N_FRAMES]

def bounds(d):
    g = geo[d]
    return np.nanmin(g["lon"]), np.nanmax(g["lon"]), np.nanmin(g["lat"]), np.nanmax(g["lat"])

fig, ax = plt.subplots(figsize=(10, 6.5))

# paint each domain directly on top of the last, at its true lon/lat location -- no zoom-box
# callouts, just higher resolution layered over the lower-resolution domain beneath it, with a
# thin white outline marking where each child's footprint actually starts
pms = {}
for i, d in enumerate(DOMAINS):
    g = geo[d]
    z = i * 2
    land = np.where(g["wet"] > 0, np.nan, 1.0)
    ax.pcolormesh(g["lon"], g["lat"], land, cmap="Greys", vmin=0, vmax=1.4, shading="auto", zorder=z)
    pms[d] = ax.pcolormesh(g["lon"], g["lat"], all_speed[d][0], cmap=cmap, norm=norm, shading="auto", zorder=z + 1)

for d in ["gibraltar_mid_child", "gibraltar_tiny_child"]:
    lon0, lon1, lat0, lat1 = bounds(d)
    ax.add_patch(mpatches.Rectangle((lon0, lat0), lon1 - lon0, lat1 - lat0,
                                     fill=False, edgecolor="white", linewidth=1.1, zorder=10))

ax.text(np.nanmean(geo["gibraltar_mid_child"]["lon"]), np.nanmax(geo["gibraltar_mid_child"]["lat"]) + 0.02,
        LABELS["gibraltar_mid_child"], fontsize=9, fontweight="bold", ha="center", va="bottom", zorder=11)
ax.text(np.nanmean(geo["gibraltar_tiny_child"]["lon"]), np.nanmin(geo["gibraltar_tiny_child"]["lat"]) - 0.015,
        LABELS["gibraltar_tiny_child"], fontsize=9, fontweight="bold", ha="center", va="top", zorder=11)
ax.text(LON_CROP[0] + 0.15, LAT_CROP[1] - 0.1, "1/12°", fontsize=9, fontweight="bold", ha="left", va="top")

ax.set_xlim(*LON_CROP)
ax.set_ylim(*LAT_CROP)
ax.set_aspect(1 / np.cos(np.deg2rad(np.mean(LAT_CROP))))
ax.set_xlabel("lon"); ax.set_ylabel("lat")

cb = fig.colorbar(pms["gibraltar_parent"], ax=ax, fraction=0.035, pad=0.02)
cb.set_label("surface current speed (m/s)")
title = ax.set_title("Strait of Gibraltar nested domains")
fig.tight_layout()

out_dir = Path.cwd()
fig.savefig(out_dir / "nesting_hierarchy_viz.png", dpi=150)
plt.show()

# animate over the run so the Gibraltar jet's day-to-day meander is visible
def update(t):
    for d in DOMAINS:
        pms[d].set_array(all_speed[d][t].ravel())
    title.set_text(f"Strait of Gibraltar nested domains, {str(dates[t])[:10]}")
    return list(pms.values()) + [title]

anim = animation.FuncAnimation(fig, update, frames=N_FRAMES, blit=False)
anim.save(out_dir / "nesting_hierarchy_movie.gif", writer=animation.PillowWriter(fps=2))
anim.save(out_dir / "nesting_hierarchy_movie.mp4", writer=animation.FFMpegWriter(fps=2, bitrate=3000))
print("saved nesting_hierarchy_viz.png, nesting_hierarchy_movie.gif, nesting_hierarchy_movie.mp4")

Summary¶

This notebook walked through a full 3-level nesting hierarchy over the Strait of Gibraltar, building all three domains together at each step, refining by 4× at each level (the max recommended child:parent ratio) so the tiny child reaches sub-kilometer resolution, and kept every domain small enough to build and run on ~10 cores or fewer:

Domain	Resolution	lon	lat	Grid size	NTASKS_OCN	OBC source
Parent (gibraltar_parent)	1/12° (~8 km)	-8.0→-2.0	34.5→37.0	~72×30	10	GLORYS
Middle child	1/48° (~2 km)	-6.0→-5.0	35.75→36.25	~48×24	8	Parent nesting output
Tiny child	1/192° (~0.5 km)	-5.9→-5.65	35.83→35.955	~48×24	6	Middle-child nesting output

Each level runs independently once its case is built. Once the parent finishes, the middle child can start;
once the middle child finishes, the tiny child can start.

See CESM_MOM_OUTPUT (CrocoDash/raw_data_access/datasets/cesm_ocean_output.py) for the native-MOM6-output readers used above, get_mom6_output_data (multi-variable-per-file history/diag_table output) and get_mom6_single_variable_data (CESM tseries-style output), and their docstrings for direct-call usage. CESM_MOM_OUTPUT is a separate product from CESM_POP_OUTPUT specifically because they need different variable/ coordinate metadata (native MOM6 names vs. CESM-POP names). Each child’s dataset_path is threaded through per-call via configure_forcings(function_overrides={"dataset_path": ...}) rather than sourced from class-level metadata, since CESM_MOM_OUTPUT has no single universal default (it’s meant to point at whichever parent run produced its OBC source). see the note in Step 6 for why that directory should hold only the 5 OBC-relevant fields (zos, thetao, so, uo, vo) rather than a parent’s full history output.

A larger, Caribbean 3-level nesting demo (loading an existing reference case for the parent instead of building it from scratch) is kept in dev/nesting_caribbean_demo/ for reference.</cell id=“cell-15”>