fin_diff_block_solve_scalar Function

  function fin_diff_block_solve_scalar(this, rhs, offset) result(solution)
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! Solves the linear(ised) system represented by this finite
    ! difference block, for a given right hand side state vector
    ! (represented by a scalar field). Optionally, the differential
    ! operator can be augmented by adding an offset, i.e. a scalar
    ! field which is added to the operator.
    !
    ! @Warning Currently this is only implemented for a 1-D field.
    !
    class(fin_diff_block), intent(inout)      :: this
    class(scalar_field), intent(in)           :: rhs
      !! The right hand side of the linear(ised) system.
    class(scalar_field), optional, intent(in) :: offset
      !! An offset to add to the differential operator
    class(scalar_field), pointer              :: solution

    real(r8), dimension(:), allocatable :: sol_vector, diag_vector
    integer                             :: flag, n, i, j
    real(r8)                            :: forward_err, &
                                           backward_err, &
                                           condition_num
    character(len=1)                    :: factor
    character(len=:), allocatable       :: msg
 
    call rhs%guard_temp()
    if (present(offset)) call offset%guard_temp()
    allocate(sol_vector(rhs%raw_size()))
    ! Allocate the arrays used to hold the factorisation of the
    ! tridiagonal matrix
    n = size(this%diagonal)
    if (.not. allocated(this%pivots)) then
      allocate(this%l_multipliers(n-1))
      allocate(this%u_diagonal(n))
      allocate(this%u_superdiagonal1(n-1))
      allocate(this%u_superdiagonal2(n-2))
      allocate(this%pivots(n))
      factor = 'N'
    else
      if (.not. present(offset) .and. .not. this%had_offset) then
        factor = 'F'
      else
        factor = 'N'
      end if
    end if

    if (present(offset)) then
#ifdef DEBUG
      if (offset%elements() /= n) then
        call logger%fatal('fin_diff_block%solve_for','Offset field has '// &
                          'different resolution than finite difference block.')
      end if
#endif
      allocate(diag_vector(n))

      where (this%diagonal == 0._r8)
        diag_vector = this%diagonal + offset%raw()
      elsewhere
        diag_vector = this%diagonal
      end where
      if (.not. any(1 == this%boundary_locs)) then
        diag_vector(1) = diag_vector(1) + offset%get_element(1)
      end if
      if (.not. any(n == this%boundary_locs)) then
        diag_vector(n) = diag_vector(n) + offset%get_element(n)
      end if
      do i=1, size(this%boundary_locs)
        j = this%boundary_locs(i)
        if ((j == 1 .or. j == n) .and. this%boundary_types(i) == free_boundary) then
          diag_vector(j) = diag_vector(j) + offset%get_element(j)
        end if
      end do
      call la_gtsvx(this%sub_diagonal, diag_vector, this%super_diagonal, &
                    rhs%raw(), sol_vector, this%l_multipliers,           &
                    this%u_diagonal, this%u_superdiagonal1,              &
                    this%u_superdiagonal2, this%pivots, factor, 'N',     &
                    forward_err, backward_err, condition_num, flag)
    else
      call la_gtsvx(this%sub_diagonal, this%diagonal, this%super_diagonal,      &
                    rhs%raw(), sol_vector, this%l_multipliers, this%u_diagonal, &
                    this%u_superdiagonal1, this%u_superdiagonal2, this%pivots,  &
                    factor, 'N', forward_err, backward_err, condition_num,      &
                    flag)
    end if
    if (flag/=0) then
      msg = 'Tridiagonal matrix solver returned with flag '//str(flag)
      call logger%error('fin_diff_block%solve_for',msg)
    else
      this%had_offset = present(offset)
#ifdef DEBUG
      msg = 'Tridiagonal matrix solver retunred with estimated condition '// &
            'number '//str(condition_num)
      call logger%debug('fin_diff_block%solve_for',msg)
#endif
    end if
    call rhs%allocate_scalar_field(solution)
    call solution%unset_temp()
    call solution%assign_meta_data(rhs)
    call solution%set_from_raw(sol_vector)
    if (present(offset)) call offset%clean_temp()
    call rhs%clean_temp(); call solution%set_temp()
  end function fin_diff_block_solve_scalar