The squirt-flow wave attenuation mechanism is implemented in Biot’s theory of poroelasticity in the form of differential equations. All the stiffnesses involved in the stress-strain relation become complex and frequency dependent, which can exactly be expressed in terms of kernels based on the Zener mechanical model. In the time domain, this approach implies time convolutions, which are circumvented by introducing memory variables. The differential equations are consistent with Gassmann’s and Mavko-Jizba equations at low and high frequencies, respectively. All the coefficients in the poro-viscoelastic differential equations have a clear physical meaning and can be obtained or estimated from independent measurements. The key additional parameters are the dry-rock bulk modulus at a confining pressure where all the compliant pores are closed, i.e., a hypothetical rock without the soft porosity, the grain-contact aspect ratio and the compliant porosity. We recasted the wave equation in the particle-velocity/stress formulation and solved it by using a time-splitting technique and the Fourier pseudospectral method to compute the spatial derivatives. The algorithm can be used to obtain synthetic wave fields in inhomogeneous media.