The advent of parallel computing has altered the landscape of turbulence simulations. Increased computer power using O(100-1000) processors permits large-eddy simulations (LESs) of atmospheric and oceanic PBLs in the presence of increasingly complex surface layers, for example turbulent flow over hills, water waves, and heterogeneous land surfaces. At the same time, advances in our ability to obtain high quality spatial turbulence measurements in the atmospheric surface layer, using arrays of sonic anemometers, has provided new insights into subfilter scale models that appear in LES codes. In this talk I will briefly describe a recently developed highly parallel algorithm for turbulence simulations. We use this code to simulate the familiar daytime convective boundary layer with grid meshes spanning 64^3 to 1024^3 using as many as 4096 CPUs. The higher resolution simulations are intriguing; coherent dust devils develop in the branches of the convective spokes and the vertical velocity skewness shows less bias and a pronounced shift towards data above mid-PBL. High resolution observational results from the "Horizontal Array Turbulence Study" field campaigns are also described. These databases are used to test subfilter scale parameterizations over a range of stratification and filter widths. We find that anisotropic production of subfilter scale flux, which is omitted in eddy-viscosity closures, is important for both momentum and scalars. For example the diagonal components of the subfilter-scale momentum flux are maintained almost entirely by anisotropic production. These results have implications for performing high Reynolds number rough-wall LES in the presence of surface layer shear.