ice_shelf_mod Module

Provides a concrete implementation of the glacier type, using a vertically integrated model of an ice shelf.


Uses

  • module~~ice_shelf_mod~~UsesGraph module~ice_shelf_mod ice_shelf_mod logger_mod logger_mod module~ice_shelf_mod->logger_mod h5lt h5lt module~ice_shelf_mod->h5lt hdf5 hdf5 module~ice_shelf_mod->hdf5 module~glacier_mod glacier_mod module~ice_shelf_mod->module~glacier_mod factual_mod factual_mod module~ice_shelf_mod->factual_mod module~jacobian_block_mod jacobian_block_mod module~ice_shelf_mod->module~jacobian_block_mod module~glacier_boundary_mod glacier_boundary_mod module~ice_shelf_mod->module~glacier_boundary_mod module~dallaston2015_glacier_boundary_mod dallaston2015_glacier_boundary_mod module~ice_shelf_mod->module~dallaston2015_glacier_boundary_mod module~newtonian_viscosity_mod newtonian_viscosity_mod module~ice_shelf_mod->module~newtonian_viscosity_mod iso_fortran_env iso_fortran_env module~ice_shelf_mod->iso_fortran_env penf penf module~ice_shelf_mod->penf module~nitsol_mod nitsol_mod module~ice_shelf_mod->module~nitsol_mod module~viscosity_mod viscosity_mod module~ice_shelf_mod->module~viscosity_mod module~boundary_types_mod boundary_types_mod module~ice_shelf_mod->module~boundary_types_mod module~glacier_mod->logger_mod module~glacier_mod->hdf5 module~glacier_mod->factual_mod module~glacier_mod->iso_fortran_env module~glacier_mod->penf module~glacier_mod->module~nitsol_mod module~jacobian_block_mod->logger_mod module~jacobian_block_mod->factual_mod module~jacobian_block_mod->iso_fortran_env module~jacobian_block_mod->penf module~jacobian_block_mod->module~boundary_types_mod f95_lapack f95_lapack module~jacobian_block_mod->f95_lapack module~glacier_boundary_mod->factual_mod module~glacier_boundary_mod->iso_fortran_env module~glacier_boundary_mod->module~boundary_types_mod module~dallaston2015_glacier_boundary_mod->factual_mod module~dallaston2015_glacier_boundary_mod->module~glacier_boundary_mod module~dallaston2015_glacier_boundary_mod->iso_fortran_env module~dallaston2015_glacier_boundary_mod->module~boundary_types_mod module~newtonian_viscosity_mod->factual_mod module~newtonian_viscosity_mod->iso_fortran_env module~newtonian_viscosity_mod->module~viscosity_mod module~nitsol_mod->iso_fortran_env module~viscosity_mod->factual_mod module~viscosity_mod->iso_fortran_env

Contents


Variables

TypeVisibility AttributesNameInitial
character(len=9), public, parameter:: hdf_type_name ='ice_shelf'
character(len=9), public, parameter:: hdf_thickness ='thickness'
character(len=8), public, parameter:: hdf_velocity ='velocity'
character(len=6), public, parameter:: hdf_lambda ='lambda'
character(len=4), public, parameter:: hdf_zeta ='zeta'
character(len=3), public, parameter:: hdf_chi ='chi'
character(len=10), public, parameter:: hdf_n_kappa ='num_kappas'
character(len=15), public, parameter:: hdf_kappa ='("kappa_",i0.4)'

Abstract Interfaces

abstract interface

  • private pure function kappa_init_func(n, location) result(kappa)

    Abstract interface for function providing the Taylor coefficients describing the distribution of internal reflectors within an ice shelf.

    Arguments

    Type IntentOptional AttributesName
    integer, intent(in) :: n

    The index of the Taylor coefficient being calculated

    real(kind=r8), intent(in), dimension(:):: location

    The position $\vec{x}$ at which to compute the coefficient

    Return Value real(kind=r8)

    The value of coefficient n at location


Derived Types

type, public, extends(glacier) :: ice_shelf

A concrete implementation of the glacier type, using a vertically integrated model of an ice shelf. This model is 1-dimensional only.

Components

