Skip to main content

Installation and running

This page covers building VILMA2, testing the build, and running a forced simulation. The build is driven by configme (one source of machine and compiler settings) and runs are staged by runme (single simulations and ensembles, locally or on SLURM). The model is built and tested on three machines: a MacBook (gfortran), the AWI cluster albedo (ifx) and the DKRZ cluster levante (ifx). Production runs are currently made on albedo.

Dependencies

Dependency Provides Source
fesm-utils, branch phase0-kryos FFTW, SHTns, the fesmutils helper library (precision, ncio, namelists) and the coords module used by the remapper cloned and built by configme
netCDF (C and Fortran) input, output and restart files system package (Homebrew, environment modules, apt)
configme, runme Makefile generation; run staging fesmc/configme, fesmc/runme (pip install)

The fesm-utils branch is pinned in .configme/manifest.toml, so a checkout records the ref it needs. VILMA2 needs the coords routine grid_init_gaussian, which is in phase0-kryos but not yet in the fesm-utils dev branch: src/vilma_remap.f90 does not compile against dev v1.4. The pin returns to dev once phase0-kryos is merged there.

Install

configme install vilma2 clones fesm-utils into the checkout, builds its serial and OpenMP variants (FFTW, SHTns, LIS and libfesmutils; about 10 minutes), and writes the top-level Makefile from config/Makefile and the machine and compiler fragments. configme show <machine> prints a fragment, configme status reports what is present and what is pending, and configme config vilma2 regenerates the Makefile of an existing checkout without cloning. To reuse an existing fesm-utils checkout instead of cloning one, add --link fesm-utils=/abs/path/to/fesm-utils.

The optimization flags differ by machine. The macbook fragment is the repository’s own (.configme/machines/macbook.mk): gfortran -O3 -mcpu=native -funroll-loops -ffast-math. The clusters use configme’s fragments: ifx -Ofast -march=znver2 -traceback on albedo (AMD EPYC 7702) and -Ofast -march=znver3 -traceback on levante (AMD EPYC 7763). The benchmark numbers in these pages were measured with the albedo build.

MacBook (gfortran)

brew install gcc netcdf netcdf-fortran
pip install git+https://github.com/fesmc/configme.git git+https://github.com/fesmc/runme.git
git clone git@github.com:fesmc/vilma2.git && cd vilma2
configme install vilma2 -m macbook -c gfortran
make -j 8 vilma
make check

AWI albedo (ifx)

The system git has been removed from the login nodes, so load the git module first. The Intel compiler and netCDF modules must be loaded before configme runs, so that it finds nf-config.

module load git intel-oneapi-compilers netcdf-c netcdf-fortran
git clone git@github.com:fesmc/vilma2.git && cd vilma2
configme install vilma2 -m awi_albedo -c ifx
make openmp=1 -j 8 vilma            # on the login node, at most -j 8

Compile on the login node and run the tests and the model through SLURM (see Testing the build). The compute nodes have no git. The build does not need it: the version string compiled into the model comes from git describe and falls back to "unknown" where git is missing.

DKRZ levante (ifx)

module load git intel-oneapi-compilers netcdf-c netcdf-fortran
git clone git@github.com:fesmc/vilma2.git && cd vilma2
configme install vilma2 -m dkrz_levante -c ifx
make openmp=1 -j 8 vilma

The build on levante has been verified (ifx 2023.2, netCDF-Fortran 4.6.2). No runs have been made there yet.

Build

make vilma        # -> bin/vilma.x        (standalone forced-run driver)
make vilma_mkref  # -> bin/vilma_mkref.x  (build a Gauss-grid reference state)
make vilma_remap  # -> bin/vilma_remap.x  (offline lon-lat -> Gauss remap)
make vilma-static # -> libvilma/include/libvilma.a and the .mod files
make bench_visc3d # -> bin/bench_visc3d.x (lateral-viscosity benchmark driver)

make vilma-static is the library that a host (kryos, CLIMBER-X) links: objects, .mod files and libvilma.a are written to libvilma/include/, and the archive is rebuilt only when an object changes. openmp=0 and openmp=1 share that directory. Changing the setting rebuilds every object, and the host must use the matching fesm-utils variant (include-serial or include-omp).

Build switches: debug = 0|1|2 (optimized, checks, profile), openmp = 0|1 (default 1) and vilma_v1 = 0|1 (the optional VILMA1 backend, see doc/vilma-v1-backend.md). The threaded degree loop is what makes production resolutions fast (see Solver). After switching to another fesm-utils checkout, run make clean: objects compiled against the old module files can fail to rebuild.

Testing the build

