Ocean mesoscale eddies strongly impact large-scale circulations via eddy-induced buoyancy, tracer, and momentum fluxes. Eddy buoyancy and momentum fluxes can be combined into the potential vorticity (PV) flux, for which parameterizations have been proposed in previous studies based on physical closure theories. However, these parameterizations have seen limited implementation in ocean models because of their idealized assumptions and numerical issues. Since modern ocean climate models can partially resolve mesoscale eddies, the parameterization for subgrid eddy mixing should account for the resolved eddy field and adapt to varying model resolutions. Here, we propose a novel approach to represent subgrid PV mixing through the implicit numerical diffusion of the PV advection scheme in isopycnal ocean models. The PV flux is reconstructed via a seventh-order weighted essentially non-oscillatory (WENO) scheme, which contains a physically constrained numerical diffusion that dissipates potential enstrophy at the grid scale while preserving the inverse energy cascade. The scheme is implemented in GFDL MOM6 and evaluated against existing energy- and/or enstrophy-conserving schemes in both idealized and realistic eddy-permitting configurations (e.g., 1/4^∘). In idealized wind-driven circulation simulations, the scheme yields more energetic eddies and large-scale circulations, with reduced viscous energy dissipation. In the 1/4^∘ global ocean simulation, it produces more realistic pathways of western boundary currents, such as Kuroshio and Gulf Stream, thereby reducing their associated sea surface temperature biases. These results suggest that the WENO PV advection scheme developed here provides physically consistent subgrid PV mixing and can complement mesoscale eddy parameterizations in climate models.