function glaze18102011_2
% set initial variables
clear 
%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 
format long

u0=200; %            eruption velocity [m/s]
r0=100;%           vent size [m]
n0=0.03; %           initial gas mass fraction
thetap0 = 1000; %           eruption temperature [K]
z0= 0; %              height of vent above ground [m]
g= 9.81; %           acceleration due to gravity [m/s^2]
omega = 0;    %          time constant for condensation [s]
al = 0.09; %            entrainment factor
l = 0; %              latent heat of water [J/kg] 2.257d06
r_v = 461; %            gas constant for volcanic gas (water vapor) [J/kg/K]
r_d = 287; %            gas constant for dry air [J/kg/K]
rho_l= 1000; %      density of liquid water [kg/m^3]
rho_s = 1000; %           density of pumice [kg/m^3]
c_d = 998; %            specific heat for dry air at constant P [J/kg/K]
c_v = 2000; %           specific heat for volcanic gas at constant P (water vapor) [J/kg/K]
c_s = 920; %            specific heat for for solids at constant P [J/kg/K]
c_l = 4200; %           specific heat for liquid water at constant P [J/kg/K]
p0 = 1e+05; %             pressure at vent level
hstep = 50; %           intergration constant for Runge Kutta [m]
%50.d00             height interval for output in file
%2                  1 read from file, 2 US Standart Atmosphere
%0                 with gas thrus =1 without = 0

  % [phi_av,phi_ad, r_ab, c_ab, rho_ad, rho_av, rho_ab] = deal(0);
  
  z_end=4000;
  nstep=ceil(z_end-z0)/hstep;

  
  
%set some addtitional values
     % heighi =  heighw ;
      ep = r_d/r_v ;
      cdcof = 0.75 ;
      u = u0   ;
      r = r0   ;
      z = z0  ;
      thetap = thetap0  ;
    %  nbh = 0  ;

     [thetaa,p]=usstand_nested(z);
     
    
%set initial conditions
      rho_d = p0 /(r_d*thetap);
      phi_d = 0;
      phi_l = 0;

     % rho_v = p*n*rho_s / ( r_v*thetap*n*rho_s-p*(1-n) );
      rho_v = p0/(r_v*thetap);  % 
      phi_s = (1 - n0)*rho_v / ( n0*rho_s + (1-n0)*rho_v );
      phi_v = 1 - phi_s - phi_d - phi_l;

% se% initial values for Runge Kutta
      %y97=zeros(5,1);
      %y97n=zeros(5,1);
      %y97=zeros(1,5)
%     
    
     % v = rho_b*u^2*r^2;
    %  t = rho_b*u*r^2*c_b*thetap ;
      
      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;
     
  % set initial condition for diff. equations, for indexes of equations see
%subroutine glaze
rho_b = rho_d*phi_d + rho_v*phi_v + rho_l*phi_l + rho_s*phi_s;
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 ;
rhob_0 = rho_b;
 m_v0 = m_v;
erupra = r^2*u*3.14159*rho_b ;
thermflux = erupra*(thetap-thetaa)*c_b;

m_d=rho_d*u*r^2*phi_d;%m_d=zeros(nstep,1); 
m_v= rho_v*u*r^2*phi_v;%m_v=zeros(nstep,1); 
m_l=rho_l*u*r^2*phi_l;%m_l=zeros(nstep,1); 
v=rho_b*u^2*r^2;%v=zeros(nstep,1); 
t= rho_b*u*r^2*c_b*thetap ;%t=zeros(nstep,1); 
    
  for istep=1:nstep    

      y97(1) = m_d;  %m_d(1)=rho_d*u*r^2*phi_d;
      y97(2) = m_v;%m_v(1)= rho_v*u*r^2*phi_v;
      y97(3) = m_l; %m_l(1)=rho_l*u*r^2*phi_l;
      y97(4) = v; %v(1)=rho_b*u^2*r^2;
      y97(5) = t; %t(1)= rho_b*u*r^2*c_b*thetap ;
      %y98(1) = m_d
      %y98(2) = v(istep)
      %y98(3) = t
      
