6.1. Model domain, grid and time¶
6.1.1. Model domain¶
6.1.1.1. Selecting a pre-defined domain¶
SICOPOLIS provides several pre-defined model domains. They can be chosen by defining a domain code in the run-specs header as follows:
#define GRL
! Simulated domain:
! ANT - Antarctica
! GRL - Greenland
! NHEM - Entire northern hemisphere
! LCIS - Laurentide and Cordilleran ice sheets
! SCAND - Fennoscandian and Eurasian ice sheets
! TIBET - Tibetan ice sheet
! ASF - Austfonna
! NPI - North Patagonian ice field
! MOCHO - Mocho-Choshuenco ice cap
! EISMINT - EISMINT (Phase 2 SGE and modifications)
! HEINO - ISMIP HEINO
! NMARS - North polar cap of Mars
! SMARS - South polar cap of Mars
! XYZ - Unspecified domain
This example would select the domain for the Greenland ice sheet. Correspondingly for the other listed domains.
6.1.1.2. Setting up a new domain¶
In addition to the pre-defined domains, there is an unspecified domain XYZ. This framework allows creating and testing new domains (e.g., some ice cap) quite easily. Create a run-specs header (by using an existing one as a template). Define the domain code as
#define XYZ
in the header. Input files (topography etc.) must be placed in sico_in/xyz and specified in the header as usual. With the parameter XYZ_SPECIAL_MODULES, it can be selected whether common or special modules will be used:
#define XYZ_SPECIAL_MODULES 0
! Only for unspecified domain XYZ:
! 0 : Common modules will be used
! 1 : Special modules 'boundary_m', 'sico_init_m'
! and 'sico_vars_m' (in 'src/subroutines/xyz/')
! will be used
In case of 0 (default), no further coding is required. In case of 1, the special modules boundary_m (climate forcing), sico_init_m (initializations) and sico_vars_m (additional global variables) must be created in src/subroutines/xyz/ (by using existing ones as templates).
6.1.2. Spatial grid¶
6.1.2.1. Horizontal coordinates¶
In principle, SICOPOLIS allows using any orthogonal coordinates on the Earth’s (or another planet’s) surface, provided that the two components \(g_{11}\) and \(g_{22}\) of the metric tensor are known (see “Orthogonal coordinates on the Earth’s surface/General considerations”). Three options are currently implemented and can be selected in the run-specs header by the parameter GRID:
0: Cartesian coordinates in the stereographic plane without distortion correction.1: Cartesian coordinates in the stereographic plane with distortion correction.2: Geographic coordinates (longitude/latitude).
For the cases 0 and 1, a polar stereoraphic projection for an ellipsoidal planet model is employed. The projection parameters are defined in the physical-parameter file:
R(\(=R_\mathrm{e}\), mean radius of planet, in m),A(\(=A\), semi-major axis of planet, in m),F_INV(\(=F_\mathrm{inv}\), inverse flattening of planet),LATD0(\(=\varphi_0\), standard parallel, in deg, +/– for N/S),LOND0(\(=\lambda_0\), central meridian, in deg, +/– for E/W).
For case 2 (geographic coordinates), only the parameters R, A and F_INV are relevant.
The semi-minor axis \(B\) can be computed from the above parameters via \(B=A(1-1/F_\mathrm{inv})\). If \(F_\mathrm{inv}=\infty\) (i.e., F_INV set to any finite value greater than 1.0e+10_dp), a spherical planet with radius \(R_\mathrm{e}\) is assumed instead of an elliptical one.
The components \(g_{11}\) and \(g_{22}\) of the metric tensor are computed for simplicity under the assumption of a spherical planet. For case 0 (distortion correction neglected), they are set to unity. For the cases 1 and 2, the derivations are found in “Orthogonal coordinates on the Earth’s surface/Polar stereographic projection” and “Orthogonal coordinates on the Earth’s surface/Geographic coordinate system”, respectively. The computations are carried out in the module metric_m.
For the most common case of Cartesian coordinates \(x\) and \(y\) in the stereographic plane (or any other projection plane), let the domain be the rectangle described by \([x_0,x_\mathrm{max}]\), \([y_0,y_\mathrm{max}]\). It is discretized by a regular (structured) grid with horizontal resolution \(\Delta{x}\), which is the same for the \(x\)- and \(y\)-directions. SICOPOLIS employs grid-line registration (Fig. 6.1). Hence, the location of the grid points (nodes) \(x_i\) and \(y_j\) is given by
where the notation \(a\,(b)\,c\) means “from \(a\) to \(c\) in steps of \(b\)”. Note that the indices \(i\) and \(j\) run from 0, so that the number of grid points is actually \(i_\mathrm{max}+1\) and \(j_\mathrm{max}+1\), respectively. In the run-specs headers, the parameters to be defined are
X0(\(=x_0\), \(x\) coordinate of the origin point in km),Y0(\(=y_0\), \(y\) coordinate of the origin point in km),DX(\(=\Delta{}x\), horizontal grid spacing in km),IMAX(\(=i_\mathrm{max}\), maximum value of the index \(i\)),JMAX(\(=j_\mathrm{max}\), maximum value of the index \(j\)).
Fig. 6.1 Two grid registration types for a regular (structured) grid. (a) Grid-line registration: the nodes (solid circles) are located at the intersections of the grid lines. (b) Pixel registration: the nodes (solid circles) are centred in the areas between grid lines. The red shades indicate the footprints of the nodes (areas represented by the values at the nodes). SICOPOLIS employs grid-line registration (a).¶
6.1.2.2. Vertical coordinate¶
For the vertical (\(z\)) direction, a terrain-following (“sigma”) transformation is employed that maps vertical columns in the physical space onto \([0,1]\) intervals. If the polythermal two-layer method (POLY, see Section “Ice thermodynamics”) is employed, this mapping is done separately for the upper cold-ice layer (\(\zeta_\mathrm{c}\) domain), the lower temperate-ice layer (\(\zeta_\mathrm{t}\) domain) and the lithosphere layer (\(\zeta_\mathrm{r}\) domain). The transformation is linear for the \(\zeta_\mathrm{t}\) and \(\zeta_\mathrm{r}\) domains. However, for the \(\zeta_\mathrm{c}\) domain, exponential stretching is used so that equidistant grid points in the transformed domain map on grid points concentrating towards the base in the physical \(z\)-coordinate:
where the geometric quantities are explained in Fig. 6.2 and \(a\) is the exponential stretch parameter for the \(\zeta_\mathrm{c}\) domain. For this parameter, \(a=2\) is a typical choice, while the limit \(a=0\) produces a linear transformation.
Fig. 6.2 Cross section through a polythermal ice sheet (vertically exaggerated).¶
The location of the grid points in the three transformed domains is given by
The numbers of grid points result as \(k_\mathrm{c,max}+1\), \(k_\mathrm{t,max}+1\) and \(k_\mathrm{r,max}+1\), respectively. The parameters in the run-specs headers are
KCMAX(\(=k_\mathrm{c,max}\), maximum value of the index \(k_\mathrm{c}\)),KTMAX(\(=k_\mathrm{t,max}\), maximum value of the index \(k_\mathrm{t}\)),KRMAX(\(=k_\mathrm{r,max}\), maximum value of the index \(k_\mathrm{r}\)),DEFORM(\(=a\), exponential stretch parameter for the \(\zeta_\mathrm{c}\) domain).
For all other thermodynamics schemes (ENTC, ENTM, COLD, ISOT; see Section “Ice thermodynamics”), the entire ice column (no matter whether cold or temperate) is mapped on the \(\zeta_\mathrm{c}\) domain. The \(\zeta_\mathrm{t}\) domain is then redundant and collapses onto the ice base:
For technical reasons, the \(\zeta_\mathrm{t}\) domain is still present and should be assigned three grid points, that is, KTMAX should be set to 2.
6.1.2.3. Staggered grid¶
A staggered Arakawa C grid is used for reasons of numerical stability (Arakawa and Lamb [1]). This means that the components of the velocity (\(v_{x}\), \(v_{y}\), \(v_{z}\)) and the volume flux (depth-integrated horizontal velocity; \(Q_{x}\), \(Q_{y}\)) are defined in between the grid points, while all other quantities are defined on the grid points.
Fig. 6.3 Arakawa C grid for the ice thickness \(H\) and the volume flux \((Q_{x},Q_{y})\) in one grid cell in the horizontal plane, centred around the grid point \((i,j)\).¶
For the example of the 2D fields ice thickness and volume flux, this is illustrated in Fig. 6.3. For the 3D fields, the principle is the same: The velocity components are defined in between the grid points:
while other fields like the temperature \(T\) are defined on the grid points: \(T_{ijk}\). Note that, depending on the layer (\(\zeta_\mathrm{c,t,r}\) domains, see above), the index \(k\) can either be \(k_\mathrm{c}\), \(k_\mathrm{t}\) or \(k_\mathrm{r}\).
Since half-integer indices are not allowed in Fortran, they are rounded down in the code. For the example of the velocity component \(v_{x}\) in the \(k_\mathrm{c}\) domain:
6.1.3. Topography¶
Gridded present-day topographies that match the horizontal grid must be provided in either NetCDF (*.nc) or ASCII (any other file extension) format. They can be specified in the run-specs header as follows (example with NetCDF files for simulation repo_grl16_bm5_ss25ka):
#define ZS_PRESENT_FILE 'grl_bm5_16_topo.nc'
! Name of the file containing the present-day
! ice-surface topography
! (if NetCDF, variable name must be 'zs')
#define ZB_PRESENT_FILE 'grl_bm5_16_topo.nc'
! Name of the file containing the present-day
! ice-base topography (only for ANF_DAT==1)
! (if NetCDF, variable name must be 'zb')
#define ZL_PRESENT_FILE 'grl_bm5_16_topo.nc'
! Name of the file containing the present-day
! lithosphere-surface topography
! (only for ANF_DAT==1)
! (if NetCDF, variable name must be 'zl')
#define ZL0_FILE 'grl_bm5_16_zl0_llra.nc'
! Name of the file containing the topography
! of the relaxed lithosphere surface
! (if NetCDF, variable name must be 'zl0')
#define MASK_PRESENT_FILE 'grl_bm5_16_topo.nc'
! Name of the file containing the present-day
! ice-land-ocean mask
! (if NetCDF, variable name must be 'mask')
#define MASK_REGION_FILE 'none'
! Name of the file containing the region mask
! ('none' if no file is to be defined)
! (if NetCDF, variable name must be 'mask_region')
6.1.4. Model time¶
Model time runs from an initial time \(t_\mathrm{init}\) until a final time \(t_\mathrm{final}\). For the numerical solution, this interval is discretized by different time steps:
\(\Delta{}t\): dynamic time step, for computing velocity and topography,
\(\Delta{}t_\mathrm{temp}\): thermodynamic time step, for computing temperature, water content and age,
\(\Delta{}t_\mathrm{wss}\): isostatic time step, for computing the isostatic steady-state displacement of the lithosphere (only if the elastic-lithosphere model is chosen).
The thermodynamic and isostatic time steps must be equal to or integer multiples of the dynamic time step. The values can be specified in the run-specs header as follows:
TIME_INIT0(\(=t_\mathrm{init}\), initial time, in a),TIME_END0(\(=t_\mathrm{final}\), final time, in a),DTIME0(\(=\Delta{}t\), dynamic time step, in a),DTIME_TEMP0(\(=\Delta{}t_\mathrm{temp}\), thermodynamic time step, in a),DTIME_WSS0(\(=\Delta{}t_\mathrm{wss}\), isostatic time step, in a).
Further, there is a parameter YEAR_ZERO that specifies the SICOPOLIS year zero in astronomical year numbering [= signed year CE (AD)]. For instance, if set to 1990, the time count of SICOPOLIS will be relative to the calendar year 1990 CE. TIME_INIT0 and TIME_END0 must be given in this SICOPOLIS calendar.