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”