A new ion transport code for planetary ionospheric studies has been developed with consideration of velocity differences among ion species involving ion-ion collision. Most of previous planetary ionosphere models assumed that ions diffuse through non-moving ion and neutral background in order to consolidate continuity and momentum equations for ions into a simple set of diffusion equations. The simplification may result in unreliable density profiles of ions at high altitudes where ion velocities are fast and their velocity differences are significant enough to cause inaccuracy when computing ion-ion collision. A new code solves explicitly one-dimensional continuity and momentum equations for ion densities and velocities by utilizing divided Jacobian matrices in matrix inversion necessary to the Newton iteration procedure. The code has been applied to Martian nightside ionosphere models, as an example computation. The computed density profiles of $O^+,;OH^+$, and $HCO^+$ differ by more than a factor of 2 at altitudes higher than 200 km from a simple diffusion model, whereas the density profile of the dominant ion, $O_2^+$, changes little. Especially, the density profile of $HCO^+$ is reduced by a factor of about 10 and its peak altitude is lowered by about 40 km relative to a simple diffusion model in which $HCO^+$ ions are assumed to diffuse through non-moving ion background, $O_2^+$. The computed effects of the new code on the Martian nightside models are explained readily in terms of ion velocities that were solved together with ion densities, which were not available from diffusion models. The new code should thus be expected as a significantly improved tool for planetary ionosphere modelling.