In this work, we propose and test the validity of a modified Phase Field Method (PFM), which was specifically developed for large scale simulations of turbulent flows with large and deformable surfactant-laden droplets. The time evolution of the phase field, φ, and of the surfactant concentration field, ψ, are obtained from two Cahn-Hilliard-like equations together with a two-order-parameter Time-Dependent Ginzburg-Landau (TDGL) free energy functional. The modifications introduced circumvent existing limitations of current approaches based on PFM and improve the well-posedness of the model. The effect of surfactant on surface tension is modeled via an Equation Of State (EOS), further improving the flexibility of the approach. This method can efficiently handle topological changes, i.e. breakup and coalescence, and describe adsorption/desorption of surfactant. The capabilities of the proposed approach are tested in this paper against previous experimental results on the effects of surfactant on the deformation of a single droplet and on the interactions between two droplets. Finally, to appreciate the performances of the model on a large scale complex simulation, a qualitative analysis of the behavior of surfactant-laden droplets in a turbulent channel flow is presented and discussed.