options = odeset('Refine', 1,'RelTol',1e-2,'AbsTol',1e-2); %'MaxStep',0.1);
[zh,y97n] = ode45(@glaze_97_nested,[z z+hstep],y97, options);
   
        m_d = y97n(end,1); %hier setze ich den letzten Eintrag der Berechnung wieder als Startwert für den nächsten Runge-Loop ein
        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;
      rcut = (rn-r)./r;
      r = rn;
      thetap = t ./ ( rho_b.*un.*r.^2.*c_b );
      phidpphiv = 1 - p.*m.*c_b./(r_v.*t.*(m_v+ep.*m_d)).*( m_l./rho_l+m_s./rho_s );
      phi_l = m_l.*rho_b ./ (m.*rho_l);
      phi_r = phidpphiv+phi_l;
      phi_s = 1 - phi_r;
      
      z = z+hstep;
      z_plot=z0:hstep:z_end;   
  end


%   
  h(1)=subplot(5,1,1);
 plot(z_plot,m_d,'k.-')%  plot(zh,y97n(:,1))
   ylabel('m_d')
 xlabel('z')
  h(2)=subplot(5,1,2);
plot(z_plot,m_v,'k.-')%  plot(zh,y97n(:,2))
  ylabel('m_v')
  xlabel('z')
  h(3)=subplot(5,1,3);
plot(z_plot,m_l,'k.-')%  plot(zh,y97n(:,3))
  ylabel('m_l')
  xlabel('z')
  h(4)=subplot(5,1,4);
plot(z_plot,v,'k.-')%  plot(zh,y97n(:,4))
  ylabel('v')
  xlabel('z')
  h(5)=subplot(5,1,5);
plot(z_plot,t,'k.-')%  loglog(zh,y97n(:,5))
  ylabel('log T')
  xlabel('log z')
 
function dy = glaze_97_nested(z,y)


dy=zeros(5,1);

m1 = y(1)+y(2)+y(3)+m_s; % gesamtmasse in column
c_b = (y(1)*c_d+y(2)*c_v+y(3)*c_l+m_s*c_s)/m1; % bulk specific heat

[thetaa,p,wa] = usstand_nested(z);


a1=r_v*y(5)/p;
b1=y(2)+ep*y(1);
c1 = b1/(m1^2*c_b);
d1 = m1/y(4);
e1=g/y(4);
f1=a1*b1/c_b;
g1=b1/m1^2*c_b;
h1=a1*b1/m1*c_b;
i1=y(3)/rho_l;
j1=m_s/rho_s;

    phi_av = wa / ( wa + ep );
      phi_ad = 1 - phi_av ;
      
      c_ab = ( c_d + wa*c_v ) / ( 1+wa );
r_ab = ( r_d + wa*r_v ) / ( 1+wa );


      rho_ad = p./( r_v.*thetaa.*(wa+ep) ) ./ phi_ad;
      rho_av = ep*rho_ad;
      rho_ab = rho_av*phi_av + rho_ad*phi_ad;
      %wa=rho_av.*phi_av./(rho_ad.*phi_ad)
      
      
dy(1)=2*al*rho_ad*phi_ad*(a1*y(4)*c1)^0.5;
dy(2)=2*al*rho_av*phi_av*(a1*y(4)*c1)^0.5-omega*y(2)*d1;
dy(3)=omega*y(2)*d1;
dy(4)=e1*(rho_ab*(f1)-m1^2);
dy(5)=2*al*rho_ab*c_ab*thetaa*(a1*y(4)*g1)^0.5+l*omega*y(2)*d1...
    -rho_ab*g*((h1)-(i1+j1));


end

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

function [te,pr,wa] = usstand_nested(z)

      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
      elseif  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 ;
%       
%       c_ab = ( c_d + wa*c_v ) / ( 1+wa );
% r_ab = ( r_d + wa*r_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;
%       %wa=rho_av.*phi_av./(rho_ad.*phi_ad)
end    
      
end