A correction for the vertical gradient of air density has been incorporated into a skewed probability density function formulation for turbulence in the convective boundary layer. The related drift term for Lagrangian stochastic dispersion modelling has been derived based on the well-mixed condition. Furthermore, the formulation has been extended to include unsteady turbulence statistics and the related additional component of the drift term obtained. These formulations are an extension of the drift formulation reported by Luhar et al. (Atmos Environ 30:1407-1418, 1996 following the well-mixed condition proposed by Thomson (J Fluid Mech 180:529-556, 1987). Comprehensive tests were carried out to validate the formulations including consistency between forward and backward simulations and preservation of a well-mixed state with unsteady conditions. The stationary state CBL drift term with density correction was incorporated into the FLEXPART and FLEXPART-WRF Lagrangian models, and included the use of an ad hoc transition function that modulates the third moment of the vertical velocity based on stability parameters. Due to the current implementation of the FLEXPART models, only a steady-state horizontally homogeneous drift term could be included. To avoid numerical instability, in the presence of non-stationary and horizontally inhomogeneous conditions, a re-initialization procedure for particle velocity was used. The criteria for re-initialization and resulting errors were assessed for the case of non-stationary conditions by comparing a reference numerical solution in simplified unsteady conditions, obtained using the non-stationary drift term, and a solution based on the steady drift with re-initialization. Two examples of "real-world" numerical simulations were performed under different convective conditions to demonstrate the effect of the vertical gradient in density on the particle dispersion in the CBL.