A new Stokesian dynamics (SD) algorithm for Brownian suspensions is presented. The implementation is based on the recently developed accelerated Stokesian dynamics (ASD) simulation method [Sierou and Brady, J. Fluid Mech. 448, 115 (2001)] for non-Brownian particles. As in ASD, the many-body long-range hydrodynamic interactions are computed using fast Fourier transforms, and the resistance matrix is inverted iteratively, in order to keep the computational cost O(N log N). A fast method for computing the Brownian forces acting on the particles is applied by splitting them into near- and far-field contributions to avoid the O(N-3) computation of the square root of the full resistance matrix. For the near- field part, representing the forces as a sum of pairwise contributions reduces the cost to O(N); and for the far-field part, a Chebyshev polynomial approximation for the inverse of the square root of the mobility matrix results in an O(N-1.25 log N) computational cost. The overall scaling of the method is thus roughly of O(N-1.25 log N) and makes possible the simulation of large systems, which are necessary for studying long-time dynamical properties and/or polydispersity effects in colloidal dispersions. In this work the method is applied to study the rheology of concentrated colloidal suspensions, and results are compared with conventional SD. Also, a faster approximate method is presented and its accuracy discussed. (C) 2003 American Institute of Physics.