Natural land surfaces are usually heterogeneous over the resolvable scales considered in atmospheric numerical models. Therefore, model surface parameterizations that assume surface homogeneity may fail to represent the surface forcing accurately. In this paper, a parameterization of the subgrid-scale forcing of heterogeneous land surfaces for atmospheric numerical models is suggested. In each surface grid element of the numerical model similar homogeneous land patches located at different places within the element are regrouped into subgrid classes. Then, for each one of the subgrid classes, a sophisticated micrometeorological model of the soil-plant-atmosphere system is applied to assess the surface temperature, humidity, and fluxes to the atmosphere. The global fluxes of energy between the grid and the atmosphere are obtained by averaging according to the distribution of the subgrid classes. In addition to the surface forcing, detailed micrometeorological conditions of the patches are assessed... Abstract Natural land surfaces are usually heterogeneous over the resolvable scales considered in atmospheric numerical models. Therefore, model surface parameterizations that assume surface homogeneity may fail to represent the surface forcing accurately. In this paper, a parameterization of the subgrid-scale forcing of heterogeneous land surfaces for atmospheric numerical models is suggested. In each surface grid element of the numerical model similar homogeneous land patches located at different places within the element are regrouped into subgrid classes. Then, for each one of the subgrid classes, a sophisticated micrometeorological model of the soil-plant-atmosphere system is applied to assess the surface temperature, humidity, and fluxes to the atmosphere. The global fluxes of energy between the grid and the atmosphere are obtained by averaging according to the distribution of the subgrid classes. In addition to the surface forcing, detailed micrometeorological conditions of the patches are assessed...