In my last post I discussed a paper that described the theory of effective static stability (O’Gorman 2011). Dr. O’Gorman posted a MATLAB script for calculating effective static stability, but since I work with NCL, I wanted convert his code into an NCL function. The original MATLAB code only worked on one profile at a time, but the function below “*should”* work for a multi-dimensional array, but I won’t guarantee that it will work for any array size.

;================================================================================ ;================================================================================ ; Calculate effective static stability - based on Paul O'Gorman's matlab code ; derived in O'Gorman, JAS, 2011, pages 75-90 (eq 8) ; inputs are 1d vertical profiles of pressure and temperature ; and the asymmetry parameter lambda (default value 0.6): ; Ta absolute temperature (K) - must be on pressure levels with a 'lev' coordinate ; lambda asymmetry parameter (defined in eq 5) - should not have a time dimension undef("calc_ess") function calc_ess (Ta,lev,lambda) local kappa,gc_ratio,P,Tc,es,L,rs,qs,exponent,theta, thetav,dtheta_dp,rho,malr,dtemp_dp_ma,dtheta_dp_ma begin cpd = 1004. ; J / (kg K) Rd = 287.04 ; J / (kg K) kappa = Rd/cpd gc_ratio = Rd/Rv Po = 100000. ; Check input if .not.iscoord(Ta,"lev") printExit("ERROR: calc_ess: Ta must have a 'lev' dimension ") end if if .not.isatt(lev,"units") then printExit("ERROR: calc_ess: the lev input must have units of 'Pa' ") end if if dimsizes(dimsizes(lev)).eq.1 then P = conform(Ta,lev,ind(getvardims(Ta).eq."lev")) else P = lev end if copy_VarAtts(lev,P) ; saturation vapor pressure [Pa] from Bolton (1980) eq 10 Tc = Ta - 273.15 es = 611.20 * exp(17.67*Tc/(Tc+243.5)) ; latent heat of condensation [J/kg] from Bolton (1980) eq 2 L = (2.501 - 0.00237*Tc)*1e6 ; 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 = Ta*(Po/P)^exponent copy_VarCoords(Ta,theta) ; derivative of potential temperature with respect to pressure dtheta_dp = calc_ddp(theta) ; calculate density from virtual potential temperature thetav = Ta*(1.0+rs/gc_ratio)/(1.0+rs) rho = P/Rd/thetav ; moist adiabatic lapse rate malr = g/cpd * (1.+rs) / (1.+cpv/cpd*rs) * (1.+L*rs/Rd/Ta) \ / ( 1. + L^2.*rs*(1.+rs/gc_ratio)/(Rv*Ta^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/Ta - exponent*theta/P ; effective static stability following equation 8 of O'Gorman, JAS, 2011 sdim = ispan(1,dimsizes(dimsizes(Ta))-1,1) tlambda = conform(Ta,lambda,sdim) dtheta_dp_eff = dtheta_dp - tlambda*dtheta_dp_ma delete(tlambda) copy_VarCoords(Ta,dtheta_dp_eff) return dtheta_dp_eff end ;======================================================================= ;=======================================================================

To use this, you simply need to copy it into a file and load it into your NCL script before the “begin” statement using a statement like this,

load “effective_static_stability.ncl”

I prefer to put files with NCL functions in my NCL root directory, so I can link to them like this,

load “$NCARG_ROOT/effective_static_stability.ncl”