% Anfangsbedingungen:

global phi_av phi_ad r_ab c_ab rho_ad rho_ab rho_av

      phi_s = (1 - n0)*rho_v / ( n0*rho_s + (1-n0)*rho_v );
      phi_v = 1 - phi_s - phi_d - phi_l;


      m_d = rho_d*u*r^2*phi_d;
      m_v = rho_v*u*r^2*phi_v;
      m_l = rho_l*u*r^2*phi_l;
      m_s = rho_s*u*r^2*phi_s;
      m = m_d+m_v+m_l+m_s; 
      c_b = ( m_d*c_d+m_v*c_v+m_l*c_l+m_s*c_s ) / m ;
      rho_b = rho_d*phi_d + rho_v*phi_v + rho_l*phi_l + rho_s*phi_s;
      rhob_0 = rho_b;
      m_v0 = m_v;

      v = rho_b*u^2*r^2;
      t = rho_b*u*r^2*c_b*thetap ;


      y97(1) = m_d;
      y97(2) = m_v;
      y97(3) = m_l;
      y97(4) = v;
      y97(5) = t;
     
    m = m_d + m_v +m_l +m_s;
    
z0=0;
z=z0;
hstep=5;  
 
  for k=0:400
 
options = odeset('Refine', 1,'RelTol',1e-2,'AbsTol',1e-2) 
[zh,y97n] = ode45(@glaze_97,[z z+k+hstep],y97, options); 
 
    
        m_d = y97n(end,1);
        m_v = y97n(end,2);
        m_l = y97n(end,3);
        v   = y97n(end,4);
        t   = y97n(end,5);
 
 
% calculate u, r, etc at new height

      m = m_d + m_v + m_l + m_s;
      un = v / m;
      c_b = ( m_d*c_d + m_v*c_v + m_l*c_l + m_s*c_s ) / m;
      rho_b = p .* m.^2 .* c_b ./( r_v .* t .* ( m_v + ep.*m_d ));
      rn = ( m.^2 ./ ( rho_b.*v ) ).^0.5;
      m_d = rho_d*u*r^2*phi_d;
      m_v = rho_v*u*r^2*phi_v;
      m_l = rho_l*u*r^2*phi_l;
      m_s = rho_s*u*r^2*phi_s;
      m = m_d+m_v+m_l+m_s; 
      c_b = ( m_d*c_d+m_v*c_v+m_l*c_l+m_s*c_s ) / m ;
      rho_b = rho_d*phi_d + rho_v*phi_v + rho_l*phi_l + rho_s*phi_s;
      rhob_0 = rho_b;
      m_v0 = m_v;
      v = rho_b*u^2*r^2;
      t = rho_b*u*r^2*c_b*thetap ;
  end