FESOM_F

FESOM_F is a clean-architecture Fortran reimplementation of the FESOM2 ocean + sea-ice model. It began as a line-for-line transcription of FESOM2's numerical kernels and reproduced a reference FESOM2 configuration bit-for-bit up to milestone M10; the surrounding "plumbing" (data types, the parallel layer, configuration, I/O) is redesigned for clarity and future flexibility. The guiding rule was "first reproduce exactly, then optimize." Since v1 the model has moved beyond FESOM2 on purpose (see below).

Main achievements beyond FESOM2:

• Bottom at vertices. The vertex column defines the bottom and each element takes the depth of its shallowest vertex, so no velocity touches topography; tracer cells are straight prisms with the same area at every level.
• Vector-invariant momentum advection (FESOM3_MOMADV_OPT=1) with an upwind-blended face vorticity, together with a Courant-dependent explicit/implicit split of the vertical velocity (FESOM3_WSPLIT).
• Quad and hybrid (triangle + quad) meshes. One set of 4-slot kernels runs triangle, quad and mixed meshes, from mesh generation and partitioning (fesom_ini) to output. FESOM_F runs on the NEMO ORCA025 quad mesh, and triangle meshes give the same results as before at roundoff level.

Status: v1 is feature-complete. Milestones M0–M10 are done and tagged (m0, m1, m2-mvp, m3 … m10): dynamical core, sea-ice EVP/mEVP, GM/Redi, KPP, zstar ALE, TKE, multi-year production runs, Zarr output, and restart/checkpoint. Through M10 the model was byte-exact against FESOM2.

Changes since v1 (September–October 2026)

The results now differ from FESOM2 on purpose. Bit-identity with the FESOM2 oracle is retired; conservation, mesh-delta, regression and self-consistency gates take its place (TESTING.md §6a).

• Bottom at vertices.
  • nlvls.out is authoritative: the vertex column defines the bottom, and each element takes the depth of its shallowest vertex (3 or 4), so velocities touching topography vanish by construction.
  • The scalar-cell area is 1-D and the same at every level, and the z* stretch covers the full column.
  • New conservation and no-leakage gates (fesom_conserve), including the Redi path.
  • See IMPLEMENTATION.md §Bottom at vertices.
• Time step: FESOM3_STEP_PER_DAY sets the number of steps per day.
• Forcing:
  • The time bracket is reset at a year rollover, which fixes the leap-year crash.
  • Three fixes synced from FESOM2 2.7.4–2.8.0.
• Vertical mixing:
  • Cubic-spline N² and shear at the interfaces (FESOM3_N2_SPLINES, FESOM3_SHEAR_SPLINES), and a ramped convective adjustment.
  • N2smth_v is removed.
• Momentum advection: a vector-invariant form (FESOM3_MOMADV_OPT=1). The face vorticity is upwind-blended (FESOM3_RVO_UPWIND), and implicit w_i advection is in advective form.
• w split (FESOM3_WSPLIT): a smooth, Courant-dependent explicit/implicit split of the vertical velocity, with implicit FCT tracer advection.
• MUSCL: up/downwind gradients are computed on the fly (reference RK3 rule), and the stored per-edge gradients are gone. The CORE2 MFCT step is about 9 % faster.
• SSH solver:
  • Jacobi preconditioner and fewer halo exchanges.
  • Tolerance and iteration limit can be set with FESOM3_SSH_SOLTOL and FESOM3_SSH_MAXITER.
  • Opt-in block-Jacobi incomplete Cholesky preconditioner (FESOM3_SSH_PRECOND=ic0): IC(0) of each rank's owned block, applied without communication. On NEMO025 (1280 ranks, 3-day runs from 1964) it needs 2.5–2.7 times fewer CG iterations and about half the solve time: AB2 at 2400 s 126 → 50 iterations (14.7 → 7.5 ms per step), RK3 at 3600 s 316 → 118 (35.5 → 16.5 ms). The iterates depend on the partition, so runs at different rank counts differ at the solver tolerance; at one rank count runs and restarts stay bitwise reproducible. FESOM2's MITgcm preconditioner did not converge on NEMO025 at 2400 s.
