  We cure the perturbative instability of the total-inclusive-photoproduction
cross sections of vector $S$-wave quarkonia observed at high
photon-proton-collision energies ($\sqrt{s_{\gamma p}}$) in Next-to-Leading
Order (NLO) Collinear-Factorisation (CF) computations. This is achieved using
High-Energy Factorisation (HEF) in the Doubly-Logarithmic Approximation (DLA),
which is a subset of the Leading-Logarithmic Approximation (LLA) of HEF which
resums higher-order QCD corrections proportional to $\alpha_s^n \ln^{n-1}
(\hat{s}/M^2)$ in the Regge limit $\hat{s}\gg M^2$ with $M^2$ being the
quarkonium mass and $\hat{s}$ is the squared partonic-center-of-mass energy.
Such a DLA is strictly consistent with the NLO and NNLO DGLAP evolutions of the
Parton Distribution Functions. By improving the treatment of the
large-$\hat{s}$ asymptotics of the CF coefficient function, the resummation
cures the unphysical results of the NLO CF calculation. The matching is
directly performed in $\hat{s}$ space using the Inverse-Error Weighting
matching procedure which avoids any possible double counting. The obtained
cross sections are in good agreement with data. In addition, the
scale-variation uncertainty of the matched result is significantly reduced
compared to the LO results. Our calculations also yield closed-form analytic
limits for $\hat{s}\gg M^2$ of the NLO partonic CF and numerical limits for
contributions to those at NNLO scaling like $\alpha_s^2 \ln(\hat{s}/M^2)$.