make check                  # the test suite: about 50 tests, ~2 min
make openmp=1 check-slow    # full-resolution benchmarks and order tests, ~10 min on 32 threads

Both targets run every test, print a summary that names the failing tests, and exit with a non-zero status if any failed. check-slow runs the four sea-level benchmark cases, the low-viscosity-zone benchmark, the off-pole rotational-invariance test at \(\ell_\max = 128\) and three convergence-order tests (strain–memory coupling, sea-level–memory coupling, adaptive step control). It raises the stack limit itself. Setting VILMA_BENCH_DUMP to a directory makes the benchmark tests write their comparison tables there (see Benchmarks).

On albedo, run both suites in a job, with the executables built beforehand on the login node (make openmp=1 -j 8 check builds them and then runs the tests; interrupt it once the build is done, or build the targets listed in TESTS and SLOW in config/Makefile_vilma.mk). A minimal job script:

#!/bin/bash
#SBATCH --account=<project> --partition=smp --qos=12h --time=01:00:00
#SBATCH --cpus-per-task=32 --mem=60G
source ~/.bash_profile
export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK
make openmp=1 check; make openmp=1 check-slow

Known state. With ifx -Ofast on albedo, all 49 tests of make check and all 9 runs of make check-slow pass, with the default stack settings. A GitHub workflow builds the model with gfortran on Ubuntu and runs make check on every push and pull request to main.

Configuration: &vilma and &ctl

A run configuration has two namelist groups:

  • &vilma: the model’s physics and numerics. This is the record a host model fills in memory.
  • &ctl: the standalone run itself: forcing, reference and output files, the time window, i_eq, online remapping and the restart path. A host model never reads it.

The complete, documented &vilma defaults are in input/vilma_defaults.nml and are loaded automatically. A run configuration such as vilma.nml (the default) or examples/deglac_lgm.nml carries the full &ctl group plus only the &vilma parameters it overrides (the yelmo defaults_file convention). A &vilma key the defaults file does not declare, such as a typo or a renamed parameter, stops the run:

./bin/vilma.x examples/deglac_lgm.nml    # default: vilma.nml
./bin/vilma.x run.nml my_defaults.nml    # optional second argument: another defaults file

vilma_remap.x and vilma_mkref.x take the same optional second argument; they too need the full &ctl group.

