Influence of finite difference schemes and subgrid-stress models on the large eddy simulation calculation of turbulent flow around a bluff body of square cylinder at a laboratory Reynolds number, has been examined. It is found that the type and the order of accuracy of finite-difference schemes and the subgrid-stress model for satisfactory results are dependent on each other, and the grid resolution and the Reynolds number. Using computational grids manageable by workstation-level computers, with which the near-wall region of the separating boundary layer cannot be resolved, central-difference schemes of realistic orders of accuracy, either fully conservative or non-conservative, suffer stability problems. The upwind-biased schemes of third order and the Smagorinsky eddy-viscosity subgrid model can give reasonable results resolving much of the energy-containing turbulent eddies in the boundary layers and in the wake and representing the subgrid stresses in most parts of the flow. Noticeable improvements can be obtained by either using higher order difference schemes, increasing the grid resolution and/or by implementing a dynamic subgrid stress model, but each at a cost of increased computational time. For further improvements, the very small-scale eddies near the upstream corners and in the laminar sublayers need to be resolved but would require a substantially larger number of grid points that are out of the range of easily accessible computers. © 2002 John Wiley and Sons, Ltd.