We describe the incorporation of polarized radiative transfer into the atmospheric radiative transfer modelling code VSTAR (Versatile Software for Transfer of Atmospheric Radiation). Using a vector discrete-ordinate radiative transfer code we are able to generate maps of radiance and polarization across the disc of a planet, and integrate over these to get the full-disc polarization. In this way we are able to obtain disc-resolved, phase-resolved and spectrally-resolved intensity and polarization for any of the wide range of atmopsheres that can be modelled with VSTAR. We have tested the code by reproducing a standard benchmark problem, as well as by comparing with classic calculations of the polarization phase curves of Venus. We apply the code to modelling the polarization phase curves of the hot Jupiter system HD 189733b. We find that the highest polarization amplitudes are produced with optically thick Rayleigh scattering clouds and these would result in a polarization amplitude of 27 ppm for the planetary signal seen in the combined light of the star and planet. A more realistic cloud model consistent with the observed transmission spectrum results is an amplitude of ~20 ppm. Decreasing the optical depth of the cloud, or making the cloud particles more absorbing, both have the effect of increasing the polarization of the reflected light but reducing the amount of reflected light and hence the observed polarization amplitude.