Creating an FVCOM nest forcing file from native NEMO output

This tutorial adapts a production workflow for creating FVCOM boundary forcing from native NEMO output. It uses a short test period before scaling to a complete month or multi-year run.

The paths and filenames are deliberately generic. Replace the configuration values with files from your own domain.

The workflow is: load an FVCOM grid, read split NEMO T/U/V/W files, interpolate physical fields, optionally add ERSEM tracers and tidal corrections, and write an FVCOM nest forcing file.

1. Imports and configuration

from datetime import datetime, timedelta
from pathlib import Path
import os
import numpy as np

from pyfvcom2.date_utils import create_datetime_array
from pyfvcom2.exceptions import PyFVCOM2ValueError
from pyfvcom2.grid import create_grid
from pyfvcom2.interpolation import NEMOInterpolator
from pyfvcom2.nemo_reader import NEMOReader, default_fvcom_to_nemo_var_names
from pyfvcom2.nest import NestManager

# This is the shared sample-data layout used by the other cookbook
# notebooks. Change only this section when moving to another domain.
data_dir = Path(os.path.expanduser('~/data/pyfvcom2_doc'))
nemo_dir = data_dir / 'AMM7_NEMO'
fvcom_dir = data_dir / 'FVCOM_tamar_estuary'

fvcom_model = "tamar_v2"
grid_file = fvcom_dir / f'{fvcom_model}_grd.dat'
obc_file = fvcom_dir / f'{fvcom_model}_obc.dat'
sigma_file = fvcom_dir / f'sigma_gen.dat'

# Use explicit lists when filenames are irregular. Glob patterns are
# convenient when there is one file per day or per month.
nemo_files = {
    'T': nemo_dir / 'amm7_1d_19820101_19820131_grid_T.nc',
    'U': nemo_dir / 'amm7_1d_19820101_19820131_grid_U.nc',
    'V': nemo_dir / 'amm7_1d_19820101_19820131_grid_V.nc',
    'W': nemo_dir / 'amm7_1d_19820101_19820131_grid_W.nc',
}
mesh_mask_file = nemo_dir / 'mesh_mask.nc'
mask_file_paths = mesh_mask_file if mesh_mask_file.exists() else None

for grid, file_path in nemo_files.items():
    print(grid, file_path, '(found)' if file_path.exists() else '(missing)')
    if not file_path.exists():
        raise FileNotFoundError(f'Missing NEMO {grid}-grid file: {file_path}')

T /users/modellers/miwa/data/pyfvcom2_doc/AMM7_NEMO/amm7_1d_19820101_19820131_grid_T.nc (found)
U /users/modellers/miwa/data/pyfvcom2_doc/AMM7_NEMO/amm7_1d_19820101_19820131_grid_U.nc (found)
V /users/modellers/miwa/data/pyfvcom2_doc/AMM7_NEMO/amm7_1d_19820101_19820131_grid_V.nc (found)
W /users/modellers/miwa/data/pyfvcom2_doc/AMM7_NEMO/amm7_1d_19820101_19820131_grid_W.nc (found)

2. Load the FVCOM grid

The NEMO interpolator uses geographic target coordinates. The FVCOM grid can still use a projected working coordinate system because longitude and latitude coordinates remain available. Replace the EPSG code with the projected CRS used by your grid.

fvcom_grid = create_grid(
    grid_file,
    mesh_type='fvcom',
    sigma_file=sigma_file,
    coordinate_system='cartesian',
    epsg_code='32630',  # replace with the CRS used by your domain
    obc_filename=obc_file,
)

