jacobian_block_solve Function

  function jacobian_block_solve(this, rhs) result(solution)
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! Solves the linear(ised) system represented by this Jacobian
    ! block, for a given right hand side state vector (represented by
    ! a scalar field).
    !
    ! @Warning Currently this is only implemented for a 1-D field. The
    ! extra derivative must be in the direction along which the field
    ! varies.
    !
    class(jacobian_block), intent(inout) :: this
    class(scalar_field), intent(in)      :: rhs
      !! The right hand side of the linear(ised) system.
    class(scalar_field), pointer         :: solution
    real(r8), dimension(:), allocatable  :: sol_vector

    integer                                   :: flag, n
    real(r8)                                  :: forward_err, &
                                                 backward_err, &
                                                 condition_num
    real(r8), dimension(:), allocatable       :: diag, subdiag, supdiag
    character(len=1)                          :: factor
    character(len=:), allocatable             :: msg

    call rhs%guard_temp()
    allocate(sol_vector(rhs%raw_size()))
    ! Construct tridiagonal matrix for this operation, if necessary
    if (.not. allocated(this%pivots)) then
      factor = 'N'
      ! Allocate the arrays used to hold the factorisation of the
      ! tridiagonal matrix
      call this%get_tridiag(diag, subdiag, supdiag)
      call move_alloc(diag, this%diagonal)
      call move_alloc(subdiag, this%sub_diagonal)
      call move_alloc(supdiag, this%super_diagonal)
      n = this%contents%raw_size()
      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))
    else
      factor = 'F'
    end if
!    print*,rhs%raw()
    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)
    if (flag/=0) then
      msg = 'Tridiagonal matrix solver returned with flag '//trim(str(flag))
!      print*,this%sub_diagonal
!      print*,
!      print*,this%diagonal
!      print*,
!      print*,this%super_diagonal
      CALL backtrace
      call logger%error('jacobian_block%solve_for',msg)
      !print*,condition_num
!      stop
#ifdef DEBUG
    else
      msg = 'Tridiagonal matrix solver returned with estimated condition '// &
            'number '//trim(str(condition_num))
      call logger%debug('jacobian_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)
    call rhs%clean_temp(); call solution%set_temp()
  end function jacobian_block_solve