This article presents a new and substantially improved finite volume procedure for simulation of incompressible flows on non-orthogonal grids. Cell-centered least-squares gradients are obtained in a robust and highly accurate way. A new discretization of the diffusive terms is employed, which is based on extension of the original cell-face gradient interpolation and is more suitable for complex grid distortions. A flexible flux-limited interpolation of dependent variables on distorted computational grids is introduced. An efficient preconditioner for Krylov method solution of linear systems is proposed, which substantially improves the solution of Poisson equation for pressure correction. The pressure-correction algorithm is adapted for efficient convergence on highly complex grids using a sequence of non-orthogonal corrector solutions and its effect on iteration convergence is analyzed. The non-orthogonalities treated by current procedure are more accustomed to numerical grids generated from a real complex terrain elevation data. The main focus is on the simulation of atmospheric micro-scale flows pertinent to wind energy application. (C) 2015 Elsevier Inc. All rights reserved.