Abstract A method for the automatic calculation of a horizontal trajectory within the boundary layer is presented. The advecting velocity is based on the surface geostrophic wind. The initial data are routine surface observations, which are objectively analysed by means of orthogonal polynomials. The trajectory equations are then integrated numerically using a variable time-step, fourth-order scheme. An estimate of the global error in the calculation is obtained by performing a number of integrations with randomly perturbed velocity fields. Comparison with a lower order, fixed time-step integration method shows that in certain cases very large errors can be produced by using a fixed time-step of 6 h, and a very small time-step is necessary to make these methods as accurate as the variable time-step method.