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
.
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 = cdms.open('nnr_2000-01-15.nc') >>> 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:
AtmQty
Using Reanalysis Data, Not Including Buoyancy FrequencyWe 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.
AtmQty
Using Reanalysis Data, Including Buoyancy FrequencyWe 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.