The evaluation of the Schur function $s_{\lambda}$ at the $n$th-roots of unity is the coefficient of $s_{\lambda}$ in the expansion in the Schur basis of the plethysm $h_k[p_n]=p_n[h_k]$ of the complete sum $h_k$ with the power sum $p_n$, when $k\cdot n=|\lambda|$ (and $0$ if $|\lambda|$ is not a multiple of $n$).
Indeed, consider the Cauchy identity,
$$
\prod_{i,j} \frac{1}{1-x_i y_j}=\sum_{\lambda} s_{\lambda}[X] s_{\lambda}[Y].
$$
Take for the $x_i$ the $n$th-roots of unity $\zeta_i$. Therefore the evaluation of $s_{\lambda}$ that you are looking for will be the coefficient of $s_{\lambda}[Y]$ in the expansion in the Schur basis of
$$
\prod_{i,j} \frac{1}{1-\zeta_i y_j}.
$$
On the other hand,
$$
\prod_{i,j} \frac{1}{1-\zeta_i y_j}=\prod_{j}\frac{1}{1-y_j^n}.
$$
This is the generating series for the complete sums $h_k$ evaluated at $y_1^n$, $y_2^n\ldots$ , that are exactly the plethysms $h_k[p_n]$:
$$
\prod_{j}\frac{1}{1-y_j^n}=\sum_k h_k[p_n[Y]].
$$
These plethysms expand in the Schur basis:
$$ h_k[p_n[Y]]=\sum_{\lambda} d(\lambda,(k),n) s_{\lambda}[Y].$$
The coefficient $d(\lambda,(k),n)$ is the evaluation you are looking for, for $k \cdot n = |\lambda|$.
For any particular computation of the $d(\lambda,(k),n)$ you may use SAGE or John Stembridge's SF package for Maple.
For a general study, you should find at least some examples of such plethysms in Macdonald's classical book. You will find more recent results by googling "plethysms power sums complete sums in Schur", for instance by William Doran, or Carbonara+Remmel+Yang.
EDIT: Actually it is easy to evaluate the Schur polynomials at the $n$-the roots of unity, from the definition of the Schur polynomials as "bialternants". One gets that the evaluation of the Schur polynomial is nonzero if and only if the classes modulo $n$ of $\lambda_n$, $\lambda_{n-1}+1, \ldots, \lambda_1+{n-1}$ are a permutation of the classes of $0,1,\ldots,n-1$. Then this evaluation is the sign of the permutation.
About the evaluation of the monomial functions, since this is what you are really interested in. A method for evaluating the monomial functions at the roots of unity is presented in Alain Lascoux and Marcel-Paul Schutzenberger's "Formulaire raisonné de fonctions symétriques", ex. 5.14, with a reference to a paper of 1881 by E. West. I think it generalizes indeed your computations by inclusion-exclusion.
The method consists in expanding the monomial functions in the power sum basis. Note that the power sum $p_r$ at the $n$-th roots of unity is $n$ if $n$ divides $k$, and $0$ else.
The expansion of monomial functions in the power sum basis is explained in Example 2.7 of the "Formulaire raisonné de fonctions symétriques". It is also explained in "On the Foundations of Combinatorial Theory. VII: Symmetric Functions through the theory of distribution and occupancy". Peter Doubilet. Studies in Applied Mathematics 51 (4), 1972.
Say you want to expand $m_{(\lambda_1,\ldots,\lambda_k)}$ in the power sum basis. Consider rather the "augmented monomial function" (which is, by the way, the function you are really interested in):
$$
M_{(\lambda_1,\ldots,\lambda_k)}=\sum x_{i_1}^{\lambda_1} x_{i_2}^{\lambda_2} \cdots x_{i_k}^{\lambda_k}
$$
where the sum is carried over all arrangements $x_{i_1},x_{i_2},\ldots,x_{i_k}$ of $k$ variables. The expansion of $M_{\lambda}$ in the power sum basis involves the set partitions $\Pi$ of $\{1,2,\ldots,k\}$. They form a lattice under refinement, whose smallest element is the partition $\hat{0}$ in $k$ singletons. This lattice admits a Möbius function. Let $B_1$, $B_2,\ldots,B_{\ell}$ be the blocks of $\Pi$. The Möbius function on the interval $[\hat{0},\Pi]$ is:
$$
\mu(\hat{0},\Pi)=(-1)^{k-\ell} \prod_i \left( \text{card} B_i -1\right)!
$$
Set $\lambda(\Pi)$ for the partition whose parts are the $\sum_{i \in B} \lambda_i$ for $B$ block of $\Pi$. Then the formula is:
$$
M_{\lambda}=\sum_{\Pi} \mu([\hat{0},\Pi]) p_{\lambda(\Pi)}
$$
where the sum is over all set partitions $\Pi$ of $\{1,\ldots,k\}$.
When evaluating at the $n$-th roots of unity you obtain:
$$
\sum_{\Pi} \mu([\hat{0},\Pi]) n^{\ell(\Pi)}
$$
where the sum is over all partitions $\Pi$ such that all parts of $\lambda(\Pi)$ are multiple of $n$, and $\ell(\Pi)$ is the number of blocks.
Douglas already commented that the asymptotics for fixed $p$ and $l\to \infty$ shoudl follow from standard methods. One gets
$$a_{\ell}^p\approx (p+1)^{2\ell+\frac{p+1}{2}}(4\pi \ell)^{-\frac{p}{2}}.$$
See theorem 4 in "Counting Abelian squares", by Richmond and Shallit. Notice that these numbers appear also in combinatorics when considering abelian squares, or more generally abelian powers, on a fixed alphabet.
For the asymptotics that you're interested in, at least in the unweighted case, one can say
$$a _{\ell} ^p=\sum _{j=0} ^{\ell} \binom{p}{j}\sum _{a _1+ \cdots +a _j = \ell \atop a _i \geq 1} \binom{\ell}{a _1,a _2,\dots,a _j}^2$$
which makes it clear that $a _{\ell}^p$ is a polynomial in $p$ of fixed degree $\ell$. The coefficient of $\binom{p}{\ell}$ is $(\ell!)^2$, and the coefficient of $\binom{p}{\ell -1}$ is $\frac{\ell-1}{4}(\ell!)^2$, so you have
$$a _{\ell}^p =\ell!p^{\ell}-\ell!\frac{\ell(\ell-1)}{4}p^{\ell-1}+O(p^{\ell-2}).$$
Best Answer
Notice that $$ y^p-x^p=\prod_{i=0}^{p-1} (y-\omega^i x)\ . $$ Therefore $$ \prod_{i=0}^{p-1} f(x\omega^i)= Res_y (y^p-x^p, f(y))\ . $$ Where $Res_y$ is the resultant of the polynomials in $y$. The above is just a particular case of the so called Poisson product formula. You can then compute the resultant using for instance Sylvester's determinant formula. See the review http://mate.dm.uba.ar/~alidick/papers/chapter1cd.pdf by Cattani and Dickenstein for a nice introduction to resultants and their properties.