I wonder what the QED Lagrangian would look like if you carefully write out all units of the terms and make sure they are consistent. I think there is something missing about Coulomb.
Can you write down the QED Lagrangian such that the physical units of all terms agree and state the units of all variables and terms?
The Wikipedia article you linked to gives the action, not the Lagrangian or Lagrangian density, so I’ll do the same.
Note that since $x^0 \equiv ct$, we have $d^4x = c \, dt \, d^3x$ and therefore
$$S \equiv \int dt \, L \equiv \int dt \int d^3x \, \mathscr L = \int \frac{d^4x}{c} \mathscr L$$
where $L$ is the Lagrangian and $\mathscr L$ is the Lagrangian density. Thus $\mathscr L$ is the quantity inside the square brackets below, and $L$ is the integral of this density over the spatial coordinates.
$S$ must have units of action, which is the same as energy times time. The Lagrangian $L$ therefore must have units of energy, and the Lagrangian density $\mathscr L$ must have units of energy density.
The units of electromagnetic quantities such as $e$, $\phi$, $\mathbf A$, $\mathbf E$, $\mathbf B$, $A^\mu$, and $F^{\mu\nu}$ are the same in QED as in classical electromagnetism. Both the classical theory and the quantum theory can be done in various unit systems. I discuss four of them below.
The Minkowski metric $\eta_{\mu\nu}$ is dimensionless in Cartesian coordinates, and few people do QED calculations in anything other than Cartesian coordinates. In any case, since $\mathscr L$ must be a Lorentz scalar, if $\mathscr L$ has the right dimensions in Cartesian coordinates, it will have the right dimensions in any coordinates.
$$S = \int\frac{d^4x}{c}\left[ -\frac{1}{4\mu_0}F^{\mu\nu}F_{\mu\nu} + \bar\psi ( i \hbar c \gamma^\mu D_\mu - m c^2) \psi \right]$$
$$D_\mu \equiv \partial_\mu + \frac{ie}{\hbar}A_\mu$$
$$x^\mu \equiv (ct, \mathbf x)$$
$$A^\mu \equiv \left( \frac1c \phi, \mathbf A \right)$$
Base units: kilograms, meters, seconds, amperes
Table of $\text{[quantity]} = \text{kg}^a \text{m}^b \text{s}^c \text{A}^d$
$$\begin{array}{c|c|c|c|} \text{quantity} & \text{kg} & \text{m} & \text{s} & \text{A} \\ \hline S & 1 & 2 & -1 & 0 \\ \hline d^4x & 0 & 4 & 0 & 0 \\ \hline c & 0 & 1 & -1 & 0 \\ \hline \mu_0 & 1 & 1 & -2 & -2 \\ \hline F^{\mu\nu} \text{ or } F_{\mu\nu} & 1 & 0 & -2 & -1 \\ \hline\ \psi \text{ or } \bar\psi & 0 & -3/2 & 0 & 0 \\ \hline \hbar & 1 & 2 & -1 & 0 \\ \hline \gamma^\mu \text{ or } \gamma_\mu & 0 & 0 & 0 & 0 \\ \hline m & 1 & 0 & 0 & 0 \\ \hline \partial_\mu \text{ or } \partial^\mu & 0 & -1 & 0 & 0 \\ \hline e & 0 & 0 & 1 & 1 \\ \hline A^\mu \text{ or } A_\mu & 1 & 1 & -2 & -1 \\ \hline g_{\mu\nu} \text{ or } g^{\mu\nu} & 0 & 0 & 0 & 0 \\ \hline \end{array}$$
$$S = \int\frac{d^4x}{c}\left[ -\frac{1}{16\pi}F^{\mu\nu}F_{\mu\nu} + \bar\psi ( i \hbar c \gamma^\mu D_\mu - m c^2) \psi \right]$$
$$D_\mu \equiv \partial_\mu + \frac{ie}{\hbar c}A_\mu$$
$$x^\mu \equiv (ct, \mathbf x)$$
$$A^\mu \equiv \left(\phi, \mathbf A \right)$$
Note: These are sometimes called "unrationalized" Gaussian units. The $F^2$ term has an extra factor of $1/4\pi$; this introduces factors of $4\pi$ into the field equations but removes them from solutions. For example, in these units the classical potential of a point charge is simply $q/r$.
In this system, $\mu_0$ and $\epsilon_0$ do not exist. This is because Gaussian units don't have a electromagnetic base unit; charge is defined in terms of the mechanical base units.
Base units: grams, centimeters, seconds
Table of $\text{[quantity]} = \text{g}^a \text{cm}^b \text{s}^c$
$$\begin{array}{c|c|c|c|} \text{quantity} & \text{g} & \text{cm} & \text{s} \\ \hline S & 1 & 2 & -1 \\ \hline d^4x & 0 & 4 & 0 \\ \hline c & 0 & 1 & -1 \\ \hline F^{\mu\nu} \text{ or } F_{\mu\nu} & 1/2 & -1/2 & -1 \\ \hline\ \psi \text{ or } \bar\psi & 0 & -3/2 & 0 \\ \hline \hbar & 1 & 2 & -1 \\ \hline \gamma^\mu \text{ or } \gamma_\mu & 0 & 0 & 0 \\ \hline m & 1 & 0 & 0 \\ \hline \partial_\mu \text{ or } \partial^\mu & 0 & -1 & 0 \\ \hline e & 1/2 & 3/2 & -1 \\ \hline A^\mu \text{ or } A_\mu & 1/2 & 1/2 & -1 \\ \hline g_{\mu\nu} \text{ or } g^{\mu\nu} & 0 & 0 & 0 \\ \hline \end{array}$$
$$S = \int\frac{d^4x}{c}\left[ -\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + \bar\psi ( i \hbar c \gamma^\mu D_\mu - m c^2) \psi \right]$$
The rest is the same as for Gaussian units.
These are sometimes called "rationalized" Gaussian units. The $F^2$ term has a prefactor of $-1/4$ instead of $-1/16\pi$. This removes factors of $4\pi$ from the field equations.
$$S = \int d^4x \left[ -\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + \bar\psi ( i \gamma^\mu D_\mu - m) \psi \right]$$
$$D_\mu \equiv \partial_\mu + ieA_\mu$$
$$x^\mu \equiv (t, \mathbf x)$$
$$A^\mu \equiv \left(\phi, \mathbf A \right)$$
Base units: GeV
Table of $\text{[quantity]} = \text{GeV}^a$
$$\begin{array}{c|c|} \text{quantity} & \text{GeV} \\ \hline c\equiv 1 & 0 \\ \hline \hbar\equiv 1 & 0 \\ \hline S & 0 \\ \hline d^4x & -4\\ \hline F^{\mu\nu} \text{ or } F_{\mu\nu} & 2 \\ \hline\ \psi \text{ or } \bar\psi & 3/2 \\ \hline \gamma^\mu \text{ or } \gamma_\mu & 0\\ \hline m & 1 \\ \hline \partial_\mu \text{ or } \partial^\mu & 1 \\ \hline e & 0 \\ \hline A^\mu \text{ or } A_\mu & 1 \\ \hline g_{\mu\nu} \text{ or } g^{\mu\nu} & 0\\ \hline \end{array}$$
With a gauge-fixing term added in: $$S = \int \left(-\frac{ε_0 c}{4} g^{μρ} g^{νσ} F_{μν} F_{ρσ} - \frac{ε_0 c}{2}\left(g^{μν}∂_μA_ν\right)^2 + \bar{ψ} \left(γ^μ (iħ \overleftrightarrow{∂_μ} + eA_μ) - mc\right) ψ\right) \sqrt{|g|} d^4 x,$$ where $$ε_0μ_0 = \frac{1}{c^2}, \quad \overleftrightarrow{∂_μ} ≡ \frac{\overrightarrow{∂_μ} - \overleftarrow{∂_μ}}{2}$$ with the arrows denoting the direction that the operator is to be applied in. This is to be understood as passing through the gammas, so it should actually be written as: $$γ^μ iħ \overleftrightarrow{∂_μ} = iħ \frac{γ^μ\overrightarrow{∂_μ} - \overleftarrow{∂_μ}γ^μ}{2}.$$
Dimensions ... set: $$[μ] = \left[x^μ\right], \quad Ω = \left[d^4x\right].$$ Then $$\left[dx^μ\right] = [μ], \quad Ω = [0][1][2][3], \quad \left[∂_μ\right] = \frac{1}{[μ]} = \left[\overleftrightarrow{∂_μ}\right].$$ For the linear-operators: $$[g(\_,\_)] = L^2, \quad \left[g^{-1}(\_,\_)\right] = \frac{1}{L^2},$$ where $L$ denotes the dimension of "length". Therefore $$\left[g_{μν} ≡ g\left(∂_μ, ∂_ν\right)\right] = \frac{L^2}{[μ][ν]}, \quad \left[\sqrt{|g|}\right] = \frac{L^4}{Ω}, \quad \left[g^{μν} ≡ g^{-1}\left(dx^μ, dx^ν\right)\right] = \frac{[μ][ν]}{L^2}.$$
And no: the metric components are not generally dimensionless! Their dimensions are coordinate-dependent. Examples: $$\left[g_{θθ} ≡ g\left(\frac{∂}{∂θ},\frac{∂}{∂θ}\right) = r^2\right] = L^2,$$ for spherical coordinates; $$\left[g_{tt} ≡ g\left(\frac{∂}{∂t},\frac{∂}{∂t}\right)\right] = \frac{L^2}{T^2},$$ since $[θ] = 1$ and $[t] = T$.
And yes: there is a $\sqrt{|g|}$ factor in the Lagrangian density, because it's a scalar density, not a scalar! And its dimension is coordinate-dependent. Example: for spherical coordinates and $x^0 = t$, $Ω = [0][1][2][3] = T·L·1·1 = TL$, and $\left[\sqrt{|g|}\right] = L^4/(TL) = L^3/T$.
For mass and (in vacuuo) light speed, Planck's constant and action, respectively: $$[m] = M, \quad [c] = \frac{L}{T}, \quad [ħ] = H ≡ \frac{ML^2}{T} = [S],$$ where $M$ denotes dimension of mass, $T$ of time, and $H$ of action.
For the charge and electromagnetic field $$[e] = Q, \quad \left[A_μ\right] = \frac{P}{[μ]}, \quad \left[F_{μν}\right] = \frac{P}{[μ][ν]}, \quad \left[ε_0\right] = \frac{QT}{PL}, \quad \left[μ_0\right] = \frac{PT}{QL}, \quad PQ = H,$$ where $Q$ denotes the dimension of electric charge and $P$ the conjugate dimension of magnetic charge/flux.
Finally, for the spinor, the gamma linear operator and its components: $$\left[γ(\_)\right] = \frac{1}{L}, \quad \left[γ^μ ≡ γ\left(dx^μ\right)\right] = \frac{[μ]}{L}, \quad \left[\bar{ψ}\right] = \sqrt{\frac{1}{L^3}} = \left[ψ\right].$$ And, again: the gamma matrices are not generally dimensionless, but are coordinate-dependent; example: $$\left[γ^t ≡ γ\left(dt\right) = \frac{γ^0}{c}\right] = \frac{[dt]}{L} = \frac{T}{L}, \quad \left[γ^θ ≡ γ\left(dθ\right) = \frac{γ^2}{r}\right] = \frac{[dθ]}{L} = \frac{1}{L},$$ since the gamma matrices that you're used to seeing are associated with orthonormal co-frames, e.g. $$γ^0 = γ\left(c dt\right), \quad γ^1 = γ\left(dr\right), \quad γ^2 = γ\left(r dθ\right), \quad γ^3 = γ\left(r \sin{θ} dφ\right),$$ when in spherical coordinates, thus: $$γ^t = \frac{γ^0}{c}, \quad γ^r = γ^1, \quad γ^θ = \frac{γ^2}{r}, \quad γ^φ = \frac{γ^3}{r \sin{θ}}.$$
Finally... the situation with the fractional dimensions for $\bar{ψ}$ and $ψ$ is a clear tell - a hallmark sign - that there's conflation going on, just like with the Gaussian units for electromagnetism, and that things are getting mixed in together and confused that should be kept separate. So, actually, we should be cutting the moorings loose here, too, like we do with the invocation of $Q$ and $P$ for electromagnetism, and set: $$\left[\bar{ψ}\right] = F = \left[ψ\right],$$ where $F$ denotes the newly-invoked dimension of "fermionity". That means modifying the Dirac part of the action to: $$S = \int \left(⋯ + \bar{ψ} \left(ϝ_0 γ^μ (iħ \overleftrightarrow{∂_μ} + eA_μ) - {ϝ'}_0 mc\right) ψ\right) \sqrt{|g|} d^4 x,$$ with two "vacuum fermittivity" coefficients with respective dimensions: $$\left[ϝ_0\right] = \frac{1}{L^3F^2} = \left[{ϝ'}_0\right].$$
Then, the renormalized Lagrangian densities, are of the very same form, except for the replacements: $$ε_0 → Z_3 ε_R, \quad \left(μ_0 → \frac{μ_R}{Z_3}\right), \quad ϝ_0 → Z_1 ϝ_R = Z_2 ϝ_R, \quad {ϝ'}_0 = Z_0 {ϝ'}_R.$$ It is not the fields and charges that are "bare" versus "dressed", but the coefficients themselves. It's the vacuum that's being regularized and renormalized, and it is what suffers the "bare" versus "dressed" distinction.
The $Z$'s, themselves, are the shadows of explicit or tacit coefficients, like $\left(ε_0, μ_0, ϝ_0, {ϝ'}_0\right)$, that is also a tell which shows that there should be coefficients there.
For reference, by the way: $$μ_0 c = «Z_0» = \sqrt{\frac{μ_0}{ε_0}} = \frac{1}{ε_0 c},$$ is the "impedance of free space" and has a more less standard symbol all its own - the $«Z_0»$ I put in scare quotes because it's a different $Z_0$ than the renormalization coefficient: a name conflict. But, I think that might be where the use of $Z$ for the renormalization coefficients themselves actually stems from. The renormalized version of it is: $$«Z_0» = \frac{«Z_R»}{Z_3},$$ and the electromagnetic part of the action is: $$S = \int \left(-\frac{g^{μρ} g^{νσ} F_{μν} F_{ρσ}}{4«Z_0»} - \frac{\left(g^{μν}∂_μA_ν\right)^2}{2«Z_0»} + ⋯\right) \sqrt{|g|} d^4 x.$$
