shelf_precondition Function

  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.
    !
    class(ice_shelf), intent(inout)          :: this
    class(glacier), dimension(:), intent(in) :: 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(r8), intent(in)                     :: water_density
      !! The density of the water below the glacier
    real(r8), dimension(:), intent(in)       :: delta_state
      !! The change to the state vector which is being preconditioned.
    real(r8), dimension(:), allocatable      :: preconditioned
      !! The result of applying the preconditioner to `delta_state`.

    type(cheb1d_scalar_field) :: delta_h
    integer :: i, sl, el, su, eu
    integer, dimension(:), allocatable :: boundary_locations
    real(r8) :: delta_t

    call melt_rate%guard_temp(); call basal_drag_parameter%guard_temp()
    allocate(preconditioned(size(delta_state)))
    select type(previous_states)
    class is(ice_shelf)
      delta_t = this%time - previous_states(1)%time
    class default
      call logger%fatal('ice_shelf%precondition','Type other than `ice_shelf` '// &
                        'passed to `ice_shelf` object as a previous state.')
      error stop
    end select
    call delta_h%assign_meta_data(this%thickness)
    call delta_h%set_from_raw(delta_state)

    associate(jac => this%thickness_jacobian, uvec => this%velocity)
      if (this%stale_jacobian) then
        sl = this%thickness_size - this%thickness_lower_bound_size + 1
        el = this%thickness_size
        su = 1
        eu = this%thickness_upper_bound_size
        boundary_locations = [(i, i=sl,el), (i, i=su,eu)]

        jac = jacobian_block(uvec%component(1), 1, boundary_locs=boundary_locations) &
            + 1._r8/delta_t

        this%stale_jacobian = .false.
      end if

      delta_h = jac%solve_for(delta_h)
    end associate
    
    preconditioned = delta_h%raw()
    call melt_rate%clean_temp(); call basal_drag_parameter%clean_temp()
#ifdef DEBUG
    call logger%debug('ice_shelf%precondition','Applied preconditioner for ice shelf.')
#endif
  end function shelf_precondition