FVCOM restart from native NEMO AMM7 output

This tutorial creates a restart file from native NEMO AMM7 physical and ERSEM output. It uses the shared sample data and a short test period.

The restart is based on a donor FVCOM restart file, which supplies the dimensions and NetCDF variable definitions. The supplied donor is physical-only, so ERSEM writing is disabled by default. Use an ERSEM-capable donor to write tracers.

1. Imports and paths

from datetime import datetime
from pathlib import Path
import os
import shutil
import numpy as np

from pyfvcom2.fvcom_reader import FVCOMReader
from pyfvcom2.interpolation import NEMOInterpolator
from pyfvcom2.nemo_reader import NEMOReader, default_fvcom_to_nemo_var_names
from pyfvcom2.restart import write_restart
from pyfvcom2.exceptions import PyFVCOM2ValueError

data_dir = Path(os.path.expanduser('~/data/pyfvcom2_doc'))
nemo_dir = data_dir / 'AMM7_NEMO'
fvcom_dir = data_dir / 'FVCOM_tamar_estuary'
donor_filepath = fvcom_dir / 'tamar_v2_donor_restart.nc'
output_dir = data_dir / 'restart_from_nemo'
output_dir.mkdir(parents=True, exist_ok=True)

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',
    'PTRC': nemo_dir / 'amm7_1d_19820101_19820131_ptrc_T.nc',
}

for path in [donor_filepath, *nemo_files.values()]:
    print(('OK' if path.exists() else 'MISSING'), path)
    if not path.exists():
        raise FileNotFoundError(path)
if shutil.which('ncks') is None:
    raise RuntimeError('write_restart requires the NCO command ncks')

2. Read the donor and NEMO data

The PTRC logical grid contains the ERSEM tracers and uses the same coordinate names as the NEMO T grid.

start_date_time = datetime(1982, 1, 1, 12)
include_biogeochemistry = False

fvcom_reader = FVCOMReader(str(donor_filepath))
nemo_reader = NEMOReader(
    nemo_files,
    grid_var_names={'PTRC': {
        'time': 'time_counter', 'lon': 'nav_lon',
        'lat': 'nav_lat', 'depth': 'deptht',
    }},
)
print('NEMO grids:', nemo_reader.grid_names)
for grid in ('T', 'U', 'V', 'PTRC'):
    print(grid, nemo_reader.time_span(grid))
    if not nemo_reader.contains_date(start_date_time, grid):
        raise ValueError(f'{start_date_time} is outside the {grid} range')

3. Select variables and validate the donor

The selected ERSEM variables are examples from ptrc_T. Benthic variables require a separate benthos source and matching variables in the donor restart.

fvcom_to_nemo = dict(default_fvcom_to_nemo_var_names)
physical_variables = ['zeta', 'u', 'v', 'temp', 'salinity']
ersem_3d_variables = ['N1_p', 'N3_n', 'O2_o', 'P1_Chl']
benthic_variables = []
donor_variables = set(fvcom_reader._metadata_dataset.variables)
missing_physical = [v for v in physical_variables if v not in donor_variables]
if missing_physical:
    raise RuntimeError(f'Donor is missing physical variables: {missing_physical}')
if include_biogeochemistry:
    requested = ersem_3d_variables + benthic_variables
    missing = [v for v in requested if v not in donor_variables]
    if missing:
        raise RuntimeError('Use an ERSEM-capable donor; missing: ' + str(missing))
    for variable in requested:
        fvcom_to_nemo[variable] = variable
print('Donor contains physical restart variables')

def source_variable_available(mapping):
    candidates = (mapping,) if isinstance(mapping, str) else mapping
    return any(
        name in dataset.variables
        for name in candidates
        for dataset in nemo_reader._metadata_datasets.values()
    )

4. Interpolate physical and optional ERSEM fields

interpolator = NEMOInterpolator(nemo_reader, fvcom_to_nemo_var_names=fvcom_to_nemo)
node_coords = fvcom_reader.get_interpolation_coordinates('node', 'layer_centre')
node_coords.dates = np.asarray([start_date_time])
element_coords = fvcom_reader.get_interpolation_coordinates('element', 'layer_centre')
element_coords.dates = np.asarray([start_date_time])
positions = {
    'temp': node_coords, 'salinity': node_coords,
    'u': element_coords, 'v': element_coords,
}
if source_variable_available(fvcom_to_nemo['zeta']):
    positions['zeta'] = node_coords
else:
    print('NEMO zeta is unavailable; donor zeta will be retained')
if include_biogeochemistry:
    positions.update({v: node_coords for v in ersem_3d_variables})
    positions.update({v: node_coords for v in benthic_variables})
var_data = {}
for variable, coordinates in positions.items():
    try:
        values = interpolator.interpolate(coordinates, variable)
    except Exception as exc:
        raise RuntimeError(f'Interpolation failed for {variable}: {exc}') from exc
    var_data[variable] = values[0] if values.ndim > 0 and values.shape[0] == 1 else values
    print(variable, values.shape, '->', var_data[variable].shape)

5. Write the restart file

write_restart preserves the donor structure and replaces the selected variables. It requires the NCO ncks command.

restart_file = output_dir / f'nsea_nemo_restart_{start_date_time:%Y%m%d}.nc'
write_restart(str(donor_filepath), str(restart_file), var_data, start_date_time)
print('Created:', restart_file)
print('Size (MB):', restart_file.stat().st_size / 1024**2)

6. Inspect and visualise the restart

from netCDF4 import Dataset
import matplotlib.pyplot as plt

with Dataset(restart_file) as ds:
    temp = np.asarray(ds.variables['temp'][0])
    u = np.asarray(ds.variables['u'][0])
    v = np.asarray(ds.variables['v'][0])
    temp_surface = temp[-1, :] if temp.ndim == 2 else temp
    u_surface = u[-1, :] if u.ndim == 2 else u
    v_surface = v[-1, :] if v.ndim == 2 else v
    speed = np.hypot(u_surface, v_surface)

fig, axes = plt.subplots(1, 2, figsize=(13, 5), constrained_layout=True)
a = axes[0].scatter(fvcom_reader.lon_nodes, fvcom_reader.lat_nodes, c=temp_surface, s=5, cmap='plasma')
axes[0].set_title('Restart surface temperature')
axes[0].set_xlabel('Longitude'); axes[0].set_ylabel('Latitude')
fig.colorbar(a, ax=axes[0], label='degC')
b = axes[1].scatter(fvcom_reader.lon_elements, fvcom_reader.lat_elements, c=speed, s=5, cmap='magma')
axes[1].set_title('Restart surface current speed')
axes[1].set_xlabel('Longitude'); axes[1].set_ylabel('Latitude')
fig.colorbar(b, ax=axes[1], label='m/s')
plt.show()

Extending to a full FVCOM-ERSEM restart

Use an ERSEM-capable donor restart, set include_biogeochemistry = True, and populate benthic_variables with fields available in both the NEMO benthos output and donor. The notebook deliberately stops when a requested variable is missing instead of silently retaining stale donor values.