Implementation of a variable stepsize variable formula method in the time-integration part of a code for treatment of long-range transport of air pollutants. A system of two partial differential equations is used for modelling of long-range transport of pollutants over Europe. The first order space derivatives (the advection terms) are discretized by a pseudospectral algorithm. A discretization of the second order spatial derivatives (the diffusion terms) is carried out by a special technique developed. The arising large systems of ordinary differential equations are integrated by a variable stepsize variable formula method, which is based on the predictor-corrector scheme. The stepsize selection strategy and the formula selection strategy are explained in detail. An accuracy control and a stability control are performed for each time step. Numerical experiments with real meteorological data prove the great efficiency of the applied method and the reliability of the obtained results.