Showing volume and surface integration is unaffected by the singularity at $\mathbf{r'}=\mathbf{r}$ 
This question is not entirely similar to the question here. Please read this question and the reader will see it is obviously not the same.

$\mathbf{M'}$ is a continuous vector field in volume $V'$ and $P$ be any point on the surface of $V'$ with position vector $\mathbf {r}$

PART I:
Consider the expression:
$\displaystyle \iiint_{V'} \mathbf{M(r')}.\nabla'  \left(
\dfrac{1}{\left|   \mathbf{r}-\mathbf{r'}  \right|}       \right) dV'$
Take the origin of our coordinate system at $P$ (see the diagram) and write $dV'$ as spherical volume element. Then the above expression can be written as:
$\displaystyle \iiint_{V'} \mathbf{M(r')}.\dfrac{\hat{r'}}{{r'}^2}
{r'}^2\ \sin\theta\ d\theta\ d\phi\ dr' = \iiint_{V'} \left[
\mathbf{M(r')}.{\hat{r'}} \right] \sin\theta\ d\theta\ d\phi\ dr' \tag{1}$
The integrand is defined everywhere except at point $P$ where $r'=0$ and the integrand is $\frac{0}{0}$.
Since $\mathbf{M(r')}$ is finite everywhere, there is no blowing up of the integrand at any point.
Therefore we can directly integrate equation $(1)$ just like ordinary integrals.
PART II:
Using the vector identity $\nabla.(\psi \mathbf{A})=\mathbf{A}.(\nabla \psi)+\psi (\nabla.\mathbf{A})$:
$\displaystyle \iiint_{V'}    \left[     \nabla' . \left(
\dfrac{\mathbf{M'}}{\left|   \mathbf{r}-\mathbf{r'}  \right|} \right)
\right]   dV' = \iiint_{V'} \mathbf{M'}.\nabla'  \left(
 \dfrac{1}{\left|   \mathbf{r}-\mathbf{r'}  \right|}       \right) dV' +
\iiint_{V'} \dfrac{\nabla' . \mathbf{M'}}{\left|   \mathbf{r}-\mathbf{r'}
\right|}  dV'\tag{2}$
Now for simplicity, let's take the origin of our coordinate system at $P$ (see the diagram). Thus equation $(2)$ becomes:
$\displaystyle \iiint_{V'}    \left[     \nabla' . \left(
\dfrac{\mathbf{M'}}{r'} \right)      \right]   dV' = \iiint_{V'}
\mathbf{M'}.\nabla'  \left(     \dfrac{1}{r'}       \right) dV' +
\iiint_{V'} \dfrac{\nabla' . \mathbf{M'}}{r'}  dV' \tag{3}$
Now by writing $dV'$ as spherical volume element, equation $(3)$ becomes:
\begin{align}
\iiint_{V'}    \left[     \nabla' . \left(
\dfrac{\mathbf{M'}}{r'} \right)      \right]   {r'}^2\ \sin\theta\ d\theta\
d\phi\ dr
&=\iiint_{V'} \mathbf{M'}.\nabla'  \left(
 \dfrac{1}{r'}  \right) {r'}^2\ \sin\theta\ d\theta\ d\phi\ dr'\\
&+\iiint_{V'} \dfrac{\nabla' . \mathbf{M'}}{r'}  {r'}^2\ \sin\theta\ d\theta\
d\phi\ dr'\\
&= \iiint_{V'} (\mathbf{M'}. \hat{r'}) \sin\theta\
d\theta\ d\phi\ dr' \\
&+ \iiint_{V'} (\nabla' . \mathbf{M'})\ r'\ \sin\theta\
d\theta\ d\phi\ dr'
\tag{4}
\end{align}
In both terms:
The integrands are defined everywhere except at point $P$ where $r'=0$ and the integrand $\frac{0}{0}$.
Since $\mathbf{M(r')}$ is finite everywhere, there is no blowing up of the integrands at any point.
Therefore we can directly integrate both terms in the $RHS$ of
equation $(4)$ just like ordinary integrals.
Therefore we can directly integrate $LHS$ of equation $(4)$
just like ordinary integrals.
Now my question is: Is Gauss divergence theorem applicable to  $LHS$ of equation $(4)$ ?
If it is applicable, then, since point $P$ lies on the surface, there would be a singularity in the equation:
$\unicode{x222F}_{S'}    \left[    \left(\dfrac{\mathbf{M'}}{\left| 
\mathbf{r}-\mathbf{r'}  \right|} \right)   . \hat{n}      \right]   dS'$
How to deal with it?
 A: I'll assume that $V'$ is bounded and that $\mathbf{M'}$ is continuously differentiable. To your primary question:

Is [the] divergence theorem applicable to of equation (4)?

The answer is, strictly speaking, no. The hypotheses of the divergence theorem include that the vector field in question is continuously differentiable, and $\frac{\mathbf{M'(r)}}{r}$ is not (in general). 
However, one can still rigorously employ the divergence theorem to yield your expression by removing a ball $B_\epsilon(P)$ from $V'$. Writing $B_\epsilon '(P) \equiv B_\epsilon(P) \cap V' $,
$$\iiint_{V'} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' = 
\iiint_{V' \backslash B_\epsilon'(P)} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right)  \right]dV' + \iiint_{B_\epsilon'(P)} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right)  \right]dV' $$
The divergence Theorem may be applied to the first term for any $\epsilon>0$, and one can show that the second term goes to $0$ in the limit $\epsilon \to 0$ via your chain of equalities in equation $(4)$, which hold $a.e.$ in $B_\epsilon'(P)$, since $\mathbf{M'}$, being continuous on the compact set $\overline{V'}$, is bounded, say by $M>0$.
With some computation, one can show that the integral over $\partial(V' \backslash B_\epsilon '(P))$ arising from the first term approaches the integral over $S' = \partial V'$ in the same limit, assuming $S'$ is smooth in a neighborhood of $P$. To be explicit, 
$$ \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'
= \left( \iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)} - \iint_{\partial B_\epsilon '(P) \backslash \partial V'} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' $$
That is, the integrals are the same up to the difference of the integrals over the piece of $\partial B_\epsilon '(P)$ from $\partial V'$ and that from $\partial B_\epsilon(P)$. Since $S'$ is smooth near $P$, it is locally the graph of a smooth function $f(x,y)$ on the tangent plane at $P$ with standard linear coordinates $(x,y)$. So, for $\epsilon$ sufficiently small, we may compute the first integral on the RHS above as an integral in the tangent plane, with area form $\sqrt{1 + f_x^2 + f_y^2}dxdy=\sqrt{1 + f_x^2 + f_y^2}sdsd\theta$ in polar coordinates $(s,\theta)$. The term $\sqrt{1 + f_x^2 + f_y^2}$ is continuous and therefore bounded, say by $C$ (independent of $\epsilon$ for $\epsilon$ sufficiently small), and further $r'(s,\theta)=\sqrt{s^2+f(s,\theta)^2} \geq s$. Hence
$$\left|\iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)}\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'\right| \leq C \int_0^\epsilon \int_0^{2\pi}|\mathbf{M}(s,\theta)| \frac{s}{\sqrt{s^2+f^2}}dsd\theta \leq 2\pi MC \epsilon$$
Note that the exact integration here may not range over all of $0 \leq s \leq \epsilon$, but the bound still holds. Finally, in the second term on the RHS above, $r'=\epsilon$, so that integral is bounded by $4\pi M\epsilon$, and we have
$$\left| \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' \right| \leq 2\pi M(2+C)\epsilon $$
which goes to $0$ as $\epsilon \to 0$. Since my first equation holds for all $\epsilon >0$, we may take the limit of both sides as $\epsilon \to 0$ to obtain
$$\iiint_{V'} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' = \unicode{x222F}_{\partial V'} \left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'$$
As a naive application of the divergence theorem would have suggested. 
Note that that "singularity" in the RHS integration does not cause the integral to diverge for essentially the same reason that it didn't in the volume integral: by passing to an appropriately defined set of coordinates in a neighborhood of $P$, we see that the area form gets small fast enough in this neighborhood for the integral to converge, just as the volume form did in your computations.