print('Nodes:', fvcom_grid.n_nodes)
print('Elements:', fvcom_grid.n_elements)
print('Sigma layers:', fvcom_grid.n_sigma_layers)
print('Open boundaries:', len(fvcom_grid.open_boundaries))

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[3], line 1
----> 1 fvcom_grid = create_grid(
      2     grid_file,
      3     mesh_type='fvcom',
      4     sigma_file=sigma_file,

File /data/proteus4/scratch/miwa/mamba_installation/miniforge3/envs/pyfvcom2/lib/python3.14/site-packages/pyfvcom2/grid.py:713, in create_grid(grid_file, mesh_type, sigma_file, coordinate_system, epsg_code, **kwargs)
    710 mesh_data = read_mesh_file(grid_file, mesh_type=mesh_type, **kwargs)
    712 # Read sigma data
--> 713 sigma_config = read_sigma_file(sigma_file)
    714 sigma_data = process_sigma_config(sigma_config, mesh_data.x3)
    716 return Grid(mesh_data, sigma_data, coordinate_system, epsg_code)

File /data/proteus4/scratch/miwa/mamba_installation/miniforge3/envs/pyfvcom2/lib/python3.14/site-packages/pyfvcom2/sigma.py:87, in read_sigma_file(filepath)
     84 zku = None
     85 zkl = None
---> 87 with open(filepath, "r") as f:
     88     lines = f.readlines()
     89     for line in lines:

FileNotFoundError: [Errno 2] No such file or directory: '/users/modellers/miwa/data/pyfvcom2_doc/FVCOM_tamar_estuary/tamar_v2_sigma_gen.dat'

3. Choose a short test period

Start with dates definitely present in the NEMO files. Three daily records are enough to validate file discovery, masks, interpolation, and output dimensions.

start_date = datetime(1982, 1, 1, 12)  # sample NEMO data
end_date = datetime(1982, 1, 3, 12)
date_times = create_datetime_array(start_date, end_date, timedelta(days=1))
print(date_times)

4. Read native NEMO output

The defaults cover common physical variables: votemper, vosaline, uo, vo, and optionally sossheig. The reader keeps logical grids separate. U and V are interpolated on their native staggered grids and rotated to east/north components by NEMOInterpolator.

nemo_reader = NEMOReader(
    nemo_files,
    mask_file_paths=mask_file_paths,
)

print('Loaded grids:', nemo_reader.grid_names)
for grid in nemo_reader.grid_names:
    print(grid, nemo_reader.time_span(grid))

for name in ('votemper', 'vosaline', 'uo', 'vo'):
    grid = nemo_reader.grid_for_variable(name)
    print(name, 'on', grid, 'with', nemo_reader.get_var_ndims(name, grid))

for date in date_times:
    for grid in ('T', 'U', 'V'):
        if not nemo_reader.contains_date(date, grid):
            raise ValueError(f'{date} is outside the {grid}-grid time range')

Optional ERSEM files

ERSEM tracer output can be registered as an additional logical grid. This sample layout provides a single ptrc_T file. Adjust the filename pattern and coordinate names if your files use separate org_T and inorg_T files.

ersem_files = {
    'PTRC': sorted(nemo_dir.glob('*ptrc_T.nc')),
}
ersem_files = {grid: files for grid, files in ersem_files.items() if files}
ersem_grid_var_names = {
    grid: {'time': 'time_counter', 'lon': 'nav_lon',
           'lat': 'nav_lat', 'depth': 'deptht'}
    for grid in ersem_files
}

# If ERSEM files are present, construct the reader with this combined map.
# Otherwise retain the physical-only reader created above.
if ersem_files:
    nemo_reader.close()
    nemo_reader = NEMOReader(
        {**nemo_files, **ersem_files},
        mask_file_paths=mask_file_paths,
        grid_var_names=ersem_grid_var_names,
    )

5. Create the nest and add physical forcing

num_grid_bands controls how far inward from the open boundary the nest extends. Two bands is a useful starting point, but use the value required by the child FVCOM configuration.

nest_manager = NestManager(
    fvcom_grid,
    num_grid_bands=2,
    weights_calculation_method='linear',
)
nest_manager.set_dates(date_times)

fvcom_to_nemo = dict(default_fvcom_to_nemo_var_names)
nemo_interpolator = NEMOInterpolator(
    nemo_reader,
    fvcom_to_nemo_var_names=fvcom_to_nemo,
)
available_ersem = []

# Add sea-surface height when sossheig, or a mapped alternative, exists.
try:
    nemo_reader.grid_for_variable(fvcom_to_nemo['zeta'])
    nest_manager.add_forcing_data(nemo_interpolator, 'zeta', 'node')
except PyFVCOM2ValueError:
    print('No sea-surface-height field found')

class ConstantZetaInterpolator:
    """Zero-elevation fallback for a test file without sossheig."""
    def interpolate(self, coordinates, fvcom_var_name):
        if fvcom_var_name != 'zeta':
            raise ValueError('This fallback only supports zeta')
        n_dates = len(np.atleast_1d(coordinates.dates))
        return np.zeros((n_dates, len(coordinates.x1)), dtype=np.float32)

try:
    has_zeta = 'zeta' in nest_manager.get_forcing_data()
except PyFVCOM2ValueError:
    has_zeta = False
if not has_zeta:
    nest_manager.add_forcing_data(ConstantZetaInterpolator(), 'zeta', 'node')

for variable, position in (
    ('u', 'element'), ('v', 'element'),
    ('temp', 'node'), ('salinity', 'node'),
):
    nest_manager.add_forcing_data(nemo_interpolator, variable, position)
    print(variable, nest_manager.get_forcing_data(variable).shape)

6. Add optional ERSEM tracers

When an ERSEM variable has the same name in NEMO and the FVCOM forcing configuration, use an identity mapping. Start with a few tracers and expand the list only after the physical fields work.

ersem_variables = ['N1_p', 'N3_n', 'O2_o', 'P1_Chl']
for variable in ersem_variables:
    try:
        nemo_reader.grid_for_variable(variable)
    except PyFVCOM2ValueError:
        print('Skipping unavailable tracer:', variable)
        continue
    fvcom_to_nemo[variable] = variable
    available_ersem.append(variable)

nemo_interpolator = NEMOInterpolator(
    nemo_reader, fvcom_to_nemo_var_names=fvcom_to_nemo
)
for variable in available_ersem:
    nest_manager.add_forcing_data(nemo_interpolator, variable, 'node')
print('ERSEM variables added:', available_ersem)

7. Optional tidal adjustment

For a tide-adjusted nest, create a TideManager from your own harmonic atlas files, register interpolators for zeta, u, and v, and call nest_manager.add_tidal_data(tide_manager). Then use:

adjust_tides = ['zeta', 'u', 'v', 'ua', 'va']

Keep this disabled until the NEMO-only forcing has been checked. Do not apply tidal corrections if the NEMO fields already include the same tidal signal.

USE_TIDES = False
# Set USE_TIDES=True only after creating and registering a TideManager.
adjust_tides = ['zeta', 'u', 'v', 'ua', 'va'] if USE_TIDES else None

8. Write and validate the forcing file

The depth-averaged ua and va fields are generated automatically when u and v are added. Extra ERSEM variables are written as node-based 3-D fields.

# create_forcing_file always requires zeta. This makes the short sample
# runnable even when the NEMO files do not contain sossheig. Replace this
# fallback with a real elevation source for production forcing.
if 'zeta' not in nest_manager.get_forcing_data():
    nest_manager.add_forcing_data(
        ConstantZetaInterpolator(), 'zeta', 'node'
    )

required = ('zeta', 'u', 'ua', 'v', 'va', 'temp', 'salinity')
missing = [name for name in required if name not in nest_manager.get_forcing_data()]
if missing:
    raise RuntimeError(f'Missing forcing variables: {missing}')

# NestManager writes the basename in the current working directory.
requested_output = data_dir / 'nest_forcing_short_test.nc'
nest_manager.create_forcing_file(
    requested_output,
    nest_type=3,
    adjust_tides=adjust_tides,
    extra_node_vars=available_ersem or None,
)
created_output = Path(requested_output).name
print('Created:', Path.cwd() / created_output)
print('Size (MB):', Path(created_output).stat().st_size / 1024**2)
nemo_reader.close()

9. Visualise the generated forcing

The following cell plots the first-time-step sea-surface elevation, surface temperature, and surface current speed over the FVCOM nest points. It uses the actual filename created by NestManager.

from netCDF4 import Dataset
import matplotlib.pyplot as plt

nest_nodes = nest_manager.get_all_nest_nodes()
nest_elements = nest_manager.get_all_nest_elements()

with Dataset(created_output) as nest_ds:
    zeta = np.asarray(nest_ds.variables['zeta'][0, :])
    temp = np.asarray(nest_ds.variables['temp'][0, 0, :])
    u_surface = np.asarray(nest_ds.variables['u'][0, 0, :])
    v_surface = np.asarray(nest_ds.variables['v'][0, 0, :])
    speed = np.hypot(u_surface, v_surface)

fig, axes = plt.subplots(1, 3, figsize=(17, 5), constrained_layout=True)

plots = [
    (axes[0], fvcom_grid.lon_nodes[nest_nodes], fvcom_grid.lat_nodes[nest_nodes], zeta, 'Sea-surface elevation (m)', 'viridis'),
    (axes[1], fvcom_grid.lon_nodes[nest_nodes], fvcom_grid.lat_nodes[nest_nodes], temp, 'Surface temperature (degC)', 'plasma'),
    (axes[2], fvcom_grid.lon_elements[nest_elements], fvcom_grid.lat_elements[nest_elements], speed, 'Surface current speed (m/s)', 'magma'),
]
for ax, lon, lat, values, title, cmap in plots:
    scatter = ax.scatter(lon, lat, c=values, s=10, cmap=cmap)
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    ax.set_title(title)
    ax.set_aspect('equal', adjustable='box')
    fig.colorbar(scatter, ax=ax)

plt.show()

A useful first check is that the plotted points cover the intended open-boundary bands, that temperature and elevation have plausible ranges, and that current speed does not contain unexpected extreme values.

Before scaling up, check that NEMO time ranges overlap the requested period, temperature and salinity are finite in the wet nest, U/V have sensible signs and magnitudes after rotation, and output dimensions match the child FVCOM setup. Also check ERSEM units and fill values, and make sure tidal corrections are not applied twice.

Once the short test is sound, expand date_times and provide the corresponding NEMO file lists. For long runs, process one output period at a time and close each NEMOReader after writing its file.