A neutral boundary layer nonhydrostatic numerical model is used to determine the characteristics of shelterbelt effects on mean wind direction and to study the processing causing wind rotation when air passes through a shelterbelt. The model uses a turbulence scheme that includes prognostic equations for turbulence kinetic energy and a master length scale proposed by Mellor and Yamada. The simulated results are in quantitative agreement with Nord's field measurements. The spatial variation of wind rotation and its dependence on incident angle and shelterbelt porosity is analysed. Dynamic processes of the wind rotation and its interactions with drag force and pressure perturbation are also discussed. It is concluded that shear of wind direction should be considered, along with shear of speed, in determining turbulent fluxes in the vicinity of a shelterbelt.