High Reynolds number boundary layers (BLs) were simulated with the help of an improved large eddy simulation model. The LES code incorporated advanced numerics such as latter modification of pressure correction method and fully conservative central difference schemes. Advanced version of dynamic mixed subfilter model (DMM) was employed as the turbulence closure. Application of the DMM was favorable to simulations within a surface layer. Simulated BLs were thoroughly examined against accurate laboratory and more scattered atmospheric measurements. It was concluded that the LES model is able to reproduce the mean quantities, second and third order turbulent statistics in high Reynolds number flows. It was also recognized that a "true" Reynolds number of simulated flow is determined both the model resolution and employed subfilter model. Skewness and flatness revealed large deviation from the Gaussian statistics in LES results. This fact suggested a presence of large scale organized structures in the simulated flow. There is no possibility to derive such structures from atmospheric observations, therefore, analysis is necessary based on LES data. Multivariate statistical analysis of the most energetic eddies was performed for simulated laboratory and atmospheric BLs. Both types of BLs consisted of essentially the same organized structures known as hairpin vortices. Earlier theoretical stability analysis had predicted existence of a larger organized structures (rolls) in the atmospheric BLs. Rolls were not recognized in this analysis. Nevertheless, many distinct properties of rolls were captured by the atmospheric BLs structures.