AbstractThe numerical treatment of an atmospheric chemical scheme, which contains 56 species, is discussed in this paper. This scheme is often used in studies of air pollution levels in different domains, as, for example, in Europe, by large-scale environmental models containing additionally two other important physical processes—transport of pollutants in the atmosphere (advection) and diffusion phenomena. We shall concentrate our attention on the efficient numerical treatment of the chemical scheme by using Implicit Runge–Kutta Methods combined with accurate and efficient advanced versions of the Richardson Extrapolation. A Variable Stepsize Variable Formula Method is developed in order to achieve high accuracy of the calculated results within a reasonable computational time. Reliable estimations of the computational errors when the proposed numerical methods are used in the treatment of the chemical scheme will be demonstrated by presenting results from several representative runs and comparing these results with “exact” concentrations obtained by applying a very small stepsize during the computations. Results related to the diurnal variations of some of the chemical species will also be presented. The approach used in this paper does not depend on the particular chemical scheme and can easily be applied when other atmospheric chemical schemes are selected.