We present a theoretical and numerical framework -- which we dub the Direct Integration Method (DIM) -- for simple, efficient and accurate evaluation of surface wave models allowing presence of a current of arbitrary depth dependence, and where bathymetry and ambient currents may vary slowly in horizontal directions. On horizontally constant water depth and shear current the DIM numerically evaluates the dispersion relation of linear surface waves to arbitrary accuracy, and we argue that for this purpose it is superior to two existing numerical procedures: the piecewise-linear approximation and a method due to textit{Dong & Kirby} [2012]. The DIM moreover yields the full linearized flow field at little extra cost. We implement the DIM numerically with iterations of standard numerical methods. The wide applicability of the DIM in an oceanographic setting in four aspects is shown. Firstly, we show how the DIM allows practical implementation of the wave action conservation equation recently derived by textit{Quinn et al.} [2017]. Secondly, we demonstrate how the DIM handles with ease cases where existing methods struggle, i.e. velocity profiles $mathbf{U}(z)$ changing direction with vertical coordinate $z$, and strongly sheared profiles. Thirdly, we use the DIM to calculate and analyse the full linear flow field beneath a 2D ring wave upon a near--surface wind--driven exponential shear current, revealing striking qualitative differences compared to no shear. Finally we demonstrate that the DIM can be a real competitor to analytical dispersion relation approximations such as that of textit{Kirby & Chen} [1989] even for wave/ocean modelling.