Suppose $D$ is a simple closed curve, and $z_0$ is away from $\partial D$. Can I apply the Cauchy–Goursat theorem and say that
$$ \oint_{\partial D} \log|z_0 - z| \, dz = 0 \ ?$$
I think that I can. $\log$ is analytic away from its singularity, so I believe the criteria for Cauchy–Goursat should be satisfied. I suspect there may be some issue with choice of a branch cut, but I think that that can be solved just by choosing an appropriate cut.
On
You can. You only have to worry about branch cuts if you are dealing with the log of a complex argument (note that $|z - z_0|$ is real). In the former case we have $\log z = \log |z| + \arg z$ from the polar form of the complex number and the only multivalued part comes from the argument function.
Fun fact: you can actually use the integral of $1/z$, defined as $\log z$ as the basis for proving the Cauchy and residue theorems and formulating the usual theory of complex integration. See Ahlfors' classic book for the details.
No, this is true for a circle by direct calculation, but not for a more general shape: for example, take $z_0=0$ (so the integrand is $\log{\lvert z \rvert}$), and $D$ as the interior of the region $0<\arg{z}<\theta$, $a<\lvert z \rvert <b$ (i.e. a sector of an annulus centred at $z=0$). Then if $z=re^{it}$, the boundary is composed of the sections $$ a<r<b, t=0, \\ r=b, 0<t<\theta, \\ a<r<b, t=\theta, \\ r=a, 0<t<\theta, $$ with anticlockwise orientation. Thus, substituting into the integral gives \begin{align} I &= \int_a^b (\log{r}) \, dr + \int_0^{\theta} (\log{b}) \, ibe^{it} \, dt + \int_b^a (\log{r}) \, e^{i\theta} \, dr + \int_{\theta}^0 (\log{a}) \, iae^{it} \, dt \\ &= (1-e^{i\theta})(b\log{b}-a\log{a}-b+a) + (b\log{b}-a\log{a})(e^{i\theta}-1) \\ &= (a-b)(1-e^{i\theta}), \end{align} which is not zero unless $\theta=0$ or $a=b$.
The problem is that $\log{\vert z_0-z \rvert}$ is the real part of a holomorphic function, namely $\log{(z_0-z)}$, and is therefore not holomorphic itself (it is real-valued, so cannot preserve angles, for example). It is true that $$ 0=\int_{\partial D} \log{(z_0-z)} \, dz=\int_{\partial D} \big( \log{\vert z_0-z \rvert} + i\arg{(z_0-z)} \big)\, dz, $$ but because, to put it crudely, $dz$ is complex, the real and imaginary parts of the integrand cannot be separated so that their individual contributions give zero: the complex path mixes them up. All you can have from Cauchy's theorem is $$ \int_{\partial D} \log{\vert z_0-z \rvert} \, dz = -i\int_{\partial D} \arg{(z_0-z)} \, dz. $$