function precip = precip_extremes_scaling(omega, temp, plev, ps)

 error(nargchk(4, 4, nargin))

 % Calculate scaling estimate for the intensity of precipitation extremes following equation 2 of O'Gorman and Schneider, PNAS, 2009 based on an approximate expression for the vertically integrated condensation rate

 % Inputs
 % omega: vertical profile of vertical velocity in Pa/s (conditioned on extremes occuring) 
 % temp: vertical profile of temperature in K (conditioned on extremes occuring)
 % plev: pressure in Pa
 % ps: surface pressure in Pa (conditioned on extremes occuring)
 
 % Output 
 % precip: precipitation estimate in kg/m^2/s (note doesn't include precipitation efficiency factor)

 % assumes that highest pressure is first element of arrays
 if plev(1)<plev(2)
  error('unexpected ordering of pressure levels')
 end
 
 % criterion for identifying tropopause
 crit_lapse_rate = 0.002; % (k/m) for tropopause
 plev_mask = 0.05e5; % (Pa) exclude levels above this as a fail-safe

 % find thermodynamic derivative of saturation specific humidity (qsat) with respect to temperature and pressure
 [dqsat_dp, dqsat_dT] = sat_deriv(plev, temp); % uses era40 formulation
 [es, qsat, rsat, latent_heat] = saturation_thermodynamics(temp, plev, 'era');

 % find moist adiabatic lapse rate
 lapse_rate = moist_adiabatic_lapse_rate(temp, plev, 'era');

 % virtual temperature
 temp_virtual = temp.*(1.0+qsat*(PARS('gas_constant_v')/PARS('gas_constant')-1.0));
 % density
 rho = plev/PARS('gas_constant')./temp_virtual; 

 dT_dp = lapse_rate/PARS('gravity')./rho;

 % find derivative of saturation specific humidity with respect to pressure along a moist adiabat at the given temperature and pressure for each level
 dqsat_dp_total = dqsat_dp+dqsat_dT.*dT_dp;

 % mask above tropopause using simple lapse rate criterion
 dT_dp_env = gradient(temp, plev);
 lapse_rate_env = dT_dp_env.*rho*PARS('gravity');

 itrop = max(find(lapse_rate_env>crit_lapse_rate));
 if itrop<length(plev)
  dqsat_dp_total(itrop+1:end)=0;
 end

 % mask above certain level as fail safe
 dqsat_dp_total(find(plev<plev_mask))=0;

 dqsat_dp_total_omega = dqsat_dp_total.*omega;

 % replaces nans with zeros as subsurface values should not contribute
 % to the column integral
 inan = find(isnan(dqsat_dp_total_omega));
 dqsat_dp_total_omega(inan)=0;

 % also use surface pressure to zero subsurface values
 kbot = max(find(plev>ps));
 if ~isempty(kbot)
  dqsat_dp_total_omega(1:kbot)=0;
 end

 % integrate in the vertical
 precip = -integrate(-dqsat_dp_total_omega',plev)/PARS('gravity');

 return

%%%%%%%%

function [es, qs, rs, L] = saturation_thermodynamics(T, p, type)
% SATURATION_THERMODYNAMICS  Saturation quantities for water.
%
% [es, qs, rs, L] = SATURATION_THERMODYNAMICS(TEMP, PRESS, type) returns
% the effective saturation vapor pressure, specific humidity, mixing
% ratio, and latent heat of vaporization/sublimation. Depending on
% type, is uses different formulations:
%    'era': Formulation used by ECMWF based on the modified
%           Tetens formula.
%    'simple': constant latent heat of condensation
  
  Rd          = PARS('gas_constant');
  Rv          = PARS('gas_constant_v');
  gc_ratio    = Rd/Rv;   % ratio of gas constants for dry air and water vapor

  if nargin < 3
    type = 'simple';
  end

  persistent been_verbose
  if isempty(been_verbose)
    disp(['SATURATION_THERMODYNAMICS: Using ', upper(type), ...
          ' formulation for computation of saturation quantities.'])
    been_verbose = 1;
  end
    
  switch lower(type)
    
   case 'era'

    % ECMWF formulation described by Simmons et al. (1999: QJRMS, 125,
    % 353--386), which uses saturation over ice for temperatures
    % less than 250 K and a quadratic interpolation between
    % saturation over ice and over liquid water for temperatures
    % between 250 K and 273 K
    
    % coefficients in Tetens saturation vapor pressure es = es0 * exp(a3 * (T-T0)/(T-a4))      
    es0       = 611.21;  % saturation vapor pressure at T0 (Pa)
    T0        = 273.16;  % (K)
    Ti        = T0 - 23; % (K)

    a3l       = 17.502;  % liquid water (Buck 1981)
    a4l       = 32.19;   % (K)
  
    a3i       = 22.587;  % ice (Alduchov and Eskridge 1996)
    a4i       = -0.7;    % (K)
  
    % saturation vapor pressure over liquid and ice
    esl       = es0 .* exp(a3l .* (T - T0)./(T - a4l));
    esi       = es0 .* exp(a3i .* (T - T0)./(T - a4i));

  
    % latent heat of sublimation 
    Ls0       = 2.834e6;  % latent heat of sublimation  (J / kg) [+- 0.01 error for 173 K < T < 273 K]
    Ls        = Ls0 * ones(size(T));

    % latent heat of vaporization
    Lv0       = 2.501e6;  % latent heat of vaporization at triple point (J / kg)
    cpl       = 4190;     % heat capacity of liquid water (J / kg / K)
    cpv       = 1870;     % heat capacity of water vapor (J / kg / K)			  
    Lv        = Lv0 - (cpl - cpv) .* (T - T0);
  
    % compute saturation vapor pressure and latent heat over liquid/ice mixture
    iice      = find(T <= Ti);
    iliquid   = find(T >= T0);
    imixed    = find(T > Ti & T < T0);
  
    es        = ones(size(T))*NaN;
    L         = ones(size(T))*NaN;
    if ~isempty(iice)
      es(iice)  = esi(iice);
      L(iice)   = Ls(iice);
    end
  
    if ~isempty(iliquid)
      es(iliquid) = esl(iliquid);
      L(iliquid)  = Lv(iliquid);
    end
    
    if ~isempty(imixed)
      a         = ( (T(imixed) - Ti)./(T0 - Ti) ).^2;
      es(imixed)= (1-a) .* esi(imixed) + a .* esl(imixed);
      L(imixed) = (1-a) .* Ls(imixed) + a .* Lv(imixed);
    end

   case 'simple'
    
    % assuming constant latent heat of vaporization and only one
    % vapor-liquid phase transition
    T0        = 273.16;
    es0       = 610.78;
    L         = PARS('latent_heat_v');   
    
    % saturation vapor pressure
    es        = es0 * exp(L/Rv*(1.0/T0 - 1.0./T));
    
   otherwise
    error(['Unknown type, ', type, 'of computing saturation quantities.'])
  end
    
  % saturation mixing ratio
  rs          = gc_ratio * es ./ (p - es);
  
  % saturation specific humidity 
  qs          = rs ./ (1 + rs); 
  
  
 return


%%%%%%%%%%%

function [dqsat_dp, dqsat_dT, dln_esat_dT] = sat_deriv(p, temp);
% Calculates derivatives of the saturation specific humidity wrt 
% temperature and pressure

  error(nargchk(1, 2, nargin))     
  
 % use finite difference approximation
 dp = 0.1;
 dT = 0.01;

 [es_p_plus, qs_p_plus] = saturation_thermodynamics(temp, p+dp, 'era');
 [es_p_minus, qs_p_minus] = saturation_thermodynamics(temp, p-dp, 'era');
 [es_T_plus, qs_T_plus] = saturation_thermodynamics(temp+dT, p, 'era');
 [es_T_minus, qs_T_minus] = saturation_thermodynamics(temp-dT, p, 'era');

 dqsat_dp   = (qs_p_plus-qs_p_minus)/(2.0*dp);

 dqsat_dT   = (qs_T_plus-qs_T_minus)/(2.0*dT);

 dln_esat_dT = (log(es_T_plus)-log(es_T_minus))/2/dT;


%%%%%%%%%%%

function F = integrate(f, x)
%INTEGRATE  Computes one-dimensional integral.
%    INTEGRATE(f, x) computes an approximation to the integral 
%    \int f(x) dx over the range of the input vector x. The input x
%    must be a vector. The input f can be a column vector or a
%    matrix. The number of rows of f must be the same as the length
%    of x. If f is a matrix, a vertical integral is computed for
%    each column.

  error(nargchk(2, 2, nargin))     
  
  x       = x(:);
  if size(f, 1) ~= length(x)
    error('Length of x must be equal to number of rows of f.')
  end
  
  dx1      = gradient(x);
  dx1(1)   = 0.5 * dx1(1);
  dx1(end) = 0.5 * dx1(end);
  dx2      = repmat(dx1, 1, size(f, 2));
  
  F        = sum(f .* dx2);

 return


%%%%%%%%%%%

function lapse_rate = moist_adiabatic_lapse_rate(T, p, type)
% MOIST_ADIABATIC_LAPSE_RATE Returns saturated moist-adiabatic lapse rate.
%  
% Units are K / m.

  g               = PARS('gravity');
  cpd             = PARS('cp');
  cpv             = PARS('cp_v');
  Rd              = PARS('gas_constant');
  Rv              = PARS('gas_constant_v');
  gc_ratio        = Rd / Rv;
  
  persistent not_first_call
  
  [es, qs, rs, L] = saturation_thermodynamics(T, p, type);

    
  switch lower(type)
   case 'simple'
    if isempty(not_first_call)
      disp(['MOIST_ADIABATIC_LAPSE_RATE: Using simplified moist adiabatic ' ...
            'lapse rate'])
      not_first_call = 1;
    end

    % cf. Holton (p. 503), obtained
    % by using rs<<1 in the more accurate expression below 
    lapse_rate      = g/cpd ...
        .* (1 + L.*rs ./ Rd ./ T)...
        ./ (1 + L.^2 .* rs ./(cpd * Rv .* T.^2));
    
   otherwise
    % cf. Emanuel (p. 131)
    lapse_rate      = g/cpd * (1 + rs) ./ (1 + cpv/cpd.*rs) ...
        .* (1 + L.*rs ./ Rd ./ T)...
        ./ (1 + L.^2 .* rs .* (1 + rs/gc_ratio)./(Rv .* T.^2 .* (cpd + rs*cpv)));
  end

 return
     
%%%%%%%%%%%

function [par] = PARS(parstr, options);
% PARS  Parameters 

  % physical parameters
  gravity 	      = 9.80665;      % in m / s^2

  gas_constant	      = 287.04;       % J / kg / K 

  kappa               = 2.0/7.0;
  cp                  = gas_constant/kappa;  % specific heat at constant pressure
  cp_v                = 1870;         % specific heat water vapor [J/kg/K]

  gas_constant_v      = 461.50;       % gas constant water vapor [J/kg/K]
  latent_heat_v       = 2.5e6;        % latent heat of evaporation [J/kg]

  eval(['par = ', parstr, ';']);
  
 return
  
