Prediction of elastic full wavefields is required for reverse time migration, full waveform inversion, borehole seismology, seismic modelling, etc. We propose a novel algorithm to solve the Navier wave equation, which is based on multi‐block methodology for high‐order finite‐difference schemes on curvilinear grids. In the current implementation, the blocks are subhorizontal layers. Smooth anisotropic heterogeneous media in each layer can have strong discontinuities at the interfaces. A curvilinear adaptive hexahedral grid in blocks is generated by mapping the original 3D physical domain onto a parametric cube with horizontal layers and interfaces. These interfaces correspond to the main curvilinear physical contrast interfaces of a subhorizontally layered formation. The top boundary of the parametric cube handles the land surface with smooth topography. Free‐surface and solid–solid transmission boundary conditions at interfaces are approximated with the second‐order accuracy. Smooth media in the layers are approximated up to sixth‐order spatial schemes. All expected properties of the developed algorithm are demonstrated in numerical tests using corresponding parallel message passing interface code.