The form of the correlation function does not flow under the coarse-graining step of the RG. Even before you do any coarse graining you would expect the correlation function to look like
$$\langle \phi(x) \phi(y) \rangle \sim a^\eta \frac{e^{-|x-y|/\xi}}{|x-y|^{d-2+\eta}}$$
(note that I have corrected a sign in the exponent). Here I have allowed for the correlation length $\xi$ to be finite. The "$\sim$" hides constant prefactors, and the $a^\eta$ is often hidden in this prefactor, which can naively make it look like dimensional analysis has broken.
If you now perform the first step of an RG analysis on this system, coarse-graining, then as long as you haven't integrated out the degrees of freedom at points $x$ and $y$ [a], the correlation function between those two points (in the original coordinates) will remain the same. The action itself acquires effective interactions to ensure ensure that this and all higher order correlation functions are unaffected.
Next, when you perform the rescaling step of the RG, setting $x \rightarrow \lambda x$ then the correlation function may look like it's changed, but that's just due to the change in scale. In particular, if $\xi \neq \infty$, then
$$\langle \phi(x) \phi(y) \rangle \sim \lambda^{-(d-2+\eta)} a^\eta \frac{e^{-|x-y|/(\xi/\lambda)}}{|x-y|^{d-2+\eta}}$$
and it appears that the correlation length has shrunk by a factor of $\lambda$. This is where the statement that finite correlation lengths decrease under RG, but it's really a property of the rescaling step (i.e., a system that is not scale-invariant will look different under a change of scale). However, if you are at the critical point (which you must tune to before doing any coarse-graining), then the exponential factor is just $1$ and you have pure power-law scaling, which gives the scaling property of the correlation function for a critical system:
$$\langle \phi(\lambda x) \phi(\lambda y) \rangle = \lambda^{-(d-2+\eta)} \langle \phi(x) \phi(y) \rangle.$$
Again, this scaling property comes from the system being scale-invariant, in principle one does not need to do RG to see it---if you tune a system to criticality and then rescale the lengths by $\lambda$, this is how the correlation function would scale. i.e., a scale invariant system will look the same under a change of scale. (See also this answer in another post.)
From this scaling property we see that if we were to also rescale the fields by $\phi \rightarrow \lambda^{-(d-2+\eta)/2}$ that this prefactor would cancel out. However, the fields $\phi$ themselves still have engineering dimension $(d-2)/2$. (The factor $\lambda$ here is dimensionless, so scaling $\phi$ by any power of it does not change the dimensions of $\phi$).
You may still be wondering where the factor of $a^\eta$ comes from in the correlation function, and in particular why it is there despite the fact that I have just claimed that the engineering dimension of $\phi$ is $(d-2)/2$ even when $\eta \neq 0$. The reasons are as follows:
The scaling of $\phi(x)$ is determined by the requirement that the action is dimensionless, which tells you that it has "engineering dimension" $(d-2)/2$ (in units of momentum), as you have noted. Dimensional analysis then demands that the correlation function must have the form
$$\langle \phi(x) \phi(y) \rangle = |x - y|^{-(d-2)} f(|x-y|/\xi,a/\xi)$$
where $f(\cdot,\cdot)$ is some function of dimensionless combinations of parameters/variables. Here we only care about variables with units of length. Since one expects that $a/\xi \ll 1$ close to criticality, the naive approximation is to assume we can set this argument to $0$ to obtain
$$\langle \phi(x) \phi(y) \rangle = |x - y|^{-(d-2)} f(|x-y|/\xi,0),$$
which gives the mean-field scaling. It turns out this is incorrect for certain dimensions (e.g., in $d < 4$ for the Ising universality class) because the function $f(\cdot,x)$ is not regular in the limit $x \rightarrow 0$, but instead behaves like $x^{\eta}$. We may make this leading order scaling behavior explicit and write
$$f(|x-y|/\xi,a/\xi) = (a/\xi)^\eta (|x-y|/\xi)^{-\eta} g(|x-y|/\xi,a/\xi) = a^\eta |x-y|^{-\eta} g(|x-y|/\xi,a/\xi),$$
introducing some new unknown function $g(\cdot,x)$ that is regular at $x = 0$, as well as the unknown power $\eta$ that cannot be determined by dimensional analysis. Now we see that for $a \ll \xi$ the correlation function must have the form
$$\langle \phi(x) \phi(y) \rangle = |x - y|^{-(d-2+\eta)} a^\eta g(|x-y|/\xi,0),$$
which is consistent with dimensional analysis but has given rise to an unexpected power $\eta$ due to the fact that the microscopic length scale $a$ cannot be neglected, even for distances $|x-y| \gg a$. This macroscopic dependence on the microscopic length scale is the "surprise" of critical systems and the origin of these "anomalous" dimensions.
Footnotes/caveats:
[a] For simplicity you can imagine we are doing a decimation-style RG in real-space, or if doing a momentum-shell type RG that the separation $|x-y|$ is larger than the integrated out length scales.