;---------------------------------------------------------------------------- ; \$Id: q1q2_y73.pro,v 1.5 2002/12/27 02:42:07 johnny Exp \$ ;+ ; NAME: ; Q1Q2_Y73 ; ; PURPOSE: ; Calculate apparent heat source and moisture sink, after Yanai et al. ; (1973). ; ; CATEGORY: ; Atmospheric Sciences. ; ; CALLING SEQUENCE: ; Result = Q1Q2_Y73(in_dims, in_tdvars) ; ; INPUTS: ; in_dims: A structure containing variables describing the space ; and time dimensions of the (thermo)dynamic variables ; in in_tdvars. Unchanged by function. The structure ; tags are: ; ; delt: Timestep, in seconds [s]. Integer or single precision ; floating scalar. ; ; lon: Longitude [deg] for each x column of (thermo)dynamic ; variables. Integer or single precision floating array ; dimensioned (NX). In deg east [0,360]. Ordered so x- ; direction increases (in deg east) with increasing index ; value, i.e. so index 1 is east of index 0, etc. ; ; lat: Latitude [deg] for each y row of (thermo)dynamic vari- ; ables. Integer or single precision floating array ; dimensioned (NY). In deg north [-90,90]. Ordered so ; y-direction increases (in deg north) with decreasing ; index value, i.e. index 0 is northern-most latitude ; band and index NY-1 is the southern-most latitude band. ; ; plev: Pressure levels [hPa] for each x-y slice of (thermo)- ; dynamic variables. Integer or single precision floating ; array dimensioned (NP). Ordered descending in magni- ; tude, i.e. with the bottom pressure level at index 0 ; and the top pressure level at index NP-1. ; ; in_tdvars: A structure containing the (thermo)dynamic variables. ; Unchanged by function. The structure tags are: ; ; u: Zonal velocity [m/s] at each (x,y,p,t) location. ; Floating array dimensioned (NX,NY,NP,NT). ; ; v: Meridional velocity [m/s] at each (x,y,p,t) location. ; Floating array dimensioned (NX,NY,NP,NT). ; ; T: Air temperature [K] at each (x,y,p,t) location. Float- ; ing array dimensioned (NX,NY,NP,NT). ; ; q: Mixing ratio of water vapor [kg/kg] at each (x,y,p,t) ; location. Floating array dimensioned (NX,NY,NP,NT). ; ; INPUT KEYWORD PARAMETERS: ; KDAY: If keyword is set true, the output Q1 and Q2 (through ; Result and keyword Q2) are given as equivalent heating ; rates from Q1 and Q2 in K/day instead of the default ; J/(kg s). ; ; XPER: If keyword is set, the input (thermo)dynamic variables ; are assumed to wrap-around in the x-direction, i.e. ; that column index NX (if there was one in the input ; arrays) is the same as column index 0. If keyword is ; not set, the input are not assumed periodic in the x- ; direction. ; ; OUTPUT KEYWORD PARAMETERS: ; Q2: Apparent moisture sink, Q2 [J/(kg s), or K/day if ; keyword Kday true]. Floating array dimensioned ; (NX,NY,NP,NT). This keyword is optional, but the ; function will always calculate it so not setting it ; will not save you on execution time. Created/over- ; written by function. ; ; OUTPUTS: ; Result: Apparent heat source, Q1 [J/(kg s), or K/day if ; keyword Kday true]. Floating array dimensioned ; (NX,NY,NP,NT). Created/overwritten by function. ; ; FILE DEVICE I/O: ; None. ; ; COMMON BLOCKS: ; None. ; ; EXAMPLE: ; None. ; ; REFERENCES: ; - Yanai, M. and R. H. Johnson (1993): "Impacts of cumulus ; convection on thermodynamic fields," The Representation of ; Cumulus Convection in Numerical Models (Boston, MA: American ; Meteorological Society), Meteorological Monographs, Vol. 24, ; No. 46, ISBN 1-878220-13-6, pp. 39-62. Abbreviated YJ93. ; - Yanai, M., S. Esbensen, and J.-H. Chu (1973): "Determination ; of bulk properties of tropical cloud clusters from large-scale ; heat and moisture budgets," J. Atmos. Sci., Vol. 30, pp. 611-627. ; Abbreviated Y73. ; ; MODIFICATION HISTORY: ; - 26 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. ; - Variable NT: This is the number of timesteps. Timesteps are ; assumed to be all the same and equal to in_dims.delt. ; - Lon-lat grid: The divergence calculation routine assumes that ; each x-y row has the same latitude value and each x-y column has ; the same longitude value and that positive directions are standard ; zonal and meridional. ; - Checks that the lon, lat, and plev vectors increase correctly are ; found in DIV_VEL_SPH and OMEGA_VEL so they aren't programmed here. ; - Calculations in this routine of dry static energy uses the pressure ; levels (in_plev) as input and computes altitude using the ICAO ; standard atmosphere (PRESS2ALT). ; - Although the original definition of Q1 and Q2 was in Y73, the ; expression used in this function for Q1 and Q2 is from YJ93, eqns. ; 4.4 and 4.5. ; - For purposes of calculating omega, the reference pressure is set ; to 1000 hPa. This is set by variable p0 in the "Error Check and ; Parameter Setting" section. ; - Specific heat at constant pressure (c_p) and gravity acceleration ; are assumed constant in calculating dry static energy. ; - Calculations in function are generally done in floating (single or ; double, depending on which is the highest precision present in ; in_tdvars) point (set by ftype). ; - All keyword parameters are optional unless otherwise stated. ; - No procedures called with _Extra keyword are invoked. ; - Non-built-in functions/procedures called: DIV_VEL_SPH, OMEGA_VEL, ; PRESS2ALT (by Dominik Brunner). ;- ; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions, ; and contact information, see http://www.johnny-lin.com/lib.html. ;---------------------------------------------------------------------------- FUNCTION Q1Q2_Y73, in_dims, in_tdvars \$ , KDAY = in_kday \$ , XPER = in_xper \$ , Q2 = out_q2 \$ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ; ; Key section input: None. ; Key section output: delt, ftype, lat, lon, p0, plev, q, T, u, v, xper ; ; Selected key variables: ; delt Timestep, scalar [sec]. ; ftype Type code to do calculations in, scalar. ; lat Latitude, dimensioned (NY) [deg]. ; lon Longitude, dimensioned (NX) [deg]. ; p0 Reference pressure, dimensioned (NX,NY) [hPa]. ; plev Pressure, dimensioned (NP) [hPa]. ; q Water vapor mixing ratio, dimensioned (NX,NY,NP,NT) [kg/kg]. ; T Air temperature, dimensioned (NX,NY,NP,NT) [K]. ; u Zonal velocity, dimensioned (NX,NY,NP,NT) [m/s]. ; v Meridional velocity, dimensioned (NX,NY,NP,NT) [m/s]. ; xper Flag if array is periodic in x-direction, scalar. ; ; Array dimensions NX,NY,NP,NT (scalars) are also calculated, set, and ; output in the section. COMPILE_OPT IDL2 ON_ERROR, 0 epsilon = (MACHAR()).eps ;- Set machine precision tolerance c_p = 1004.0 ;- Atmospheric specific heat at constant ; pressure [J/(kg K)] d2r = !PI / 180.0 ;- Factor to mult. to go from deg to rad g = 9.81 ;- Gravity acceleration [m/s^2] L = 2501e3 ;- Latent heat of vaporization [J/kg] R_e = 6.371e6 ;- Earth's mean radius [m] lon = FLOAT(in_dims.lon) ;- Protect input parameters lat = FLOAT(in_dims.lat) plev = FLOAT(in_dims.plev) delt = FLOAT(in_dims.delt) u = in_tdvars.u v = in_tdvars.v T = in_tdvars.T q = in_tdvars.q ;- Set number of points in lon, lat directions, pressure levels, and ; number of timesteps: NX = N_ELEMENTS(lon) NY = N_ELEMENTS(lat) NP = N_ELEMENTS(plev) NT = (SIZE(u, /Dimensions))[3] ;- Make sure input are all floating or double and set the type to do ; calculations in (ftype): inp_types = [ SIZE(u, /Type), SIZE(v, /Type) \$ , SIZE(T, /Type), SIZE(q, /Type) ] if ((MAX(inp_types) gt 5) or (MIN(inp_types) lt 4)) then \$ MESSAGE, 'error--input not single or double floating type' ftype = MAX(inp_types) ;- Check all input (thermo)dynamic variables have the same ; dimensions and that they are the same as lon, lat, plev, NT: tmp = [ (SIZE(u))[0:4] ne (SIZE(v))[0:4] \$ , (SIZE(u))[0:4] ne (SIZE(T))[0:4] \$ , (SIZE(u))[0:4] ne (SIZE(q))[0:4] \$ , (SIZE(q))[0:4] ne [4,NX,NY,NP,NT] ] if (TOTAL(tmp) gt epsilon) then MESSAGE, 'error--dim(s) not consistent' ;- Supplementary checks for realism of values: if (delt le 0.0) then MESSAGE, 'error--neg. delt' if (delt gt 60.0*60.0*24.0*14.0) then \$ MESSAGE, 'error--routine delt probably inaccurately large for d/dt' if (MIN(T) le 50.0) then MESSAGE, 'error--T unlikely to be air temp. in K' ;- Set flag for whether x is periodic or not (xper=1 means is periodic): if (KEYWORD_SET(in_xper) eq 1) then \$ xper = 1 \$ else \$ xper = 0 ;- Set press. at reference level [hPa] as scalar or (NX,NY) array: p0 = 1000.0 ; ----------------------- Calculate Dry Static Energy ----------------------- ; ; Key section input: plev, T ; Key section output: s (dimensioned NX,NY,NP,NT) ; ; Algorithm: s = c_p*T + g*z, where z is the altitude (see YJ93, p. 39). ; First we calculate plev as (NX,NY,NP) array (plev_xyp) and calculate ; altitude z_xyp (NX,NY,NP). ; ; Units: s [J/kg], plev [hPa], T [K]. plev_xyp = REBIN(REFORM(plev, 1,1,NP), NX,NY,NP) ;- plev as (NX,NY,NP) array z_xyp = PRESS2ALT(TEMPORARY(plev_xyp)) ;- Alt. as (NX,NY,NP) array s = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) ;- Init. divergence array ;- Calculate s at each time point (done as a loop in time to help ; conserve memory): for it=0,NT-1 do begin s[*,*,*,it] = ( c_p * T[*,*,*,it] ) + ( g * z_xyp ) endfor tmp = TEMPORARY(z_xyp) tmp = 0 ; ----------- Calculate Divergence and Vertical Pressure Velocity ----------- ; ; Key section input: ftype, lat, lon, p0, plev, u, v, xper ; Key section output: omega (dimensioned NX,NY,NP,NT) ; ; Units: div [1/s], lat [deg], lon [deg], omega [hPa/s], p0 [hPa], ; plev [hPa], u [m/s], v [m/s]. div = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) ;- Init. divergence array omega = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) ;- Init. omega array for it=0,NT-1 do begin for ik=0,NP-1 do begin div[*,*,ik,it] = DIV_VEL_SPH( u[*,*,ik,it], v[*,*,ik,it] \$ , lon, lat, Xper=xper ) endfor endfor for it=0,NT-1 do begin omega[*,*,*,it] = OMEGA_VEL( plev, div[*,*,*,it], P0=p0 ) endfor tmp = TEMPORARY(div) ;- Deallocate to save memory tmp = 0 ; ------- Calculate Spatial, Pressure, and Time Derivatives of s and q ------ ; ; Key section input: ftype, lon, lat, plev, q, s, u, v ; Key section input: dsdt, dsdp, dsdx, dsdy, dqdt, dqdp, dqdx, dqdy ; ; Selected key variables: ; dqdp Pressure derivative of q, dimensioned (NX,NY,NP,NT) ; dqdt Time derivative of q, dimensioned (NX,NY,NP,NT) ; dqdx Derivative of q wrt x, dimensioned (NX,NY,NP,NT) ; dqdy Derivative of q wrt y, dimensioned (NX,NY,NP,NT) ; dsdp Pressure derivative of s, dimensioned (NX,NY,NP,NT) ; dsdt Time derivative of s, dimensioned (NX,NY,NP,NT) ; dsdx Derivative of s wrt x, dimensioned (NX,NY,NP,NT) ; dsdy Derivative of s wrt y, dimensioned (NX,NY,NP,NT) ; ; Units: omega [hPa/s], dsdp [J/(kg hPa)], dsdt [J/(kg s)], plev [hPa], ; s [J/kg], u [m/s], v [m/s]. ;- Initialize output arrays: dqdp = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dqdt = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dqdx = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dqdy = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dsdp = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dsdt = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dsdx = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) dsdy = MAKE_ARRAY(NX,NY,NP,NT, Type=ftype) ;- Compute d/dp: for it=0,NT-1 do begin for iy=0,NY-1 do begin for ix=0,NX-1 do begin dqdp[ix,iy,*,it] = DERIV( plev, q[ix,iy,*,it] ) dsdp[ix,iy,*,it] = DERIV( plev, s[ix,iy,*,it] ) endfor endfor endfor ;- Compute d/dt: for ik=0,NP-1 do begin ;+ d/dt with dt units of delt for iy=0,NY-1 do begin for ix=0,NX-1 do begin dqdt[ix,iy,ik,*] = DERIV( q[ix,iy,ik,*] ) dsdt[ix,iy,ik,*] = DERIV( s[ix,iy,ik,*] ) endfor endfor endfor dqdt = TEMPORARY(dqdt) / delt ;+ dq/dt with dt units of sec dsdt = TEMPORARY(dsdt) / delt ;+ ds/dt with dt units of sec ;- Compute d/dx and d/dy: lon_iy_m is longitude in lon at a given ; latitude value in meters E from 0E, lat_m is all latitudes in lat ; (for all longitudes cases) in meters N from equator: lat_m = R_e * lat * d2r for it=0,NT-1 do begin for ik=0,NP-1 do begin for iy=0,NY-1 do begin lon_iy_m = R_e * COS(lat[iy] * d2r) * lon * d2r dqdx[*,iy,ik,it] = DERIV( lon_iy_m, q[*,iy,ik,it] ) dsdx[*,iy,ik,it] = DERIV( lon_iy_m, s[*,iy,ik,it] ) endfor for ix=0,NX-1 do begin dqdy[ix,*,ik,it] = DERIV( lat_m, q[ix,*,ik,it] ) dsdy[ix,*,ik,it] = DERIV( lat_m, s[ix,*,ik,it] ) endfor endfor endfor ; --------------------------- Calculate Q1 and Q2 --------------------------- ; ; Key section input: dsdt, dsdp, dsdx, dsdy, dqdt, dqdp, dqdx, dqdy, ; omega, u, v ; Key section output: q1, q2 (both dimensioned NX,NY,NP,NT) ; ; If keyword Kday is true, q1 and q2 are given in K/day. All dq* and ; ds* variables deleted on section output. q1 = TEMPORARY(dsdt) \$ + ( (u * TEMPORARY(dsdx)) + (v * TEMPORARY(dsdy)) ) \$ + (omega * TEMPORARY(dsdp)) q2 = -L * ( TEMPORARY(dqdt) \$ + ( (u * TEMPORARY(dqdx)) + (v * TEMPORARY(dqdy)) ) \$ + (omega * TEMPORARY(dqdp)) ) if (KEYWORD_SET(in_kday) eq 1) then begin q1 = TEMPORARY(q1) / c_p * 86400.0 q2 = TEMPORARY(q2) / c_p * 86400.0 endif ; ------------------------------ Prepare Output ----------------------------- out_q2 = TEMPORARY(q2) ;- Output Q2 (optional) Result = TEMPORARY(q1) ;- Output Q1 RETURN, Result END ;=== end of function === ; ========== end file ==========