jacobian_block_multiply Function

  recursive function jacobian_block_multiply(this, rhs) result(product)
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! Provides a matrix multiplication operator between a Jacobian
    ! block and a scalar field (which corresponds to a state vector).
    !
    class(jacobian_block), intent(in) :: this
    class(scalar_field), intent(in)   :: rhs
      !! A field corresponding to a state vector being multiplied by
      !! the Jacobian block.
    class(scalar_field), pointer :: product
    class(scalar_field), pointer :: tmp
    real(r8), dimension(:), allocatable :: bounds
    integer :: i
    call rhs%guard_temp()
    call this%contents%allocate_scalar_field(product)
    call product%unset_temp()
    if (this%extra_derivative==no_extra_derivative) then
      select case(this%has_increment)
      case(0)
        product = this%coef*this%derivative * rhs + this%contents * rhs%d_dx(this%direction)
      case(1)
        product = (this%coef*this%derivative + this%real_increment) * rhs &
                + this%contents * rhs%d_dx(this%direction)
      case(2)
        product = (this%coef*this%derivative + this%field_increment) * rhs &
                + this%contents * rhs%d_dx(this%direction)
      case default
        error stop ('Invalid increment has been added.')
      end select
    else
      call rhs%allocate_scalar_field(tmp)
      call tmp%guard_temp()
      tmp = rhs%d_dx(this%extra_derivative)
      select case(this%has_increment)
      case(0)
        product = this%coef*this%derivative * tmp + this%contents * tmp%d_dx(this%direction)
      case(1)
        product = this%coef*this%derivative * tmp + this%contents * tmp%d_dx(this%direction) &
                + this%real_increment*rhs
      case(2)
        product = this%coef*this%derivative * tmp + this%contents * tmp%d_dx(this%direction) &
                + this%field_increment*rhs
      case default
        error stop ('Invalid increment has been added.')
      end select
      call tmp%clean_temp()
    end if
    if (associated(this%block_increment)) then
      product = product + this%block_increment*rhs
    end if
    call this%get_boundaries(this%contents,this%derivative,rhs,      &
                            this%boundary_locs,this%boundary_types, &
                            bounds)
    do i = 1, size(this%boundary_locs)
      call product%set_element(this%boundary_locs(i), bounds(i))
    end do
    call rhs%clean_temp()
#ifdef DEBUG
    call logger%debug('jacobian_block%multiply','Multiplied vector by '// &
                      'Jacobian block.')
#endif
    call product%set_temp()
  end function jacobian_block_multiply