tensor_trace Function

real(dp) pure function tensor_trace(A, ax, ay, az, order)

    real(dp), dimension(:), intent(in) :: A
    integer, intent(in) :: ax
    integer, intent(in) :: ay
    integer, intent(in) :: az
    integer, intent(in), optional :: order

    integer :: ti
    integer :: bx, by, bz
    integer :: trace_order

    if (present(order)) then
        trace_order = order
    else
        trace_order = 1
    end if

    tensor_trace = 0.0_dp

    do bx = trace_order, 0, -1
        do by = trace_order - bx, 0, -1
            bz = trace_order - bx - by
            ti = convert_index(ax+2*bx, ay+2*by, az+2*bz)
            tensor_trace = tensor_trace + trinomial(bx, by, bz) * A(ti)
        end do
    end do

end function tensor_trace