Procedures

Safely assigns the value of one coriolis block to another.

Returns the number of elements in the plume's state vector. This is the size of the vector returned by state_vector and taken as an argument by update.

Computes and returns a quantity which may be necessary to determine the frictional drag the plume exerts on the bottom of the ice shelf. The plume would actually tend to exert no drag on the bottom of the ice shelf, but this method is present so that there is a consistent interface with the ground data type.

Instantiates an asym_plume object with initial coniditions provided by the arguments.At present only a 1D model is supported. If information is provided for higher dimensions then it will be ignored.

Computes and returns the melt rate at the bottom of the ice shelf due to interaction with the plume.

Reads the state of the plume object from an HDF file in the specified group. This sets the thickness, velocity, temperature, salinity dataset, and parameter values.

Solves the state of the plume for the specified ice properties, at the specified time. This is done using the a quasilinearisation method and a GMRES iterative linear solver.

Returns the state vector for the current state of the plume. This takes the form of a 1D array.

Updates the state of the plume from its state vector. The state vector is a real array containing the value of each of the plume's properties at each of the locations on the grid used in discretization.

Computes and returns the density of the plume water beneath the ice shelf. The density of this water would vary depending on how much saline ambient water has been entrained into the plume versus how much fresh water has been released due to melting. However, the Boussinesq approximation is used here and only a single reference density is returned.

Writes the state of the plume object to an HDF file in the specified group. This will consist of a thickness, a velocity, a temperature, and a salinity dataset.

A wraper for the nitsol implementation of the biconjugate gradient stabilized method (BiCGSTAB) for iteratively solving linear systems. This provides a more general interface not specifically intended for use with Newton iteration. It also uses Fortran 90 features to provide a more convenient call signature.

A root finder using the hybrid bisection-secant algorithm. Returns a computed value of the root within xleft and xright once error < maxerr or has run maximum number of iterations.

Default implementation of the method getting lower and upper boundary information, which is then passed to the methods for getting and setting raw representations of fields. It returns a 1D array of length 2, indicating free boundaries (the raw data should represent all cells contained in the field, not excluding any near the boundaries).

Provides information about the type of boundary, and the number of layers of data-points needed to describe it.

Default implementation of the methods getting the boundary types for a glacier. It returns an array which indicates free boundaries.

Default implementation of the boundary_residuals method. It returns a zero-length array, effectively indicating free boundaries.

************80

Creates a new scalar field with a uniform value across all of space but a non-zero gradient.

Create a preconditioner object with the desired tolerance and maximum number of iterations.

Builds a Chebyshsev pseudospectral differentiation matrix block which can be used to solve the inverse problem. The result can only be used with fields having the same grid as the template.

Build a block in a Jacobian matrix, with the form $$\frac{\partial F}{\partial x_i} + F\Delta_i,$$ where $F$ is a scalar field and $\Delta_i$ is the differentiation operator in the $i$-direction. Additionally, a further differentiation operator may be added to the right hand side of this matrix block. Optional arguments allow for handling of boundary conditions. See the end of the documentation of the jacobian_block type for a description of how boundary conditions are treated.

Builds a Coriolis block which can be used to solve the inverse problem for the linear components of the plume momentum equations. The result can only be used with fields having the same grid as the template.

Instantiates a ground object.

Creates a new ice_sheet object with initial conditions provided by the arguments. At present only a 1D model is supported. If information is provided for higher dimensions then it will be ignored.

Build a tridiagonal matrix block for finite differences. See the end of the documentation of the fin_diff_block type for a description of how boundary conditions are treated.

Produces an ambient object which will return the specified salinity and temeprature values.

Constructs a boundary condition object for an ice shelf based on the conditions used in Dallaston et al. (2015).

Constructs a boundary condition object for an ice shelf based on the conditions used in Dallaston et al. (2015), but with seasonal variations in subglacial discharge.

