A new low-order scaling algorithm for second-order Moller-Plesset perturbation theory (MP2) is described. The method employs localized orbitals. However, in contrast to other local MP2 schemes reported in the literature, all atomic orbitals with significant contribution to the correlation energy for a given pair or orbitals are retained. Our approach introduces no significant approximation if all significant pairs are taken into account. However, the method is intended for the description of correlation between spatially close ("strong") pairs; for distant pairs, more economical alternative methods can be used. The algorithm, including a new implementation of the simultaneous two-index integral transformation scheme, is described in detail. Test calculations demonstrate asymptotic linear scaling for large systems. Strong-pair MP2 calculations with more that 1800 contracted basis functions performed on a single PC are reported. As a consequence of its lower scaling, the calculation of the MP2 energy is less expensive than the calculations of the Hartree-Fock energy for large systems. This has been demonstrated for strong pair correlation but is expected to hold also for full MP2 calculations.