Irrelevant! This is a two dimensional surface and so can be written in terms of two parameters. Here, they have chosen to use x and y as those parameters- z is then a function of x and y. In fact, here's how you could have calculated that:

from x^2/a^2+ y^2/b^2+ z^2= 1, we can solve for z: z= (1- x^2/a^2- y^2/b^2)^{1/2} and write the "vector equation" of the surface- r(x,y)= <x, y, z>= <x, y, (1- x^2/a^2- y^2/b^2)^{1/2}> in terms of the two "parameters" x and y.

Now, differentiate that with respect to x and y:

r_x= <1, 0, -x/a^2(1- x^2/a^2- y^2/b^2)^{-1/2}>= <1, 0, -x/(a^2z)>

r_y= <0, 1, -y/b^2(1- x^2/a^2- y^2/b^2)^{-1/2}>= <0, 1, -y/(b^2z)>

Take the cross product of those two vectors., the "fundamental vector product" for this surface

<x/(a^2z), y/(a^2z), 1>

That, multiplied by the differentials of the parameters, x and y, gives

<x/(a^2z), y(a^2z), 1> dxdy as the problem says.

You are completely misunderstanding the problem and making it much harder than it is. You are NOT to integrate over the ellipsoid or even the top part of it. You are to integrate over the region S with is the part of the ellipsoid directly above the square bounded by x= 0, y= 0, x= 1, y= 1. That does NOT make a closed surface even adding the "base" (by which I take it you mean the ellipse where the ellipsoid passes through the xy plane) so the divergence theorem does not apply. And since the are of integration is a square, with no circular symmetry, cylindrical coordinates would make no sense at all.I have a feeling I may need to use cylindrical polar coordinates to solve this problem? Or maybe I should use the divergence theorem to calculate F.ds over the closed surface and then evaluate it over the base and take that away?

Just do what they say:

integrate F dot dS as given, with limits of integration, x= 0 to 1, y= 0 to 1.

Many thanks!