function res=fall_speed(T,p,pden,r)
% calculate aerosol fall velocity using Stokes Law 
% dynamic viscosity of air is only a function of temperature,  
% Rogers and Yau (R+Y) pg 103.
% The fall speed (v1) is from Pruppacher + Klett (P+K), pg 362
% Another approach from R+Y for fall speed (pg 97)
% v2    = pden2*(2.0*r*1.e-6)^2 * g / (18.0 * dvisc)  ; m/s
% input:
%       T - temperature (K)
%       p - pressure (mb)
%       pden - particle density (g/cm3)
%       r - particle radius (microns)
% output:
%       v - aerosol fall speed (m/s)

  pden2=pden*1000;                            % particle density kg/m3
  g=9.81;                                     % gravity (m/s2)
  Rd=287.04;                                  % gas constant (J/kg/K)
  aden=100.0*p/(Rd*T);                        % air density (kg/m3) 
  mfp=6.6e-2*(1013.25/p)*(T/293.15);          % P+K pg 323, mean free path (um)
  KN=mfp/r;                                   % Knudson number
  a=1.257+0.4*exp(-1.1/KN);                   % P+K pg. 261, eqn 12-16
  
  dvisc=1.72e-5*(393.0/(T+120.0))*(T/273.0)^(1.5);         % dynamic viscosity air (kg/m/s)
  res=2.*(1.+a*KN)*(pden2-aden)*g*(r*1.e-6)^2 /(9.*dvisc); % fall velocity, m/s