Constructs a boundary condition object for an ice shelf based on the conditions used in Dallaston et al. (2015), but with the ice flux at the grounding line varying in time.

Constructs a boundary condition object which integrates an IVP from actual boundary values to calculate the staet of the plume a little upstream. This can be used to avoid boundary layers.

Constructs a boundary condition object for an ice shelf based on the conditions used in Dallaston et al. (2015).

Instantiates an instance of a viscosity object implementing Glen's flow law.

Returns the current date and time in the same format as the compile_time function.

Indicates that one layer of cells at the lower boundary in the first dimension should be omitted. This is appropriate for Dirichlet boundary conditions.

Specifies that the lower boundary in the first dimension has Dirichlet boundary conditions.

Returns the difference between the glacier conditions of the plume and the Dirichlet conditions prescribed in the model of Dallaston et al. (2015)

Indicates that one layer of cells at the upper boundary in the first dimension should be omitted for thickness. This is appropriate for Dirichlet boundary conditions.

Specifies that the upper boundary in the first dimension has Neumann boundary conditions.

A dummy subroutine which does not calculate the function.

A dummy subroutine which does not apply a preconditioner or calculate an analytic Jacobian. This can be passed to nitsol for the argument jacv.

************80

Solves the linear(ised) system represented by this finite difference block, for a given right hand side state vector (represented by a scalar field). Optionally, the differential operator can be augmented by adding an offset, i.e. a scalar field which is added to the operator.

Solves the linear(ised) system represented by this finite difference block, for a given right hand side state vector (represented by a vector field). Optionally, the differential operator can be augmented by adding an offset, i.e. a vector field which is added to the operator.

Returns the current time of the cryosphere system.

Integrates the glacier's state to time. This is done using the NITSOL package of iterative Krylov solvers. If a different algorithm for the integration is desired, then this method may be overridden in the concrete implementations of the glacier type.

Dummy routine which does nothing.

Calculates the viscosity of ice using Glen's flow law. See the documentation of the glens_law_viscosity object for a description of this parameterisation.

A subroutine which finds the values of all roots of a function (or as many as will fit into the provided arrays) within a given range of values. This subroutine uses the hybrid bisection-secant root-finding algorithm.

A global bracket finder. For a given function it finds values on each side of each of the function's roots within a given range.

A wraper for the nitsol implementation of the generalised minimal residual method (GMRES) for iteratively solving linear systems. This provides a more general interface not specifically intended for use with Newton iteration. It also uses Fortran 90 features to provide a more convenient call signature.

Returns the number of elements in the ground's state vector. This is the size of the vector returned by state_vector and taken as an argument by update.

Computes and returns a quantity which may be necessary to determine the frictional drag the ground exerts on the bottom of the ice sheet. An example would be the coefficient of friction. The description of this method is left deliberately vague so that as not to constrain how the drag is parameterized.

Computes and returns the melt rate at the bottom of the ice sheet due to interaction with the ground.

Reads the state of the ground object from the specified group in an HDF file.

Solves the state of the ground for the specified ice properties, at the specified time.

Returns the state vector for the current state of the ground. This takes the form of a 1D array.

Updates the state of the ground from its state vector. The state vector is a real array containing the value of each of the ground's properties at each of the locations on the grid used in discretization.

Computes and returns the density of the water beneath the ice sheet. This water would be subglacial discharge and would tend to lubricate the motion of the ice sheet. The density probably won't be important in the case of an ice sheet, but is included so that the ground data type can have the same interface as the plume data type.

Writes the state of the ground object to an HDF file in the specified group.

Integrate the Taylor coefficients representing the vertical structure of internal reflectors forward to the specified time. This is done using an implicit method, with the resulting linear system solved using GMRES.

Increases the time step by a factor of 2, unless doing so would take it above the maximum.

Initialise a cryosphere object from the provided components. This object will model the evolution of a glacier/ice shelf/ice sheet and its surroundings.

Integrates the cryosphere forward until the specified time is reached.

