The Reynolds-averaged equations for turbulent flow over a deep-water sinusoidal gravity wave, z = a cos kx = h(0)(x), are formulated in the wave-following coordinates xi, eta, where, x = xi, z = eta + h(xi, eta), h(xi, 0) = h(0)(xi) and h is exponentially small for k eta much greater than 1, and closed by a viscoelastic constitutive equation (a mixing-length model with relaxation). This closure is derived from Townsend's boundary-layer-evolution equation on the assumptions that: the basic velocity profile is logarithmic in eta + z(0), where z(0) is a roughness length determined by Charnock's similarity relation; the lateral transport of turbulent energy in the perturbed flow is negligible; the dissipation length is proportional to eta + z(0). A counterpart of the Orr-Sommerfeld equation for the complex amplitude of the perturbation stream function is derived and used to construct a quadratic functional for the energy transfer to the wave. A corresponding Galerkin approximation that is based on independent variational approximations for outer (quasi-laminar) and inner (shear-stress) domains yields an energy-transfer parameter beta that is comparable in magnitude with that of the quasi-laminar model (Miles 1957) and those calculated by Townsend (1972) and Gent & Taylor (1976) through numerical integration of the Reynolds-averaged equations. The calculated limiting values of beta for very slow waves, with Charnock's relation replaced by kz(0) = constant, are close to those inferred from observation but about three times the limiting values obtained through extrapolation of Townsend's results.