Values for angular momentum in energy levels of isotropic Harmonic oscillator Consider the Harmonic oscillator.
Some considerations:
$H_{xyz}=H$  is not a C.S.C.O. in $\mathcal{E}=\mathcal{E}_{xyz}$  by itself: even if we know the eigenvalue (the energy) of $H$ , we may not not the state of the system: e.g if $E=1$
 , all we know is that the state is a linear combination of $\vert100\rangle,\vert010\rangle$  and $\vert001\rangle$ .
However, $\{H_{x},H_{y},H_{z}\}$  is a C.S.C.O. This is easy to see.
Now, $[H,L_{i}]=0$ since $H$  is rotationally invariant. 
One can prove that $\{H,L^{2},L_{z}\}$ is a C.S.C.O.
Thus, we get two basis for $\mathcal{E}$ : $\{\vert n_{x}\,n_{y}\,n_{z}\rangle\}$
  and $\{\vert n,l,m\rangle\}$, one for each of the C.S.C.O.s found.
We can of course go from one basis to the other.
A common question asked is: for a given n, what are the possible values of l? 
I will try to answer this using a method which I haven't found anywhere, so I would like someone to confirm it.
For example, for $n=1$, I know that $\mathcal{E}(n=1)$ is three-dimensional, since it is spanned by $\{\vert100\rangle,\vert010\rangle,\vert001\rangle\}$. Thus, one of the following is true:
$$\mathcal{E}(n=1)=3\mathcal{E}(n=1,l=0)$$
or 
$$\mathcal{E}(n=1)=\mathcal{E}(n=1,l=1)$$
So, if I find one state of $\mathcal{E}(n=1)$
  with $m=1$, then that state cannot belong to $\mathcal{E}(l=0)$
  and so $\mathcal{E}(n=1)=\mathcal{E}(n=1,l=1)$
  holds.
Another example: $n=2$
 . I know that $\mathcal{E}(n=2)$
  is six-dimensional, since it is spanned by $\{\vert200\rangle,\vert020\rangle,\vert002\rangle,\vert110\rangle,\vert101\rangle,\vert011\rangle\}$
 . Thus, one of the following is true:
$$\mathcal{E}(n=2)=6\mathcal{E}(n=2,l=0)$$
or 
$$\mathcal{E}(n=2)=\mathcal{E}(n=2,l=0)\oplus\mathcal{E}(n=2,l=2)$$
So, if I find one state of $\mathcal{E}(n=2)$
  with $m=-1$, for instance, then that state cannot belong to $\mathcal{E}(l=0)$  and so $\mathcal{E}(n=2)=\mathcal{E}(n=2,l=0)\oplus\mathcal{E}(n=2,l=2)$  holds.
Is this argument correct? The results are, but I'm not sure the reasoning is: I assumed that, if I have a state with a certain l
  in $\mathcal{E}(n)$
 , then all the $2l+1$
  states $\vert n,l,m\rangle$
 are in $\mathcal{E}(n)$. I do not know how to justify this.
If the method or the question are still unclear, please tell me and I'll try to clarify them.

After reading what Emilio and ZeroTheHero and some other things, I tried another method which seems more general. Here it goes.
Let
$$
\begin{cases}
V_{1}=-\frac{a_{x}+ia_{y}}{\sqrt{2}}\\
V_{0}=a_{z}\\
V_{-1}=\frac{a_{x}-ia_{y}}{\sqrt{2}}
\end{cases}\mbox{ and }\begin{cases}
V_{1}^{\dagger}=-\frac{a_{x}^{\dagger}-ia_{y}^{\dagger}}{\sqrt{2}}\\
V_{0}^{\dagger}=a_{z}^{\dagger}\\
V_{-1}^{\dagger}=\frac{a_{x}^{\dagger}+ia_{y}^{\dagger}}{\sqrt{2}}
\end{cases}
$$
We have 
$$
[J_{z},V_{q}^{\dagger}]=q\hbar V_{q}^{\dagger}
$$
$$
[J_{z},V_{q}]=q\hbar V_{q}
$$
$$
[N,V_{q}^{\dagger}]=V_{q}^{\dagger}
$$
$$
[N,V_{q}]=-V_{q}
$$
So the “daggered” operators increase the eigenvalue $n$
  of $N$
  by one (and thus the energy by $\hbar\omega$
 ), while the “undaggered” operators decrease it by one. Also, the $V_{q}^{(\dagger)}$
  change the eigenvalue of J_{z}
  by $q\hbar$
 .
