function dtheta_dp_eff = eff_stat_stab(p, temp, lambda) % calculates the effective static stability derived in O'Gorman, JAS, 2011, pages 75-90 according to equation 8 % % inputs are 1d vertical profiles of pressure and temperature and the asymmetry parameter lambda (default value 0.6): % p pressure (Pa) % temp temperature (K) % lambda asymmetry parameter defined by equation 5 of O'Gorman, JAS, 2011 error(nargchk(2, 3, nargin)) if nargin<3 lambda = 0.6; % default value of lambda (appropriate for midlatitudes) end if ~isequal(size(temp), size(p)) | min(size(p))~=1 error('temp and p must be one-dimensional and of same size') end % constants Rd = 287.04; % gas constant for dry air [J/kg/K] Rv = 461.5; % gas constant for water vapor [J/kg/K] cpd = 1005.7; % specific heat dry air [J/kg/K] cpv = 1870; % specific heat water vapor [J/kg/K] g = 9.80665; % gravitational acceleration [m/s^2] p0 = 1e5; % reference pressure [Pa] kappa = Rd/cpd; gc_ratio = Rd/Rv; % saturation vapor pressure [Pa] Tc = temp-273.15; es = 611.20*exp(17.67*Tc./(Tc+243.5)); % Bolton 1980 equation 10 % latent heat of condensation [J/kg] L = (2.501-0.00237*Tc)*1e6; % Bolton 1980 equation 2 % saturation mixing ratio rs = gc_ratio*es./(p-es); % saturation specific humidity qs = rs./(1+rs); % potential temperature exponent = kappa.*(1+rs./gc_ratio)./(1+rs.*cpv./cpd); theta = temp.*(p0./p).^exponent; % derivative of potential temperature with respect to pressure dtheta_dp = gradient(theta, p); % density temp_virtual = temp.*(1.0+rs/gc_ratio)./(1.0+rs); rho = p/Rd./temp_virtual; % moist adiabatic lapse rate malr = g/cpd*(1+rs)./(1+cpv/cpd.*rs) ... .*(1+L.*rs./Rd./temp)... ./(1+L.^2.*rs.*(1+rs/gc_ratio)./(Rv.*temp.^2.*(cpd+rs*cpv))); % derivative of potential temperature wrt pressure along a moist adiabat % (neglects small contribution from vertical variations of exponent) dtemp_dp_ma = malr/g./rho; dtheta_dp_ma = dtemp_dp_ma.*theta./temp-exponent.*theta./p; % effective static stability following equation 8 of O'Gorman, JAS, 2011 dtheta_dp_eff = dtheta_dp-lambda.*dtheta_dp_ma;