shelf_initialise Subroutine

  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.
    !
    class(ice_shelf), intent(out)        :: this
    real(r8), dimension(:,:), intent(in) :: 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, dimension(:), intent(in)    :: 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(r8), intent(in), optional       :: temperature
      !! The temperature of the ice in the ice shelf.
    class(abstract_viscosity), allocatable, optional, &
                           intent(inout) :: 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), allocatable, optional, &
                         intent(inout)   :: boundaries
      !! An object specifying the boundary conditions for the ice
      !! shelf. Will be unallocated on return.
    real(r8), intent(in), optional       :: lambda
      !! The dimensionless ratio 
      !! $\lambda \equiv \frac{\rho_0m_0x_0}{\rho_ih_0u_0}$.
    real(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(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(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(r8), intent(in), optional       :: max_dt
      !! The maximum allowable time step. This defaults to \(1\times
      !! 10^{99}\) (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, optional, intent(in)        :: 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.

    integer :: n
    integer, dimension(:), allocatable :: lower, upper

    this%thickness = cheb1d_scalar_field(resolution(1),thickness,domain(1,1), &
                                         domain(1,2))
    this%velocity = cheb1d_vector_field(resolution(1),velocity,domain(1,1), &
                                        domain(1,2))
    this%eta = cheb1d_scalar_field(resolution(1), lower_bound=domain(1,1), &
                                   upper_bound=domain(1,2))
    this%thickness_size = this%thickness%raw_size()
    this%velocity_size = this%velocity%raw_size()

    if (present(kappa) .and. present(n_kappa)) then
      allocate(this%kappa(n_kappa))
      do n = 1, n_kappa
        this%kappa(n) = cheb1d_scalar_field(resolution(1),kappa_func,domain(1,1), &
                                            domain(1,2))
      end do
    end if

    if (present(temperature)) then
      continue ! This doesn't do anything at the moment
    else
      continue ! Again, doesn't do anything at the moment
    end if

    if (present(viscosity_law)) then
      call move_alloc(viscosity_law, this%viscosity_law)
    else
      allocate(newtonian_viscosity:: this%viscosity_law)
    end if

    if (present(boundaries)) then
      call move_alloc(boundaries, this%boundaries)
    else
      allocate(dallaston2015_glacier_boundary :: this%boundaries)
    end if

    if (present(lambda)) then
      this%lambda = lambda
    else
      this%lambda = 0.37_r8
    end if

    if (present(chi)) then
      this%chi = chi
    else
      this%chi = 4.0_r8
    end if

    if (present(zeta)) then
      this%zeta = zeta
    else
      this%zeta = 1.3e-11_r8
    end if

    if (present(courant)) then
      this%courant = courant
    else
      this%courant = 1e2_r8
    end if

    if (present(max_dt)) then
      this%max_dt = max_dt
    else
      this%max_dt = 1e99_r8
    end if

    lower = this%boundaries%thickness_lower_bound()
    upper = this%boundaries%thickness_upper_bound()
    this%thickness_lower_bound_size = this%thickness_size &
                                    - this%thickness%raw_size(lower)
    this%thickness_upper_bound_size = this%thickness_size &
                                    - this%thickness%raw_size(upper)

    lower = this%boundaries%velocity_lower_bound()
    upper = this%boundaries%velocity_upper_bound()
    this%velocity_lower_bound_size = this%velocity_size &
                                   - this%velocity%raw_size(lower)
    this%velocity_upper_bound_size = this%velocity_size &
                                   - this%velocity%raw_size(upper)

    this%boundary_start = this%thickness_size + this%velocity_size + 1 &
                        - this%thickness_lower_bound_size &
                        - this%thickness_upper_bound_size &
                        - this%velocity_lower_bound_size &
                        - this%velocity_upper_bound_size

    this%time = 0.0_r8
    this%stale_eta = .true.
    this%stale_jacobian = .true.
#ifdef DEBUG
    call logger%debug('ice_shelf','Initialised new ice shelf object')
#endif

  contains

    pure function kappa_func(location) result(thickness)
      !* Author: Chris MacMackin
      !  Date: April 2016
      !
      ! Wrapper around user-provided `kappa` routine, converting it to
      ! the form needed to initialise field-types.
      !
      real(r8), dimension(:), intent(in) :: location
        !! The position $\vec{x}$ at which to compute the thickness
      real(r8) :: thickness
        !! The thickness of the glacier at `location`
      thickness = kappa(n, location)
    end function kappa_func

  end subroutine shelf_initialise