In this work we present a novel method for the solution of gear contact problems in flexible multibody. These problems are characterized by significant variation in the location and size of the contact area, typically requiring a high number of degrees of freedom to correctly capture deformation and stress fields. Therefore fully dynamic simulation is computationally prohibitive. To overcome these limitations, we exploit a combined analytic-numerical contact model within a Parametric Model Order Reduction (PMOR) scheme. The reduction space consists of a truncated set of eigenvectors augmented with a parameter dependent set of residual static shape vectors. Each static shape is computed by interpolating among a set of displacement modes of the interacting bodies, obtained from a series of precomputed static contact analyses. During the contact analyses, an analytic model based on the Hertz theory describes the teeth local deformation. We implement the proposed method in an in-house code and we apply it to spur and helical gears dynamic contact analyses. We compare the results with classical PMOR schemes highlighting how the combined use of the semianalytic contact model allows to decrease further the model complexity as well as the computational burden, for both static and dynamic cases. Finally, we validate the methodology by means of a comparison with experimental data found in literature, showing that the numerical method is able to capture quantitatively the static transmission error measurements in case of both helical and spur geared transmission for different torque levels.