#!/usr/bin/python -tt #======================================================================= # General Documentation """Single-function module. See function docstring for description. """ #----------------------------------------------------------------------- # Additional Documentation # # RCS Revision Code: # $Id: static_stability.py,v 1.4 2004/03/16 02:23:05 jlin Exp $ # # Modification History: # - 15 Mar 2004: Original by Johnny Lin, Computation Institute, # University of Chicago. Passed reasonable tests. # # Notes: # - Written for Python 2.2.2. # - Module docstrings can be tested using the doctest module. To # test, execute "python static_stability.py". # - Non-"built-in" packages and modules required: atmqty, gemath. # See also import statements throughout for more information. # # Copyright (c) 2004 by Johnny Lin. For licensing, distribution # conditions, contact information, and additional documentation see # the URL http://www.johnny-lin.com/py_pkgs/atmqty/doc/. #======================================================================= #---------------- Module General Import and Declarations --------------- #- Set module version number etc. to package version number etc.: import package_version __version__ = package_version.version __author__ = package_version.author __date__ = package_version.date __credits__ = package_version.credits del package_version #-------------------------- General Function --------------------------- def static_stability(T_in, z_in, missing=1e+20): """Calculate static stability. The static stability is also known as the square of the buoyancy frequency, or the square of the Brunt-Vaisala frequency. If the static stability is N**2, the period of a buoyancy oscillation is 2*pi/N. Method Positional Arguments: * T_in: Air temperature [K]. 1-D Numeric floating point vector. * z_in: The elevation [m] of each element of T_in. 1-D Numeric floating point vector of same size as T_in. * missing: If T_in and/or z_in has missing values, this is the missing value value. Floating point scalar. Default is 1e+20. Output: * Static stability [1/s**2] at each element of T_in. Numeric floating point array of the same size and shape as input T_in. If there are any missing values in output, those values are set to the value in argument missing from the input. For instance, if T_in is only one element, a one-element vector is returned as the derivative with the value of argument missing. If there are missing values in the output due to math errors and missing is set to None, output will fill those missing values with the MA default value of 1e+20. References: * Iredell, Mark: "How does NCEP/NCAR reanalaysis compute the potential vorticity on isentropic surfaces?" NCEP/NCAR Reanaly- sis FAQ: http://dss.ucar.edu/pub/reanalysis/FAQ.html. * Wallace, J. M., and P. V. Hobbs (1977): Atmospheric Science: An Introductory Survey. San Diego, CA: Academic Press, ISBN 0-12-732950-1, p. 237. Example without missing values: >>> from static_stability import static_stability >>> import Numeric as N >>> T = N.array([283.1, 280.2, 278.4, 277.1]) >>> z = N.array([ 0.0, 100.1, 230.9, 350.5]) >>> N2 = static_stability(T, z) >>> ['%.5g' % N2[i] for i in range(4)] ['-0.00066521', '-0.00037055', '-9.2029e-05', '-3.9e-05'] Example with missing values: >>> T = N.array([283.1, 280.2, 278.4, 277.1]) >>> z = N.array([ 0.0, 100.1, 1e+20, 350.5]) >>> N2 = static_stability(T, z, missing=1e+20) >>> ['%.5g' % N2[i] for i in range(4)] ['-0.00066521', '1e+20', '-9.2029e-05', '1e+20'] Note how N2[2] is not missing, even though z[2] is, because the algorithm for static stability at N2[2] does not use z[2], al- though it does use T[2]. """ import gemath from gemath.deriv import deriv import MA import Numeric as N from atmconst import AtmConst from is_numeric_float import is_numeric_float #- Check input is of the correct type: if (is_numeric_float(T_in) != 1) or (is_numeric_float(z_in) != 1): raise TypeError, "static_stability: Arg(s) not Numeric floating" #- Import constants: const = AtmConst() #- Calculate vertical temperature gradient: dTdz_N = deriv(z_in, T_in, missing=missing, algorithm='default') #- Prep arrays for masked array calculation: if missing == None: dTdz = MA.masked_array(dTdz_N) T = MA.masked_array(T_in) else: dTdz = MA.masked_values(dTdz_N, missing, copy=0) T = MA.masked_values(T_in, missing, copy=0) #- Calculate and return static stability: return MA.filled( const.g / T * (dTdz + (const.g/const.c_p)) \ , missing ) #-------------------------- Main: Test Module ------------------------- #- Define additional examples for doctest to use: __test__ = {'Additional Examples': """ >>> from static_stability import static_stability >>> import Numeric as N >>> T = N.array([283.1, 280.2]) >>> z = N.array([ 0.0, 100.1]) >>> N2 = static_stability(T, z) >>> ['%.5g' % N2[i] for i in range(2)] ['-0.00066521', '-0.0006721'] >>> T = [283.1, 280.2] >>> z = N.array([ 0.0, 100.1]) >>> N2 = static_stability(T, z) Traceback (most recent call last): ... TypeError: static_stability: Arg(s) not Numeric floating >>> T = N.array([283.1, 280.2]) >>> z = N.array([ 0.0]) >>> N2 = static_stability(T, z) Traceback (most recent call last): ... ValueError: array dimensions must agree >>> T = N.array([283.1]) >>> z = N.array([ 0.0]) >>> N2 = static_stability(T, z, missing=1e+25) >>> ['%.5g' % N2[i] for i in range(1)] ['1e+25'] >>> T = N.array(283.1) >>> z = N.array( 0.0) >>> N2 = static_stability(T, z, missing=1e+25) Traceback (most recent call last): ... ValueError: deriv: Inputs not a vector """} #- Execute doctest if module is run from command line: if __name__ == "__main__": """Test the module. Tests the examples in all the module documentation strings, plus __test__. Note: To help ensure that module testing of this file works, the parent directory to the current directory is added to sys.path. """ import doctest, sys, os sys.path.append(os.pardir) doctest.testmod(sys.modules[__name__]) # ===== end file =====