Quarkonium production in high-energy proton (deuteron)–nucleus collisions is investigated in the color glass condensate framework. We employ the color evaporation model assuming that the quark pair produced from dense small-x gluons in the nuclear target bounds into a quarkonium outside the target. The unintegrated gluon distribution at small Bjorken x in the nuclear target is treated with the Balitsky–Kovchegov equation with running coupling corrections. For the gluons in the proton, we examine two possible descriptions, unintegrated gluon distribution and ordinary collinear gluon distribution. We present the transverse momentum spectrum and nuclear modification factor for J/ψ production at RHIC and LHC energies, and those for ϒ(1S) at LHC energy, and discuss the nuclear modification factor and the momentum broadening by changing the rapidity and the initial saturation scale.