Implementation and attribution
FastEarth3D is a modern-Fortran implementation of published GIA physics, built on a small set of well-established numerical libraries. This page records what the model is built on, and what informed its development and validation.
Physics and references
FastEarth3D reproduces the physics of VILMA: an incompressible, self-gravitating, Maxwell-viscoelastic sphere with rotational feedback and a full migrating-coastline sea-level equation. VILMA’s character is documented in Albrecht et al. (2024) and Bagge et al. (2021). The numerical method follows Martinec (2000), with rotation from Martinec and Hagedoorn (2014) and the sea-level equation from Kendall et al. (2005) and Martinec et al. (2018).
The implementation is modern and minimal: a small, dependency-light Fortran codebase organized as single-purpose modules (fe_*), validated rung-by-rung against the community benchmarks. It is 3D-ready from the start — laterally varying viscosity costs nothing extra in the solver — while validated first in 1D.
Libraries we build on
| Library | Role | Reference |
|---|---|---|
| SHTns | fast spherical-harmonic transforms (forward/inverse, quadrature, point & gradient evaluation) | Schaeffer (2013) |
| FFTW3 | the longitudinal FFTs underlying SHTns | Frigo and Johnson (2005) |
| netCDF | restart and output I/O | — |
| fesm-utils | shared build system + helper library (precision, ncio, namelists) | — |
The FFT and spherical-harmonic transform are central, not incidental. A pseudo-spectral method moves between grid and spectral space many times per step: the SLE forms its ocean-function product on the grid and its load–response convolution in spectral space. Using SHTns and FFTW for that kernel, rather than a hand-rolled transform, keeps the model fast at the resolutions of interest. The transform convention (orthonormal, no Condon–Shortley phase, Gauss grid) is pinned across the whole code. See Spectral reduction.
The linear solver, by contrast, is hand-rolled: a dependency-free, re-entrant, pivoted banded LU (fe_band). This choice avoids an external dependency (LIS, LAPACK) that would complicate linking into an OpenMP host. See Solver.
Codebases and datasets that informed the validation
Several existing codebases provided independent reference solutions and datasets used during validation of FastEarth3D. Each implements the physics in a somewhat different way, which was helpful for understanding the method, and each supplied reference data used in the benchmark pages:
- giapy (S. B. Kachuck, MIT) — source of the vendored M3–L70–V01 Love-number table, the SBK sea-level reference curves, and the analytic load/basin specifications for the Martinec et al. (2018) cases.
- TABOO (D. Melini, GPLv3) — a normal-mode code used as an independent reference to regenerate the M3–L70–V01 Love numbers (NV=3, CODE=7).
- isostasy_data / FastIsostasy.jl (J. Jereczek) — source of the Spada et al. (2011) disc-load spatial-response reference.
Full provenance and licensing of every vendored dataset is in the repository’s data/benchmarks/PROVENANCE.md.
Build and host coupling
The model uses a configme build: a machine fragment inserted into a Makefile template. Dependencies are provided through a fesm-utils checkout (FFTW, SHTns) plus system netCDF, and runs are staged with runme. Both are covered in Install & run.
The host-coupling contract, designed for CLIMBER-X, is intentionally narrow. The host passes ice thickness in and receives relative sea level and bedrock out, on its own grid. The model owns its Gauss grid and the host↔︎Gauss remap (conservative in, bilinear out), built from the parameter record, and keeps all spherical-harmonic transforms inside.
License and citation
See the repository for the current license. When citing the method, cite the primary sources — above all Martinec (2000) for the spectral–finite-element approach, and Spada et al. (2011) and Martinec et al. (2018) for the benchmarks against which FastEarth3D is validated.