average_kp_matrices Subroutine

    subroutine average_kp_matrices(nrepeats, matrices, mean, se)

        integer, intent(in) :: nrepeats
        real(dp), intent(in) :: matrices(:, :, :)
        real(dp), intent(out) :: mean(:, :), se(:, :)

        integer :: irepeat

        mean = 0.0_dp
        se = 0.0_dp

        do irepeat = 1, nrepeats
            mean = mean + matrices(:, :, irepeat)
        end do
        mean = mean / nrepeats

        if (nrepeats > 1) then
            do irepeat = 1, nrepeats
                se = se + (matrices(:, :, irepeat) - mean)**2
            end do
            se = se / ((nrepeats - 1) * nrepeats)
            se = sqrt(se)
        end if

    end subroutine average_kp_matrices