The majority of the published two- and three-dimensional simulations of the interstellar medium, including the disk-halo interaction, are set up under the assumption of collisional ionization equilibrium (CIE). However, even disregarding the dynamics of the medium, that is, letting the plasma cool isochorically or isobarically, the rates of collisional ionization and radiative recombination do not balance, so that as time goes by, the ionization structure is severely driven out of equilibrium. Consequently, the cooling rates, which depend on the presence of certain ions at a given time, vary as a function of time. Taking into account the dynamics of the plasma, the imbalance becomes even more severe, because adiabatic compression or expansion can change the kinetic temperature of the electrons significantly, while leaving the ionization structure at first unchanged. Hence, ionization and recombination can be delayed, if for example a shock or
a rarefaction wave passes over a certain fluid element, respectively. It is important to realize that the dynamical and thermal history of the plasma are inextricably intertwined. For instance the dynamics depends on the heating and cooling history and vice versa.
Recently, owing to the development of the Atomic+Molecular Plasma Emission Code (EA+MPEC) and its coupling to a PPM based AMR code, it has been possible to carry out multi-fluid calculations of the ISM tracing both the thermal and dynamical evolutions of the gas self-consistently. A far reaching consequence of these investigations is the insight that cooling is both spatially and temporally dependent, and can differ by more than an order of magnitude from CIE. Other results from these NEI simulations are: (i) Hitopology derived from the genus similar to that determined for the IGALFA survey data, (ii) an enhanced electron density, especially at lower temperatures, which can very well explain observed pulsar dispersion measures up to 4 kpc in the Galactic plane, (iii) an increased concentration of Ovi (up to 70% of the total mass) in
thermally unstable temperature regimes, and (iv) a more realistic distribution of Li-like ions in the local ISM, allowing a better constraint of the age of the Local Bubble. In this talk a review of these developments (with an eye towards consistency of atomic and molecular data, validation of simulations against observations, acceleration of calculations with GPUs, and simulations results) is presented.