With this in mind, we can face the problem.
Assumption: All states of the system can be obtained by applying these operators to the ground state $\vert0\rangle$
 . I don't know why!-Please explain it if you do.
For $n=1$
 , the biggest value for $m$
  one can get for a state of $\mathcal{E}(n=1)$
  is $1$, since $(V_{1}^{\dagger})^{1}\vert0\rangle$
 is the only way to get the maximum $m$
  while (increasing $n$
  to $1$) by application of the operators defined above, and indeed $L_{z}V_{1}^{\dagger}\vert0\rangle=\hbar V_{1}^{\dagger}\vert0\rangle$
 . Now, this state, having $m=1$
 , must have $l\geq1$
 . It cannot be $l>1$
  since this would imply that I had states with $m>1$
  (obtained by applying $L_{+}$
  to $V_{1}^{\dagger}\vert0\rangle$
 ) in $\mathcal{E}(n=1)$
  (the Hamiltonian commutes with $L_{+}$
  because of rotational invariance and so does N
 ), which is absurd. For $m=0$
 , we easily see that there is only one state: $V_{0}^{\dagger}\vert0\rangle$
 . Since there is only one, it must be the one of $\mathcal{E}(l=1)$
 . There is no need to analyze other values for m
  (negative values will bring nothing new). In conclusion, the only value for $l$
  for states with $n=1$
  is $1$
 .
Let us try something harder.
For $n=3$
 , the biggest value for $m$
  for a state of $\mathcal{E}(n=3)$
  is three, corresponding to the state $(V_{1}^{\dagger})^{3}\vert0\rangle$
 . Clearly, there are no other states with $m=3$
  in $\mathcal{E}(n=3)$
 . Just like before, this means that the space $\mathcal{E}(l=3)$
  is in the direct sum decomposition of $\mathcal{E}(n=3)$
 . Now, for $m=2$
  we only have one state $V_{0}^{\dagger}(V_{1}^{\dagger})^{2}\vert0\rangle$
 , so there is no $l=2$
  representation in the direct product decomposition. For $m=1$
 , we have $V_{-1}^{\dagger}(V_{1}^{\dagger})^{2}\vert0\rangle$
  but also $(V_{0}^{\dagger})^{2}V_{1}^{\dagger}\vert0\rangle$
 , so there is one more representation (besides $l=2$
 ) with states with $m=1$
 . It can't be $l=2$
  because there are no $m=2$
  states. It is therefore an $l=1$
  representation. Finally, for $m=0$
  I only have two states: $V_{-1}^{\dagger}V_{0}^{\dagger}(V_{1}^{\dagger})^{1}\vert0\rangle$
  and $(V_{0}^{\dagger})^{3}
 \vert0\rangle$
 , which must belong to the two representations already found. Therefore, $l\in\{1,3\}$
 .
 A: Your method is laborious and will not practical except for small dimensions.  
First, on a side note, you are attempting to redo or rediscover much of the work that was done in the nuclear physics community in the 60s and 70s.  The original paper that got the ball rolling is by Phil Elliott (this one: J.P. Elliott, Proc. Roy. Soc. A 245 (1958) 562).  A good source the specialized kind of tensor work you want to do is: 


*

*Polynomial bases and Wigner coefficients for SU(3) ⊃ R3, by 
R.T. Sharp, Hans C. Von Baeyer and S.C. Pieper, Nuclear Physics A
Volume 127, Issue 3, 21 April 1969, Pages 513–524. 


