  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.
