detrace_tensor Subroutine

pure subroutine detrace_tensor(A)

    real(dp), dimension(:), intent(inout) :: A

    integer :: rank
    integer :: norm
    integer :: ti
    integer :: ax, ay, az
    integer :: bx, by, bz
    real(dp) :: xn, yn, zn
    real(dp) :: sgn
    real(dp) :: divisor
    real(dp), dimension(:), allocatable :: trace_tensor

    rank = tensor_rank(A)

    divisor = double_factorial(2 * rank - 1)

    allocate(trace_tensor(size(A)))

    trace_tensor = 0.0_dp

    do ax = rank, 0, -1
        do ay = rank - ax, 0, -1
            az = rank - ax - ay
            ti = convert_index(ax, ay, az)
            do bx = 0, floor(real(ax)/2.0)
                xn = some_nomial(ax, bx)
                do by = 0, floor(real(ay)/2.0)
                    yn = some_nomial(ay, by)
                    do bz = 0, floor(real(az)/2.0)
                        zn = some_nomial(az, bz)
                        norm = bx + by + bz
                        if (norm /= 1) then
                            cycle
                        end if
                        sgn = (-1.0_dp)**(norm)
                        trace_tensor(ti) = trace_tensor(ti) &
                                           + sgn * double_factorial(2*rank-2*norm-1) * xn * yn * zn &
                                             * tensor_trace(A, ax-2*bx, ay-2*by, az-2*bz) / divisor
                    end do
                end do
            end do
        end do
    end do

    A = A + trace_tensor

    deallocate(trace_tensor)

    contains

    real(dp) pure function some_nomial(ai, bi)

        integer, intent(in) :: ai
        integer, intent(in) :: bi

        some_nomial = factorial(ai) / (2.0_dp**(bi) * factorial(bi) * factorial(ai-2*bi))

    end function some_nomial

end subroutine detrace_tensor