A methodology termed ``velocity-scalar filtered mass density function'' (VSFMDF) is developed and implemented for large eddy simulation (LES) of variable density, turbulent reacting flows. In this methodology, the effects of unresolved subgrid scales (SGS) are taken into account by considering the joint probability density function (PDF) of the velocity-scalar field [1]. This is the most general form of the filtered density function for reacting flow simulations. In the VSFMDF transport equation, the effects of SGS convection and chemical reaction are in closed forms. The unclosed terms are modeled in a fashion similar to PDF methods [1]. The modeled VSFMDF transport equation is solved by a Lagrangian Monte Carlo method. The methodology is employed to simulate turbulent shear flows. The predicted results are assessed via comparison with data obtained by direct numerical simulation and laboratory experiments. [1] Pope, S. B., Turbulent flows, Cambridge University Press, Cambridge, UK (2000).