  We study the effects of including a running coupling constant in high-density
QCD evolution. For fixed coupling constant, QCD evolution preserves the initial
dependence of the saturation momentum $Q_s$ on the nuclear size $A$ and results
in an exponential dependence on rapidity $Y$, $Q^2_s(Y) = Q^2_s(Y_0) \exp{[
\bar\alpha_s d (Y-Y_0) ]}$. For the running coupling case, we re-derive
analytical estimates for the $A$- and $Y$-dependences of the saturation scale
and test them numerically. The $A$-dependence of $Q_s$ vanishes $\propto 1/
\sqrt{Y}$ for large $A$ and $Y$. The $Y$-dependence is reduced to $Q_s^2(Y)
\propto \exp{(\Delta^\prime\sqrt{Y+X})}$ where we find numerically
$\Delta^\prime\simeq 3.2$. We study the behaviour of the gluon distribution at
large transverse momentum, characterizing it by an anomalous dimension
$1-\gamma$ which we define in a fixed region of small dipole sizes. In contrast
to previous analytical work, we find a marked difference between the fixed
coupling ($\gamma \simeq 0.65$) and running coupling ($\gamma \sim 0.85$)
results. Our numerical findings show that both a scaling function depending
only on the variable $r Q_s$ and the perturbative double-leading-logarithmic
expression, provide equally good descriptions of the numerical solutions for
very small $r$-values below the so-called scaling window.
