  We propose a novel family of asymptotically stable, implicit-explicit,
adaptive, time integration method (denoted with the $\theta$-method) for the
solution of the fractional advection-diffusion-reaction (FADR) equations. This
family of time integration method generalized the computationally explicit
$L_1$-method adopted by Brunner (J. Comput. Phys. {\bf 229} 6613-6622 (2010))
as well as the fully implicit method proposed by Jannelli (Comm. Nonlin. Sci.
Num. Sim., {\bf 105}, 106073 (2022)). The spectral analysis of the method
(involving the group velocity and the phase speed) indicates a region of
favorable dispersion for a limited range of Peclet number. The numerical
inversion of the coefficient matrix is avoided by exploiting the sparse
structure of the matrix in the iterative solver for the Poisson equation. The
accuracy and the efficacy of the method is benchmarked using (a) the
two-dimensional (2D) fractional diffusion equation, originally proposed by
Brunner, and (b) the incompressible, subdiffusive dynamics of a planar
viscoelastic channel flow of the Rouse chain melts (FADR equation with
fractional time-derivative of order $\alpha=\nicefrac{1}{2}$) and the Zimm
chain solution ($\alpha=\nicefrac{2}{3}$). Numerical simulations of the
viscoelastic channel flow effectively capture the non-homogeneous regions of
high viscosity at low fluid inertia (or the so-called `spatiotemporal
macrostructures'), experimentally observed in the flow-instability transition
of subdiffusive flows.
