We develop a numerical algorithm for simulation of wave propagation in linear nonisothermal poroelastic media, based on Biot theory and a generalized Fourier law of heat transport in analogy with Maxwell model of viscoelasticity. A plane wave analysis indicates the presence of the classical P and S waves and two slow waves, namely, the Biot and the thermal slow modes of propagation, which present diffusive behavior under certain conditions, depending on viscosity, frequency, and the thermoelastic constants. The wavefield is computed with a direct meshing method using the Fourier differential operator to calculate the spatial derivatives. We propose two alternative time‐stepping algorithms, namely, a first‐order explicit Crank‐Nicolson method and a second‐order splitting method. The Fourier differential operator provides spectral accuracy in the calculation of the spatial derivatives. Modeling the thermal diffusive mode is relevant for high‐temperature high‐pressure fields and since it leads to mesoscopic attenuation by mode conversion of the fast waves to the thermal waves.