Abstract. We present SELEN4 (SealEveL EquatioN solver), an open-source program written in Fortran 90 that simulates the glacial isostatic adjustment (GIA) process in response to the melting of the Late Pleistocene ice sheets. Using a pseudo-spectral approach complemented by a spatial discretization on an icosahedron-based spherical geodesic grid, SELEN4 solves a generalized sea-level equation (SLE) for a spherically symmetric Earth with linear viscoelastic rheology, taking the migration of the shorelines and the rotational feedback on sea level into account. The approach is gravitationally and topographically self-consistent, since it considers the gravitational interactions between the solid Earth, the cryosphere, and the oceans, and it accounts for the evolution of the Earth's topography in response to changes in sea level. The SELEN4 program can be employed to study a broad range of geophysical effects of GIA, including past relative sea-level variations induced by the melting of the Late Pleistocene ice sheets, the time evolution of paleogeography and of the ocean function since the Last Glacial Maximum, the history of the Earth's rotational variations, present-day geodetic signals observed by Global Navigation Satellite Systems, and gravity field variations detected by satellite gravity missions like GRACE (the Gravity Recovery and Climate Experiment). The “GIA fingerprints” constitute a standard output of SELEN4. Along with the source code, we provide a supplementary document with a full account of the theory, some numerical results obtained from a standard run, and a user guide. Originally, the SELEN program was conceived by Giorgio Spada (GS) in 2005 as a tool for students eager to learn about GIA, and it has been the first SLE solver made available to the community.