The elastic Helmholtz equation is capable of readily simulating attenuation and dispersion behaviors of the elastic wave and performing full-wavefield modeling in wave-equation-based elastic inversions and migrations. However, solving the elastic Helmholtz equation using a finite-difference frequency-domain (FDFD) method is computationally prohibitive especially in heterogeneous media with fine-scale heterogeneities. The FDFD method usually leads to a large discrete linear system of the elastic Helmholtz equation. We develop a multiscale method of FDFD to solve the elastic Helmholtz equation in isotropic media based on the general framework of heterogeneous multiscale method (HMM). The HMM framework decomposes the elastic Helmholtz problem into a series of microscale problems and a macroscale problem. The idea of multiscale basis functions is introduced to decouple the coupled microscale and macroscale problems and to capture fine-scale heterogeneity in medium properties. A reconstruction-based downscaling coupling and a flux-based upscaling coupling are used to convey the fine-scale medium heterogeneity to a coarse scale. The dimension of the resulting linear system is much smaller than those of linear systems generated with the conventional FDFD methods. We use a homogeneous model and a heterogeneous model to investigate the effects of the size of local sampling domains and the coarse-element number per S-wave wavelength on the accuracy of our new method, and employ two highly heterogeneous models to demonstrate the superiority in terms of the efficiency and memory consumption of our method based on the optimal local sampling-domain size and stable coarse-mesh discretization. The results demonstrate that our new method can approximate the fine-scale reference FDFD solutions with a significant decrease in computational complexity.