In this paper, the Approximate Deconvolution Method (ADM) is applied for modeling of sub-grid terms in Large Eddy Simulation of turbulent multi-species flows. ADM is based on a repetitive filtering according to van Cittert method, which is used to approximate an inverse filter operator that enables the computation of unfiltered variables. Regarding the scalar field, ADM has some drawbacks as it does not ensure the essential conservation and boundedness conditions. We overcome this problem by applying normalization and regularization methods. ADM is applied for 2D vorticity-stream function formulation of the Navier-Stokes equations and scalar transport equations. The test cases include a double jet and mixing layer configuration characterized by significant shear stresses causing the development of small-scale phenomena. The computations are performed applying various explicit and compact discretization schemes with the approximation order up to 20th for the vorticity-stream function equation and Total Variation Diminishing method for the convective terms of scalar transport equations. The inverse filter is based on a top-hat filter and we analyze how a particular discretization scheme and the inverse filter approximation accuracy affect the solutions. The results obtained using ADM are compared with the results of the classical LES with the Smagorinsky model and the DNS data. It is shown that ADM provides accurate solutions both in terms of instantaneous local quantities as well as spatially computed time-varying statistical moments. It is found that the accuracy of ADM is significantly dependent on the method of calculation of the inverse filter and the discretization method only when small turbulent scales are present in the flow field.