Motivated by the recent work of Brzeminski, Motyka, Sadzikowski and Stebel in arXiv:1611.04449, where forward Drell--Yan production is studied in proton-proton collisions at the LHC, we improve their calculation by introducing an unintegrated gluon density obtained in arXiv:1209.1353 and arXiv:1301.5283 from a fit to combined HERA data at small values of Bjorken $x$. This gluon density was calculated within the BFKL formalism at next-to-leading order with collinear corrections. We show that it generates a good description of the forward Drell--Yan cross section dependence on the invariant mass of the lepton pair both for LHCb and ATLAS data.