To permit understanding and improvement of subgrade load bearing performance, this paper presents the analysis of a three-dimensional (3D), half-space, elastic subgrade coupling of multi-fields using normal mode analysis (NMA) under Ezzat's fractional order generalized theory of thermoelasticity. Using this method, the analytical equation can be divided into two parts, and the analysis is done with neither integral transformation nor inverse transformation, thereby increasing the speed of decoupling and eliminating the limitation of numerical inverse transformation. The influence of fractional order parameter, load frequency, and permeability coefficient on the thermo-hydro-mechanical dynamic (THMD) coupling media is analyzed for two loading protocols, including thermal impact and normal load. The distributions of the dimensionless vertical displacement, temperature, stress, excess pore water pressure, and strain in the porous structure are analyzed and discussed in detail. The results of the analysis show that the load frequency plays an important role in determining the result values of all dimensionless physical variables. The fractional order parameter only has a significant effect on the calculations of physical variables of the upper surface subjected to thermal impact, and the permeability coefficient mainly affects the stress, strain, and excess pore water pressure. The developed computational method can improve the decoupling speed without any transformation, and it eliminates the limitation of numerical inverse transformation. The method has wide potential application in many engineering fields.