Quick answer:
For a critical spin chain of length $L$, the formula for the entanglement entropy of an interval $S(l) = \frac{c}{3}\ln \left(\frac{L}{\pi} \sin{\frac{\pi l}{L} }\right) +C$ for an interval of length $l$ is only valid for large intervals. It's simply inaccurate* (though see below) at small $l$, even when you include the constant. And at large $l$ where the formula holds, the constraint $S(l) \leq l \log(d)$ is easily satisfied, so there's no problem.
General correspondence between lattice model and CFT
Stepping back, it's good to remember that when we say a critical spin chain is described by a CFT, usually we just mean that the long-distance properties of the spin chain match those of the CFT.
For instance, when we say a spin chain correlation function $\langle O(x) O(y)\rangle$ (for a lattice operator $O$) matches a CFT correlation function $\langle \phi(x) \phi(y)\rangle$ (for some CFT field $\phi$), this relation only holds at long distance $|x-y|$. Here's one reason we can't expect a match at short distances: the CFT correlation function $\langle \phi(x) \phi(y)\rangle$ will generally blow up as $x \to y$, whereas on a spin chain, the correlation function $\langle O(x) O(y)\rangle$ is bounded even as $x \to y$. Instead, they match at long distance, where they have the same power-law decay. One way to express this is $\langle \phi(0) \phi(r)\rangle \propto \lim_{\lambda \to \infty} \langle O(0) O(\lambda r) \rangle \lambda^{2\Delta}$ for some scaling dimension $\Delta>0$. On the other hand, when $|x-y|$ is just a few lattice spacings, the spin chain correlation function $\langle O(x) O(y)\rangle$ and the CFT correlation function $\langle \phi(x) \phi(y)\rangle$ may bear no particular relation.
Just like a correlation function, the function $S(l)$ at small $l$ will depend on details of the lattice model. For instance, you could have two different critical spin chain Hamiltonians described by the same CFT in the IR, and their correlation functions and entropies at small length scales could be very different, even though they have the same universal behavior at long distance.
Addendum 1: Accuracy of CFT formula?
[*] I said the formula $S(l) = \frac{c}{3}\ln \left(\frac{L}{\pi} \sin{\frac{\pi l}{L} }\right) +C$ is simply wrong for small $l$. That might seem surprising if you look at the figure in the question you linked, where the entropies are calculated for a spin chain of length 12, and even for regions of size $l=1,2,\dots$ the fit looks perfect. Apparently, for this particular critical spin chain Hamiltonian, the long-distance scaling of entropy already provides a very accurate approximation at short distances. See Fig. 1 of Calabrese et al. where they plot the error between an exact spin chain calculation of $S(l)$ and the CFT formula above; in their particular case they also find a very good match at short distances. Maybe that's somehow related to the exact solvability of these particular lattice Hamiltonians?
If you chose a more generic lattice Hamiltonian that exhibited the same behavior in the IR, you might not be so lucky to have a high-accuracy match with the $S(l)$ formula at short distances.
Addendum 2: Upper bound on central charge?
Thanks to Brandon Rayhaun for discussion; we have wondered about this question several times actually.
There was a suggestion in the comments for another possible resolution: maybe there's a universal upper bound for the central charge $c$ of any CFT that can be realized by range-$r$ Hamiltonians on spin chains of local dimension $d$, with the upper bound depending only on ($r,d$). Well, such an upper bound is no longer needed to resolve the original question, but I will try to comment anyway.
First, any such bound would need to depend on both $r$ and $d$, since you can always coarse-grain the qudits so that they only have $r=2$ nearest-neighbor interactions, at the price of increasing $d$. Second, even if you restricted to e.g. nearest-neighbor qubit interactions, I think you would need to impose translation-invariance (under single-qubit translations). Otherwise, I think there are tricks to "simulate" more general 1D Hamiltonians (with longer-range interactions, or higher-dimensional qudits) using only nearest-neighbor qubit spin chains. (Such tricks are discussed by Piddock et al. in the 2D case here.) A more complicated CFT might then appear in the IR of a simple nearest-neighbor qubit chain.
However, restricting finally to nearest-neighbor translation-invariant qubit chains, there may well be an upper bound on the possible central charge of CFTs in the IR. The parameter space of such Hamiltonians (up to additive and multiplicative scaling, and some other redundancies, see e.g. below Eq. 5.8 here) is just a compact 9-dimensional real manifold. Imagine the set of critical spin chains (say we mean: gapless spin chains corresponding to some CFT in the IR) as a subset of this compact parameter space. If there were indeed nearest-neighbor qubit chains corresponding to CFTs with arbitrarily large central charge, then there would be a sequence of critical spin chains with central charges $c_1,c_2,\ldots$ diverging to infinity, and the corresponding sequence in the compact parameter space would need to have an accumulation point.
So if there's no upper bound on the central charge of a CFT in the IR of a translation-invariant nearest-neighbor qubit chain, there would need to be some nearest-neighbor chain [the above accumulation point] such that, with arbitrarily small perturbations, you could tune the central charge of the IR CFT to vary by arbitrary amounts. Maybe that's possible, I don't know. But since it seems strange, it seems to me good evidence for having an upper bound instead. I don't know that anyone has shown such an upper bound, or better yet, classified the space of CFTs that can be realized in the IR of translation-invariant spin chain Hamiltonians of fixed range and local dimension.
