upstream_calculate Subroutine

  subroutine upstream_calculate(this, t, func, b)
    !* Author: Chris MacMackin
    !  Date: July 2017
    !
    ! Calculates the boundary values to use at the current time with
    ! the current ice thickness.
    !
    class(upstream_plume_boundary), intent(inout) :: this
    real(r8), intent(in)                          :: t
      !! The time at which to calculate the boundary values.
    procedure(non_diff)                           :: func
      !! A function which returns the non-diffusive, non-inertial
      !! components of the ODEs describing the plume.
    class(scalar_field), intent(in)               :: b
      !! The depth of the ice shelf base.

    integer :: n
    real(r8) :: D_0, U_0, S_0, T_0, DU_0, DUS_0, DUT_0
    real(r8), dimension(:), allocatable :: Uvec_0, DUvecU_0
    class(scalar_field), pointer :: b_x
    real(r8), dimension(:), allocatable :: y0, ygot, yderiv_got
    real(r8) :: xgot
    integer :: flag
    type(rk_comm_real_1d) :: comm

    call b%guard_temp()
    call this%get_boundaries(t, D_0, Uvec_0, T_0, S_0)

    U_0 = Uvec_0(1)
    DU_0 = D_0*U_0
    DUvecU_0 = DU_0*Uvec_0
    DUT_0 = DU_0*T_0
    DUS_0 = DU_0*S_0
    y0 = [DU_0, DUvecU_0, DUT_0, DUS_0]
    n = size(y0)
    allocate(ygot(n), yderiv_got(n))
    if (.not. allocated(this%thresholds)) then
      allocate(this%thresholds(n))
      this%thresholds = 1._r8
    end if
    b_x => b%d_dx(1)
    call b_x%guard_temp()

    call setup(comm, 0._r8, y0, this%distance, 1e-6_r8, this%thresholds, &
               method='h', h_start=2e-3_r8*this%distance)
    call range_integrate(comm, integrand, this%distance, xgot, ygot, &
                         yderiv_got, flag)
    if (flag /= 1) then
      call logger%error('upstream_plume_boundary%calculate',            &
                        'Could only integrate plume to x = '//str(xgot))
    end if
    !this%thresholds = abs(ygot) + 1e-10
    this%thickness = ygot(1)**2/ygot(2)
    this%velocity = ygot(2:n-2)/ygot(1)
    this%temperature = ygot(n-1)/ygot(1)
    this%salinity = ygot(n)/ygot(1)
    call logger%trivia('upstream_plume_boundary%calculate',             &
         'Calculated boundary values D='//str(this%thickness)//', U='// &
         str(this%velocity)//', T='//str(this%temperature)//', S='//    &
         str(this%salinity))
    if (this%thickness < 0._r8) error stop ('Negative plume thickness found')
    call collect_garbage(comm)
    call b%clean_temp(); call b_x%clean_temp()

  contains

    function integrand(x, y) result(f)
      real(r8), intent(in)               :: x
        !! The independent variable
      real(r8), dimension(:), intent(in) :: y
        !! The dependent variable
      real(r8), dimension(size(y))       :: f
        !! The derivatives

      type(uniform_scalar_field)   :: D, T, S, DU_x, DUT_x, DUS_x
      type(uniform_vector_field)   :: Uvec, DUU_x
      type(uniform_gradient_field) :: b_loc
      D = uniform_scalar_field(y(1)**2/y(2))
      Uvec = uniform_vector_field(y(2:n-2)/y(1))
      T = uniform_scalar_field(y(n-1)/y(1))
      S = uniform_scalar_field(y(n)/y(1))
      b_loc = uniform_gradient_field(b%interpolate([x]), [b_x%interpolate([x])])
      call func(D, Uvec, T, S, b_loc, DU_x, DUU_x, DUT_x, DUS_x)
      f(1) = DU_x%get_value()
      f(2:n-2) = DUU_x%get_value()
      f(n-1) = DUT_x%get_value()
      f(n) = DUS_x%get_value()
    end function integrand
    
  end subroutine upstream_calculate