How would you calculate the flux of a vector field $\displaystyle F(x,y,z)=(y,x,z)$ downwards through the surface $\displaystyle z=16-x^2-y^2$?
The first thing I did is write the parametrization of the surface as $r(x,y)=(x,y,16-x^2-y^2)$ and calculated the $\partial_xz $ and $\partial_yz$. From there I got $r_x \times r_y=(2x,2y,1)$ and tried to calculate $$\int_PF \cdot (r_x\times r_y) dA$$ where P is projection of given surface onto $xy$-plane.
But i got $z$ coordinate in the integral, and I'm integrating in 2D, so what's the mistake?
As in my comment, since you have parametrized the position of each coordinate in the $x,y,z$ space with two variables, it is no problem that you have a third coordinate, since $z$ of course will be expressed as the other two ($x,y$), which means you will have a double integral (as you would expect for a surface integral)
Your integral (assuming your normal vector is correct (this is your vector after you took the cross product)) will expand to
$$\iint_{S} \begin{pmatrix} y \\ x \\ 16-x^2-y^2 \end{pmatrix}\cdot \begin{pmatrix} 2x \\ 2y \\ 1 \end{pmatrix} \mathrm{d}x\mathrm{d}y \qquad $$ and you can then put in the appropriate bounds for $x$ and $y$ for the limits.