This paper presents a three-dimensional stochastic linear model of the atmospheric flow induced by the variability of heat flux over land surface. The primitive equations relating perturbation terms of wind field, geopotential and buoyancy are formulated as a system of stochastic partial differential equations and solved analytically. The solution is based on spectral representations of the homogeneous random fields. The flow intensity is found to be proportional to the standard deviation of the heat flux into the atmosphere. The intensity of the vertical motion becomes more sensitive to the differential heating with a larger length scale as altitude goes higher. Stability and synoptic wind inhibit the development of the flow. The proposed theory improves the understanding of the role that heterogeneous land surface plays in atmospheric circulations at the mesoscale.