A coupled mode model is presented for the propagation of acoustic-gravity waves in layered ocean waveguides. The analysis extends previous work for acoustic waves in inhomogeneous environment. The coupled mode system is derived by means of a variational principle in conjunction with local mode series expansion, obtained by utilizing eigenfunction systems defined in the vertical section. These are obtained through the solution of vertical eigenvalue problems formulated along the waveguide. A crucial factor is the inclusion of additional modes accounting for the effects of spatialy varying boundaries and interfaces. This enhancement provides an implicit summation for the slowly convergent part of the localmode series, rendering the series rapidly convergent, increasing substantialy the efficiency of the method. Particular aspects of the method include high order Lagrange Finite Element Methods for the solution of local vertical eigenvalue problems in the case of multilayered waveguides, and Gauss-type quadrature for the computation of the coupled-mode system coefficients. The above aspects make the present method quite efficient for long range propagation in extended waveguides, such as the ones found in geophysical applications, e.g. ocean basins, as only few modes are needed for the accurate representation of the wave field.