Note there is also the review paper by Slansky (Group theory for model building or something like this) that will provide table of branching rules for $SU(3) ⊃ SO(3)$.
You are generally right in your approach but counting states become impractical.  Of course, I take it you realize that, if you solve the Schr\"{o}dinger equation in spherical coordinates, it follows immediately that the possible values of $\ell$ are $N,N-2,N-4$ all the way down to $0$ or $1$.  
The most intuitive way to see this is to look at $(a_z^\dagger)^N\vert 0\rangle$.  This state certainly is in your representation, and in terms of the position basis, the polynomial $z^N$ will be in your basis.  $(a_z^\dagger)^N$ is a reducible tensor with $M=0$, and its use is that it must contain all possible values of $\ell$ since, if you write $(a_z^\dagger)^N\vert 0\rangle$ in a position basis, it will contain the polynomial $z^N$.  If you expand $z^N$ in terms of Legendre polynomial (this expansion can be found in many places), you will also get the angular momentum contents of $z^N$ and thus all states in the $N$ shell.  
The disadvantage of this methods is that the expansion coefficients of $z^N$ in terms of Legendre polymomials does not give the correct normalization of the $\vert \ell,M=0\rangle$ states.  This normalization problem has plagued the Elliott formalism since its inception, but if you do not care about the states and just want to get the angular momentum content this will work.
The paper by Sharp actually constructs the angular momentum states directly using creation and destruction operators if you are interested in directly constructing states.
As additional information, your problem is limited to $SU(3)$ irreps of the type $(N,0)$ (in the Dynkin notation).  For general representations of the $(\lambda,\mu)$ type, the branching rules were very very laboriously worked out by Elliott, hence the appeal of the shortcut to be found in Slansky (some sort of matrix projection method - I haven't read this trick in a while).  In particular, some $SO(3)$ irreps will occur more than once, and there is no canonical way of lifting this degeneracy, although the closest thing to a solution to this problem was developed by Rowe in 
 1. D. J. Rowe, Resolution of missing label problems; a new perspective on K‐matrix theory, Journal of Mathematical Physics, Volume 36, Issue 3, November 1994,
 2. D J. Rowe, R. Le Blanc and J. Repka, A rotor expansion of the $su(3)$ Lie algebra, Journal of Physics A: Math. Gen. volume 22 L309-316 (1989)
Pre $K$-matrix theory, most of the work was entirely numerical since the multiple occurrences of a given $\ell$ value were resolved using Gram-Schmidt process.  There is considerable work by Jerry Draayer on this topic.

Edit: in answer to your comment.  Your assumption is 99% correct. The more correct statement is that any polynomial of degree $N$ in the raising operators will give you an state in the shell, and any state with exactly $N$ total excitations constructed as you suggest will be in the irreducible representation $(N,0)$. The reason it works is that the $u(3)$ Lie algebra (or more correctly its complexification) is spanned by the $9$ number-conserving operators $a^\dagger_i a_k$.  The $su(3)$ algebra is obtained by removing from those $9$ the total number operator $\hat N=\sum_k a^\dagger_k a_k$, which, up to a proportionality factor, is just the scalar you constructed out of $V$ and $V^\dagger$.   $\hat N$ commutes with everything else so it's usually taken out of $u(3)$ to get $su(3)$.
Your method is in general sound in the sense that if you start with any polynomial $P(V_k^\dagger), k=1,2,3$ in the creation operators you will get a vector in the irrep, and repeatedly acting with all possible elements of $su(3)$ on this vector will generate the entire irrep.  This is because an irreducible representation does not contain invariant subspace, so (in a simplified interpretation) it must be possible to reach any vector from any other vector by the repeated action of something in the algebra.  The state counting that you do it largely correct but you would have to show, in the case $N=3$, that you cannot find a polynomial of degree $3$ in the $V^\dagger_k$ to give you something with $\ell=2$.  This is not automatic but follows from accidental zeroes of the CG coefficients.  Otherwise you are on the right track.
A: I don't think your method is sound. Dimension-counting is a nice idea, but it doesn't really work. For example, nothing in your arguments prohibits the decomposition
$$\mathcal{E}(n=2) = 3\mathcal{E}(n=2,l=0) \oplus \mathcal{E}(n=2,l=1),$$
where $\mathcal{E}(n=2)$ means the eigenspace of $H$ with $n=2$, and the informal notation $3\mathcal{E}(n=2,l=0)$ means the direct sum of three independent eigenspaces of $H$ and $L^2$ with $n=2$ and $l=0$ respectively. Moreover, this gets worse and worse at higher dimensions. You are arriving at your conclusions because you know what the results are to begin with.

This much is true, though:

I assumed that, if I have a state with a certain $l$  in $\mathcal{E}(n)$, then all the $2l+1$   states $\vert n,l,m\rangle$ are in $\mathcal{E}(n)$.

This is because $H$ commutes with all the components of $\mathbf L$, and you can get from any $|l,m\rangle$ to any other $|l,m'\rangle$ by applying suitable combinations of the components of $\mathbf L$. (This is most easily done, for those states, by applying the ladder operators $L_\pm = L_x \pm i L_y$.) More fundamentally, this is because the phrase "a subspace of definite $l$" really denotes an irreducible representation of the rotation group, and those are an everything-or-nothing kind of deal.

If you want an elegant method to resolve which $l$s contribute to the eigenspace at any given $n$, I would recommend the methods in this question or these notes.
