Carleson measures are a useful tool in harmonic analysis. They are connected with the nontangential maximal function and the so called square function. However, they have been mostly studied in the Euclidean setting. We prove characterization of Carleson measures in Heisenberg group. Moreover, for a subelliptic harmonic function we give a bound on a norm of its square function in terms of norm of its boundary data and prove that square of norm of gradient of a subelliptic harmonic function gives rise to a Carleson measure. We mostly work with nontangentially accessible domains. The talk is based on joint work with Tomasz Adamowicz.