Some mathematical models of interest, e.g., for Meteorology, can be formulated in terms of diffusion equations with time and/or space fractional derivatives. The usual time derivative can be replaced, for instance, by the so-called Caputo fractional derivative (of order gamma epsilon (0, 1)), while the space derivatives can be written as a Riemann-Liouville fractional derivatives (of order alpha epsilon (1, 2)). In this paper, we implement third-order accurate in time numerical algorithms to solve two-dimensional fractional diffusion equations. These are new finite difference schemes, based on the Grunwald-Letnikov difference operator and some ADI methods, combined with an optimized extrapolation strategy. Numerical examples, concerning model-problems as well as real-world applications, are given.