atmqty Examples: Isentropic Potential Vorticity


This web page provides an example of calculating isentropic potential vorticity (IPV) using the ipv method in the AtmQty class. The data used in this example comes from the NCEP/NCAR reanalysis. These plots operate as a "unit test" of the function.

A link to the full source code needed and the data file used in these examples is at the bottom of this page. For more information on the ipv method see the pydoc description for ipv.

Import Statements and Data

All examples below use the following import statements and data:

>>> import atmqty
>>> import cdms
>>> import MV
>>> import Numeric as N

The data is from a netCDF data file of fields from the NCEP/NCAR reanalysis for 15 Jan 2000 (0Z) on isentropic levels. We read in the data using netCDF utilities from cdms, convert the arrays to Numeric arrays, and convert units as needed to AtmQty standard units. Vectors lon and lat are longitude and latitude, respectively, in degrees, and lev is isentropic levels, in K. We transpose the arrays so that the first index cycles through latitude, the second through longitude, and the third through levels. We also reverse the latitude order so that latitude increases with increasing index (in the dataset latitudes start from 90°N and progress south):

>>> f ='')

>>> T_obj = f['TMP']
>>> lon = T_obj.getLongitude()[:]
>>> lat = T_obj.getLatitude()[::-1]
>>> lev = T_obj.getLevel()[:].astype(N.Float)

>>> T = MV.filled( T_obj )
>>> u = MV.filled( f['UGRD'] )
>>> v = MV.filled( f['VGRD'] )
>>> rh = MV.filled( f['RH']/100.0 )
>>> n2 = MV.filled( f['BVF2'] )

>>> T = N.transpose(T, axes=(1,2,0))[::-1,:,:]
>>> u = N.transpose(u, axes=(1,2,0))[::-1,:,:]
>>> v = N.transpose(v, axes=(1,2,0))[::-1,:,:]
>>> rh = N.transpose(rh, axes=(1,2,0))[::-1,:,:]
>>> n2 = N.transpose(n2, axes=(1,2,0))[::-1,:,:]

Temperature is in K, winds are in m/s, relative humidity is a fraction, and buoyancy frequency (n2) is in 1/s2. The levels are at 300, 315, 330, and 350 K. Longitude and latitude values are separated every 2.5 degrees. The missing value in these arrays are specified by miss_val = T_obj.missing_value[0].

Finally, we read in IPV (and convert units to SI) from the NCEP/NCAR reanalysis by:

>>> pv = MV.transpose( f['PV___'] \
...                  , axes=(1,2,0) )[::-1,:,:] * 1e-6

This IPV at isentropic level 315 K (units K m2 s-1 kg-1) looks like:


Relative Vorticity Calculated By AtmQty Using Reanalysis Data, Not Including Buoyancy Frequency

We create an AtmQty object (x1) for the above ECMWF domain and fill it with the above atmospheric quantities, except buoyancy frequency; we want the AtmQty object to calculate buoyancy frequency. Using those quantities, we compute isentropic potential vorticity:

>>> x1 = atmqty.AtmQty( lat, lon, lev, lev_type="isentropic" \
...                   , missing=miss_val )
>>> x1.add_qty({'T':T, 'u':u, 'v':v, 'RH':rh})
>>> x1.ipv()

Note that using the method ipv in the fashion above calculates altitude automatically, using default settings. If the default settings (e.g. surface pressure) do not apply, you need to manually call method alt.

The 315 K level is plotted below. The contour intervals are the same as in the reanalysis IPV above (units K m2 s-1 kg-1):

Because the level below 315 K (i.e. 300 K) has a number of missing values in temperature in the Tropics, when you calculate the vertical temperature gradient at 315 K (and buoyancy frequency), you get a missing value in those areas. Thus, IPV is undefined too. Outside those regions, however, you get a decent match with IPV from reanalysis.

Relative Vorticity Calculated By AtmQty Using Reanalysis Data, Including Buoyancy Frequency

We create and fill the following AtmQty object:

>>> x2 = atmqty.AtmQty( lat, lon, lev, lev_type="isentropic" \
...                   , missing=miss_val )
>>> x2.add_qty({'T':T, 'u':u, 'v':v, 'RH':rh, 'N2':n2})
>>> x2.ipv()

The IPV at 315 K is:

By using the buoyancy frequency provided by the reanalysis (array n2), many of the IPV missing value areas in the Tropics are filled in with data.

The contour intervals are the same as in the reanalysis IPV above. Again, we get good agreement with the reanalysis IPV.

Full Code and Data File For Examples

Return to atmqty Home Page.
Return to Johnny Lin's Python Library.

Valid HTML 4.01! Valid CSS! Updated: May 26, 2004 by Johnny Lin <email address>. License.