We study quasi-spherical, supersonic accretion flows around black holes using high-accuracy numerical simulations. We describe a code, the Lagrangian Total Variation Diminishing (TVD), and a remap routine to address a specific issue in the Advection Dominated Accretion Flow (ADAF) that is, appropriately handling the angular momentum even near the inner boundary. The Lagrangian TVD code is based on an explicit finite difference scheme on mass-volume grids to track fluid particles with time. The consequences are remapped on fixed grids using the explicit Eulerian finitedifference algorithm with a third-order accuracy. Test results show that one can successfully handle flows and resolve shocks within two to three computational cells. Especially, the calculation of a hydrodynamical accretion disk without viscosity around a black hole shows that one can conserve nearly 100% of specific a ngular momentum in one-and twodimensional cylindrical coordinates. Thus, we apply this code to obtain a numerically similar ADAF solution. We perform simulations, including viscosity terms in one-dimensional spherical geometry on the non-uniform grids, to obtain greater quantitative consequences and to save computational time. The error of specific angular momentum in Newtonian potential is less than 1% between r~10rs and r~10 4 rs, where rs is sink size. As Narayan et al. (1997) suggested, the ADAFs in pseudo-Newtonian potential become supersonic flows near the black hole, and the sonic point is rsonic~5.3rg for flow with α =0.3 and  =1 .5. Such simulations indicate that even the ADAF with  =5/3 is differentially rotating, as Ogilvie (1999) indicated. Hence, we conclude that the Lagrangian TVD and remap code treat the role of viscosity more precisely than the other scheme, even near the inner boundary in a rotating accretion flow around a nonrotating black hole.