The wide use of satellite-based instruments provides measurements in climatology on a global scale, which often have nonstationary covariance structure. In this paper we address the issue of modeling axially symmetric spatial random fields on sphere with a kernel convolution approach. The observed random field is generated by convolving a latent uncorrelated random field with a class of Matern type kernel functions. By allowing the parameters in the kernel functions to vary with locations, we are able to generate a flexible class of covariance functions and capture the nonstationary properties. Since the corresponding covariance functions generally do not have a closed form, numerical evaluations are necessary and a pre-computation table is used to speed up the computation. For regular grid data on sphere, the circulant block property of the covariance matrix enables us to use Fast Fourier Transform (FFT) to get its determinant and inverse matrix efficiently. We apply this approach to the Total Ozone Mapping Spectrometer (TOMS) ozone data and compare it with other existing models.