A few conventions worth knowing:

  • Time fields (dt_*, time_*, equil_time_max) are given in years in the namelist and converted to SI seconds on load; equil_rate_tol is in m/yr. Everything else is SI. A host that fills se%par in memory sets the time fields in seconds.
  • Earth structure is chosen by earth: a named built-in ("M3-L70-V01", or "PREM", which fills \(\rho,\mu\) from shell averages of incompressible PREM on the given layer geometry) or "custom" to assemble from the surface-first layer arrays. The radial mesh follows the depth of each element (see Radial finite elements).
  • Response solver is chosen by earth_response: "ve" (full viscoelastic, default), "elastic", or "null".
  • Memory scheme is chosen by scheme: "fe" (explicit, default) or "trap" (2nd-order, adaptive); see Rheology.
  • Degree-1 frame is chosen by deg1_frame: "cm" (default) or "cf"; see Sea-level equation.
  • Reference state is chosen by i_eq (&ctl). Relative sea level is measured against this relaxed reference.
    • 0: the state at the start of the run window is the equilibrium: the forcing’s ice at that slice and the bed name_z_bed_start, read from file_forcing (lon, lat, time) or, with remap_input = .false., from the 2-D file_ref.
    • 1 (default): a present-day reference from z_bed_ref_file and h_ice_ref_file, so that rsl is about zero at the present day.
    • 2: an equilibrium state read from z_bed_eq_file and h_ice_eq_file (e.g. an LGM state).
    • 3: the present-day reference bed plus the rsl field of rsl_restart_file.
    The variable names in these files are set by name_z_bed_ref, name_h_ice_ref and name_rsl.
  • Spin-up is controlled by equil_time_max, equil_rate_tol (the memory is relaxed under the start-of-run load until the bed is stationary) and pre_spinup_1d (a cheap 1-D pre-equilibration first); see Spin-up.
  • 3-D viscosity is enabled with l_visc_3d = .true. and a lateral log10(eta) field (visc_3d_file). With visc_3d_structure = n (n ≥ 1) the file is in the FesmData layout (log10_visc on lon, lat, depth, structure) and structure n is read. The default is the FesmData product fesmdata/Global/GLOBAL-T170/GLOBAL-T170_SolidEarth-Bagge2021.nc (Bagge et al. (2021); structure 6, v_1.0_s16). The path is relative to the run folder, which needs a fesmdata link to the FesmData store (on levante /work/ba1442/fesmdata). The model does not ship a viscosity file. The field can be shifted by f_visc_sd standard deviations (σ from the variable name_visc_sd, or f_visc_rel |log10 η|) and by region and depth with visc_shift_file. It is then clamped to visc_log10_min–visc_log10_max, and visc_3d_lid_depth and visc_3d_lid_log10max make its top effectively elastic. visc_3d_tol sets which elements are treated as genuinely 3-D. The order of these steps is given in Rheology.
  • Toroidal flow (l_toroidal) is on by default: with lateral viscosity it changes local relative sea level near former ice sheets by up to about 90 m in a deglaciation, though global means only by centimetres. l_toroidal = .false. makes such runs 1.5–1.8 times faster, for quick runs that do not depend on local sea level (see Toroidal flow).
  • Rotational feedback (l_rotation) is on by default for real-Earth runs. Set l_rotation = .false. for the non-rotating community benchmarks. The secular Love number in the Liouville equation defaults to the model’s own fluid limit; rotation_k_s (> 0) prescribes it and rotation_c_minus_a sets C − A (default 2.63e35 kg m²).
  • Forcing. The ice-thickness forcing is a lon-lat field (global or regional) on a time axis in years (kyr is converted; ages and other units stop the run) that must be strictly ascending. The ice ramps linearly between slices. The axes name_lon, name_lat are cell centres; longitudes must ascend, latitudes may ascend or descend. A regional source is remapped with its mass conserved and no ice outside it. Missing, non-finite or out-of-range values stop the run (put zero where there is no ice). The remap weights are cached in maps/ in the working directory, keyed by the grid dimensions and a hash of the axis values.
  • Output. file_out receives the diagnostic surface fields at every forcing slice, including the convergence of the sea-level solve (sle_resid, the worst relative inner residual of the interval, and sle_n_cap, the number of solves that ended at an iteration cap; the model also warns once per interval when a cap is hit) and the uplift rate dz_bed_rate [m/yr], the change of the bed dz_bed over the interval that ends at the slice divided by its length: the mean rate of that interval, so a present-day slice after a short interval (a century, say) gives the rate a GNSS station measures (in the frame of deg1_frame). The folder of file_out is created if needed. file_hor (optional) receives the horizontal surface displacement. The variables and their output precision are listed in input/vilma-variables.md. Diagnostics are rounded to that precision. write_lossless = .true. (&ctl) writes them lossless, as needed for bit-for-bit comparisons. Restart files are always lossless.
  • Restart. The driver writes the full model state to spinup/vilma_restart.nc after the spin-up and to final/vilma_restart.nc at the end of the run, in the run directory. restart_in_file (&ctl) resumes from such a file: the model is initialized at the reference and the saved state is restored, so a resumed run continues the original one bit for bit. If a spin-up is enabled, it then runs from the restored state. A restart from a lower degree is padded: its memory fills the shared low degrees and the higher degrees start at zero (not bit for bit). The run resumes at the restart’s time: a time_init that selects another start stops it, unless restart_time_reset = .true. (&ctl), which takes the restart as the state at the start of the window. A restart from another radial mesh is refused. A restart from another viscosity (radial profile or 3-D field, both stamped in the file) is refused unless restart_other_visc = .true. (&ctl), which uses its memory as a seed. Restarts of the default PREM model written before the depth-based mesh rule (220 instead of 219 elements) are therefore refused; mesh_depth = 80e3, 670e3 rebuilds the old mesh exactly. What the restart carries and what is checked on read is described in Time stepping.
  • Backend is chosen by solver: "v2" (this model, default) or "v1", the VILMA1 library behind the same driver, forcing and output for like-for-like comparison (requires a build with vilma_v1=1).

Parameter reference

Every &vilma parameter and its default is documented in input/vilma_defaults.nml, and every &ctl parameter in src/vilma_control.f90. The physics and numerics pages describe most of them in context. The remaining ones are:

