function fall_speed_pmc, P,T,rad,AR,shape

;------------------------------------------------------------------------------
; Compute the fall sped of a PMC ice particle according to Turco et al. [1982]
;
; Input:
;   P.........pressure (hPa)
;   T.........temperature (K)
;   rad.......particle radius (nm)
;   AR........axial ratio,  eccentricity or length / diameter, i.e. rotational over equatorial axis 
;   shape.....0 = sphere, 1 = cylinder, 2 = hexagon
;
; Output:
;   v_fall....fall speed (m/s)
;
; Source: Mark Hervig
;-------------------------------------------------------------------------------

   rad_m = rad * 1e-9  ; convert radius in nm to m
   
;- Constants

   Di   = 930.0           ; ice density, kg/m3
   g    = 9.57            ; gravity, 84 km altitude, m/s2
   Ma   = 4.845d-26       ; air molecule mass, kg
   Rd   = 287.04          ; gas constant (J/kg/K)
   Na   = 6.02252e23      ; Avagadros #, molec/mol
   Md   = 28.964          ; molec wt dry air, g/mol
   Bc   = 1.380662e-23    ; Boltzmans constant, J/K
      
;- air density

   aden1 = 100.0 * p / (Rd * t)      ; air density, kg/m3    
   aden  = 1e3 * aden1 * Na / Md     ; air density, #/m3
   
;- Use the expressions in Turco et al. [1982]

   if (shape eq 0) then phi_ls = 1.0                                      ; sphere
   if (shape eq 1) then phi_ls = (2/3.)^(1/3.) *AR^(1/6.)                 ; Table 1, cylinder      
   if (shape eq 2) then phi_ls = (!pi/(9.*tan(!pi/6.)))^(1/3.) *AR^(1/6.) ; Table 1, hexagon
   
   beta   = sqrt(0.75 + 1/(4*AR) + !pi/8.) * sqrt(0.5 + AR/2. + !pi/8.)  ; eqn 10
   
   phi_ld = phi_ls / beta
   
   v_fall = Di * g * rad_m / (2. * aden) * sqrt( !pi / (2.*Ma *Bc * T)) * phi_ld  ; m/s, eqn 11
   
;- done

   return,v_fall
   
   end