Parameterization of gravity waves due to subgrid-scale orography is now included in most existing large-scale models of the atmosphere. Parameterization schemes, however, have so far been evaluated mainly in view of the overall performance of the large-scale models. This may lead to an inappropriate assessment of the schemes since errors from various sources may interact with one another. To avoid this situation, an approach is taken in which a numerical model that explicitly resolves gravity waves is used to evaluate the performance of the schemes. For this purpose, a mesoscale two-dimensional nonlinear anelastic nonhydrostatic model is developed and used to numerically simulate gravity waves for a variety of orographic conditions. Regarding a subdomain bf the mesoscale model as the horizontal grid interval of a large-scale model; two vertical profiles of gravity wave drag are compared-one for the subdomain-averaged values of the drag simulated by the mesoscale model and the other for the drag calculated by a parameterization scheme applied to the subdomain-averaged variables. A test parameterization scheme is constructed by adopting the essential features of the existing schemes. An extensive evaluation of the test parameterization scheme with the aid of the dataset obtained from the mountain wave simulations shows that the scheme does not properly treat the enhancement of drag due to low-level wave breaking through the resonant amplification of nonhydrostatic waves. The authors show that the standard deviation of orography and the tuning coefficient in the scheme alone are not sufficient fdr properly representing this effect. The authors discuss the approach taken to overcome this deficiency by including additional statistical information on subgrid-scale orography in the input to the parameterization. A revised parameterization scheme constructed following this approach is presented.