Numerical precision

Working precision

Yelmo and its Fortran dependencies are compiled in single precision, except elsa. The working precision wp is hard-coded in each package:

Package	File	Setting
yelmo	`src/yelmo_defs.f90`	`wp = sp`
fesm-utils	`src/precision.f90`	`wp = sp`
FastHydrology	`src/fast_hydrology.f90`, `src/bucket.f90`, `src/closures.f90`	`wp = sp`
tracer	`src/tracer_precision.f90`	`wp = sp`, `prec_time = sp`
elsa	`src/elsa_precision.f90`	`wp = dp`

There is no build switch to change this. Interfaces between packages pass wp arrays, so yelmo, fesm-utils, FastHydrology and tracer must agree on wp. elsa works in double precision; Yelmo converts its fields at the elsa interface (src/yelmo_tracers.f90).

Model time

The model time is real(dp) in Yelmo (the public routines, e.g. yelmo_init, yelmo_init_state, yelmo_update, yelmo_restart_write, and the clocks of the components), in the drivers’ clock (fesm-utils tstep_class, ts%time, and timeout), in FastHydrology (hydro_update), in elsa (all entry points, also those for sp fields) and in tracer (prec_time = dp: the time, the deposition times and its archive). In single precision the spacing of floating-point numbers at |time| = 150 kyr is 0.016 yr, so steps would not be exact.

Some solvers work internally in double precision:

Lis (src/physics/solver_linear.F90) solves in double precision; the matrix is assembled in wp and converted.
FastHydrology’s K24 transport copies its wp inputs to dp, solves, and converts the result back (hydro_update, TRANSPORT_K24 branch).

Symmetry regression check

The EISMINT-moving symmetry check (symmetry_check in the &ctrl group, src/yelmo_symmetry.f90) compares the solution with its mirror images, normalised by the maximum ice thickness, against the tolerance symmetry_tol. EISMINT-moving and EXPA check the final state at 1e-3. EXPF’s symmetric state is thermomechanically unstable and shows transient symmetry-breaking bursts with either thermal solver, so it checks the mean over the last 20 kyr (symmetry_avg_time) at 2e-2; persistent asymmetry (“spokes”, ~0.1) still fails.

With the SIA solver the check passes in single precision. With ydyn.solver = "diva" it fails (Linf/Hmax ≈ 1.7e-3), although every DIVA operation is symmetric: with exactly symmetric inputs the solver output is symmetric to 4e-7 (sp) and 2e-13 (dp). The cause is dynamical. Thin grounded frontal cells switch on and off every 100-300 yr (grounded cells have f_ice = 1 whenever H_ice > 0), and the eight mirror-image cells slowly drift out of phase with an e-folding time of ~1.7 kyr. Single-precision round-off (~1e-7) grows to a phase slip after ~18 kyr. In double precision the asymmetry at 25 kyr is ~4e-12.

The symmetry check is therefore meaningful for DIVA only in double precision.

Building in double precision (for regression checks)

Build a separate copy; do not modify the shared checkouts.

Copy fesm-utils, FastHydrology, tracer and yelmo into a scratch tree and link them to each other there (fesm-utils, FastHydrology, tracer symlinks inside each package, as in a normal setup). elsa can be linked as is.
Set wp = dp in yelmo, fesm-utils, FastHydrology and tracer (table above), and also prec_time = dp in tracer.
Build fesm-utils (make fesmutils-static openmp=0 in fesm-utils), then make yelmo-static benchmarks openmp=0 in yelmo, which also builds FastHydrology, elsa and tracer. Start from a tree without object or module files. No extra compiler flags are needed; -fdefault-real-8 makes sp = kind(1.0) equal to dp and breaks the sp/dp generic interfaces (e.g. is_finite in yelmo_tools.f90).

Run the check with, for example:

runme -r -e benchmarks -o output/sym-diva-dp -n par/yelmo_EISMINT_moving.nml -p ctrl.time_end=25e3 ydyn.solver=diva