  In order to describe forward hadron productions in high-energy nuclear
collisions, we propose a Monte-Carlo implementation of
Dumitru-Hayashigaki-Jalilian-Marian formula with the unintegrated gluon
distribution obtained numerically from the running-coupling BK equation. We
discuss influence of initial conditions for the BK equation by comparing a
model constrained by global fit of small-x HERA data and a newly proposed one
from the running coupling MV model.
