In order to do regrid the stellar parameters and find the best-fit model, I did the polynomial fitting of a combined P(Teff, log(g)). The original steps of P(Teff, log(g)) before the regridding is 50 K in Teff and 0.25 in log(g/(cm/s^2)), respectively.
 
Example:
Simulation spectrum: Teff=7245K, log(g)=3.62, Z=-0.2, [a/Fe]=0.25, r_mag = 15.0
 
Marginal probability distribution functions (PDF) along Teff and log(g). The original PDFs (blue) and polynomial fit (orange). The fitting was done only around the peak of the original PDF.

 
 
The PDF contours of the original (black) and the polynomial fit (blue). The polynomial fitting of P(Teff, log(g)) finds the best-fit model at its peak (a blue star), although this stellar typing itself has a systematic offset from the input Teff and log(g).