  We present the scale setting for a new set of gauge configurations generated
with $N_f=3+1$ Wilson quarks with a non-perturbatively determined clover
coefficient in a massive O($a$) improvement scheme. The three light quarks are
degenerate, with the sum of their masses being equal to its value in nature and
the charm quark has its physical mass. We use open boundary conditions in time
direction to avoid the problem of topological freezing at small lattice
spacings and twisted-mass reweighting for improved stability of the
simulations. The decoupling of charm at low energy allows us to set the scale
by measuring the value of the low-energy quantity $t_0^\star/a^2$, which is the
flow scale $t_0$ at our mass point, and comparing it to an $N_f=2+1$ result in
physical units. We present the details of the algorithmic setup and tuning
procedure and give the bare parameters of ensembles with two lattice spacings
a=0.054 fm and a=0.043 fm. We discuss finite volume effects and lattice
artifacts and present physical results for the charmonium spectrum. In
particular the hyperfine splitting between the $\eta_c$ and $J/\psi$ mesons
agrees very well with its physical value.
