We developed a new model for the X-ray spectral fitting \xspec package which takes into account the effects of both thermal and dynamical (i.e. bulk) Comptonization. The model consists of two components: one is the direct blackbody-like emission due to seed photons which are not subjected to effective Compton scattering, while the other one is a convolution of the Green's function of the energy operator with a blackbody-like seed photon spectrum. When combined thermal and bulk effects are considered, the analytic form of the Green's function may be obtained as a solution of the diffusion Comptonization equation. Using data from the BeppoSAX, INTEGRAL and RXTE satellites, we test our model on the spectra of a sample of six persistently low magnetic field bright neutron star Low Mass X-ray Binaries, covering three different spectral states. Particular attention is given to the transient powerlaw-like hard X-ray (> 30 keV) tails that we interpret in the framework of the bulk motion Comptonization process. We show that the values of the best-fit delta-parameter, which represents the importance of bulk with respect to thermal Comptonization, can be physically meaningful and can at least qualitatively describe the physical conditions of the environment in the innermost part of the system. Moreover, we show that in fitting the thermal Comptonization spectra to the X-ray spectra of these systems, the best-fit parameters of our model are in excellent agreement with those of COMPTT, a broadly used and well established XSPEC model.