In this paper we develop a theoretical framework for the phase-space Gaussian beamlet decomposition of the beam-like wavefields in plasmas. It is based on the discrete Gabor expansion of the wavefield distribution along a pre-selected equiphase surface. In particular, we consider the situations where the Gaussian beamlets can be chosen much narrower than the parent beam, though much wider than the wavelength of the radiation. Eventual summation of the independently computed beamlets will thus give the net wavefield that properly accounts for the effects of smaller-scale inhomogeneities. While using the standard beam-tracing (BT) method for the computation of each propagating beamlet, this approach is able to improve the accuracy of the stand-alone BT technique applied to the wide parent beam in cases of strong refraction or oblique absorption. We demonstrate the potency and relatively low computational cost of the introduced method in numerical examples, which reproduce different patterns of microwave propagation in the TJ-II stellarator plasma.