This paper investigates the control volume method with multipoint flux approximation (MPFA) applied to the discetization of –div(K(x)grad u) = f(x), the equation describing fluid flow through anisotropic porous media. The permeability tensor K(x) is allowed to have discontinuities. Transmissibility coefficients are obtained from local numerical flow experiments (transmissibility upscaling) for each cell face. Monotonicity of the solution matrix is discussed and a version of the method that provides an M-matrix is described. Convergence and numerical flux consistency of the scheme for uniform permeability tensor are studied. This discretization scheme is applied to reservoir simulation on 3D structured grids with distorted geometry, highly anisotropic media and discontinuities in the permeability tensor. Simulation results are presented and compared with results from other two-point and multipoint flux approximations.