In this paper we report upon our computational methodology for numerically integrating the quantum Hamilton-Jacobi equations using hydrodynamic trajectories. Our method builds upon the moving least squares method developed by Lopreore and Wyatt [Phys. Rev. Lett. 82, 5190 (1999)] in which Lagrangian fluid elements representing probability volume elements of the wave function evolve under Newtonian equations of motion which include a nonlocal quantum force. This quantum force, which depends upon the third derivative of the quantum density, rho, can vary rapidly in x and become singular in the presence of nodal points. Here, we present a new approach for performing quantum trajectory calculations which does not involve calculating the quantum force directly, but uses the wavefront to calculate the velocity field using mv=delS, where S/(h) over bar is the argument of the wave function psi. Additional numerical stability is gained by performing local gauge transformations to remove oscillatory components of the wave function. Finally, we use a dynamical Rayleigh-Ritz approach to derive ancillary equations-of-motion for the spatial derivatives of rho, S, and upsilon. The methodologies described herein dramatically improve the long time stability and accuracy of the quantum trajectory approach even in the presence of nodes. The method is applied to both barrier crossing and tunneling systems. We also compare our results to semiclassical based descriptions of barrier tunneling.