seasonal_residuals Function

  function seasonal_residuals(this, thickness, velocity, viscosity, t) &
                                   result(residuals)
    !* Author: Chris MacMackin
    !  Date: October 2017
    !
    ! Returns the difference between the glacier conditions of the
    ! plume and the Dirichlet conditions prescribed in the model of
    ! Dallaston et al. (2015)
    !
    class(seasonal_glacier_boundary), intent(in) :: this
    class(scalar_field), intent(in)     :: thickness
      !! A field containing the thickness of the glacier
    class(vector_field), intent(in)     :: velocity
      !! A field containing the flow velocity of the glacier
    class(scalar_field), intent(in)     :: viscosity
      !! A field containing the viscosity of the ice in the glacier.
    real(r8), intent(in)                :: t
      !! The time at which the boundary conditions are to be
      !! calculated.
    real(r8), allocatable, dimension(:) :: residuals
      !! An array containing the difference between the required
      !! boundary values and those which are actually present. They
      !! are stored in the order: lower thickness boundary, upper
      !! thickness boundary, lower velocity boundary, and upper
      !! velocity boundary.
    real(r8) :: vel, tm
    class(scalar_field), pointer :: thickness_bound,      &
                                    velocity_bound_upper, &
                                    velocity_deriv
    class(vector_field), pointer :: velocity_bound_lower
    call thickness%guard_temp(); call velocity%guard_temp(); call viscosity%guard_temp()
    call thickness%allocate_scalar_field(thickness_bound)
    thickness_bound = thickness%get_boundary(-1,1) - this%thickness
    call velocity%allocate_vector_field(velocity_bound_lower)
    if (this%square) then
      vel = this%mean + sign(this%amplitude, sin(this%frequency*t))
    else
      vel = this%mean + this%amplitude*sin(this%frequency*t)
    end if
    velocity_bound_lower = velocity%get_boundary(-1,1) - [vel]
    call velocity%allocate_scalar_field(velocity_deriv)
    velocity_deriv = velocity%component_d_dx(1,1)
    call velocity%allocate_scalar_field(velocity_bound_upper)
    velocity_bound_upper = velocity_deriv%get_boundary(1,1)  &
                         * thickness%get_boundary(1,1)       &
                         - (0.25_r8*this%chi)*thickness%get_boundary(1,1)**2 &
                         / viscosity%get_boundary(1,1)
    residuals = [thickness_bound%raw(), velocity_bound_lower%raw(), &
                 velocity_bound_upper%raw()]
    call thickness%clean_temp(); call velocity%clean_temp(); call viscosity%clean_temp()
  end function seasonal_residuals