gemath Examples: Calculating Curl (Additional)

Description and Import Statements

This web page provides additional examples of using curl_2d to calculate curl on winds from some synthetic data.

The basic examples page was kept limited to avoid there being too much clutter. That page, plus this one, when compared to function output, operate as a "unit test" of the function.

A link to the full source code and the data file used in these examples is at the bottom of this page. For more information on the curl_2d command, type help(curl_2d).

All examples below use the following import statements:

>>> import cdms, MV
>>> import gemath
>>> from gemath.curl_2d import curl_2d

Additional Examples Using CCM3 Output

Data

The data is the same CCM3 output as used in the basic examples page; see that page for plots of the data. The data is from a netCDF data file of zonal (u) and meridional (v) winds at a pressure level of 501.275 hPa for a moment of time in a CCM3 run (the moment in time is at 2970.0 days since 0000-09-01 00:00:00). Winds are in m/s. We read in the winds using netCDF utilities from cdms and convert the arrays to Numeric. Vectors x and y are longitude and latitude, respectively, in degrees:

>>> f = cdms.open('uv_501hPa_ccm3.nc')
>>> u = f('U')
>>> v = f('V')
>>> x = MV.filled( u.getAxis(-1)[:] )
>>> y = MV.filled( u.getAxis(-2)[:] )
>>> u = MV.filled( u[:,:] )
>>> v = MV.filled( v[:,:] )

Default Spherical Algorithm With Edge Missing Values

Same as the example with missing values on the basic examples page, but we also make some edge points "missing." This is to check how the order 1 (spherical) algorithm functions with missing edge points, which is chosen since spectral methods cannot work in those cases. We set a block of u in the "middle" of the domain to the default missing value 1e+20.

>>> u[30:41,20:37] = 1e+20
>>> u[0:9,50:67]   = 1e+20
>>> curl = curl_2d(x, y, u, v, algorithm='default_spherical')

Default Algorithm

For this domain, setting algorithm to 'default' will select order 1 Cartesian. Here x and y, although the values do not change, are assumed to be in meters, in which case curl is still in 1/s. Missing values are as above.

>>> u[30:41,20:37] = 1e+20
>>> u[0:9,50:67]   = 1e+20
>>> curl = curl_2d(x, y, u, v, algorithm='default')

Order 1 Cartesian

Here x and y, although the values do not change, are assumed to be in meters, in which case curl is still in 1/s. Missing values are as above. This plot should be identical to the default case above.

>>> u[30:41,20:37] = 1e+20
>>> u[0:9,50:67]   = 1e+20
>>> curl = curl_2d(x, y, u, v, algorithm='order1_cartesian')

Additional Examples Using Synthetic Data

Data

The velocities are generally non-zero, except poleward of 60N/S, in which they are set to zero. Assume velocities are in m/s. Vectors x and y are longitude and latitude, respectively, in degrees, unless a Cartesian algorithm is used in which x and y are in m:

>>> nx = 181
>>> ny = 91
>>> x = N.arange(nx) * 2.0 - 180.0
>>> y = N.arange(ny) * 2.0 - 90.0

>>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )
>>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )
>>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)

>>> N.putmask( u, N.where(N.absolute(u) < 1e-10, 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(v) < 1e-10, 1, 0), 0. )
>>> N.putmask( u, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )

The zonal (u) wind looks like:

and the meridional (v) wind looks like:

Order 1 Spherical

The algorithm used is order 1 differencing in spherical coordinates. Curl is in 1/s, if the winds are in m/s and the domain are in degrees longitude and latitude.

>>> curl = curl_2d(x, y, u, v, algorithm='order1_spherical')

SPHEREPACK

NCAR SPHEREPACK, which solves the curl in spectral space, is used as the algorithm. Since the grid is regular in latitude and longitude and include the poles, SPHEREPACK can be used.

>>> curl = curl_2d(x, y, u, v, algorithm='spherepack')

Order 1 Spherical (On Mars)

Curl is in 1/s, if the winds are in m/s and the domain are in degrees longitude and latitude. We assume we're on the planet Mars.

>>> curl = curl_2d( x, y, u, v, algorithm='order1_spherical' \
...               , R_sphere=3390e+3 )

SPHEREPACK (On Mars)

We assume we're on the planet Mars.

>>> curl = curl_2d( x, y, u, v, algorithm='spherepack' \
...               , R_sphere=3390e+3 )

Order 1 Cartesian

Here x and y, although the values do not change, are assumed to be in m, in which case curl is still in 1/s. Missing values are as above. We assume we're on the Earth.

>>> curl = curl_2d(x, y, u, v, algorithm='order1_cartesian')

Order 1 Spherical (On Mars), With Missing Values

Here x and y, although the values do not change, are assumed to be in meters, in which case curl is still in 1/s. Missing values are as above.

>>> u[30:46,90:114] = 1e+20
>>> curl = curl_2d( x, y, u, v, algorithm='order1_spherical' \
...               , R_sphere=3390e+3 )

Full Code and Data File For Examples

• test_curl_2d_add.py
• uv_501hPa_ccm3.nc