shelf_residual Function

  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.
    !
    class(ice_shelf), intent(in)             :: 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 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(r8), intent(in)                     :: water_density
      !! The density of the water below the glacier.
    real(r8), dimension(:), allocatable      :: residual
      !! The residual of the system of equations describing the glacier.
    type(cheb1d_scalar_field) :: scalar_tmp
    integer :: start, finish, bounds_start, bounds_finish
    integer, dimension(:), allocatable :: lower, upper
    real(r8), dimension(:), allocatable :: bounds
    logical :: success

    call melt_rate%guard_temp(); call basal_drag_parameter%guard_temp()
    allocate(residual(this%data_size()))
    start = 1

    ! Use same or similar notation for variables as in equations
    select type(previous_states)
    class is(ice_shelf)
      associate(h => this%thickness, h_old => previous_states(1)%thickness, &
                uvec => this%velocity, m => melt_rate, eta => this%eta,     &
                lambda => this%lambda, t_old => previous_states(1)%time)
        ! Boundary conditions
        if (this%stale_eta) then
          bounds = this%boundaries%boundary_residuals(h, uvec, &
               this%viscosity_law%ice_viscosity(uvec, this%ice_temperature(), &
                                               this%time), this%time)
        else
          bounds = this%boundaries%boundary_residuals(h, uvec, eta, this%time)
        end if

        ! Continuity equation
        scalar_tmp = (h - h_old)/(this%time - t_old) + .div.(h*uvec) + lambda*m
        
        lower = this%boundaries%thickness_lower_bound()
        upper = this%boundaries%thickness_upper_bound()
        ! TODO: Figure out how to make this independent of order which
        ! values are stored in the field
        finish = start + this%thickness_upper_bound_size - 1
        bounds_start = this%thickness_lower_bound_size + 1
        bounds_finish = this%thickness_lower_bound_size &
                      + this%thickness_upper_bound_size
        residual(start:finish) = bounds(bounds_start:bounds_finish)
        start = finish + 1
        finish = start + scalar_tmp%raw_size(lower,upper) - 1
        residual(start:finish) = scalar_tmp%raw(lower,upper)
        start = finish + 1
        finish = start + this%thickness_lower_bound_size - 1
        bounds_start = 1
        bounds_finish = this%thickness_lower_bound_size
        residual(start:finish) = bounds(bounds_start:bounds_finish)
        start = finish + 1

      end associate
    class default
      call logger%fatal('ice_shelf%residual','Type other than `ice_shelf` '// &
                        'passed to `ice_shelf` object as a previous state.')
      error stop
    end select

    call melt_rate%clean_temp(); call basal_drag_parameter%clean_temp()
#ifdef DEBUG
    call logger%debug('ice_shelf%residual','Calculated residual of ice shelf.')
#endif
  end function shelf_residual