Inversion of magnetic data is an important step towards interpretation of the practical data. Smooth inversion is a common technique for the inversion of data. Physical bound constraint can improve the solution to the magnetic inverse problem. However, how to introduce the bound constraint into the inversion procedure is important. Imposing bound constraint makes the magnetic data inversion a non-linear inverse problem. In this work, a new algorithm is developed for the 3D inversion of magnetic data, which uses an efficient penalization function for imposing the bound constraint and Gauss Newton method to achieve the solution. An adaptive regularization method is used in order to choose the regularization parameter in this inversion approach. The inversion results of synthetic data show that the new method can produce models that adequately match the real location and shape of the synthetic bodies. The test carried out on the field data from Mt. Milligan copper-gold porphyry deposit shows that the new inversion approach can produce the magnetic susceptibility models consistent with the true structures.