At the bottom of the cylinder, $z=0$. Since $F_z(x,y,0)=0$ we find
$$\begin{align}
\int_{S_1} \vec F \cdot \hat n\, dS&=\int_{S_1}\vec F(x,y,0)\cdot \hat z\,dS\\\\
&=0
\end{align}$$
and there is no flux contributed.
At the top of the cylinder, $z=5$. Since $F_z(x,y,5)=15$ we find
$$\begin{align}
\int_{S_2} \vec F \cdot \hat n\, dS&=\int_{S_2}\vec F(x,y,5)\cdot \hat z\,dS\\\\
&=\int_0^{2\pi}\int_0^2 (15)\,\rho\,d\rho\,d\phi\\\\
&=60\pi
\end{align}$$
Putting it all together, we find
$$\begin{align}
\int_S \vec F\cdot \hat n \,dS&=\int_V \nabla \cdot \vec F\,dV-\int_{S_2} \vec F \cdot \hat n\, dS-\int_{S_1} \vec F \cdot \hat n\, dS\\\\
&=120\pi-60\pi-0\\\\
&=60\pi
\end{align}$$
Note that
\begin{align}
\nabla \cdot v = \dfrac{\partial v_x}{\partial x} + \dfrac{\partial v_y}{\partial y} + \dfrac{\partial v_z}{\partial z}
\end{align}
And in each partial derivative, since $x,y,z$ have units of length, the units of the divergence of $v$ is $\dfrac{1}{\text{length}}$. Multiplying by the volume element $d\tau$ implies that $(\nabla \cdot v)\, d\tau$ has units of $(\text{length})^2 = \text{area}$. So, there is no contradiction anywhere, even the LHS has units of area.
So, the key part you were missing is that you didn't take into account the units of divergence. You incorrectly assumed that just because $v$ is dimensionless, $\nabla \cdot v$ is also dimensionless. This is incorrect, because the differentiation process of $\dfrac{\partial}{\partial x}$ etc. introduces an extra factor of $\dfrac{1}{\text{length}}$.
Now you should convince yourself that in general, both sides of the equation
\begin{align}
\iiint_{\Omega} (\nabla \cdot v) \, d\tau = \iint_{\partial \Omega} (v \cdot n) \, da
\end{align}
have units equal to $(\text{(units of $v$)} \cdot \text{area})$.
Best Answer
The divergence theorem is about a three-dimensional body $B$, its boundary surface $\partial B$, and a vector field ${\bf v}$ with domain $\Omega\supset B$. The theorem says that $$\int_{\partial B}{\bf v}\cdot{\bf n}\>dS=\int_B{\rm div}\,{\bf v}({\bf x})\> dV\ .\tag{1}$$ In your case $B$ is a half ball, and ${\bf v}={\rm curl}(\vec F)$. It follows that ${\rm div}\,{\bf v}({\bf x})\equiv 0$ in $B$. But $\partial B$ is not just the spherical part of $\partial B$. The boundary of $B$ also contains a plane circular disc in the $(x,y)$-plane. It follows that the right hand side of $(1)$ consists of two terms, one of which you have already computed. Computing the flow of ${\bf v}$ through the disc is simpler, since ${\bf n}\equiv(0,0,-1)$ and $dS=dA$, where $dA$ is the standard area element in the $(x,y)$-plane.