Does using the Goertzel algorithm actually give better frequency resolution? I am reading this article, and I am getting a little confused by the author's liberal use of 'frequency resolution' regarding the Goertzel Algorithm.

Basic question: Does using the Goertzel algorithm actually give you more frequency resolution over a specific band of interest, or does it simply efficiently compute the FFT over only the specified band of interest, but at the same frequency resolution specified by sampling frequency divided by number of samples? 

For example, lets say the $F_s$ is 100 KHz, (fixed) and the number of data samples $N$ is 10000. (Also fixed). If I compute a normal FFT, where FFT length is also $N$, my frequency resolution is $\frac{F_s}{N}$ as to be expected, and it will be equal to 10 Hz. This means my bins are separated by 10 Hz, from -50,000 Hz to 50,000 Hz.

Now let us say I want to use the Geortzel Algorithm to only look at frequencies in the range of say, 20,000-21,000 Hz. If I use the same $N$ for number of samples, and use the same $N$ for my FFT size, then what is my frequency resolution? Still 10 Hz? Or is it $\frac{21,000-20,000}{10000} = 0.1$Hz? 

I have a feeling that I am not really increasing my frequency resolution, as much as simply interpolating points on the main lobe, by using the same $N$ in to evaluate the frequencies from 21,000 to 20,000 as I did from 0 to 50,000. 

Is this a correct understanding?
Your understanding is correct.

Goertzel algorithms give almost exactly the same result as 1 bin of a DFT or FFT of the same length or number of samples (and where the FFT twiddle factors are generated by a trig recursion), when used for frequencies that are exactly integer periodic in the Goertzel length.  But many forms of the Goertzel algorithm provide only the magnitude, not the complex result or the phase of the FFT 1 bin result.  Also, computationally, a generic Goertzel may be numerically slightly less stable than a generic FFT.  For non-integer-periodic-in-aperture frequencies, the result is equivalent to a sinc interpolation between bins of a DFT or FFT of the same length (which can be slightly more accurate of an interpolation than the more typical parabolic interpolation of FFT results).

The interpolation might be said to increase resolution in a graphical sense (more plot points), or making it easier to visually spot the maxima, but not in the information theoretic sense, nor for better separating 2 close spaced spectral lines as 2 separate peaks.