ice_sheet – ISOFT

type, public, extends(glacier) :: ice_sheet

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

Inherits

Help

Components

Type	Visibility	Attributes		Name
type(cheb1d_scalar_field),	private		::	thickness	Thickness of ice sheet, $h$
type(cheb1d_vector_field),	private		::	velocity	Flow velocity of ice sheet, $\vec{u}$
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 $\chi \equiv \frac{\rho_igh_0x_x}{2\eta_0u_0}$
class(abstract_viscosity),	private,	allocatable	::	viscosity_law	An object representing the model used for ice viscosity.
real(kind=r8),	private		::	time	The time at which the ice sheet is in this state

Constructor

public interface ice_sheet

private function constructor(domain, resolution, thickness, velocity, viscosity_law, lambda, chi) result(this)

Author: Christopher MacMackin
Date: April 2016

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.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=r8),	intent(in),		dimension(:,:)	::	domain	An array containing the upper and lower limits of the domain for the ice sheet. 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 sheet 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 sheet.
class(abstract_viscosity),	intent(in),	optional		::	viscosity_law	An object which calculates the viscosity of the ice.
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}$

Return Value type(ice_sheet)

An ice sheet object with its domain and initial conditions set according to the arguments of the constructor function.

Type-Bound Procedures

generic, public :: assignment(=) => assign

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.

Arguments

Type	Intent	Optional	Attributes		Name
class(ice_shelf),	intent(out)			::	this
class(glacier),	intent(in)			::	rhs	The ice shelf to be assigned to this one.

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.

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

Author: Chris MacMackin
Date: November 2016

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.

Arguments

Type	Intent	Attributes		Name
class(glacier),	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

procedure, public :: integrate_layers => glacier_integrate_layers

Dummy routine which can be over-ridden to integrate internal layers of the glacier to the specified time.

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

Author: Chris MacMackin
Date: September 2018

Dummy routine which does nothing.

Arguments

Type	Intent	Attributes		Name
class(glacier),	intent(inout)		::	this
class(glacier),	intent(in),	dimension(:)	::	old_states	Previous states of the glacier, with the most recent one first.
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

procedure, public :: ice_thickness => sheet_thickness

private pure function sheet_thickness(this) result(thickness)

Author
Christopher MacMackin

Date
April 2016

Returns the thickness of the ice sheet across its domain.

Arguments

Type Intent Optional Attributes Name

class(ice_sheet), intent(in) :: this

Return Value class(scalar_field), pointer

The ice thickness.

procedure, public :: ice_density => sheet_density

private pure function sheet_density(this) result(density)

Author
Christopher MacMackin

Date
April 2016

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

Arguments

Type Intent Optional Attributes Name

class(ice_sheet), intent(in) :: this

Return Value real(kind=r8)

The ice density.

procedure, public :: ice_temperature => sheet_temperature

private pure function sheet_temperature(this) result(temperature)

Author
Christopher MacMackin

Date
April 2016

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

Arguments

Type Intent Optional Attributes Name

class(ice_sheet), intent(in) :: this

Return Value real(kind=r8)

The ice density.

procedure, public :: residual => sheet_residual

private function sheet_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	Intent	Attributes		Name
class(ice_sheet),	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 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.

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

The residual of the system of equations describing the glacier.

procedure, public :: update => sheet_update

private subroutine sheet_update(this, state_vector)

Author: Christopher MacMackin
Date: April 2016

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.

Arguments

Type	Intent	Optional	Attributes		Name
class(ice_sheet),	intent(inout)			::	this
real(kind=r8),	intent(in),		dimension(:)	::	state_vector	A real array containing the data describing the state of the glacier.

procedure, public :: precondition => sheet_precondition

private function sheet_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	Intent	Attributes		Name
class(ice_sheet),	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.

procedure, public :: set_time => sheet_set_time

private subroutine sheet_set_time(this, time)

Author: Christopher MacMackin
Date: November 2016

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

Arguments

Type	Intent	Optional	Attributes		Name
class(ice_sheet),	intent(inout)			::	this
real(kind=r8),	intent(in)			::	time	The time at which the glacier is in the present state.

procedure, public :: data_size => sheet_data_size

private pure function sheet_data_size(this)

Author
Christopher MacMackin

Date
August 2016

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.

Arguments

Type Intent Optional Attributes Name

class(ice_sheet), intent(in) :: this

Return Value integer

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

procedure, public :: state_vector => sheet_state_vector

private pure function sheet_state_vector(this) result(state_vector)

Author
Christopher MacMackin

Date
April 2016

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

Arguments

Type Intent Optional Attributes Name

class(ice_sheet), intent(in) :: this

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

The state vector describing the ice sheet.

procedure, public :: solve_velocity => sheet_solve_velocity

private subroutine sheet_solve_velocity(this, basal_drag, success)

Author: Chris MacMackin
Date: May 2017

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

Arguments

Type	Intent		Name
class(ice_sheet),	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

procedure, public :: read_data => sheet_read_data

private subroutine sheet_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	Intent		Name
class(ice_sheet),	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.

procedure, public :: write_data => sheet_write_data

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

Author: Chris MacMackin
Date: November 2016

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.

Arguments

