In this paper, a modified lattice Bhatnagar-Gross-Krook (LBGK) model for axisymmetric thermal flow is proposed. In the model, two tunable parameters related to the relaxation times are available by introducing additional gradient terms into the equilibrium functions. By doing so, the values of the dimensionless relaxation times in the LBGK model can be kept in a proper range, and thus the stability of the LBGK model can be improved. Local schemes for the calculation of the additional gradient terms are presented to keep the intrinsic merit of LB. Also, some correction terms are introduced into the evolution function for the temperature field to eliminate the deviation terms exist in the previous LBGK model. Chapman-Enskog analysis demonstrates that the macro equations in the cylindrical coordinates system can be exactly recovered. Numerical tests, including thermal flow in a pipe, thermal Womersley flow and natural convection in a vertical annulus, have been carried out, and the results predicted by the present LBGK model agree well with the existing numerical data. Also, numerical results demonstrate that the present model is superior to the existing LBGK model for axisymmetric thermal flow in term of numerical stability. (C) 2016 Elsevier Ltd. All rights reserved.