The computation of NCEP Potential Vorticity

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 gridlengths.  (The vorticity is computed in the T36 spectral domain
for the 2.5 x 2.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.

Mark Iredell
EMC/NCEP