Group Parameter Meaning
&vilma grid lmax maximum spherical-harmonic degree (64 in the defaults file; the in-memory record has none, so a host must set it)
nlat, nphi Gauss latitudes and longitudes; 0 gives the de-aliased defaults \(2\ell_\max + 2\) and \(4\ell_\max\)
mmax, mres maximum order (< 0: lmax) and order stride
eps_polar SHTns polar-optimization threshold (< 0: library default)
&vilma custom earth (earth = "custom" or "PREM") n_layer, r_earth, r_core number of layers (surface first), surface and core radius [m]
r_bot, r_top inner and outer radius of each layer [m]
rho, mu, eta density, shear modulus and viscosity of each layer (rho, mu filled from PREM for "PREM")
rheology per layer: 0 elastic, 1 Maxwell, 2 fluid
&vilma VILMA1 backend (solver = "v1" only) vilma_v1_jmax, vilma_v1_nsub VILMA1 degree and sub-steps per coupling interval
vilma_v1_input_dir, vilma_v1_out_dir, vilma_v1_grid_file VILMA1 input directory, scratch directory and grid
vilma_v1_visc_1d_file, vilma_v1_visc_3d_file, vilma_v1_l_prem VILMA1 viscosity files and earth-structure choice
&ctl forcing file_forcing, name_h_ice, name_time ice-thickness forcing (lon, lat, time [yr]) and its variable names
name_lon, name_lat longitude and latitude axes of every 2-D input file [deg], in both modes; with remap_input = .false. they must be the Gauss grid’s (longitudes from 0)
time_init, time_end run window [yr]; default the first and the last forcing slice
remap_input .true. (default): remap lon-lat input online; .false.: input already on the Gauss grid
&ctl start-slice bed (i_eq = 0) name_z_bed_start, file_ref bed variable, read from file_forcing (lon, lat, time), or from the 2-D Gauss-grid file_ref with remap_input = .false.
&ctl restart restart_in_file, restart_time_reset, restart_other_visc see above
&ctl output file_out, file_hor, write_lossless see above

The drivers

Executable Role
bin/vilma.x reads a reference state and an ice-thickness forcing h_ice(lon,lat,time), marches the model across the forcing (online lon-lat → Gauss remap by default), and writes the diagnostic surface fields (rsl, z_bed, h_ice, ocean function) to file_out, and optionally the horizontal surface displacement to file_hor
bin/vilma_mkref.x generates the canonical Gauss-grid present-day reference (e.g. RTopo) used by i_eq = 1; other resolutions remap it online (cached in maps/)
bin/vilma_remap.x offline conservative lon-lat → Gauss remap, to preprocess a forcing onto the model grid (remap_input = .false.)

A complete worked configuration, an LGM to present-day deglaciation, is in examples/deglac_lgm.nml.

Running with runme

runme stages a clean run directory (linking in data and input, copying the namelist, recording the parameters), then runs or submits it. The executable aliases match the drivers: main → vilma.x, mkref → vilma_mkref.x, remap → vilma_remap.x, visc3d → bench_visc3d.x.

# single local run, 8 OpenMP threads, overriding a few parameters
runme -o runs/deglac -e main --omp 8 -r \
      -p vilma.lmax=128 vilma.earth_response=ve ctl.file_out=out.nc

How a run is launched is set by two flags:

  • -r runs the executable;
  • -s prepares a SLURM submit script (and submits it, with -r);
  • -s -r prepares and submits on HPC; -s alone stages the script to inspect; omitting both stages the run directory only.

A comma list, range, or distribution in -p turns a parameter into an ensemble dimension (-p vilma.lmax=64,128), with -a naming the member directories from their parameter values.

SLURM settings live in .runme/config.toml, which is not tracked. Copy .runme/config.default.toml or run runme config init, then set the cluster (hpc, a key of runme queues --all) and the account; runme info shows what runme reads. On albedo:

hpc     = "awi_albedo"
account = "<project>"        # e.g. the p_forclima project account
omp     = 32
runme -o runs/deglac_lgm -e main -n examples/deglac_lgm.nml --omp 32 \
      -q 12h -w 06:00:00 -s -r -p vilma.lmax=128

On levante (hpc = "dkrz_levante", account = "<project>"):

runme -o runs/deglac_lgm -e main -n examples/deglac_lgm.nml --omp 64 \
      -q compute -w 06:00:00 -s -r -p vilma.lmax=128

Embedding in a host model

The CLIMBER-X coupling path uses the same code behind a single use vilma. The contract is 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 (built from the parameter record), and keeps all spherical-harmonic work inside (see Implementation):

use vilma
type(solid_earth) :: se
call vilma_par_load(se%par, "vilma.nml", defaults_file=DEFAULTS_FILE)  ! configuration lives in se%par
call solid_earth_init(se, z_bed_eq, h_ice_eq, grid=host_grid, time=t0)  ! Gauss grid + remap; clock at t0 [yr]
call solid_earth_spinup(se, h_ice_lgm)                ! optional LGM-memory spin-up (holds the clock)
do
   call solid_earth_update(se, h_ice, dt_yr)          ! advance dt_yr [years]
   ! ... read se%rsl, se%z_bed, se%dz_bed, se%dz_ss, se%C on the host grid ...
   ! ... rsl = se%dz_ss - se%dz_bed: the sea-surface change N + dphi and the bed change u ...
   ! ... se%dz_bed_rate = (u(t) - u(t - dt_yr))/dt_yr [m/yr]: the uplift rate over the interval ...
