API
The public routines, the state object and the kernels
A host model uses the module elsa and nothing else. This module re-exports every type and procedure that a host needs. The modules behind it are described at the end of this page, since the benchmarks call them directly.
Module layout
| Module | File | Role |
|---|---|---|
elsa |
src/elsa.f90 |
the public interface and the update sequence |
elsa_defs |
src/elsa_defs.f90 |
the derived types, allocation and deallocation |
elsa_physics |
src/elsa_physics.f90 |
the layer kernels: advection, normalization, mass balance |
elsa_interp |
src/elsa_interp.f90 |
the horizontal maps and the vertical velocity average |
elsa_io |
src/elsa_io.f90 |
NetCDF output and restart |
elsa_precision |
src/elsa_precision.f90 |
the kinds sp, dp, wp and the constants |
The kernels in elsa_physics and elsa_interp take and return plain arrays. They hold no state, open no file and know no derived type, which makes them directly testable and safe to call inside an OpenMP region.
Public routines
elsa_init
call elsa_init(els,filename,group,time,time_end,x,y,zeta,H_ice,stagger [,restart])| Argument | Type | Content |
|---|---|---|
els |
type(elsa_class), inout |
the elsa object |
filename, group |
character | the namelist file and the group to read |
time |
real, yr | the initial time |
time_end |
real, yr | the end of the run, later than time, which sizes the layer stack |
x, y |
real (nx), (ny), m |
the host’s cell-centre axes |
zeta |
real (nz) |
the host’s sigma levels, 0 at the bed to 1 at the surface |
H_ice |
real (nx,ny), m |
the initial ice thickness |
stagger |
character | "acx_acy" or "aa", the location of the host velocities |
restart |
character, optional | a restart file; "None" or empty means a cold start |
The routine reads the parameters, builds the isochrone schedule and the grid maps, allocates the state and fills the initial layers. It ends by printing a summary of the configuration. All real arguments are either single or double precision. elsa_init may be called again on an object after elsa_end.
elsa_update
call elsa_update(els,time,H_ice,ux,uy,smb,bmb)| Argument | Type | Content |
|---|---|---|
time |
real, yr | the current absolute model time |
H_ice |
real (nx,ny), m |
ice thickness |
ux, uy |
real (nx,ny,nz), m yr\(^{-1}\) |
horizontal velocity on the host’s sigma levels |
smb, bmb |
real (nx,ny), m yr\(^{-1}\) |
surface and basal mass balance, positive for gain |
Every call adds the interval since the previous call to the time integrals of ux, uy, smb and bmb, without converting whole arrays in either precision. The routine then returns unless dt_coupling has elapsed since the last update. Otherwise it performs the update sequence given under Design. It stops the program if an isochrone is due and the layer stack is exhausted, which cannot happen with a schedule built at init.
elsa_end
call elsa_end(els)Deallocates the state, the maps and the isochrone schedule.
elsa_write_init and elsa_write_step
call elsa_write_init(els,filename,time)
call elsa_write_step(els,filename,time,n)elsa_write_init creates the output file and its axes. elsa_write_step writes the state as slice n (1-based) along the time axis. Both take time in the working precision wp, which is double. See Output.
elsa_restart_write
call elsa_restart_write(els,filename)Writes the complete state to a new file. Reading is done through the restart argument of elsa_init.
Other public names
elsa_version is a character constant holding the version string. The kinds sp, dp and wp are re-exported, as are the types below.
The elsa object
type(elsa_class) has three components: par (the parameters), now (the state) and map (the grid and the interpolation weights). A host may read any of them. It should not modify them.
els%par
The seven namelist parameters (see Parameters), and two derived quantities:
| Field | Content |
|---|---|
n_layers |
total layers allocated: n_layers_init + size(time_add) + 1 |
time_add(:) |
the times at which a new isochrone is laid down, in years |
els%now
| Field | Shape | Units | Content |
|---|---|---|---|
time |
scalar | yr | time of the last completed update |
time_acc |
scalar | yr | time of the last call, up to which the forcing is integrated |
n_top |
scalar | 1 | index of the topmost active layer |
i_add |
scalar | 1 | next entry of par%time_add to be applied |
n_reseed |
scalar | 1 | columns reseeded at the bed in the last update |
n_reseed_total |
scalar | 1 | columns reseeded since init |
t_dep |
(n_layers) |
yr | deposition time of each layer; MV if undated |
d_iso |
(nx,ny,n_layers) |
m | layer thickness, cell centres |
dsum_iso |
(nx,ny,n_layers) |
m | height of each layer top above the bed |
ux_iso |
(nx,ny,n_layers) |
m yr\(^{-1}\) | layer-mean velocity, acx faces |
uy_iso |
(nx,ny,n_layers) |
m yr\(^{-1}\) | layer-mean velocity, acy faces |
H_ice |
(nx,ny) |
m | host ice thickness at this update |
H_ice_prev |
(nx,ny) |
m | host ice thickness at the previous update |
smb, bmb |
(nx,ny) |
m yr\(^{-1}\) | host mass balance, mean over the last coupling period |
ux_lev, uy_lev |
(nx,ny,nz) |
m yr\(^{-1}\) | period-mean host velocity on elsa’s faces and the host’s sigma levels |
smb_acc, bmb_acc |
host (nx,ny) |
m | host mass balance integrated since the last update |
ux_acc, uy_acc |
host (nx,ny,nz) |
m | host velocity integrated since the last update |
Here nx and ny are the dimensions of elsa’s grid, which equal those of the host only for grid_factor = 1. The layer convention is described under Output.
els%map
| Field | Content |
|---|---|
nx, ny, nz |
dimensions of elsa’s grid and the number of host levels |
dx, dy |
grid spacing of elsa’s grid, in m |
x(:), y(:) |
cell-centre axes of elsa’s grid, in m |
zeta(:) |
the host’s sigma levels |
nx_src, ny_src, x_src(:), y_src(:) |
the host grid |
cons_x, cons_y |
conservative weights, host cells to elsa cells |
ux_x, ux_y, uy_x, uy_y |
bilinear weights, host velocity nodes to elsa faces |
Kernels
The following routines are not exported by elsa, but are public in their own modules and are exercised by the physics and interp benchmarks.
elsa_physics
| Routine | Action |
|---|---|
calc_n_substeps(ux,uy,dx,dy,dt,cfl) |
number of substeps that keep the explicit update positive |
advect_layer(d,ux,uy,dx,dy,dt,cfl,fx,fy) |
advects one layer over dt, sub-stepped; fx, fy are scratch arrays |
normalize_layers(d,H_ice,n_top) |
rescales each column to sum to H_ice |
reseed_empty_columns(d,H_ice,n_top,n_reseed) |
returns the host’s ice to emptied columns, in layer 1 |
calc_dsum(dsum,d,n_top) |
cumulative layer heights |
apply_smb(d,dm,n_top) |
adds accumulation to the top layer, or removes ablation from the top downward |
apply_bmb(d,dm,n_top,allow_pos_bmb) |
removes basal melt from the bottom upward, or adds freeze-on to layer 1 |
In these routines ux(i,j) is the velocity on the face between cells (i,j) and (i+1,j), and uy(i,j) on the face between (i,j) and (i,j+1). The last face in each direction is not read. dm is a thickness change in m, i.e., the mass balance multiplied by the time step.
elsa_interp
| Routine | Action |
|---|---|
elsa_map_init(map,x_src,y_src,zeta,stagger,grid_factor) |
defines elsa’s grid and builds all weights |
elsa_map_end(map) |
deallocates the map |
map_scalar(map,f_src,f) |
conservative remap of a cell-centred field |
map_velocity(map,ux_src,uy_src,ux,uy) |
bilinear map of the velocity onto elsa’s faces, level by level |
interp_u_column(u_layer,u_lev,zeta,H,dsum) |
thickness-average of a velocity profile over each layer of one column |
Error handling
elsa stops the program with error stop 1 and a message naming the routine for every violation of its contract: invalid parameters, a time_end that is not later than the start time, a non-uniform axis, an invalid zeta, an unknown stagger, a missing file, an isochrone time that is not later than the initial time, or a restart file that does not match the grid. It issues warnings, without stopping, for isochrones that are closer together than the coupling period and for a restart time that differs from the host time.