• Hybrid meshes (triangles and quads):
  • Python mesh tools (tools/mesh/): NEMO quad mesh, checks, splitting, cut-outs.
  • fesom_ini, a port of FESOM2's mesh_part with METIS 5.1.0 that handles mixed meshes.
  • Unified 4-slot kernels with kite weights for ocean dynamics, tracers, mixing and GM, forcing and ice coupling, EVP/mEVP and the sea-ice FCT.
  • Triangle-only meshes give the same results at roundoff level. See §5.
• Shared-element velocity sync:
  • A border element belongs to every rank that owns one of its vertices. Each step, update_vel now sends one designated copy to the other owners (FESOM3_SHARED_SYNC=0 turns this off).
  • This fixes the volume leak of the earlier multi-rank runs.
  • A run resumed from a restart now matches a straight run bit for bit at np > 1.
• Zarr I/O: the chunk writer and reader no longer copy through transfer(). This fixes the end-of-run restart crash on large meshes.
• Split-explicit barotropic solver (FESOM3_SSH_SE=1, branch split-explicit):
  • Shchepetkin–McWilliams AB3–AM4 subcycling of the 2-D system after Sergey Danilov's notes, in place of the CG solve; linfs and zstar, any rank count, exact restarts. Off by default, and with it off the results are byte-identical to the build before it.
  • The barotropic substeps carry their own bottom drag (on by default, since 2026-10-18, after Sergey Danilov's answers). The Coriolis term stays out of the barotropic forcing (removed with stored transport sums, as Banerjee et al. 2024 omit it); the first-substep form tried on 2026-10-18 blew up and was reverted.
  • The dissipation parameter χ is mesh-dependent: use the smallest value that is stable on the mesh. The default is 0.1, the smallest value with a clean zstar year on CORE2 with the drag. χ = 0.05 passes 3000-step screens (linfs and zstar) and a linfs year, but its zstar year blew up on day 229 (a grid-scale mode on one 30 m Kara-shelf node), and 0.07 is marginal there. Without the drag χ = 0.1 blows up a cold-started CORE2 run after about 7 days (a shelf mode fed through the sea ice).
  • An SE run can resume a checkpoint written with the CG solve (AB2 or RK3), e.g. to branch an SE run from a CG run. The barotropic rings then start cold from the restored state, as in a cold start; under AB2 the momentum also takes one cold (Euler) first step.
  • See IMPLEMENTATION.md §14.
• RK3 time stepping (FESOM3_RK3=1, branch rk3, opt-in):
  • The three-stage Wicker–Skamarock Runge–Kutta step of Sergey Danilov's quad FV prototype in place of AB2. Stages 1–2 are cheap (tracer advection; momentum advection, Coriolis, pressure gradient), stage 3 does the rest (viscosity, mixing, diffusion, GM/Redi, surface fluxes, the FCT limiter, the w split).
  • Both barotropic solvers: the semi-implicit CG solve (default) or the split-explicit substeps (FESOM3_SSH_SE=1); linfs and zstar; any rank count; exact restarts. Off by default, and with it off the results are byte-identical to the build before it.
  • Stages 1–2 advect with the full w of the stage, so a uniform tracer stays uniform; the w split is made in stage 3 only. Momentum and non-FCT tracers are then stable up to a vertical Courant number of about 9 at FESOM3_WSPLIT_MAXCFL=1 (the RK3_CHECK line reports it as cflz=).
  • With the CG solver the stage-1/2 trims go to the depth-mean velocity extrapolated in time (AB2's weights for the stage-3 barotropic Coriolis; rk3.tex's "AB-way"). The trim to Ūⁿ of rk3.tex made the barotropic Coriolis forward Euler and grew a spurious Weddell cyclone on NEMO025 (−3.3 m ssh vs −1.9 m AB2; 2400 s runs blew up in 10–13 days). The checkpoint carries Ūⁿ⁻¹ (rk3_ubold_u/_v, optional on read). Plan Addendum A.
  • Meant to run at about twice the AB2 time step; the longer-Δt screens and the timing on compute nodes are open (see the env block in §4). See IMPLEMENTATION.md §15.
• Conserving sea-ice cut_off (FESOM3_ICE_CUTOFF_CONSERVE=1, opt-in). FESOM2's cut_off deletes the ice and snow where a_ice or m_ice falls below 1e-9 without any flux. The ice FCT limits a_ice and m_ice separately and leaves a_ice ~ 0 with m_ice 1e-5 to 2e-4 m at open-water nodes ahead of a moving ice edge. On CORE2 this deletes about 65 km³ of ice a year, which under zstar is -2.3e11 psu m³ of ocean + ice salt a year (-4.8e-9 of the content). With the switch, cut_off records the deleted ice and snow per node and the thermodynamics of the same step hands them to the ocean as melt: their fresh water, the ice salt (zstar: real_salt_flux; linfs: the virtual-salt form of therm_ice) and negative growth in thdgr/thdgrsn, which the zstar freshwater balance uses. Ocean + ice salt then closes to roundoff under zstar. The latent heat is not taken from the ocean (FESOM2's cut_off never accounted heat): about 3.5e-6 K a year of ocean warming on CORE2. Off (the default) the model is bit-identical to before. Gate: tools/run_cutoff_conserve_pi.sh (TESTING.md §6d).
• Exact zstar restart (opt-in): with FESOM3_RESTART_ZSTAR_EXACT=1 a zstar checkpoint also holds the SSH stiffness values and dhe, and a resumed run matches the straight run bit for bit. Without it a zstar resume continues from the cold stiffness matrix and geometry (not bitwise).
• Redi diffusivity in the tapered zones scales with resolution. Where the slope taper acts (mixed layer, N² ≤ 0, isopycnal slopes above ODM95_Scr), the Redi tensor turns horizontal and its diffusivity blends to Redi_Kmin. FESOM2 keeps that at 100 m²/s, a 1° value, while K_GM and the interior Redi diffusivity scale with √(2·area)/100 km. With scaling_resolution (set in every driver's GM/Redi setup) it now follows the same factor, clamped at K_GM_min: on NEMO025 about 30 m²/s in the Gulf Stream beside an interior 300 m²/s; unchanged on cells with √(2·area) ≥ 100 km. Every run with FESOM3_REDI changes (CORE2 in its refined regions). Test: test_redi_kmin.

Production runs on Levante with these changes conserve volume to roundoff: CORE2 on 512 ranks (about 95 model years per day) and the NEMO ORCA025 quad mesh on 1280 ranks (about 29 model years per day).

See IMPLEMENTATION.md for the architecture and design decisions, and TESTING.md for how it is validated against FESOM2 and how it scales.

A note on the name. The project is called FESOM_F in this documentation, but it lives in a repository still named fesom3, and you will see that working name in real identifiers: the environment-variable prefix FESOM3_, build directories like build_intel_dp, driver binaries fesom_*, and git tags m0…m10. Those are the actual names — use them verbatim. Only the prose uses "FESOM_F".

---

1. Requirements

FESOM_F currently builds and runs on the Levante supercomputer (DKRZ). Everything is provided through environment modules, loaded automatically by env.sh:

	Anchor (bit-identity)	Portability
Compiler	intel-oneapi-compilers/2022.0.1	gcc/11.2.0
MPI	openmpi/4.1.2-intel-2021.5.0	openmpi/4.1.2-gcc-11.2.0
netCDF	netcdf-c/4.8.1 + netcdf-fortran/4.5.3	same

• netCDF-Fortran is required (forcing and initial-condition reads); the build aborts if nf-config is not on the PATH.
• liblz4 is optional — it enables the LZ4 output compressor; without it the writer is uncompressed-only.
• A C compiler is required for the mesh partitioner fesom_ini (mesh_part/, with METIS 5.1.0 vendored in lib/metis-5.1.0/); env.sh sets CC=mpicc CXX=mpicxx and configure.sh passes them to CMake. Builds that bypass configure.sh must configure a C compiler themselves.
• The Python mesh tools (tools/mesh/*.py, tools/make_dist1.py) need numpy/scipy: module load python3 (the system python3 has neither).
• The Intel + double-precision build is the anchor: its compiler flags are copied verbatim from FESOM2 v2.7.3, which is what makes bit-identical results possible. The GNU build is for portability testing.

---

2. Build

configure.sh sources env.sh for you (loading the modules and the Levante MPI workaround), runs CMake, and builds:

./configure.sh --compiler intel --precision dp --clean --build   # anchor build

Binaries land in build_intel_dp/bin/ — the model drivers and the mesh partitioner fesom_ini (§5) — and Fortran modules in build_intel_dp/module/. A binary built before the hybrid-mesh work cannot read quad meshes and has no fesom_ini: rebuild before using a NEMO/quad mesh. The portability build is the same with --compiler gnu (→ build_gnu_dp/):

./configure.sh --compiler gnu --precision dp --clean --build

configure.sh options: --compiler intel|gnu, --precision dp|sp|hp, --debug (builds into a separate build_<compiler>_<precision>_debug/ so it never overwrites the anchor), --clean, --build, --jobs N. The build directory name encodes the compiler and precision, so several builds coexist.

Build and run on a Levante node — env.sh only knows how to set up Levante and will stop on any other host.

---

3. Quick check

Run the self-test suite on a login node:

cd build_intel_dp && ctest --output-on-failure

For a stronger check, run the conservation gate — an unforced 20-step pi run that verifies total heat and salt content hold to round-off under zstar, that no velocity survives in a partly-land prism, and that the layer thicknesses stay consistent:

bash tools/run_conserve_pi.sh      # expect: GATE OK
bash tools/run_bottom_audit.sh     # mesh-delta audit on pi + core2

The FESOM2 oracle byte-gates (run_step_gate.sh and friends) are retired — see TESTING.md §6.

What these checks mean — and the full validation story — is in TESTING.md.

---

4. Running a simulation

The production driver is fesom_lifecycle_native_mr (multi-rank, computes its own air–sea fluxes from atmospheric forcing). A run is configured through FESOM3_* environment variables; the example below mirrors tools/run_lifecycle_2yr_freerun_f3_dist864.sbatch:

source env.sh intel                         # modules + Levante MPI workaround

POOL=/pool/data/AWICM/FESOM2
export FESOM3_MESH_DIR="$POOL/MESHES_FESOM2.1/core2"
export FESOM3_IC_FILE="$POOL/INITIAL/phc3.0/phc3.0_winter.nc"
export FESOM3_FORCING_DIR="$POOL/FORCING/JRA55-do-v1.4.0"
export FESOM3_FORCING=JRA55
export FESOM3_RUNOFF_FILE="$POOL/FORCING/JRA55-do-v1.4.0/CORE2_runoff.nc"
export FESOM3_SSS_FILE="$POOL/FORCING/JRA55-do-v1.4.0/PHC2_salx.nc"
export FESOM3_START_CLOCK="0 1 1958"        # "timeofday day year"
export FESOM3_NSTEPS=35040                   # 2 model years at 48 steps/day

# physics (the default FESOM2 production configuration)
export FESOM3_WHICH_ALE=zstar
export FESOM3_FER_GM=1 FESOM3_REDI=1 FESOM3_MIX_TKE=1 FESOM3_SW_PENE=1 FESOM3_CHL_SWEENEY=1

srun -l -n 864 ./build_intel_dp/bin/fesom_lifecycle_native_mr > run.log 2>&1

The physics options added since v1 are off by default (the shared-element sync is on). The current CORE2 and NEMO production runs use FESOM3_WHICH_ALE=linfs and add:

export FESOM3_MOMADV_OPT=1 FESOM3_RVO_UPWIND=1.0     # vector-invariant momentum advection
export FESOM3_WSPLIT=1                               # Courant-dependent explicit/implicit w split
export FESOM3_SHEAR_SPLINES=1 FESOM3_N2_SPLINES=1    # spline shear and N² (set both)
export FESOM3_STEP_PER_DAY=48                        # 1800 s on CORE2; 72 (1200 s) on NEMO 0.25°

On/off switches are presence-based. FESOM3_FER_GM, FESOM3_REDI, FESOM3_MIX_TKE, FESOM3_MIX_KPP, FESOM3_SW_PENE, FESOM3_CHL_SWEENEY, FESOM3_WSPLIT, FESOM3_SHEAR_SPLINES, FESOM3_N2_SPLINES and FESOM3_STEP_DIAG are on whenever the variable is set and not empty, also for =0; leave it unset to turn the option off. FESOM3_MOMADV_OPT, FESOM3_RVO_UPWIND, FESOM3_WSPLIT_MAXCFL/_MINCFL, FESOM3_SHARED_SYNC, FESOM3_SSH_SE, the FESOM3_SE_* knobs, FESOM3_RK3 and the FESOM3_RK3_* knobs take values.

Where the knobs are defined. The authoritative, always-current list of run variables is the environment block at the top of src/drivers/fesom_lifecycle_native_mr.F90 (each is read with get_environment_variable); the SSH solver knobs, FESOM3_SHARED_SYNC, FESOM3_LINFS_CONSERVE and FESOM3_RESTART_ZSTAR_EXACT in the table below are read where they act (the solver, src/infra/mod_shared_elem.F90, src/oce/oce_ale_tracer.F90, src/io/mod_io_restart.F90). The tools/*.sbatch and tools/run_lifecycle_*.sh scripts are working examples. The most common variables:

Variable	Meaning
FESOM3_MESH_DIR	mesh directory (must contain dist_<NP>/)
FESOM3_IC_FILE	initial T/S netCDF
FESOM3_IC_EXTRAP_DET, FESOM3_IC_EXTRAP_TOL	cold start: fill the land-masked points of the initial T/S with FESOM2's deterministic fill (ic_extrap_det, PR #979): the same initial state on every dist_<NP>. Off by default (the historical fill, whose result depends on the partition and the node numbering); the tolerance of its relaxation defaults to 1e-3. Read in src/oce/oce_initial_state.F90
FESOM3_FORCING_DIR, FESOM3_FORCING	atmospheric forcing directory and set (JRA55/CORE2)
FESOM3_RUNOFF_FILE, FESOM3_SSS_FILE	runoff and sea-surface-salinity restoring climatologies
FESOM3_START_CLOCK	cold-start clock, "timeofday day year"
FESOM3_NSTEPS	number of time steps (or use FESOM3_RUN_LENGTH + FESOM3_RUN_UNIT)
FESOM3_STEP_PER_DAY	time steps per day (default 48 = 1800 s; CORE2 uses 48, finer meshes need more)
FESOM3_WHICHEVP	sea-ice solver: 0 = EVP, 1 = mEVP
FESOM3_ICE_CUTOFF_CONSERVE	sea-ice cut_off: 1 = conserving (the ice and snow it deletes melt into the ocean in the same step, with their fresh water and ice salt), 0 or unset = FESOM2's (deleted without a flux); any other value stops. See Changes since v1 above
FESOM3_WHICH_ALE	vertical coordinate: linfs (default) or zstar. zlevel is not ported, and the driver does not reject it: do not use it
FESOM3_LINFS_CONSERVE	S, T or TS (unset = FESOM2): with linfs the layers are fixed while w at the surface is not zero, so the vertical advection moves tracer through the surface and the ocean's salt and heat content drift (CORE2 JRA55-do, Dmitry's setup: whole-ocean S -5.5e-5 in the cold-start year, ~-2e-6 a year later). This returns that flux, summed over the ocean, as a uniform surface flux of S and/or T every step (MITgcm's linFSConserveTr). linfs and the FCT limiter only. Read in src/oce/oce_ale_tracer.F90
FESOM3_FER_GM, FESOM3_REDI	Gent–McWilliams bolus advection, Redi isopycnal diffusion. FESOM3_REDI also switches GM on
FESOM3_MIX_KPP, FESOM3_MIX_TKE, FESOM3_SW_PENE	vertical mixing and shortwave penetration
FESOM3_SHEAR_SPLINES, FESOM3_N2_SPLINES	cubic-spline shear and N² at the interfaces (set both)
FESOM3_MOMADV_OPT	momentum advection: 2 (default) = flux form, 1 = vector invariant
FESOM3_RVO_UPWIND	upwind blend of the face vorticity in [0, 1] (default 0; vector invariant only)
FESOM3_WSPLIT, FESOM3_WSPLIT_MAXCFL, FESOM3_WSPLIT_MINCFL	Courant-dependent explicit/implicit w split; explicit CFL cap (default 1.0) and onset (default half the cap)
FESOM3_OUTPUT, FESOM3_OUTPUT_EVERY	output directory (enables output) and step cadence
FESOM3_RESTART, FESOM3_RESTART_LENGTH, FESOM3_RESTART_UNIT	checkpoint directory and cadence
FESOM3_RESTART_ZSTAR_EXACT	1 = exact zstar resume: the checkpoint adds the stores ssh_stiff (SSH stiffness values) and dhe, and the resume rebuilds the zstar geometry from hnode; set it for the run that writes and the run that reads. 0 or unset = off (default; a zstar resume then differs from the straight run). Ignored with one log line under linfs, whose resume is exact anyway; resuming a checkpoint without the two stores with it on stops. Read in src/io/mod_io_restart.F90
FESOM3_SHARED_SYNC	shared border elements: 1 (default) = every step the other owners take the designated owner's uv (no volume leak, exact np>1 restarts); 0 = off, the old unsynced copies (A/B comparisons only). Read in src/infra/mod_shared_elem.F90
FESOM3_STEP_DIAG	per-step global diagnostic line (max|eta|, max|uv|, ice, Tmax, Smax, and the shared-element numbers divsum=, dUV=, dAB=, dsig=); two extra collectives per step, off by default. Under zstar a second line per step, salt step: the ocean + sea-ice salt content, the step's residual (change minus the balanced SSS restoring; roundoff when salt is conserved) and the ice salt the conserving cut_off returned (one more collective)
FESOM3_SSH_SOLTOL, FESOM3_SSH_MAXITER	SSH solver tolerance and iteration cap — test/diagnostic only, read in src/oce/oce_ssh_solve.F90 (defaults 1e-5, 2000)
FESOM3_SSH_PRECOND	SSH CG preconditioner: jacobi (default) or ic0 (block-Jacobi incomplete Cholesky, src/oce/oce_ssh_ic0.F90); any other value stops
FESOM3_SSH_SE	free surface: 0 (default) = semi-implicit CG solve, 1 = split-explicit barotropic subcycling (production driver and fesom_conserve only; AB2, no cavities). A checkpoint written with the CG solve resumes with cold barotropic rings (one log line). Read in src/oce/oce_ssh_se.F90, like the knobs below
FESOM3_SE_M, FESOM3_SE_M_FORCE	substeps per baroclinic step (default 50); the run stops if M is below the CFL-guard minimum printed at startup (CORE2 at 1800 s: 40 at the default χ 0.1, 38 at χ 0.05, 48 at χ 0.25) unless FESOM3_SE_M_FORCE=1
FESOM3_SE_CHI	dissipation parameter χ in [0, 0.25] (default 0.1). Mesh-dependent: use the smallest value stable on the mesh; the default is that value for CORE2 with the barotropic drag (0.05 blew up a zstar year, 0.07 is marginal; see "Changes since v1")
FESOM3_SE_VISC_AH	barotropic viscosity A_h = value × edge length, the value in m/s (default 0.01; 0 = off)
FESOM3_SE_VISC_FLOW	1 = the flow-dependent barotropic viscosity of the Zenodo SE code in place of A_h (FESOM3_SE_VISC_GAMMA0/_GAMMA1, default 10 / 2750); default 0
FESOM3_SE_BOTTDRAG	barotropic bottom drag in the substeps, r = C_d·|u_bottom| frozen over the step: 1 (default) = on, 0 = off (the 3-D drag stays in the forcing either way)
FESOM3_SE_CHECK	1 = one SE_CHECK invariant line per step (η compatibility, trim, mean sea level, shared-copy spreads; zstar adds SE_CHECK_TAIL); default 0
FESOM3_RK3	time stepping: 0 (default) = AB2, 1 = RK3 (Wicker–Skamarock; production driver and fesom_conserve only; linfs or zstar, double precision, no cavities, α = θ = 1). Combines with FESOM3_SSH_SE. Read in src/oce/oce_rk3.F90, like the knobs below; a value other than 0/1 stops
FESOM3_RK3_LO_NSUB	substeps of the FCT low order in stage 3: 1 (default, the AB2 low order with the w split), 2 or 3 (only for a horizontal low-order Courant number near 1, RK3_CHECK cflh=)
FESOM3_RK3_CHECK, FESOM3_RK3_CHECK_TOL	1 = one RK3_CHECK line per step: the stage trims trim1/2/3 (and trimM with SE), the η compatibility etacomp, maxu, nbad, the horizontal low-order Courant number cflh and max CFL_z cflz; with _TOL > 0 a trim or etacomp above it stops the run (a non-finite value always does). Default 0 / 0

Configuration today vs. intent. The production drivers are configured by these FESOM3_* variables plus in-code defaults. A write-once configuration module (mod_config, reading a single namelist.config) already exists, but the lifecycle drivers do not read it yet; consolidating run configuration into one namelist file is the intended future interface.

RK3 runs. The env block for a CORE2 RK3 run on top of the production physics (the twice-longer time step is the point of RK3; compare with the AB2 run at 48 steps per day first):

export FESOM3_RK3=1                     # RK3 in place of AB2 (unset = AB2)
export FESOM3_STEP_PER_DAY=24           # 3600 s; 48 (1800 s) for the AB2-vs-RK3 twin
# optional
export FESOM3_RK3_CHECK=1               # one RK3_CHECK line per step (trims, etacomp, cflh, cflz)
export FESOM3_SSH_SE=1 FESOM3_SE_M=80   # split-explicit; M must be >= the CFL-guard minimum printed at
                                        # startup, which scales with the step (40 at 1800 s on CORE2)

An RK3 checkpoint carries "time_scheme": "rk3" in checkpoint.json. An AB2 run that resumes it starts AB2 cold (an Euler first step, one log line); an RK3 run resumes any checkpoint.

Restart. Set FESOM3_RESTART=<dir> to write checkpoints (cadence via FESOM3_RESTART_LENGTH / FESOM3_RESTART_UNIT). On the next launch, if a complete checkpoint is found under that directory (a restart.latest pointer to a folder containing checkpoint.json), the run resumes from it byte-exactly; otherwise it cold-starts. Under FESOM3_WHICH_ALE=zstar the resume is byte-exact only with FESOM3_RESTART_ZSTAR_EXACT=1 (table above). See IMPLEMENTATION.md §I/O.

---

5. Using different meshes

A mesh is a directory holding the global mesh files plus one sub-directory of partition files per MPI rank count:

core2/
  nod2d.out  elem2d.out  aux3d.out  nlvls.out  edges.out  edge_tri.out  edgenum.out
  elvls.out                                       # present but NOT read (see below)
  dist_1/  dist_2/  dist_8/  dist_864/   ...      # one per rank count (partition files)

Triangle, quad and mixed meshes. elem2d.out may have 3 or 4 node ids per line (the column count is taken from the first line; a triangle in a 4-column file repeats its 1st vertex). Meshes made from NEMO grids carry a depth.out instead of aux3d.out; the model reads the levels from depth.out when aux3d.out is absent.

elvls.out is no longer read. FESOM_F places the bottom at vertices: nlvls.out (the per-vertex level count) is authoritative and each element's bottom is derived as the shallowest of its 3 or 4 vertices, so an element is wet only where all of its corners are. The file can stay in the mesh directory; nothing opens it. See IMPLEMENTATION.md §Bottom at vertices.

Point FESOM3_MESH_DIR at this directory. The driver auto-selects dist_<NP>/ from the number of MPI ranks (srun -n <NP>); a single-rank run reads the global nod2d.out/elem2d.out directly. So switching meshes — or switching rank counts on the same mesh — is just changing FESOM3_MESH_DIR and -n, provided the matching dist_<NP>/ exists.

These meshes are used routinely:

• pi — a small test mesh (~3k nodes), shipped with the FESOM2 oracle at /home/a/a270088/port2/fesom2/tests/data/MESHES/pi (with dist_1/, dist_2/, dist_8/). This is the default when FESOM3_MESH_DIR is unset, and what the byte-gates and ctest use — it runs in seconds on a login node.
• CORE2 — the global ~1° triangle production mesh (126,858 nodes, 244,659 triangles, 47 levels) at /pool/data/AWICM/FESOM2/MESHES_FESOM2.1/core2, with partitions up to dist_864.
• NEMO 0.25° — the ORCA025 grid as a quad mesh, nemo025q_fesom (915,866 nodes, 902,408 quads, 46 levels), and the mixed test mesh nemo025mix (685,452 quads and 433,912 triangles on the same nodes), at /work/ab0995/a270029/meshes/. Production runs use 1280 ranks and a 1200 s step.

To run on a new mesh, create its dist_<NP>/ partition for your target rank count with the in-tree partitioner build_*/bin/fesom_ini (FESOM2's mesh_part, ported): run it in a work directory whose namelist.config sets MeshPath, n_part=<NP> and, for a NEMO mesh with depth.out, use_depthfile="nemo". It writes the edge and level files and dist_<NP>/ into the mesh directory, including dist_<NP>/format.out, which announces the appended halo lists needed by exchange_edge/exchange_nod_full. Then point FESOM3_MESH_DIR at the mesh directory.

A mesh made from a NEMO grid goes through tools/mesh/nemo_meshq.py → remove_deadends.py → [split_quads.py for a mixed mesh, cut_region.py for a regional cut-out] → check_mesh.py → fesom_ini → check_dist.py. The exact commands, the namelist.config and the pitfalls are in docs/HANDOFF.md, section "Hybrid (mixed triangle + quad) meshes" → "Mixed-mesh setup workflow".

---

6. Output and reading results

Set FESOM3_OUTPUT=<dir> to enable output. Variables and writer options are read from namelist.io in the run directory (see config/namelist.io); without it, a default variable set is written at the FESOM3_OUTPUT_EVERY step cadence. Output is written as Zarr stores — one per variable per period:

<dir>/ssh.fesom.1958.zarr
<dir>/temp.fesom.1958.zarr
<dir>/fesom.mesh.diag.zarr        # mesh geometry / connectivity (UGRID-style)

The stores are partition-independent (the same regardless of rank count) and open directly with xarray:

import xarray as xr
ds = xr.open_zarr("ssh.fesom.1958.zarr")   # CF time axis, embedded lon/lat

Format change (hybrid meshes). fesom.mesh.diag.zarr now holds 4-column per-face arrays (dim n4: face_nodes, face_edges, face_links, gradient_sca_x/y). A triangle's 4th slot is the _FillValue (−1; −999 in face_links; coefficient 0 in gradient_sca_x/y), so xarray returns the integer arrays as float with NaN unless opened with mask_and_scale=False. Element lon/lat are now the cyclic-wrapped centroid (elements across the rotated dateline moved). Details: docs/HANDOFF.md → "Gate migration for hybrid meshes" → "Output format changes".

---

7. Repository layout

src/            Fortran source
  params/         precision, constants, configuration            (6 modules)
  types/          mesh / partition / dynamics / tracers / ice    (5)
  infra/          MPI halo layer, partitioning, clock, timers    (10)
  mesh/           mesh read, geometry, rotation                  (4)
  oce/            ocean physics kernels                          (28)
  ice/            sea-ice dynamics, thermodynamics, coupling     (6)
  forcing/        atmospheric forcing read + bulk fluxes         (3)
  io/             Zarr output, decomposition, restart            (8)
  step/           timestep orchestration                         (2)
  drivers/        program entry points + byte-gate drivers       (33)
mesh_part/      fesom_ini sources (FESOM2's mesh partitioner, mixed-mesh capable)
lib/metis-5.1.0/  vendored METIS used by fesom_ini
cmake/          compiler-flag and helper modules
config/         namelist.io (output configuration)
env/            per-host module/setup files (levante.dkrz.de)
test/           ctest self-tests
tools/          byte-gate scripts, oracle runner, diff tools, sbatch jobs
  mesh/           Python mesh tools (nemo_meshq, remove_deadends, split_quads,
                  cut_region, check_mesh, check_dist) + test_*.sh
docs/           HANDOFF, LESSONS, plans (internal working logs)
configure.sh    build helper          env.sh   environment loader

---

8. Documentation map

• IMPLEMENTATION.md — architecture, design decisions, differences from FESOM2.
• TESTING.md — validation approach, byte-gate methodology, results, scalability.
• docs/HANDOFF.md — the detailed internal development log (milestone-by-milestone state).
• docs/LESSONS.md — accumulated gotchas and subtle reproducibility traps.
• docs/plans/ — per-milestone implementation plans (and completed/).