;----------------------------------------------------------------------------
; $Id$
;
;                           function reduce_to_slp
;
; Purpose:
;   Given a value of air pressure, air temp., and water vapor mixing raio,
;   at a given elevation, reduces the pressure value to sea level.
;
; File I/O:
;   None.
;
; Entry and Exit States:
;   N/A.
;
; Input Parameter:
;   in_ps          Air pressure (in hPa).  Can be scalar, or any size
;                  array.  Unchanged by procedure.
;   in_Ts          Surf. temperature (in K).  Can be scalar, or any size
;                  array.  Unchanged by procedure.
;   in_qs          Water vapor mixing ratio (in kg/kg).  Same size as
;                  in_ps.  Unchanged by procedure.
;   in_zs          Elevation (in m) above mean sea level.  Same size as 
;                  in_ps.  Unchanged by procedure.
;
; Output Parameter:
;   out_slp        Pressure reduced to sea level (hPa).  Created.  Same 
;                  size as in_ps.  Double precision floating.
;
; Keywords (optional unless otherwise noted):
;   None.
;
; Revision History:
; - 13 Apr 2001:  Orig. ver. by Johnny Lin, CIRES/Univ. of Colo.  Email:
;   air_jlin@yahoo.com.  Passed moderately reasonable tests.
;
; Notes:
; - Written for IDL 5.2.1.
; - Copyright (c) 2001 Johnny Lin.  For licensing and contact information
;   see http://www.johnny-lin.com/lib.html.  
; - Calculations are all done in DOUBLE.
; - No procedures called with _EXTRA keyword invoked.
; - No user-written procedures called.
; - No common blocks are used in this program.
; - References:  See code comments below.
;----------------------------------------------------------------------------

FUNCTION REDUCE_TO_SLP, in_ps, in_Ts, in_qs, in_zs  $
                      , _EXTRA = extra



; -------------------- Error Check and Parameter Setting -------------------- 

ON_ERROR, 0

if (N_PARAMS() ne 4) then MESSAGE, 'error-bad param list'

ps = DOUBLE(in_ps)      ;- Protect input
Ts = DOUBLE(in_Ts)
qs = DOUBLE(in_qs)
zs = DOUBLE(in_zs)

nps = N_ELEMENTS(ps)

if (nps ne N_ELEMENTS(Ts)) then MESSAGE, 'error--arrays different sizes'
if (nps ne N_ELEMENTS(qs)) then MESSAGE, 'error--arrays different sizes'
if (nps ne N_ELEMENTS(zs)) then MESSAGE, 'error--arrays different sizes'



; ---------------------- Calculate Virtual Temperature ----------------------
;
; - See Wallace, J.M. and P.V. Hobbs (1977):  Atmospheric Science:  An
;   Introductory Survey, San Diego, Academic Press, pp. 50, 72.
; - The mean virtual temperature is computed assuming the water vapor 
;   mixing ratio is const. from surf. to sea level, that the temp. from 
;   the surf. to sea level has a lapse rate of 6.50 K/km, and that the
;   mean is the arithmetic average of virtual temp. at sea level and the
;   surf.
; - Rd_Rv = 0.622    Ratio of the gas constants for dry air:water vapor

econst  = 0.607717041d0                   ;- Equals (1 - Rd_Rv) / Rd_Rv
T0      = (0.0065d0 * zs) + Ts            ;- Sea level temp. [K]
q0      = qs                              ;- Sea level water vapor mix ratio

T_vs    = Ts * ( 1.d0 + (econst * qs) )   ;- Surf. virtual temp. [K]
T_v0    = T0 * ( 1.d0 + (econst * q0) )   ;- Sea level virtual temp. [K]
T_vmean = 0.5d0 * (T_v0 + T_vs)           ;- Layer mean virtual temp. [K]



; ---------------------- Calculate Sea Level Pressure -----------------------
;
; - Use hypsometric equation.  See Wallace, J.M. and P.V. Hobbs (1977):
;   Atmospheric Science:  An Introductory Survey, San Diego, Academic 
;   Press, pp. 59.
; - Set surf. geopotential height to geometric height (i.e. assume g does
;   not change appreciably with height.

p_g = ps * 100.d0    ;- Pressure at ground surface [Pa]
g0  = 9.81d0         ;- Gravitational constant at earth's surf. [m/s^2]
R_d = 287.0d0        ;- Gas constant for dry air [J/(K kg)]
Z_g = zs             ;- Surface geopotential height [m]

out_slp = p_g * EXP( (g0*Z_g) / (R_d*T_vmean) )   ;- Sea level press. [Pa]
out_slp = TEMPORARY(out_slp) / 100.d0             ;- Convert to hPa



; ---------------------------- Clean-Up and Output --------------------------

RETURN, out_slp



END         ; ===== end of function =====

; ========== end file ==========