  We present a formulation of the initial state of heavy ion collisions that
generalizes the 2+1D boost invariant IP-Glasma \cite{Schenke:2012wb} to 3+1D
through JIMWLK rapidity evolution of the pre-collision Wilson lines. The
rapidity dependence introduced by the JIMWLK evolution leads us to modify the
initial condition for the gauge fields, and to solve Gauss' law iteratively in
order to allow for temporal evolution on a 3-dimensional lattice.
  While the transverse physics of QGP has been studied nearly exhaustively, the
effect of longitudinal fluctuations introduced by the JIMWLK evolution has yet
to be studied in detail phenomenologically. Hence, we couple our 3+1D IP-Glasma
model to MUSIC+UrQMD, for completely 3+1D simulations of heavy ion collisions.
Specifically, we consider Pb-Pb collisions at $\sqrt{s} = 2.76\, {\rm TeV}$ and
study the rapidity dependence of the charged hadron $v_n(\eta)$ via the
$\eta$-dependent flow factorization ratios $r_n(\eta_a,\eta_b)$ as measured by
CMS \cite{Khachatryan:2015oea}, as well as the charged hadron multiplicity
$dN_{ch}/d\eta$.
