Shelterbelts are used world-wide for such purposes as reduction of soil errosion, control of snow drift, and provision of an effective agrometeorological method of field microclimate management and yield enhancement. Whether performing a wind tunnel test, conducting a field observation, or implementing a numerical simulation to investigate shelterbelt effects, researchers are more interested in an optimum reduction in a thin air lasier near the ground on the leeside of the shelterbelt rather than total wind-speed reduction in the whole flowfield. The purpose of this study is to formulate a Navier-Stokes based scheme to simulate the turbulent aerodynamic characteristics of a shelterbelt. Qualitative results from field observation of a living-tree shelterbelt under real atmospheric flow conditions and a wind-tunnel flow visualization of scale-model fences were used to explore the fundamental phenomena of the shelterbelt flow to help in the numerical modeling. A modified higher-order numerical scheme using the Lagrange interpolation to represent the interface convection terms is developed and applied to better simulate the turbulent shelterbelt flowfield. It is shown that this new scheme not only can enhance accuracy during computation but also is capable of retaining the numerical stability and good convergence characteristics which are lost in most higher-order numerical schemes. The flow retardation and porosity of shelterbelts are modelled via momentum sources with the help of the aerodynamic parameters, normal pressure drag and skin friction drag. The results obtained from this newly developed numerical scheme show satisfactory agreement with both field experiments and other numerical simulations. In addition, this procedure offers a generalized technique for simulating more complicated shelterbelt configurations.