We use a new approach to color evaporation model (CEM) for quarkonium production. The production of $c\bar c$ pairs is performed within $k_T$-factorization approach using different unintegrated gluon distribution functions (UGDF) from the literature. We include all recent improvements to color evaporation model. We get poor description of the large transverse momentum distributions of $J/\psi$ with the JH-2013 CCFM-based UGDF. Here explicit inclusion of $2 \to 3$ processes considerably improves the situation. Similar effects are discussed in the context of the KMR UGDF.