To solve this ODE system, begin by determining its behavior near r == 0.
Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]
(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0 *)
In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)
eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1,
a[eps] == 1, g[eps] == eps gp,
WhenEvent[g[r] > 12/10, {bool = 1, "StopIntegration"}],
WhenEvent[{g[r] < 8/10, a[r] < 0}, {bool = 0, "StopIntegration"}]},
{a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 60;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp];
rm = (First[st] /. r -> "Domain")[[1, 2]];
If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i;
If[bool == -1, Return[]], {i, imax}];
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]
The solutions takes only a few seconds.