end do
call solid_earth_finalize(se)                         ! frees the grid too

use vilma is the only module a host needs: it exports the parameter record, solid_earth and its life cycle, the checks and the restart routines, and nothing internal. h_ice is the whole ice column; the model applies its own flotation rule to it. The model keeps its own clock: solid_earth_init sets it (time, default 0), each update advances it, and solid_earth_set_time relabels it without changing the state. A host never writes se%time. A negative dt_yr stops the run; dt_yr = 0 re-solves the surface fields under the ice passed in.

grid (a coords lon-lat grid) is optional: omit it when the fields are already on the model’s Gauss grid in the SHTns layout (rows north first). A grid that is the Gauss grid passes the fields through exactly, in its own row order (latitude ascending, as coords builds it, or descending). Any other global grid is remapped; a regional grid is refused. The host-grid fields are always in the order the host supplied; se%gg is the internal Gauss-grid state. The full prognostic state is persisted on restart (vilma_restart_write, vilma_restart_read), so a coupled run resumes bit for bit; the host-grid outputs are re-derived on read.

vilma_par_load overlays the user file on defaults_file; without it the user file must hold every &vilma parameter. A host that fills se%par in memory instead starts from the in-code defaults, which leave lmax unset. The rest of the interface:

Procedure Use
solid_earth_check_solver(par) validate solver (and that a "v1" build has VILMA1) before any work
solid_earth_set_time(se, time) relabel the clock [yr] without changing the state
solid_earth_sync_host(se) refresh the host-grid outputs from the Gauss-grid state
solid_earth_enable_visc_3d(se, sht) switch to the 3-D memory advance, keeping the memory (e.g. after a 1-D spin-up); init calls it when l_visc_3d is on
solid_earth_host_check(w, I_g, C, F, msg, tol) admissibility of host-resolved terms (tol default 1e-9)
vilma_restart_write(se, time, filename, folder, init) write the full state; defaults ./vilma_restart.nc, init = .false. appends a snapshot
vilma_restart_read(se, filename, time, other_visc) restore the state at time (default the last snapshot) into an initialized se
vilma_restart_time(filename) time [yr] of the last snapshot
vilma_restart_compatible(se, filename, msg, other_visc) the checks vilma_restart_read applies, without reading
vilma_write_step, vilma_write_horizontal diagnostic output, as the standalone driver writes it
vilma_io_set_table(filename) path of the variable table (default input/vilma-variables.md)
vilma_version the model version string

Two files are read relative to the working directory: the variable table input/vilma-variables.md, which every output and restart write uses (set another path with vilma_io_set_table), and the remap-weight cache maps/, keyed by the grid dimensions and a hash of the axis values.

Host-resolved terms (optional)

A host that resolves the ice margin on its own grid, inside regions of weight \(w\), can pass on that grid, per unit cell area, the grounded-ice column \(I_g = w H f_g\), the ocean fraction \(C = w (1 - f_g)\) and the flooded-area fill \(F = w \zeta^{(0)} (C - C_\text{ref})\):

call solid_earth_init(se, z_bed_eq, h_ice_eq, grid, w=w, I_g=I_g, C=C)    ! reference terms
call solid_earth_spinup(se, h_ice_lgm, w, I_g, C, F)                      ! held over every pass
call solid_earth_update(se, h_ice, dt_yr, w, I_g, C, F)                   ! ramped like the ice
ok = solid_earth_host_check(w, I_g, C, F, msg)                            ! admissibility

Inside \(w\) the sea-level solve uses these terms. On \(1 - w\) it applies its own flotation and coastline rule, and h_ice, h_ice_eq are the ice of that part (cells with \(w = 1\) ignore them). The terms are all or none and fixed for the run. They cannot be combined with solver = "v1" or sle_fixed_ocean, and \(F\) enters only through the sub-grid term (sle_subgrid). Without them, or with \(w \equiv 0\), a run is bit-identical to one without. The terms are part of the restart. Details are in the sea-level equation page and in doc/design.md §15.

Back to top

References

Bagge, M., V. Klemann, B. Steinberger, M. Latinović, and M. Thomas. 2021. “Glacial-Isostatic Adjustment Models Using Geodynamically Constrained 3D Earth Structures.” Geochemistry, Geophysics, Geosystems 22 (11): e2021GC009853. https://doi.org/10.1029/2021GC009853.