pro pmc_0d_model2a,z,p,t,h2o,rad,vpop,sfac, t_ice,p_ice,s_ice,h2o_sat,v_ice,m_ice,h2o_ice,Ztop,Zmax,Zbot,iwc,rc

;-------------------------------------------------------
;  purpose:
;     Routine computes various properties of ice assuming bulk thermodynamic equilibrium. 
;     The user inputs profile data and gets back various profiles plus some scalars.     
;
;  input:
;     z.......altitude (km), fltarr(nz)
;     t.......temp (K), fltarr(nz)
;     p.......pressure (mb), fltarr(nz)
;     h2o.....h2o mixing ratio (ppmv), fltarr(nz)
;
;     rad.....ice particle radius (nm) for the Kelvin effect, fltarr(nz), if rad=0 the Kelvin effect is ignored.
;              The Kelvin effect here is the vapor pressure ratio according to eqn 6-38 
;              in Pruppacher & Klett, "Microphysics of Clouds and Precipitation"
; 
;     vpop....saturation vapor pressure option:
;              1 = use Murphy & Koop [2005], recommended by Rapp & Thomas [2006]  (very close to Marti & Mauersberger [1993])
;              2 = use Mauersberger & Krankowsky [2003], here the expression for T < 169K is used
;              3 = use Marti & Mauersberger [1993]
;
;     sfac....a scale factor to apply to the calculated mass densities (affects h2o_ice,m_ice,v_ice,iwc),  
;             set sfac = 1 to do no scaling.  
;
;     Note: the input Z,T,P,H2O,rad,ARo must all have the same dimension
;
;  output:
;     t_ice......frost point temperature (K), fltarr(nz)
;     p_ice......saturation vapor pressure over ice (mb), fltarr(nz)
;     s_ice......saturation ratio WRT ice (unitless), fltarr(nz)
;     h2o_sat....saturation H2O mixing ratio (ppmv), fltarr(nz)
;     v_ice......ice volume density (um3 / cm3), fltarr(nz)
;     m_ice......ice mass density (ng / m3), fltarr(nz)
;     h2o_ice....gas phase equivalent H2O in ice (ppmv), fltarr(nz)
;     ztop.......ice layer top altitude (km)
;     zmax.......ice layer mass density peak altitude (km) (same as SOFIE Zmax, which are based on IR extinction)
;     zbot.......ice layer bottom altitude (km)
;     iwc........column ice abundance (ug/m2) 
;     rc.........critical radius (nm)
;
;  Source:    Mark Hervig, GATS Inc.
;
;  Revision:  4/3/09
;             9/20/2011,  added rc
;             3/21/12: MH removed the capacitance calculation,  added "sfac"
;
;-------------------------------------------------------------   
   
;- Constants

   Mww = 18.0        ; molec wt. h2o, g/mol
   di  = 0.93        ; density of ice, g/cm3
   R   = 8.314       ; J/mol/K
   Sti = 0.12        ; surface tension of ice in the presence of air, J/m2

;- Housekeeping
   
   nz = n_elements(z)
   
   if (vpop lt 1 or vpop gt 3) then stop,'vpop not valid in pmc_0D_model.pro'

   t_ice   = dblarr(nz) 
   p_ice   = dblarr(nz)
   s_ice   = dblarr(nz)
   v_ice   = dblarr(nz)
   m_ice   = dblarr(nz)
   h2o_sat = dblarr(nz) + 999
   h2o_ice = dblarr(nz)
   
   kelv     = dblarr(nz) + 1d0   
   rc       = dblarr(nz)
   
   Zmax = 0.
   Ztop = 0.
   Zbot = 0.
   iwc  = 0.
   
;- Loop over altitude

   for i = 0,nz-1 do begin

    if (t(i) lt 180 and t(i) gt 0) then begin