TypeVisibility AttributesNameInitial
type(cheb1d_scalar_field), private :: thickness

Thickness of ice shelf,

type(cheb1d_vector_field), private :: velocity

Flow velocity of ice shelf,

type(cheb1d_scalar_field), private :: eta

Viscosity of the ice,

type(cheb1d_scalar_field), private, dimension(:), allocatable:: kappa

Taylor coefficients for the vertical structure of a Lagrangian tracer representing englacial layers/internal reflectors.

real(kind=r8), private :: lambda

The dimensionless ratio $\lambda \equiv \frac{\rho_0m_0x_0}{\rho_iH-0u_0}$

real(kind=r8), private :: chi

The dimensionless ratio

real(kind=r8), private :: zeta

The dimensionless ratio , corresponding to the Reynolds number. Currently unused.

real(kind=r8), private :: courant

The Courant number to use when calculating the time step.

class(abstract_viscosity), private, allocatable:: viscosity_law

An object representing the model used for ice viscosity.

class(glacier_boundary), private, allocatable:: boundaries

An object specifying the boundary conditions for the ice shelf.

real(kind=r8), private :: max_dt

The maximu allowable time step

real(kind=r8), private :: time

The time at which the ice shelf is in this state.

integer, private :: thickness_size

The number of data values in the thickness field.

integer, private :: velocity_size

The number of data values in the velocity field.

integer, private :: boundary_start

The number of data values needed to represent the boundary conditions.

integer, private :: thickness_lower_bound_size

The number of data values needed to represent the lower boundary conditions for thickness.

integer, private :: thickness_upper_bound_size

The number of data values needed to represent the upper boundary conditions for thickness.

integer, private :: velocity_lower_bound_size

The number of data values needed to represent the lower boundary conditions for velocity.

integer, private :: velocity_upper_bound_size

The number of data values needed to represent the upper boundary conditions for thickness.

type(jacobian_block), private :: thickness_jacobian

A representation of the Jacobian for the ice shelf thickness.

type(jacobian_block), private :: velocity_jacobian

A representation of the Jacobian for the ice shelf velocity.

logical, private :: stale_eta

Indicates whether the viscosity needs updating.

logical, private :: stale_jacobian

Indicates if the Jacobians are stale and in need of updating.

Type-Bound Procedures

generic, public :: assignment(=) => assign
procedure, public :: integrate => glacier_integrate

Performs a time-step of the integration, taking the state of the glacier to the specified time using the provided melt-rate data.

procedure, public :: initialise => shelf_initialise
procedure, public :: ice_thickness => shelf_thickness
procedure, public :: ice_density => shelf_density
procedure, public :: ice_temperature => shelf_temperature
procedure, public :: residual => shelf_residual
procedure, public :: update => shelf_update
procedure, public :: precondition => shelf_precondition
procedure, public :: set_time => shelf_set_time
procedure, public :: data_size => shelf_data_size
procedure, public :: state_vector => shelf_state_vector
procedure, public :: kappa_vector => shelf_kappa_vector
procedure, public :: read_data => shelf_read_data
procedure, public :: write_data => shelf_write_data
procedure, public :: time_step => shelf_time_step
procedure, public :: solve_velocity => shelf_solve_velocity
procedure, private :: assign => shelf_assign
procedure, public :: integrate_layers => ice_shelf_integrate_layers

Functions

private function shelf_thickness(this) result(thickness)

Author
Christopher MacMackin
Date
April 2016

Returns the thickness of the ice shelf across its domain.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value class(scalar_field), pointer

The ice thickness.

private function shelf_velocity(this) result(velocity)

Author
Christopher MacMackin
Date
July 2016

Returns the velocity of the ice shelf across its domain.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value class(vector_field), pointer

The ice velocity.

private pure function shelf_density(this) result(density)

Author
Christopher MacMackin
Date
April 2016

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

Read more…

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value real(kind=r8)

The ice density.

private pure function shelf_temperature(this) result(temperature)

Author
Christopher MacMackin
Date
April 2016

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

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value real(kind=r8)

The ice density.

private function shelf_residual(this, previous_states, melt_rate, basal_drag_parameter, water_density) result(residual)

