You lack sufficient information to construct a meaningful surface, since you have only 2 pieces of information in terms of Reynolds number.
First, plot what you have. A 3-d plot is effectively useless here, and makes it far more difficult to understand what you do have.
plot(Cl1,Cd1,'-or',Cl2,Cd2,'-*b')
legend('RE 400K','RE 500K')
So we see two curves that are not terribly similar in shape. There is a bit of noise seen. And the abscissa for these two curves are not the same, so you could not just use a tool like interp2 to interpolate. Extrapolating those curves, or for Reynold's numbers that fall (significantly) outside the interval [4e5,5e5] would seem to be a foolish (risky) task, due to the high amount of curvature of the respective curves, and due to the noise in the data.
If you merely wish to predict a point for some other Reynolds number however, you could simply use a spline interpolant, a linear interpolant, or perhaps a smoothing spline for each curve. Then, since you have only two distinct Reynolds numbers and thus only two distinct curves, just use linear interpolation between the two. For example, if you have the curve fitting toolbox, I might do this:
CF1 = fit(Cl1,Cd1,'smoothingspline');
CF2 = fit(Cl2,Cd2,'smoothingspline');
CFinterp = @(Cl,Re) (5e5 - Re)/(5e5-4e5)*CF1(Cl) + (Re - 4e5)/(5e5-4e5)*CF2(Cl);
plot(Cl,CFinterp(Cl,Reinterp),'g-')
legend('RE 400K','RE 500K','RE 450K')
I hard coded the Reynolds numbers in the interpolant there, but it is easy enough to fix that.
