#!/usr/bin/env python
# hypsometric.py
import numpy as num

def hypsometric(pressure,temperature,mixing_ratio,surface_altitude):
  """altitude = hypsometric(pressure,temperature,mixing_ratio,surface_altitude)
 
NAME:
      geopotential_altitude   - Paul van Delst's name of the routine
                      (DT renamed it to this).

PURPOSE:
 This function calculates geopotential altitudes using the hypsometric equation.
 It is like the "hypsometric" function I have implemented, but this one
 also accounts for water vapor in the atmosphere.


CATEGORY:
      Meteorology

CALLING SEQUENCE:
      altitude = hypsometric( pressure, $          ; Input
                              temperature, $       ; Input
                              mixing_ratio, $      ; Input
                              surface_altitude, $  ; Input

INPUTS:
      pressure:        Array of pressure in mb
      temperature:     Array of temperature in K
      mixing_ratio:    Array of water vapor mixing ratio in g/kg
      surface_height:  Height of surface in km. 

KEYWORD PARAMETERS:
      None.

OUTPUTS:
      altitude:        Geopotential altitudes in km in the same order as the
                       input data.

CALLS:
      None.


COMMON BLOCKS:
      None

SIDE EFFECTS:
      None

RESTRICTIONS:
   - Input pressure, temperature and mixing_ratio MUST be arrays with at
     LEAST two elements where the element with the highest pressure corresponds
     to the surface altitude passed.
   - Input surface height must be > or = to 0.0km and a SCALAR.
   - No allowance is made, yet, for the change in the acceleration due to 
     gravity with altitude.

PROCEDURE:
      Geopotential heights are calculated using the hypsometric equation:

                        -
                   Rd * Tv    [  p1  ]
        z2 - z1 = --------- ln[ ---- ]
                      g       [  p2  ]

      where Rd    = gas constant for dry air (286.9968933 J/K/kg),
            g     = acceleration due to gravity (9.80616m/s^2),
            Tv    = mean virtual temperature for an atmospheric layer,
            p1,p2 = layer boundary pressures, and
            z1,z2 = layer boundary heights.

      The virtual temperature, the temperature that dry air must have in
      order to have the same density as moist air at the same pressure, is
      calculated using:

                 [      1 - eps      ]
        Tv = T * [ 1 + --------- * w ]
                 [        eps        ]

  where T   = temperature,
        w   = water vapor mixing ratio, and
        eps = ratio of the molecular weights of water and dry air (0.621970585).

EXAMPLE:
      Given arrays of pressure, temperature, and mixing ratio:

        IDL> PRINT, p, t, mr
              1015.42      958.240
              297.180      291.060
              7.83735      5.71762

      the geopotential altitudes can be found by typing:

        IDL> result = geopotential_altitude( p, t, mr, 0.0, alt )
        IDL> PRINT, result, alt
               1      0.00000     0.500970

MODIFICATION HISTORY:
      Written by:     Paul van Delst, CIMSS/SSEC, 08-Dec-1997

Log: hypsometric.pro,m  
converted to matlab code - leslie moy 02-02-02

$Log: hypsometric2.pro,v $
Revision 1.1  2002/01/16 13:28:15  dturner
Initial revision

Revision 1.1  1999/03/03 16:24:55  paulv