In this work a nonlinear CW ultrasound field propagation model based on a second-order operator splitting approach is studied and a number of significant enhancements are introduced and implemented. In this model the ultrasound field is calculated and propagated plane by plane and the effects of diffraction, nonlinearity and absorption are applied independently over incremental steps. This work completes the preceding works (Christopher and Parker 1991, Tavakkoli et al. 1998, Zemp et al. 2003, Williams et al. 2006) by introducing an arbitrary source geometry and excitation definition, full diffraction solution, enhanced pressure, enhanced power deposition rate and temperature prediction capabilities. The result is a particularly useful tool in carrying out simulations of high intensity focused ultrasound (HIFU) that includes temperature rise predictions. Comparisons are made with other codes in both linear and nonlinear regimes. Different dynamics of lesion formation are obtained in linear versus nonlinear models, specially at the onset of lesion creation during HIFU exposure.