Type	Intent		Name
class(ice_sheet),	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 sheet's data.
integer,	intent(out)	::	error	Flag indicating whether routine ran without error. If no error occurs then has value 0.

procedure, public :: time_step => sheet_time_step

private function sheet_time_step(this) result(dt)

Author
Chris MacMackin

Date
December 2016

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

Arguments

Type Intent Optional Attributes Name

class(ice_sheet), intent(in) :: this

Return Value real(kind=r8)

The time-step to use

procedure, private :: assign => sheet_assign

private subroutine sheet_assign(this, rhs)

Author: Chris MacMackin
Date: February 2017

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.

Arguments

Type	Intent	Optional	Attributes		Name
class(ice_sheet),	intent(out)			::	this
class(glacier),	intent(in)			::	rhs	The ice sheet to be assigned to this one.

Source Code

  type, extends(glacier), public :: ice_sheet
    !* Author: Chris MacMackin
    !  Date: April 2016
    !
    ! A concrete implementation of the [[glacier]] type, using a vertically
    ! integrated model of an ice sheet. This model is 1-dimensional only.
    !
    private
    type(cheb1d_scalar_field) :: thickness
      !! Thickness of ice sheet, $h$
    type(cheb1d_vector_field) :: velocity
      !! Flow velocity of ice sheet, $\vec{u}$
    real(r8)                  :: lambda
      !! The dimensionless ratio 
      !! $\lambda \equiv \frac{\rho_0m_0x_0}{\rho_iH-0u_0}$
    real(r8)                  :: chi
      !! The dimensionless ratio
      !! $\chi \equiv \frac{\rho_igh_0x_x}{2\eta_0u_0}$
    class(abstract_viscosity), allocatable :: viscosity_law
      !! An object representing the model used for ice viscosity.
    real(r8)                  :: time
      !! The time at which the ice sheet is in this state
  contains
!$    procedure            :: t => sheet_dt
!$    procedure            :: local_error => sheet_local_error
!$    procedure            :: integrand_multiply_integrand => sheet_m_sheet
!$    procedure            :: integrand_multiply_real => sheet_m_real
!$    procedure, pass(rhs) :: real_multiply_integrand => real_m_sheet
!$    procedure            :: add => sheet_add
!$    procedure            :: sub => sheet_sub
!$    procedure            :: assign_integrand => sheet_assign
    procedure :: ice_thickness => sheet_thickness
!$    procedure :: ice_velocity => sheet_velocity
    procedure :: ice_density => sheet_density
    procedure :: ice_temperature => sheet_temperature
    procedure :: residual => sheet_residual
    procedure :: update => sheet_update
    procedure :: precondition =>  sheet_precondition
    procedure :: set_time => sheet_set_time
    procedure :: data_size => sheet_data_size
    procedure :: state_vector => sheet_state_vector
    procedure :: solve_velocity => sheet_solve_velocity
    procedure :: read_data => sheet_read_data
    procedure :: write_data => sheet_write_data
    procedure :: time_step => sheet_time_step
    procedure, private :: assign => sheet_assign
  end type ice_sheet

ice_sheet Derived Type

Contents

Variables

Constructor

Type-Bound Procedures

Source Code

type, public, extends(glacier) :: ice_sheet

Inherits

Contents

Variables

Constructor

Type-Bound Procedures

Source Code

Components

Constructor

public interface ice_sheet

private function constructor(domain, resolution, thickness, velocity, viscosity_law, lambda, chi) result(this)

Arguments

Return Value type(ice_sheet)

Type-Bound Procedures

generic, public :: assignment(=) => assign

private subroutine shelf_assign(this, rhs)

Arguments

procedure, public :: integrate => glacier_integrate

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

Arguments

procedure, public :: integrate_layers => glacier_integrate_layers

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

Arguments

procedure, public :: ice_thickness => sheet_thickness

private pure function sheet_thickness(this) result(thickness)

Arguments

Return Value class(scalar_field), pointer

procedure, public :: ice_density => sheet_density

private pure function sheet_density(this) result(density)

Arguments

Return Value real(kind=r8)

procedure, public :: ice_temperature => sheet_temperature

private pure function sheet_temperature(this) result(temperature)

Arguments

Return Value real(kind=r8)

procedure, public :: residual => sheet_residual

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

Arguments

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

procedure, public :: update => sheet_update

private subroutine sheet_update(this, state_vector)

Arguments

procedure, public :: precondition => sheet_precondition

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

Arguments

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

procedure, public :: set_time => sheet_set_time

private subroutine sheet_set_time(this, time)

Arguments

procedure, public :: data_size => sheet_data_size

private pure function sheet_data_size(this)

Arguments

Return Value integer

procedure, public :: state_vector => sheet_state_vector

private pure function sheet_state_vector(this) result(state_vector)

Arguments

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

procedure, public :: solve_velocity => sheet_solve_velocity

private subroutine sheet_solve_velocity(this, basal_drag, success)

Arguments

procedure, public :: read_data => sheet_read_data

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

Arguments

procedure, public :: write_data => sheet_write_data

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

Arguments

procedure, public :: time_step => sheet_time_step

private function sheet_time_step(this) result(dt)

Arguments

Return Value real(kind=r8)

procedure, private :: assign => sheet_assign

private subroutine sheet_assign(this, rhs)

Arguments

Source Code