A conservative discrete velocity method (DVM) is developed for the ellipsoidal Fokker-Planck (ES-FP) equation in prediction of nonequilibrium neutral gas flows in this paper. The ES-FP collision operator is solved in discrete velocity space in a concise and quick finite difference framework. The conservation problem of the discrete ES-FP collision operator is solved by multiplying each term in it by extra conservative coefficients whose values are very close to unity. Their differences to unity are in the same order of the numerical error in approximating the ES-FP operator in discrete velocity space. All the macroscopic conservative variables (mass, momentum, and energy) are conserved in the present modified discrete ES-FP collision operator. Since the conservation property in a discrete element of physical space is very important for the numerical scheme when discontinuity and a large gradient exist in the flow field, a finite volume framework is adopted for the transport term of the ES-FP equation. For nD-3V (n<3) cases, a nD-quasi nV reduction is specifically proposed for the ES-FP equation and the corresponding FP-DVM method, which can greatly reduce the computational cost. The validity and accuracy of both the ES-FP equation and FP-DVM method are examined using a series of 0D-3V homogenous relaxation cases and 1D-3V shock structure cases with different Mach numbers, in which 1D-3V cases are reduced to 1D-quasi 1V cases. Both the predictions of 0D-3V and 1D-3V cases match well with the benchmark results such as the analytical Boltzmann solution, direct full-Boltzmann numerical solution, and DSMC result. Especially, the FP-DVM predictions match well with the DSMC results in the Mach 8.0 shock structure case, which is in high nonequilibrium, and is a challenge case of the model Boltzmann equation and the corresponding numerical methods.