The Helmholtz Equation (- Δ - K2n2)u = 0 with a variable index of refraction, n, and a suitable radiation condition at infinity serves as a model for a wide variety of wave propagation problems. A numerical algorithm has been developed and a computer code implemented that can effectively solve this equation in the intermediate frequency range. The equation is discretized by using the finite element method, thus allowing for the modeling of complicated geometries (including interfaces) and complicated boundary conditions. A global radiation boundary condition is imposed at the far-field boundary that is exact for an arbitrary number of propagating modes. The resulting large, nonselfadjoint system of linear equations with indefinite symmetric part is solved by using the preconditioned conjugate-gradient method applied to the normal equations. A new preconditioner based on the multigrid method is developed. This preconditioner is vectorizable and is extremely effective over a wide range of frequencies, provided the number of grid levels is reduced for large frequencies. A heuristic argument is given that indicates the superior convergence properties of this preconditioner. The relevant limit to analyze convergence is for K increasing and a fixed prescribed accuracy level. The efficiency and robustness of the numerical algorithm are confirmed for large acoustic models, including interfaces with strong velocity contrasts.