First numerical results from the newly developed pseudo-spectral code P-FLARE (Penalised FLux-driven Algorithm for REduced models) are presented. This flux-driven turbulence/transport code uses a pseudo-spectral formulation with the penalisation method to impose radial boundary conditions. Its concise, flexible structure allows implementing various quasi-two-dimensional reduced fluid models in flux-driven formulation. Here, results from simulations of the modified Hasegawa–Wakatani system are discussed, where particle transport and zonal flow formation, together with profile relaxation, are studied. It is shown that coupled spreading/profile relaxation that one obtains for this system is consistent with a simple one-dimensional model of coupled spreading/transport equations. Then, the effect of a particle source is investigated, which results in the observation of sandpile-like critical behaviour. The model displays profile stiffness for certain parameters, with very different input fluxes resulting in very similar mean density gradients. This is due to different zonal flow levels around the critical value for the control parameter (i.e. the ratio of the adiabaticity parameter to the mean gradient) and the existence for this system of a hysteresis loop for the transition from two-dimensional turbulence to a zonal flow dominated state.