;-   Compute the Kelvin term

     if (rad(i) gt 0) then kelv(i) = exp(1d3 * 2. *Mww *Sti / ( R *t(i) *di *rad(i)) )

;-   Saturation vapor pressure at T,  Kelvin term if r > 0
     
     if (vpop eq 1) then p_ice(i) = kelv(i) *0.01*exp(9.550426-5723.265/t(i)+3.53068*alog(t(i))-0.00728332*t(i))    
     if (vpop eq 2) then p_ice(i) = kelv(i) *0.01*10^(14.88-3059.0/t(i))
     if (vpop eq 3) then p_ice(i) = kelv(i) *0.01*exp(28.868 - 6132.935 / t(i))      
   
;-   Saturation mixing ratio, etc...

     h2o_sat(i)  = 1d6 *p_ice(i) / p(i)  ; saturation mixing ratio, ppmv
     
     s_ice(i) = h2o(i) / h2o_sat(i)      ; saturation ratio

;-   Critical radius

     if (s_ice(i) gt 1.2) then begin
       rc(i) = critical_radius_ice( t(i),s_ice(i) )   ; critical radius (nm)
     endif
        
;-   Equilibrium ice properties

     q_xs = h2o(i) - h2o_sat(i)    ; excess h2o mix ratio, ppmv      

     q_xs = q_xs * sfac   ; scale it by sfac
          
     if (q_xs gt 0.0) then begin     ; if saturated then go on
     
       n_xs = p(i)*1d2 *q_xs*1d-6/(R*t(i))  ; excess mols h2o per m3 air
       
         v_ice(i) = 1d6 * n_xs * Mww / di   ; ice volume density, microns3 / cm3     
         m_ice(i) = 1e9 * n_xs * Mww        ; ice mass density, ng/m3       
       h2o_ice(i) = q_xs                    ; H2O(ice), ppmv
       
     endif

;-   Find the frost point temperature, t_ice.  Use the fact that
;    a linear relationship exists between 1/T and alog10(p_h2o).
;    Do a linear interpolation between some hi and low temp.

     p_h2o = p(i) * h2o(i) *1d-6  ; h2o partial pressure, mb

     t1 = t(i) - 20.  
     t2 = t(i) + 20.  
          
     if (vpop eq 1) then begin
       p1 = 0.01*exp(9.550426-5723.265/t1+3.53068*alog(t1)-0.00728332*t1) 
       p2 = 0.01*exp(9.550426-5723.265/t2+3.53068*alog(t2)-0.00728332*t2) 
     endif
      
     if (vpop eq 2) then begin
       p1 = 0.01*10^(14.88-3059.0/t1) 
       p2 = 0.01*10^(14.88-3059.0/t2) 
     endif
     
     if (vpop eq 3) then begin
       p1 = 0.01*exp(28.868 - 6132.935 / t1)
       p2 = 0.01*exp(28.868 - 6132.935 / t2)
     endif
              
     dsds     = (alog10(p2) - alog10(p_h2o)) / (alog10(p2) - alog10(p1))
     t_ice(i) = 1. / ( 1/t2 - dsds * (1/t2 - 1/t1) )
          
    endif  ; if cold enough to bother

   endfor  ; loop over altitude

;- Find Zmax, Ztop, Zbot, IWC

   k = where(z gt 80 and z lt 95 and m_ice gt 0., nk)
   
   if (nk gt 0) then begin 
     
     l = where(m_ice(k) eq max(m_ice(k)) ) & l = k(l(0))
     
     zmax = z(l)       
     ztop = max(z(k))
     zbot = min(z(k))
     
     for i = 0,nk-1 do begin
      
       if ( z(k(i)) gt zmax-2.) then begin           
         dz  = abs( z(k(i)+1) - z(k(i)-1) )*0.5  ; layer spacing (km)       
         iwc = iwc + m_ice(k(i)) * dz            ; micro-g/m2    
       endif
        
     endfor
       
   endif
    
;- done

   return
   end