A two-dimensional in space fractional diffusion equation with functional delay of a general form is considered. For this problem, the Crank-Nicolson method is constructed, based on shifted Grunwald-Letnikov formulas for approximating fractional derivatives with respect to each spatial variable and using piecewise linear interpolation of discrete history with continuation extrapolation to take into account the delay effect. The Douglas scheme is used to reduce the emerging high-dimensional system to tridiagonal systems. The residual of the method is investigated. To obtain the order of the method, we reduce the systems to constructions of the general difference scheme with heredity. A theorem on the second order of convergence of the method in time and space steps is proved. The results of numerical experiments are presented.