;---------------------------------------------------------------------------- ; $Id: omega_vel.pro,v 1.6 2002/11/23 21:57:35 johnny Exp $ ;+ ; NAME: ; OMEGA_VEL ; ; PURPOSE: ; Calculate vertical pressure velocity given divergence at pressure ; levels. ; ; CATEGORY: ; Atmospheric Sciences. ; ; CALLING SEQUENCE: ; Result = OMEGA_VEL(in_plev, in_div) ; ; INPUTS: ; in_plev: Pressure levels [hPa] at which input divergence is ; given, dimensioned (NP). Ordered descending in ; magnitude, i.e. with the bottom-most pressure level ; at index 0 and the top-most pressure level at index ; NP-1. Unchanged by function. ; ; in_div: Velocity divergence [1/s], dimensioned (NX,NY,NP), ; where (NX,NY) define a 2-D horizontal grid, given at ; the NP pressure levels given in in_plev. If the ; horizontal grid is only a single point, in_div should ; be dimensioned (1,1,NP). Unchanged by function. ; ; INPUT KEYWORD PARAMETERS: ; P0: Pressure at reference level [hPa]. Scalar or array ; dimensioned (NX,NY). If undefined, defaults to 1000. ; Vertical pressure velocity (and thus horizontal ; divergence) at reference level is assumed to be 0. ; Unchanged by function. ; ; OUTPUT KEYWORD PARAMETERS: ; None. ; ; OUTPUTS: ; Result: Vertical pressure velocity [hPa/s], dimensioned ; (NX,NY,NP). Floating or double, depending on input ; type. Created/overwritten by function. ; ; FILE DEVICE I/O: ; None. ; ; COMMON BLOCKS: ; None. ; ; EXAMPLE: ; Given the following values of divergence at these levels: ; lev = [950, 900, 800, 700] ; div = ([-4, -4, -3, -2]) * 1e-6 ; div = REFORM(div,1,1,4) ; result = OMEGA_VEL(lev,div) ; PRINT, lev, REFORM(result) * 3600. ; ; IDL prints out: ; 950 900 800 700 ; -0.360000 -1.08000 -2.34000 -3.24000 ; ; These values of divergence are similar to values given in Fig. 6 ; of Yanai et al. (1973, J. Atmos. Sci., Vol. 30, pp. 611-627). ; ; REFERENCES: ; - Holton, J. R. (1992): An Introduction to Dynamic Meteorology ; (San Diego, CA: Academic Press), ISBN 0-12-354355-X, Section ; 3.5. ; ; MODIFICATION HISTORY: ; - 23 Nov 2002: Orig. ver. Johnny Lin, CIRES/University of Colorado. ; Email: air_jlin@yahoo.com. Passed passably adequate tests. ; ; NOTES: ; - Written for IDL 5.3. ; - Vertical pressure velocity is calculated by integrating the con- ; tinuity equation (e.g. eqn. 3.38 in Holton 1992). ; - Computations are done with the highest floating point precision ; of the in_plev and in_div inputs. If neither are float or double, ; function defaults to single precision float. ; - For synoptic-scale motions, it is approximately true that pressure ; velocity omega = - rho * g * w, where w is the vertical velocity ; (eqn. 3.37 in Holton 1992). ; - All keyword parameters are optional unless otherwise stated. ; - No procedures called with _Extra keyword are invoked. ; - Only built-in functions and procedures called. ;- ; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions, ; and contact information, see http://www.johnny-lin.com/lib.html. ;---------------------------------------------------------------------------- FUNCTION OMEGA_VEL, in_plev, in_div $ , P0 = in_p0 $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ; ; Key section input: None. ; Key section output: div, NX, NY, NP, plev, p0 ; ; Selected key variables: ; div Divergence, dimensioned (NX,NY,NP). ; ftype Type code that all computations will be done in. ; plev Pressure, dimensioned (NP). ; p0 Reference pressure, dimensioned (NX,NY). ; ; On section exit, div, plev, and p0 are all set to the type specifed by ; ftype. COMPILE_OPT IDL2 ON_ERROR, 0 epsilon = (MACHAR()).eps ;- Set machine precision tolerance div = in_div ;- Protect input parameters plev = in_plev NP = N_ELEMENTS(plev) ;- Set number of pressure levels dims = SIZE(div) ;- Set NX,NY, check div dimensions, NX = dims[1] NY = dims[2] if (dims[0] ne 3) then MESSAGE, 'error--div bad num of dims' if (dims[3] ne NP) then MESSAGE, 'error--div/plev mismatch' ;- Set calculation type code (ftype) and convert input to that type: ftype = MAX([SIZE(div, /Type), SIZE(plev, /Type)]) ;+ Initial setting if (ftype ne 5) then ftype = 4 ;+ Final setting case ftype of 4: div = FLOAT(div) 5: div = DOUBLE(div) else: MESSAGE, 'error--bad ftype' endcase case ftype of 4: plev = FLOAT(plev) 5: plev = DOUBLE(plev) else: MESSAGE, 'error--bad ftype' endcase ;- Set array for press. at reference level (set to keyword P0 if ; defined, to 1000 if not defined) and convert type if needed: if (N_ELEMENTS(in_p0) eq 0) then begin p0 = REPLICATE(1000.0, NX,NY) endif else begin case N_ELEMENTS(in_p0) of 1: p0 = REPLICATE( in_p0, NX,NY ) NX*NY: p0 = REFORM( in_p0, NX,NY ) else: MESSAGE, 'error--bad P0' endcase endelse case ftype of 4: p0 = FLOAT(p0) 5: p0 = DOUBLE(p0) else: MESSAGE, 'error--bad ftype' endcase ;- Check input plev vector increases correctly: tmp = SORT(-plev) ne LINDGEN(NP) if (TOTAL(tmp) gt epsilon) then MESSAGE, 'error--dim(s) in bad order' ;- Check minimum of p0 is larger than the maximum of plev: if (MAX(plev) ge MIN(p0)) then $ MESSAGE, 'error--p0 value must always be larger than plev value' ; ---------------- Vertically Integrate Horizontal Divergence --------------- ; ; Key section input: div, p0, plev ; Key section output: omega ; ; Algorithm: Integration is by the trapezoidal rule, the piecewise sum ; of the area under straight-line pieces of the divergence curve, from p0 ; to p. ; ; Selected variable key: ; plev_xyp Array plev as an array dimensioned (NX,NY,NP). ; plev_m1_xyp Array of pressure levels one element below/before the ; current one, as an array dimensioned (NX,NY,NP). ; dp Pressure weighting between current level and the one ; below/before it, dimensioned (NX,NY,NP), in hPa. ;- Create *_xyp pressure arrays and calc. press. weightings for each ; level: plev_xyp = REBIN(REFORM(plev, 1,1,NP), NX,NY,NP) plev_m1_xyp = SHIFT(plev_xyp, 0,0,1) ;+ Shift to get below/before plev plev_m1_xyp[*,*,0] = p0[*,*] ;+ Set 1st below/before level to p0 dp = TEMPORARY(plev_xyp) - TEMPORARY(plev_m1_xyp) ;- Integrate horizontal divergence from p0 to p level: div_m1 = SHIFT(div, 0,0,1) ;+ Shift to get below/before div div_m1[*,*,0] = 0.0 ;+ Set divergence at level p0 is 0 omega = -TOTAL( (0.5 * (div + div_m1) * dp), 3, /Cumulative ) ; ------------------------------ Prepare Output ----------------------------- Result = TEMPORARY(omega) ;- Output vertical pressure velocity RETURN, Result END ;=== end of function === ; ========== end file ==========