An efficient, numerically robust algorithm for calculating acoustic normal modes in a fluid, layered ocean environment is described. Each layer is assumed to have a sound-speed profile for which the mode functions may be expressed analytically in terms of Airy functions. Attenuation is included as a perturbation. The propagator matrix method is used to satisfy the boundary conditions at layer interfaces to obtain a characteristic equation, the zeros of which correspond to the eigenvalues of the trapped normal modes. The algorithms for finding eigenvalues and for computing mode functions avoid the numerical instabilities associated with evanescent fields by removing potentially large exponential terms from the propagator matrices. The likelihood of missing modes in multichannel environments is minimized by finding eigenvalues accurately and sampling the wave-number spectrum finely. The CPU time of the model is shown to increase linearly with frequency and waveguide thickness, whereas the CPU time of models that solve the wave equation numerically in layers a fraction of a wavelength thick increases quadratically. Comparisons with two existing mode models show that the new algorithm can be significantly faster, especially at moderate to high frequencies and for moderate- to deep-water environments.© 1995, Acoustical Society of America. All rights reserved.