% calculate the atmospheric temperature and presssure for a
% US Standart atmosphere

function [te,pr] = usstand(z)
global phi_av phi_ad r_ab c_ab rho_ad rho_ab rho_av wa m_s c_d c_l c_v c_s al r_v ep  omega g  c_b  l rho_l rho_s r_d 
%global g omega al l r_v r_d rho_l rho_s c_d c_v c_s c_l hstep

      t0 = 288.86;
      p0 = 1.01325e+5;
      h1= 11;
      h2= 20;
      xmu1 = 6.5;
      xmu2 = -2.0;

      zkm = z / 1000 ;
      exp1 = g*1000  / (r_d*xmu1);
      exp2 = g*1000 / (r_d*xmu2);

      if  zkm < h1 
        te = t0-xmu1*zkm;
        pr = p0*((te/t0)^(exp1));
       elseif ( (zkm >= h1) && (zkm <= h2) )
        te = t0-xmu1*h1;
        pr = p0*((te/t0)^(exp1))*exp(g*(h1-zkm)*1000.d00/(r_d*te)); % dexp
       else  zkm > h2 ; 
        te = t0-xmu1*h1-xmu2*(zkm-h2);
        c1 = exp(g*(h1-h2)*1000.d00/(r_d*(t0-xmu1*h1))); %dexp
        c2 = (te/(t0-xmu1*h1))^(exp2);
        pr = p0*((t0-xmu1*h1)/t0)^(exp1)*c1*c2;
      end
      
      if zkm > 50
          print 'Plume to high'
      end 
          wa=0;
     
 %calculate volume fractions and densities
      phi_av = wa / ( wa + ep );
      phi_ad = 1 - phi_av ;
      r_ab = ( r_d + wa*r_v ) / ( 1+wa );
      c_ab = ( c_d + wa*c_v ) / ( 1+wa );

      rho_ad = pr / ( r_v*te*(wa+ep) ) / phi_ad;
      rho_av = ep*rho_ad;
      rho_ab = rho_av*phi_av + rho_ad*phi_ad;

   
end