Produces a Jacobian block which is the sum of two existing blocks. Boundary conditions are set by the first operand (this).

Produces a Jacobian block which has been offset by a scalar field.

Produces a Jacobian block which has been offset by some constant increment.

Copies the contents of the rhs Jacobian block into this one. It will safely deallocate any data necessary.

A default implementation of the get_boundaries procedure pointer for the jacobian_block type. It corresponds to setting all boundaries to 0 when multiplying a field by the Jacobian block.

Computes the tridiagonal matrix used to solve for this Jacobian block.

Provides a matrix multiplication operator between a Jacobian block and a scalar field (which corresponds to a state vector).

Solves the linear(ised) system represented by this Jacobian block, for a given right hand side state vector (represented by a scalar field).

$$e = E_0 |\vec{U}\sin(\theta) \simeq E_0|\vec{U}||\nabla b|$$ Here, $E_0$ is a coefficient typically taken to be 0.036 (the default value), $\vec{U}$ is the velocity of the plume, $\theta$ is the angle of slope of the ice shelf base, and $b$ is the basal depth of the ice shelf.

$$e = \frac{c_L^2}{S_m}\sqrt{|\vec{U}|^2+\frac{g'D}{S_m}}.$$ Here, $c_L$ is an entrainment coefficient, $\vec{U}$ is the velocity of the plume, $g'$ is the reduced gravity, and $S_m$ is the turbulent Schmidt number. The Schmidt number is a function of the Richardson number $Ri = g'D/|\vec{U}|^2$: $$S_m = \frac{Ri}{0.0725(Ri + 0.186 - \sqrt{Ri^2 - 0.316Ri + 0.0346})}.$$

Returns the haline contraction coefficient.

Returns the haline contraction coefficient.

Returns the thermal contraction coefficient.

Returns the thermal contraction coefficient.

Calculates the density of the water from the temperature and salinity, using a linear equatino of state, $$\rho = \rho_0[1-\beta_T(T-T_0) + \beta_S(S-S_0)].$$

Calculates the density of the water from the temperature and salinity, using a linear equation of state, $$\rho(x,y) = \rho_0[1-\beta_T(T(x,y)-T_0) + \beta_S(S(x,y)-S_0)].$$

Calculates one form of the horizontally-averaged density of the water from the temperature and salinity, using a linear equation of state, $$\bar{rho}(x) = \rho_0[1-\beta_T(\alpha_{DT}T(x)-T_0) + \beta_S(\alpha_{DS}S(x)-S_0)].$$

Calculates another form of the horizontally-averaged density of the water from the temperature and salinity, using a linear equation of state, $$\tilde{\rho}(x) = \rho_0[1-\beta_T(\tilde{\alpha}_{DT}T(x)-T_0) + \beta_S(\tilde{\alpha}_{DS}S(x)-S_0)].$$

Calculates the derivative of the water density from the temperature and salinity, using a linear equatino of state, $$\rho = \rho_0[1-\beta_T(T-T_0) + \beta_S(S-S_0)].$$

Calculates the derivative of the average water density from the temperature and salinity, using a linear equation of state with the second type of averaging, $$\tilde{\rho} = \rho_0[1-\beta_T(\tilde{\alpha}_{DT}T-T_0) + \beta_S(\tilde{\alpha}_{DS}S-S_0)].$$

FCNRM is no longer an argument, as it needs to be calculated from the residual, which isn't known until a little way into this routine. It could be calculated by the calling routine, but that would mean doing the calculation twice.

FCNRM is no longer an argument, as it needs to be calculated from the residual, which isn't known until a little way into this routine. It could be calculated by the calling routine, but that would mean doing the calculation twice.

FCNRM is no longer an argument, as it needs to be calculated from the residual, which isn't known until a little way into this routine. It could be calculated by the calling routine, but that would mean doing the calculation twice.

Returns the number of elements in the plume's state vector. This is the size of the vector returned by state_vector and taken as an argument by update.

Computes and returns a quantity which may be necessary to determine the frictional drag the plume exerts on the bottom of the ice shelf. The plume would actually tend to exert no drag on the bottom of the ice shelf, but this method is present so that there is a consistent interface with the ground data type.

Instantiates a plume object with initial coniditions provided by the arguments.At present only a 1D model is supported. If information is provided for higher dimensions then it will be ignored.

Computes and returns the melt rate at the bottom of the ice shelf due to interaction with the plume.

Reads the state of the plume object from an HDF file in the specified group. This sets the thickness, velocity, temperature, salinity dataset, and parameter values.

Solves the state of the plume for the specified ice properties, at the specified time. This is done using the a quasilinearisation method and a GMRES iterative linear solver.

Returns the state vector for the current state of the plume. This takes the form of a 1D array.

Updates the state of the plume from its state vector. The state vector is a real array containing the value of each of the plume's properties at each of the locations on the grid used in discretization.

Computes and returns the density of the plume water beneath the ice shelf. The density of this water would vary depending on how much saline ambient water has been entrained into the plume versus how much fresh water has been released due to melting. However, the Boussinesq approximation is used here and only a single reference density is returned.

Writes the state of the plume object to an HDF file in the specified group. This will consist of a thickness, a velocity, a temperature, and a salinity dataset.

Use Picard iteration to approximately multiply the state vector by the inverse Jacobian. The details for this procedure are in the documentation of the preconditioner type.

Assigns this object to another equation of state object, allowing the definided assignment for the precalculated density field to work correctly.

Returns the haline contraction coefficient.

Returns the thermal contraction coefficient.

Returns the density corresponding to the prescribed salinity, as calculated in the constructor.

Calculates the derivative of the water density.

Solves the linear(ised) system represented by this finite difference block, for a given right hand side state vector (represented by a scalar field).

Solves the linear(ised) system represented by this finite difference block, for a given right hand side state vector (represented by a vector field).

This is an iterative method to solve nonlinear systems of ODEs. Consider $$L^{(n)}\vec{u}(x) = \vec{f}(\vec{u}(x), \vec{u}^{(1)}(x), \ldots, \vec{u}^{n-1}(x), x), \qquad \vec{u} \in \mathbb{R}^{m}.$$ For a domain from $(0,b)$, boundary conditions are specified by $$g_k(\vec{u}(0) = 0, \vec{u}^{(1)}(0), \ldots, \vec{u}^{n-1}(0)), \qquad k = 1, \ldots, l$$ $$g_k(\vec{u}(b), \vec{u}^{(1)}(b), \ldots, \vec{u}^{n-1}(b)) = 0, \qquad k = l+1, \ldots, mn.$$ Here, $L^{(n)}$ is an $n$th order ordinary differential operator, and $\vec{f}, g_1, g_2, \ldots, g_{nm}$ are nonlinear functions of $\vec{u}(x)$ and its first $n-1$ derivatives.

Reads the data describing the cryosphere from an HDF5 file on the disc. h5open_f must have been called once prior to using this method. After the method has been used, h5close_f must be called once before the end of the program.

Reads the data describing the ice component of the cryosphere from an HDF5 file on the disc. Data on anything below the ice is ignored. h5open_f must have been called once prior to using this method. After the method has been used, h5close_f must be called once before the end of the program.

Reads the data describing the part of the cryosphere beneath the ice from an HDF5 file on the disc. Data on the ice itself is ignored. h5open_f must have been called once prior to using this method. After the method has been used, h5close_f must be called once before the end of the program.

Reuces the time step by a factor of 2, unless doing so would take it below the minimum value.

Returns a field containing the boundary values for the specified boundary location.

Indicates that the lower boundary is Dirichlet and the upper boundary is Neumann.

Indicates that one layer of cells at the lower boundary in the first dimension should be omitted. This is appropriate for Dirichlet boundary conditions.

Specifies that the lower boundary in the first dimension has Dirichlet boundary conditions.

Returns the difference between the glacier conditions of the plume and the Dirichlet conditions prescribed in the model of Dallaston et al. (2015)

Returns a field containing the salinity boundary values for the specified boundary location.

Sets the time at which boundary conditions are to be calculated.

Returns a field containing the temperature boundary values for the specified boundary location.

Returns a field containing the thickness boundary values for the specified boundary location.

Indicates that the lower boundary is Dirichlet and the upper boundary is free.

Indicates that one layer of cells at the upper boundary in the first dimension should be omitted for thickness. This is appropriate for Dirichlet boundary conditions.

Specifies that the upper boundary in the first dimension has Neumann boundary conditions.

Returns a field containing the velocity boundary values for the specified boundary location.

Sets the time at which boundary conditions are to be calculated.

Copies the data from one ice sheet into another. This is only needed due to a bug in gfortran which means that the intrinsic assignment for glacier types is not using the appropriate defined assignment for the field components.

Returns the number of elements in the ice sheet's state vector. This is the size of the vector returned by residual and state_vector and taken as an argument by update.

Returns the density of the ice in the sheet, which is assumed to be uniform across its domain.

Provides a preconditioner for the nonlinear solver trying to bring the residual to zero. The Jacobian is approximated as a block matrix, where each block is a tridiagonal matrix using a finite difference method for differentiation.

Reads the state of the ice shelf object from the specified group in an HDF5 file. This sets the thickness, the velocity, and parameter values.

Returns the residual when the current state of the glacier is run through the system of equations describing it. The residual takes the form of a 1D array, with each element respresenting the residual for one of the equations in the system.

Sets the time information held by the ice sheet object. This is the time at which the ice sheet is in its current state.

Computes the ice sheet velocity at the current time with the current ice thickness.

Returns the state vector for the current state of the ice sheet. This takes the form of a 1D array.

Returns the density of the ice in the sheet, which is assumed to be uniform across its domain.

Returns the thickness of the ice sheet across its domain.

Calculates the time step for integrating the ice sheet, using the CFL condition.

Updates the state of the ice sheet from its state vector. The state vector is a real array containing the value of each of the ice sheet's properties at each of the locations on the grid used in descretization.

Returns the velocity of the ice sheet across its domain.

Writes the state of the ice sheet object to an HDF file in the specified group. This will consist of a thickness and a velocity dataset.

Copies the data from one ice shelf into another. This is only needed due to a bug in gfortran which means that the intrinsic assignment for glacier types is not using the appropriate defined assignment for the field components.

Returns the number of elements in the ice shelf's state vector. This is the size of the vector returned by residual and state_vector and taken as an argument by update.

Returns the density of the ice in the shelf, which is assumed to be uniform across its domain.

Initialises an ice_shelf object with initial conditions provided by the arguments. At present only a 1D model is supported. If information is provided for higher dimensions then it will be ignored.

Integrates the glacier's state forward to time. This is done using an explicit method for the thickness and a Newton's solver for velocity.

Returns the a vector representing the current state of the internal reflectors in the ice shelf. This takes the form of a 1D array. The routien is only used for debugging purposes.

Provides a preconditioner for the nonlinear solver trying to bring the residual to zero. The Jacobian is approximated as a block matrix, where each block is a tridiagonal matrix using a finite difference method for differentiation.

Reads the state of the ice shelf object from the specified group in an HDF5 file. This sets the thickness, the velocity, and parameter values.

Returns the residual when the current state of the glacier is run through the system of equations describing it. The residual takes the form of a 1D array, with each element respresenting the residual for one of the equations in the system.

Sets the time information held by the ice shelf object. This is the time at which the ice sheet is in its current state.

Computes the ice shelf velocity at the current time with the current ice thickness.

Returns the state vector for the current state of the ice shelf. This takes the form of a 1D array.

Returns the density of the ice in the shelf, which is assumed to be uniform across its domain.

Returns the thickness of the ice shelf across its domain.

Calculates the time step for integrating the ice shelf, using the CFL condition.

Updates the state of the ice shelf from its state vector. The state vector is a real array containing the value of each of the ice shelf's properties at each of the locations on the grid used in discretization.

Returns the velocity of the ice shelf across its domain.

Writes the state of the ice shelf object to an HDF file in the specified group. This will consist of a thickness and a velocity dataset.

Indicates that the lower boundary is Dirichlet and the upper boundary is Neumann.

Returns a field containing the salinity boundary values for the specified boundary location.

Returns a field containing the temperature boundary values for the specified boundary location.

Returns a field containing the thickness boundary values for the specified boundary location.

Indicates that the lower boundary is Dirichlet and the upper boundary is free.

Returns a field containing the velocity boundary values for the specified boundary location.

Inverts the linear portions of the plume momentum equation with the provided data. This is done by solving the linear ODE described in the documentation for the coriolis_block type. The block object must first have been initialised using the constructor.

Returns the state vector for the current state of the cryosphere. This takes the form of a 1D array. This routine is mainly useful for unit-testing.

Returns the number of elements in the plume's state vector. This is the size of the vector returned by state_vector and taken as an argument by update.

Computes and returns a quantity which may be necessary to determine the frictional drag the plume exerts on the bottom of the ice shelf. The plume would actually tend to exert no drag on the bottom of the ice shelf, but this method is present so that there is a consistent interface with the ground data type.

Instantiates a plume object with initial coniditions provided by the arguments.At present only a 1D model is supported. If information is provided for higher dimensions then it will be ignored.

Computes and returns the melt rate at the bottom of the ice shelf due to interaction with the plume.

Reads the state of the plume object from an HDF file in the specified group. This sets the thickness, velocity, temperature, salinity dataset, and parameter values.

Would normally solve, but for this static implementation of the plume it does nothing.

Returns the state vector for the current state of the plume. This takes the form of a 1D array.

Updates the state of the plume from its state vector. The state vector is a real array containing the value of each of the plume's properties at each of the locations on the grid used in discretization.

Computes and returns the density of the plume water beneath the ice shelf. The density of this water would vary depending on how much saline ambient water has been entrained into the plume versus how much fresh water has been released due to melting. However, the Boussinesq approximation is used here and only a single reference density is returned.

Writes the state of the plume object to an HDF file in the specified group. This will consist of a thickness, a velocity, a temperature, and a salinity dataset.

A wraper for the nitsol implementation of the transpose-free quasi-minimal residual method (TFQMR) for iteratively solving linear systems. This provides a more general interface not specifically intended for use with Newton iteration. It also uses Fortran 90 features to provide a more convenient call signature.

Calculates an appropriate time step with which to integrate the cryosphere so as not to cause numerical instability.

${\rm field} = {\rm field}$

$\frac{\partial^{\rm order}}{\partial x_{\rm dir}^{\rm order}}{\rm field}$

$\nabla{\rm field}$

Evaluates whether two scalar fields are equal within a tolerance, specified by set_tol.

Returns the ambient ocean salinity.

Returns the ambient ocean temperature.

Calculates the boundary values to use at the current time with the current ice thickness.

Indicates that the lower boundary is Dirichlet and the upper boundary is Neumann.

Returns a field containing the salinity boundary values for the specified boundary location.

Returns a field containing the temperature boundary values for the specified boundary location.

Returns a field containing the thickness boundary values for the specified boundary location.

Indicates that the lower boundary is Dirichlet and the upper boundary is free.

Returns a field containing the velocity boundary values for the specified boundary location.

Returns a field containing the boundary values for the specified boundary location.

Writes the data describing the cryosphere to the disc as an HDF5 file. h5open_f must have been called once prior to using this method. After the method has been used, h5close_f must be called once before the end of the program.