Using Goertzel to find phase shift
If the complex frequency bin outputs of the two Goertzel filters are $a$ and $b$, you can calculate the phases as $\alpha = \text{atan2}\left(\text{Im}(a),\ \text{Re}(a)\right)$ and $\beta = \text{atan2}\left(\text{Im}(b),\ \text{Re}(b)\right)$, where $\text{Re}$ returns the real component and $\text{Im}$ returns the imaginary component of the argument. Most programming languages will have the two-argument $\text{atan2}$ function. The phase difference can be calculated as $\alpha - \beta = \text{atan2}\left(\text{Im}\left(a\ \text{conj}(b)\right),\ \text{Re}\left(a\ \text{conj}(b)\right)\right)$, where $\text{conj}$ denotes taking the complex conjugate, saving you one $\text{atan2}$ evaluation and wrapping of the phase, at the cost of one complex multiplication. 

There's one more trick available to calculate $a\ \text{conj}(b)$. Multiply the Goertzel output $a$ by the second input signal and lowpass filter the result. In frequency domain it looks like this:



The lowpass filtering isolates the circled spectral peak that corresponds to $a\ \text{conj}(b)$. This saves you one Goertzel filter at the cost of a lowpass filter. If you use this trick you must also take into account the phase shift due to the lowpass filter.

For fast $\text{atan2}$ see Methods of computing fixed point atan2 on FPGA.