function glaze_programm
% set initial variables
clear 

format long

u0=200; %            eruption velocity [m/s]
r0=100; %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 =100; %           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

  
  z_end=10000;
  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;

% set initial values for Runge Kutta   
    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); 
    m_s = rho_s*u*r^2*phi_s;
    m = m_d+m_v+m_l+m_s; 
     
  % 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;
    c_b = ( m_d*c_d+m_v*c_v+m_l*c_l+m_s*c_s ) / m ;
    erupra = r^2*u*3.14159*rho_b ;
    thermflux = erupra*(thetap-thetaa)*c_b;

    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;
      
      if v <= 0
      v=0;
      end
      
      y97(5) = t; %t(1)= rho_b*u*r^2*c_b*thetap ;
      
      if t<=0
      t=0;
      end
      
      %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);
             if m_v <= 0
             m_v = 0;
             end       
             
        m_l = y97n(end,3);
        
        v = y97n(end,4);
            if v <= 0
            v=0;
            end
            
         t = y97n(end,5);
            if t <= 0
            t=0;
            end
        
        
% calculate u, r, etc at new height
      m = m_d+ m_v + m_l + m_s;
            if  m<=0
            m=0;
            end
            
      u = v/ m;
           if u <= 5
           u=0;
           end
             if m <= 0
             u=0;
             end
      
      
      c_b = ( m_d*c_d + m_v*c_v + m_l*c_l + m_s*c_s ) / m;
          if m <= 0
          c_b =0;
          end
          
      %rho_b = rho_d*phi_d + rho_v*phi_v + rho_l*phi_l + rho_s*phi_s;
      rho_b = p /(r_v.*t).* (m.^2 .* c_b ./( m_v + ep.*m_d ));
      rn = ( m.^2 ./ ( rho_b.*v) ).^0.5;
          if rho_b<=0
          rn=0;
          end
          
      if v <=0
      rn=0;
      end
      
          if rn <= 0
          rn=0;
          end
      
      rcut = (rn-r)./r;
          if rcut <=0
          rcut=0;
          end
      
      r = rn;
         if r<=0
          r = 0;
         end
      
      thetap = t ./ ( rho_b.*u.*r.^2.*c_b ); % da stand un anstatt u??????????????????????
         if rho_b==0
         thetap = 0;
         end
         
      if u== 0
      thetap=0;
      end
      
         if r == 0
         thetap = 0;
         end
         
      if c_b==0
      theta_p=0;
      end
      
         if thetap <=0
         thetap=0;
         end
         
      phidpphiv = 1 - p./(r_v.*t).*(m.*c_b)./(m_v+ep.*m_d).*( m_l./rho_l+m_s./rho_s );
      phi_l= 1/rho_l.*(p/r_v.*t).*(m_l.*m.*c_b)./(m_v+ep.*m_d);
      phi_s = 1-(phidpphiv)-phi_l ;

      z = z+hstep;
     
%plot erster versuch mit temperatur   

 %aa=-z:z; %imagesc(aa,r,t)
% plot(aa,r)
% hold on;
% plot (z,r)
% plot (-z,r)
% %set(gca,'YDir','Reverse')
% hold on
% 
h(1)=subplot(5,1,1);
   plot(zh,m_d)%  
   ylabel('m_d')
 xlabel('z')
 hold on
  h(2)=subplot(5,1,2);
plot(zh,m_v)%  
  ylabel('m_v')
  xlabel('z')
  hold on
 
  h(3)=subplot(5,1,3);
plot(zh,m_l)% 
  ylabel('m_l')
  xlabel('z')
  hold on
  h(4)=subplot(5,1,4);
plot(zh,v)%  
  ylabel('v')
  xlabel('z')
  hold on
  h(5)=subplot(5,1,5);
plot(zh,t);%  
  ylabel('T')
  xlabel('z')
  hold on

  if r==0
      stop glaze_programm2;
  end


  end
  
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);

    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;
      
    a1=r_v*y(5)/p;
    b1=y(2)+ep*y(1);
    m11=m1^2;
    c1 = b1/(m11*c_b); 
    
         if c_b<=0
         c1 = 0;
         end
         
    d1 = m1/y(4);
         if y(4)<=0
         d1 = 0;
         end
        
    e1=g/y(4);
         if y(4)<=0
         e1 = 0;
         end
         
    f1=a1*b1/c_b;
    g1=b1/(m1*c_b); 
    
    h11=(a1*b1)/c_b;
         if c_b <= 0
         h11 = 0;
         end
         
    i1=y(3)/rho_l;
    j1=m_s/rho_s;
% 
    dy(1)=2*al*rho_ad*phi_ad*(sqrt (a1 *y(4)* c1));
        if y(4)<=0
        dy(1) = 0;
        end
        
    dy(2)=2*al*rho_av*phi_av*(sqrt(a1 *y(4)*c1)) - omega*y(2)* d1;
   
        if dy(2) <= 0
        dy(2) = 0;
        end
    
    dy(3)=omega*y(2)*d1;
        if y(4)<=0
        dy(3) = 0;
        end
        
    dy(4)=e1* (rho_ab* h11-m11);
        if y(4)<=0
        dy(4) = 0;
        end
   if y(5)<=0
   dy(4) = 0;
   end
     
    dy(5)= 2*al*rho_ab*c_ab*thetaa*(sqrt(a1*y(4)*c1))+ l*omega*y(2)*d1 -rho_ab*g*((a1*g1)-(i1+j1)); 
        if y(4)<=0
        dy(5) = 0;
        end
   if y(5)<=0
   dy(5) = 0;
   end
    
end

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

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

      ta0 = 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 = ta0-xmu1*zkm;
        t1=te/ta0;
        pr = p0*((t1)^(exp1));
       elseif ( (zkm >= h1) && (zkm <= h2) )
        te = ta0-xmu1*h1;
        t2=te/ta0;
        ab1=1000/(r_d*te);
        pr = p0*(t2^(exp1))*exp(g*(h1-zkm)*ab1); % dexp
      elseif  zkm > h2 ; 
        te = ta0-xmu1*h1-xmu2*(zkm-h2);
        ab2=1000/(r_d*(ta0-xmu1*h1));
        c1 = exp(g*(h1-h2)*ab2); %dexp
        ab3=te/(ta0-xmu1*h1);
        c2 = (ab3)^(exp2);
        ab_1=(ta0-xmu1*h1)/ta0;
        pr = p0*(ab_1)^(exp1)*c1*c2;
      end
      
      if zkm > 50
          print 'Plume to high'
      end 
     wa=0;
   

end    
      
end