Abstract: With the increasing popularity of Numerical Weather Prediction (NWP) models and the wide applications of Potential Vorticity (PV) at various weather and climate scales, further improvement on the accuracy of the Ertel PV calculation with model output is imperative. The main source of computational errors can be attributed to interpolation (interpolation method) in terrain-following coordinates, which are commonly used in mesoscale models. Variables used for the calculation of Ertel PV need to be interpolated from the z or p coordinates to terrain-following coordinate first. In order to improve the accuracy of Ertel PV calculation, the algorithm for the Ertel PV defined in z coordinate is converted to a terrain-following coordinate form using the transformation relation between terrain-following coordinate and z(p) coordinate. In this way, the Ertel PV defined in z coordinate can be calculated directly on model grids (the direct method) without interpolation. Meanwhile, the physical insight of Ertel PV in the z coordinate is kept to make it convenient for analyses and applications. To further reveal the effect of the direct method in reducing computational errors, an extratropical cyclone outbreak occurred in the complex terrain area in North China was simulated using the Weather Research Forecasting (WRF) model, and detailed analysis of the computational errors in Ertel PV using the interpolation method was conducted. The result shows that by applying the direct method, the RMS errors in the middle-lower and middle-upper levels can be reduced by 0.5 PVU and 0.3 PVU, respectively, which will benefit the analysis of fine structures of convective storms and the study of climate statistics.