What is the formula for monthly isentropic potential vorticity in the Reanalysis dataset?

I found the following explanation on the FAQ page on the NCAR Reanalysis site.

"How does NCEP/NCAR reanalaysis compute the potential vorticity on isentropic surfaces?" [12/06/96]

First the static stability N**2=g/T*(dT/dz+g/cp) is computed on model levels. Then the winds, temperature and static stability are interpolated to isentropic surfaces linearly in log(theta). (Outside the model domain, the fields are held constant for now and will later be compressed out of the final product.) Then the absolute vorticity zeta is spectrally computed on the isentropic surfaces, where the shortest wavelength in the spectral domain is about 4 grid lengths. (The vorticity is computed in the T36 spectral domain for the 2.5x2.5 degree grid.) The density rho=(T/theta)**(cp/R)*p0/(R*T) is also computed at this time directly on the isentropic surfaces. Finally, the NCEP potential vorticity is computed as zeta*N**2/(g*rho). Thus the units of NCEP PV are m**2/s/kg, which is different from the usual units of K*m**2/s/kg; one must multiply the NCEP PV by theta to compare them. The NCEP PV is still a form of the Ertel potential vorticity, since log(theta) is as well conserved as theta. Packing the PV in these units is simpler. The NCEP PV is currently rounded to the nearest 1e-10 m**2/s/kg for packing. The physical constants used are g=9.8, R=287.05, cp=1004.6 and omega=7.2921e-5. (by Dr.Mark Iredell)

-- Michael Bell