function glaze%(u0,r0,n0,thetap0)
% set initial variables
%clear 

format long

 u0=200; %            eruption velocity [m/s]
 r0=50; %100;%           vent size [m]
 n0=0.02; %           initial gas mass fraction
 thetap0 = 800; %           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=20000;
  nstep=ceil(z_end-z0)/hstep;
  r_vec=zeros(nstep,1);
%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; 
    t= rho_b*u*r^2*c_b*thetap ;
    
  for istep=1:nstep    

      y97(1) = m_d;    
      y97(2) = m_v;    
      y97(3) = m_l;    
      y97(4) = v;      
                if v <= 0
                v=0;
                end
      y97(5) = t;      
                if t<=0
                t=0;
                end
      
      
      
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
 
         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 u <= 0
                        v=0;
                    end
                         if m <= 0
                         u=0;
                         end
                         if v<=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 = p /(r_v.*t).* (m.^2 .* c_b ./( m_v + ep.*m_d ));
         if r_v <= 0
             rho_b =0;
         end
         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
                    if r<=0
                        rcut=0;
                    end
          r = rn;
          r_vec(istep) = r;
          
                    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 ;     

%         hold on;
%         plot (z,r)
%         plot (-z,r)

 
  
      y98(1) = m_d;
      y98(2) = v;
                if v <= 0
                v=0;
                end
      y98(3) = t;      
                if t<=0
                t=0;
                end
           
        
options = odeset('Refine', 1,'RelTol',1e-2,'AbsTol',1e-2); %'MaxStep',0.1);
[zh,y98n] = ode45(@glaze_98_nested,[z z+hstep],y98, options);

              
        m_d  = y98n(end,1); %hier setze ich den letzten Eintrag der Berechnung wieder als Startwert für den nächsten Runge-Loop ein
        
        v  = y98n(end,2);
                    if v <= 0
                    v=0;
                    end
        t = y98n(end,3);
                    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 = p /(r_v.*t).* (m.^2 .* c_b ./( m_v + ep.*m_d ));
         if r_v<=0
             rho_b=0;
         end
         if t<=0
             rho_b=0;
         end
         
         
         rn = ( m.^2 ./ ( rho_b.*v) ).^0.5;
                    
                                if rn <= 0
                                rn=0;
                                end
                                if r_v <=0
                                    rn=0;
                                end
                               if m_v<=0
                                rn=0;
                               end
                               if rho_b<=0
                    rn=0;
                    end
                        if v <=0
                        rn=0;
                        end
  
          rcut = (rn-r)./r;
          if r<=0
              rcut=0;
          end
                    if rcut <=0
                    rcut=0;
                    end
          r = rn;
                    if r<=0
                        disp('break');
                        return
                    r = 0;
                    end
          thetap = t ./ ( rho_b.*u.*r.^2.*c_b ); % da stand un anstatt u??????????????????????
          if c_b<=0
              thetap =0;
          end
          if r <= 0
              thetap = 0;
                                     end        
          if rho_b <=0
                    thetap = 0;
                    end
                            if u<= 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 ;

  if y97n(4)- y98n(2)>0 & y97n(5)- y98n(3)>0  
      
      
        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
    
  else
    
   m_d  = y98n(end,1); %hier setze ich den letzten Eintrag der Berechnung wieder als Startwert für den nächsten Runge-Loop ein
        
        v  = y98n(end,2);
                    if v <= 0
                    v=0;
                    end
        t = y98n(end,3);
                    if t <= 0
                    t=0;
                    end       
  end
       z = z+hstep;
     


  end
  hold on;
  z_vec = (hstep:hstep:z)';
  plot (r_vec,z_vec,'r-')
  plot (-r_vec,z_vec,'r-') 
  
  z2 = [z_vec; z_vec(end:-1:1)];
  r2 = [r_vec; -r_vec(end:-1:1)];
  patch(r2,z2,ones(size(r2)),'r');
  h_vulkan = 5000;
  

hold on
h_vulkan = 5000;
  

x =[-h_vulkan              h_vulkan       -h_vulkan           h_vulkan ;          
   -h_vulkan               h_vulkan          -200             200;
   h_vulkan               -h_vulkan          -200             200];

y =[-h_vulkan             -h_vulkan      -h_vulkan+2200     -h_vulkan+2200;
    -h_vulkan+2200        -h_vulkan+2200     0                 0;
    -h_vulkan             -h_vulkan+2200  -h_vulkan+2200     -h_vulkan+2200];



z = [1 1 1 1;
    1 1  1 1;
    1 1  1 1];
    

patch(x,y,z,'FaceColor','black')

hold on
x = [-200                  -200;
     +200                   200;
     -200                   200];
  
y = [ -h_vulkan+2200         0;
      -h_vulkan+2200        -h_vulkan+2200;
         0                   0];
 
 z= [1 1;
    1 1;
    1 1];
 
 patch(x,y,z,'FaceColor','red')


  axis equal

% 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

  % -----------------  ohne Gasschub  ------------------------------
  
  
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 
 
   %  --------------   mit gasschubregion  ----------------------
   
   function dy = glaze_98_nested(z,y)
       
    dy=zeros(3,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;
    m = y(1)+m_v+m_l+m_s; % gesamtmasse in column
  
    c_b = (y(1)*c_d+m_v*c_v+m_l*c_l+m_s*c_s)/m; % 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;
      
  
% 
    dy(1)= 0.125*(rho_ab*y(2))^0.5;
             
    dy(2)= g/y(2)*(rho_ab*r_v*y(3)/p*(m_v+ep*y(1))/c_b-m^2);
    
    dy(3) =0.125*(rho_ab*y(2))^0.5*c_ab*thetaa-rho_ab*g*((m_v+ep*y(1))/(m*c_b)-m_s/rho_s);
   
     
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