Author
Christopher MacMackin
Date
April 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this
class(glacier), intent(in), dimension(:):: previous_states

The states of the glacier in the previous time steps. The first element of the array should be the most recent. The default implementation will only make use of the most recent state, but the fact that this is an array allows potential other implementations to use older states for higher-order integration methods.

class(scalar_field), intent(in) :: melt_rate

Thickness of the ice above the glacier.

class(scalar_field), intent(in) :: basal_drag_parameter

A paramter, e.g. coefficient of friction, needed to calculate the drag on basal surface of the glacier.

real(kind=r8), intent(in) :: water_density

The density of the water below the glacier.

Return Value real(kind=r8), dimension(:), allocatable

The residual of the system of equations describing the glacier.

private function shelf_precondition(this, previous_states, melt_rate, basal_drag_parameter, water_density, delta_state) result(preconditioned)

Author
Chris MacMackin
Date
January 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
class(glacier), intent(in), dimension(:):: previous_states

The states of the glacier in the previous time steps. The first element of the array should be the most recent. The default implementation will only make use of the most recent state, but the fact that this is an array allows overriding methods to use older states for higher-order integration methods.

class(scalar_field), intent(in) :: melt_rate

Thickness of the ice above the glacier

class(scalar_field), intent(in) :: basal_drag_parameter

A paramter, e.g. coefficient of friction, needed to calculate the drag on basal surface of the glacier.

real(kind=r8), intent(in) :: water_density

The density of the water below the glacier

real(kind=r8), intent(in), dimension(:):: delta_state

The change to the state vector which is being preconditioned.

Return Value real(kind=r8), dimension(:), allocatable

The result of applying the preconditioner to delta_state.

private pure function shelf_data_size(this)

Author
Christopher MacMackin
Date
August 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value integer

The number of elements in the ice shelf's state vector.

private function shelf_state_vector(this) result(state_vector)

Author
Christopher MacMackin
Date
April 2016

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

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value real(kind=r8), dimension(:), allocatable

The state vector describing the ice shelf.

private function shelf_kappa_vector(this) result(kappa_vector)

Author
Christopher MacMackin
Date
April 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value real(kind=r8), dimension(:), allocatable

The state vector describing the ice shelf.

private function shelf_time_step(this) result(dt)

Author
Chris MacMackin
Date
December 2016

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

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this

Return Value real(kind=r8)

The time-step to use


Subroutines

private subroutine shelf_initialise(this, domain, resolution, thickness, velocity, temperature, viscosity_law, boundaries, lambda, chi, zeta, courant, max_dt, kappa, n_kappa)

Author
Christopher MacMackin
Date
April 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(out) :: this
real(kind=r8), intent(in), dimension(:,:):: domain

An array containing the upper and lower limits of the domain for the ice shelf. The first index represents the dimension for which the boundaries apply. If the second index is 1 then it corresponds to the lower bound. If the second index is 2 then it corresponds to the upper bound.

integer, intent(in), dimension(:):: resolution

The number of data points in each dimension.

procedure(thickness_func) :: thickness

A function which calculates the initial value of the thickness of the ice shelf at a given location.

procedure(velocity_func) :: velocity

A function which calculates the initial value of the velocity (vector) of the ice at a given location in an ice shelf.

real(kind=r8), intent(in), optional :: temperature

The temperature of the ice in the ice shelf.

class(abstract_viscosity), intent(inout), optional allocatable:: viscosity_law

An object which calculates the viscosity of the ice. If not specified, then Glen's law will be used with $n=3$. Will be unallocated on return.

class(glacier_boundary), intent(inout), optional allocatable:: boundaries

An object specifying the boundary conditions for the ice shelf. Will be unallocated on return.

real(kind=r8), intent(in), optional :: lambda

The dimensionless ratio $\lambda \equiv \frac{\rho_0m_0x_0}{\rho_ih_0u_0}$.

real(kind=r8), intent(in), optional :: chi

The dimensionless ratio $\chi \equiv \frac{\rho_igh_0x_x}{2\eta_0u_0}\left(1 - \frac{\rho_i}{\rho_0}\right)$.

real(kind=r8), intent(in), optional :: zeta

