Photoacoustic (PA) tomography (PAT) is a rapidly developing imaging modality that can provide high contrast and spatial-resolution images of light-absorption distribution in tissue. However, reconstruction of the absorption distribution is affected by nonuniform light fluence. This paper introduces a reconstruction method for reducing amplification of noise and artifacts in low-fluence regions. In this method, fluence compensation is integrated into model-based reconstruction, and the absorption distribution is iteratively updated. At each iteration, we calculate the residual between detected PA signals and the signals computed by a forward model using the initial pressure, which is the product of estimated voxel value and light fluence. By minimizing the residual, the reconstructed values converge to the true absorption distribution. In addition, we developed a matrix compression method for reducing memory requirements and accelerating reconstruction speed. The results of simulation and phantom experiments indicate that the proposed method provides a better contrast-to-noise ratio (CNR) in low-fluence regions. We expect that the capability of increasing imaging depth will broaden the clinical applications of PAT.