The dimensionless ratio $\zeta \equiv \frac{\rho_iu_0x_0}{\eta_0}$, corresponding to the Reynolds number. Currently this is unused and always treated as 0.

real(kind=r8), intent(in), optional :: courant

The Courant number to use when calculating the time step. Defaults to 100. Too large a value will pose difficulties for the nonlinear solver, while too small a value can be numerically unstable. Typically, smaller values are needed for lower resolution.

real(kind=r8), intent(in), optional :: max_dt

The maximum allowable time step. This defaults to (effectively no maximum).

procedure(kappa_init_func), optional :: kappa

A function which specifies the initial values of the Taylor coefficients describing the vertical distribution of internal reflectors within the ice. The initial conditions at the grounding line will provide the boundary conditions there throughout the simulation. If this parameter is not provided then these layers will not be included in the integration. Both this parameter and n_kappa must be specified for the calculation to take place.

integer, intent(in), optional :: n_kappa

The number of Taylor coefficients used to describe internal reflectors. If not provided then these reflectors will not be included in the integration. Both this parameter and kappa must be specified for the calculation to take place.

private subroutine shelf_update(this, state_vector)

Author
Christopher MacMackin
Date
April 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
real(kind=r8), intent(in), dimension(:):: state_vector

A real array containing the data describing the state of the glacier.

private subroutine shelf_set_time(this, time)

Author
Christopher MacMackin
Date
November 2016

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

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
real(kind=r8), intent(in) :: time

The time at which the glacier is in the present state.

private subroutine shelf_read_data(this, file_id, group_name, error)

Author
Chris MacMackin
Date
April 2017

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
integer(kind=hid_t), intent(in) :: file_id

The identifier for the HDF5 file/group in which this data is meant to be written.

character(len=*), intent(in) :: group_name

The name to give the group in the HDF5 file storing the ice shelf's data.

integer, intent(out) :: error

Flag indicating whether routine ran without error. If no error occurs then has value 0.

private subroutine shelf_write_data(this, file_id, group_name, error)

Author
Chris MacMackin
Date
November 2016

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(in) :: this
integer(kind=hid_t), intent(in) :: file_id

The identifier for the HDF5 file/group in which this data is meant to be written.

character(len=*), intent(in) :: group_name

The name to give the group in the HDF5 file storing the ice shelf's data.

integer, intent(out) :: error

Flag indicating whether routine ran without error. If no error occurs then has value 0.

private subroutine shelf_solve_velocity(this, basal_drag, success)

Author
Chris MacMackin
Date
May 2017

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

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
class(scalar_field), intent(in) :: basal_drag

A paramter, e.g. coefficient of friction, needed to calculate the drag on basal surface of the glacier.

logical, intent(out) :: success

True if the integration is successful, false otherwise

private subroutine shelf_integrate(this, old_states, basal_melt, basal_drag, water_density, time, success)

Author
Chris MacMackin
Date
May 2017

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
class(glacier), intent(in), dimension(:):: old_states

Previous states of the glacier, with the most recent one first.

class(scalar_field), intent(in) :: basal_melt

The melt rate that the bottom of the glacier experiences during this time step.

class(scalar_field), intent(in) :: basal_drag

A paramter, e.g. coefficient of friction, needed to calculate the drag on basal surface of the glacier.

real(kind=r8), intent(in) :: water_density

The density of the water below the glacier.

real(kind=r8), intent(in) :: time

The time to which the glacier should be integrated

logical, intent(out) :: success

True if the integration is successful, false otherwise

private subroutine shelf_assign(this, rhs)

Author
Chris MacMackin
Date
February 2017

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.

Read more…

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(out) :: this
class(glacier), intent(in) :: rhs

The ice shelf to be assigned to this one.

private subroutine ice_shelf_integrate_layers(this, old_states, time, success)

Author
Chris MacMackin
Date
September 2018

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.

Arguments

Type IntentOptional AttributesName
class(ice_shelf), intent(inout) :: this
class(glacier), intent(in), dimension(:):: old_states

Previous states of the ice_shelf, with the most recent one first.

real(kind=r8), intent(in) :: time

The time to which the ice_shelf should be integrated

logical, intent(out) :: success

True if the integration is successful, false otherwise