function glaze%(u0,r0,n0,thetap0)
% set initial variables
%clear 
 
format long
% % % 
  u0=350;        %eruption velocity [m/s]
  r0=150; %100;%           vent size [m]
  n0=0.045; %           initial gas mass fraction
  thetap0 = 1200; %           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 =500; %           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=40000;
  nstep=ceil(z_end-z0)/hstep;
  r_vec=zeros(nstep,1);
  t_vec=zeros(nstep,1);
%set some addtitional values
 
     % heighi =  heighw ;
      ep = r_d/r_v ;%0.622
      cdcof = 0.75 ;
      u = u0   ;
      r = r0   ;
      z = z0  ;
      thetap = thetap0  ;
 
 
     [thetaa,p]=usstand_nested(z);
 
%set initial conditions
      rho_d = p0 /(r_d*thetap);
      phi_d = 0;
      phi_l = 0;
      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_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; 
     
  % 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;
    rho_b0=rho_b;
    m_v0=m_v;
    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
    
    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 ;
    
igasthrust = 1   ; 
 
  
      y97(1) = m_d;    
      y97(2) = m_v;    
      y97(3) = m_l;    
      y97(4) = v;              
      y97(5) = t;      
 
      y98(1) = m_d;
      y98(2) = v;       
      y98(3) = t;      
                    
       
                
  for istep=1:nstep    
 
      
options = odeset('Refine', 1,'RelTol',1e-1,'AbsTol',1e-1); %'MaxStep',0.1);
[zh,y97n] = ode45(@glaze_97_nested,[z z+hstep],y97,options);
  
options = odeset('Refine', 1,'RelTol',1e-1,'AbsTol',1e-1); %'MaxStep',0.1);
[zh,y98n] = ode45(@glaze_98_nested,[z z+hstep],y98, options);
 
 
    if igasthrust ==1
        
               
        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);        
        t    = y98n(end,3);         
                    
   if y97n(4)- y98n(2)>0 && y97n(5)- y98n(3)>0  
        igasthrust =0;
   end          
                    
    else
 
        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);                
    end
    
  % calculate u, r, etc at new height
 
         m = m_d+ m_v + m_l + m_s;   
         u = v/ m;
                    if u <= 5
                    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 t<=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
        rcut = (rn-r)/r;
        if r<=0;
            rcut=0;
        end
        r = rn;
        r_vec(istep) = r;
        thetap = t / ( rho_b*u*r^2*c_b ); 
        if rho_b<=0;
            thetap=0;
        end
        if u<=5;
            thetap=0;
        end
        if r<=0;
            thetap=0;
        end
        if c_b<=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 );   
        if r_v <=0;
            phidpphiv =0 ;
        end
        if t <=0;
            phidpphiv =0 ;
        end
     
            
        
        phi_l= 1/rho_l*(p/r_v*t)*(m_l*m*c_b)/(m_v+ep*m_d);
        if r_v<=0;
            phi_l=0;
        end
       if t<=0;
            phi_l=0;
        end 
        phi_s = 1-(phidpphiv)-phi_l ;
 
        t_vec(istep)=t;
        
        
      
        if igasthrust ==1
    
    y98(1)=y98n(end,1);
    y98(2)=y98n(end,2);
    y98(3)=y98n(end,3);
    y97(1)=y98n(end,1);
    y97(2)= m_v;
    y97(3)= 0;
    y97(4)=y98n(end,2);
    y97(5)=y98n(end,3) ;
    
else
    
       y97(1)=y97n(end,1);
       y97(2)=y97n(end,2);
       y97(3)=y97n(end,3);
       y97(4)=y97n(end,4);
       y97(5)=y97n(end,5); 
       
        end
 
 
       z = z+hstep;
            
  end
 
%   %figure(1);clf
% 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  
%  colorbar('peer',axes,...
hold on
%colorbar('location','east')
%caxis([0 20])
cbh=colorbar('peer',gca,...
    [0.852213541666667 0.234071093226023 0.0260416666666666 0.385647216633132]);
label=get(cbh,'title');
set(label,'string','Temperatur');
%set(label, 'title', 'String', 'Strahlung'); 
set(cbh,'location','east');
%set(cbh,[35 40 15 -10]);

 hold on


%[0.885515873015873 0.127801621363853 0.0158730158730159 0.400095374344301]);
% % %cb=colorbar
 hold all

backgroundImage = importdata('vulkan_Bild.jpg');
%select the axes
%axes(handles.Image);
%place image onto the axes
image([-40 40 ],[35 -15],backgroundImage);
hold on
% 
if max(r_vec) > 8000
brennend_image = importdata('brennendes_haus.jpg');
%select the axes
%axes(handles.Image);
%place image onto the axes
image([-15 -23 ],[3 -3],brennend_image)

else max(r_vec) <= 8000
 nichtbrennend_image = importdata('nicht_brennend1.jpg');
%select the axes
%axes(handles.Image);
%place image onto the axes
image([-15 -23 ],[3 -3],nichtbrennend_image)

end
    

hold on

%  Auto-generated by MATLAB on 16-Jan-2012 14:42:03
hold all
% Create colorbar

 
  z_vec = (hstep:hstep:z)';
 


 r_vec=r_vec/1000;
 z_vec=z_vec/1000;
 t_vec=t_vec/1000;
 
  plot (r_vec(r_vec~=0),z_vec(r_vec~=0),'r-')
  plot (-r_vec(r_vec~=0),z_vec(r_vec~=0),'r-')
  z2 =[z_vec; z_vec(end:-1:1)];
  r2= [r_vec;-r_vec(end:-1:1)];
 
 t3=[t_vec;t_vec(end:-1:1)];

  hold all
 % axis ([-35000 35000 -5000 35000]);
  axis([-40000/1000 40000/1000 -15000/1000 35000/1000])
  %set(gca,'xTickLabel',{'30','20','10','0','10','20','30'} )
  hold on
  xlabel('Breite (km)')
  ylabel('Höhe über Kraterradius (km)')
  set(gca,'xtick',[-40,-30,-20,-10,0,10,20,30,40],'xLim',[-40 40],...
      'XTickLabel',{'40','30','20','10',0,'10','20','30','40'});
  
  % ------------   Der Elipsenabschluss beim Plume  -----------------
  
  indx=find(r_vec);  % bringt mir nur die positiven Werte von r_vec hervor
  z_indx = z_vec(indx);
  z_plot = z_indx(end);
  t=linspace(0,2*pi,500); %von 0 bis 2*pi mit 500 schritten
  x= 0;
  y=z_plot;
  a=r_vec;
  b=1/10*r_vec; 
  plot(x-a*cos(t),y-b*sin(t),'Color',[.85 0 0]);  %dunkelrot
    indxnull= r_vec==0;  % sucht mir alle werte für r_vec=0 raus
    
%plot(z_vec(indxnull),t_vec(indxnull),'color','white' )
%t4=t3(r_vec~=0);
r3=r2(r2>0);
z4=z2(r2~=0);
r4=r2(r2~=0);



 hold on
%  -------------  Färbung des Kegels ------------------------------
 patch(r2(r2~=0),z2(r2~=0),t3(r2~=0),'edgecolor','none');
 hold on % t3(r2~=0)
% plot(z_vec(indxnull),'color','white' );
% z_vec(r_vec==0)=0;
%  pl=plot(z_vec(indxnull))
%  set(pl,'color','white')
%plot(z2(indxnull),t3(indxnull),'color','white' );% ich habe eigentlich einen langen...
  % strich bei r=0, von z=-50000 bis z_ende. Um den wegzubekommen plotte
  % ich alle z_vec-werte, die bei r_vec=0 sind weiss. so ist der strich
  % weg.
%   % ----------------------  Färbung des Dreiecks --------------------------
% h_vulkan = 5.000;
%  
% hold on
%  
% 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];
%    
%   h_vulkan = 5.000;
%  
% hold on
%  
% 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+ 2.200       -h_vulkan+ 2.200;
%     -h_vulkan+ 2.200      -h_vulkan+ 2.200     0                     0;
%     -h_vulkan            -h_vulkan+ 2.200  -h_vulkan+ 2.200       -h_vulkan+ 2.200];
%  
%  z = [0.001 0.001 0.001 0.001;
%       0.001 0.001 0.001 0.001;
%       0.001 0.001 0.001 0.001];
%    
%   patch(x,y,z,'FaceColor','black')
%   hold on
  
  
  % --------------  Färbung des schwarzen Balkens ------------------------
  
%   x =[-40000                  -500             40000                     500 ;          
%       -40000                  -500             40000                     500;
%       -500                   -40000              500                  40000];
%  
% y =  [-h_vulkan            -h_vulkan           -h_vulkan         -h_vulkan;
%       -h_vulkan+ 2200      -h_vulkan+ 2200     -h_vulkan+2200    -h_vulkan+2200;
%       -h_vulkan            -h_vulkan+ 2200     -h_vulkan         -h_vulkan+ 2200];
%  
% z=[1 1 1 1
%   1 1 1 1
%   1 1 1 1];
% 
%  
%   x =[-40.000                  -.500             40.000                     .500 ;          
%       -40.000                  -.500             40.000                     .500;
%       -.500                   -40.000              .500                  40.000];
%  
% y =  [-h_vulkan            -h_vulkan           -h_vulkan         -h_vulkan;
%       -h_vulkan+ 2.200      -h_vulkan+ 2.200     -h_vulkan+2.200    -h_vulkan+2.200;
%       -h_vulkan            -h_vulkan+ 2.200     -h_vulkan         -h_vulkan+ 2.200];
%  
%   z = [0.001 0.001 0.001 0.001;
%       0.001 0.001 0.001 0.001;
%       0.001 0.001 0.001 0.001];
% 
%     
% patch(x,y,z,'FaceColor','black')
%  
%  
% hold on
% h_vulkan = 5.000;
%   
%  
% 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+2.200      -h_vulkan+2.200;
%     -h_vulkan+2.200        -h_vulkan+2.200       0                       0;
%     -h_vulkan             -h_vulkan+2.200  -h_vulkan+2.200      -h_vulkan+2.200];
%  
%  z = [0.001 0.001 0.001 0.001;
%       0.001 0.001 0.001 0.001;
%       0.001 0.001 0.001 0.001];
% 
%  
%  
% patch(x,y,z,'FaceColor','black')
%  
% hold on
% x = [-.300                  -.300;
%      +.300                   .300;
%      -.300                   .300];
%   
% y = [ -h_vulkan+2.200         .500;
%       -h_vulkan+2.200        -h_vulkan+2.200;
%          .500                    .500];
% 
%  
%   z = [.001 .001 ;
%       .001 .001 ;
%       .001 .001 ];
%   
% %   x = [-300                  -300;
% %      +300                   300;
% %      -300                   300];
% %   
% % y = [ -h_vulkan+2200         500;
% %       -h_vulkan+2200        -h_vulkan+2200;
% %          500                    500];
% %  
% %   z = [1 1 ;
% %       1 1 ;
% %       1 1 ];
% 
%  patch(x,y,z,'FaceColor','red','edgecolor','none')
%  
%  G1=imread('brennendes haus.jpg')
%  image(G1);% axis image
% set(gca,'color','none')
 
 
 
%   --------------- Die Magmakammer ----------------------
%  f=inline('3*sin(2*x.^2)');
% A=-pi/2;
% B=pi/2;
% x=linspace(A,B,200);
% y=f(x);
% 
% area(x,y);
% 
% x=linspace(0,2*pi);
% y=-5 +sin(x)%cos(x);
% area(x,y);

%  t=linspace(0,2*pi,'r'); 
%   x= 0;
%   y=-5.000+1.400;
%   a=1.900;
%   b=1.100; 
%   area(x-a*cos(t),y-b*sin(t));
%   h=area(x-a*cos(t),y-b*sin(t));
%   set(h,'FaceColor',[.7 0 0]); %dunkelrot

%     
%-------------------------  Das Haus  mit Dach -------------------------------
% x = [14.000          15.700;      
%      15.700          15.700;        
%      14.000          14.000];
%  
%  
%  
% y = [-h_vulkan+2.200    -h_vulkan+2.200;       
%      -h_vulkan+2.200     -h_vulkan+4.500;
%      -h_vulkan+4.500     -h_vulkan+4.500];
%  
%  z= [.001 .001;
%     .001 .001;
%     .001 .001];

% x = [14000          15700;      
%      15700          15700;        
%      14000          14000];
%  
%  
%  
% y = [-h_vulkan+2200    -h_vulkan+2200;       
%      -h_vulkan+2200     -h_vulkan+4500;
%      -h_vulkan+4500     -h_vulkan+4500];
%  
%  z= [1 1;
% %     1 1;
% %     1 1];
% % if  r_vec > 1500
%  patch(x,y,z,'FaceColor','black','edgecolor','none')
% %  else
% %      patch(x,y,z,'FaceColor','none','edgecolor','none')
% % end
%  
% x=[ 14.000
%     15.700
%     14.850];
%     
% y=[-h_vulkan+4.500
%    -h_vulkan+4.500
%    -h_vulkan+6.000];
%  
% z =[ .001
%      .001
%      .001];
% %  
% %  x=[ 14000
% %     15700
% %     14850];
% %     
% % y=[-h_vulkan+4500
% %    -h_vulkan+4500
% %    -h_vulkan+6000];
% %  
% % z =[ 1
% %      1
% %      1];
% %  if  r_vec > 1500
%  patch(x,y,z,'FaceColor','green','edgecolor','none')
% %  else
%      patch(x,y,z,'FaceColor','none','edgecolor','none')
%  end
 
% img=imread('nicht_brennend.jpg');
% imagesc(img)
 
 
 set(gca,'YDir','normal')  % beim Anzeigen im Gui wird das Bild umgedreht;
 % um das zu verhindern muss die Y-Achse hier angepasst werden fürs
 % GUI-Bild
   
 if v<=0
    v=0;
end
 
if isnan(v)||isinf(v) 
    disp('break')
     return
end
 
 
 if isnan(u)||isinf(u) 
    disp('break')
     return
 end
 
 
   if u == 0
   disp('break')
    return
   end
   
  if r<=0
    r=0;
  end
 
 if isnan(r)|| isinf(r)
     disp('break')
     return
 end
 
 
 
  if isnan(rn) || isinf(rn)
     disp('break')    
     return
  end
  
   if m<=0
    m=0;
   end
 
if isnan(m) || isinf(m)
     disp('break')
     return
end
 
 
  
  % -----------------  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
    if m1<=0;
        c_b=0;
    end
 
    [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))*(1/(wa+ep) ) )/ phi_ad;
    if r_v<=0;
        rho_ad=0;
    end
    if phi_ad<=0;
        rho_ad=0;
    end
    rho_av = ep*rho_ad; %( (p/( r_v*thetaa))*(wa/(wa+ep)) )/ phi_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
    if m11<=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;
    if c_b<=0;
        f1=0;
    end
    g1=b1/(m1*c_b);
    if m1<=0;
        g1=0;
    end
    if c_b<=0;
        g1=0;
    end
    h11=(a1*b1)/c_b; 
    if c_b<=0;
        h11=0;
    end
    i1=y(3)/rho_l;
    if rho_l<=0;
        i1=0;
    end
    j1=m_s/rho_s;
    if rho_s<=0;
        j1=0;
    end
   
 
    dy(1)=2*al*rho_ad*phi_ad*(sqrt (a1 *y(4)* c1));
    w1=(a1*y(4)*c1);
    if w1 <=0
        dy(1)=0;
    end
    dy(2)=2*al*rho_av*phi_av*(sqrt(a1 *y(4)*c1)) - omega*y(2)* d1;
     if w1 <=0
        dy(2)=0;
    end
    dy(3)=omega*y(2)*d1;   
    dy(4)=e1*(rho_ab* h11-m11);
    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 w1 <=0
        dy(5)=0;
    end
    
   
   
end 
 
   %  --------------   mit gasschubregion  ----------------------
   
   function dy = glaze_98_nested(z,y)
       
    dy=zeros(3,1);
 
    m2 = y(1)+m_v+m_l+m_s;
   
  
    c_b = (y(1)*c_d + m_v*c_v + m_l*c_l + m_s*c_s)/m2; 
 
    [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))*(1/(wa+ep) ) )/ phi_ad;
    if r_v<=0;
       rho_ad=0;
    end
    if thetaa<=0;
        rho_ad=0;
    end
    rho_av = ep*rho_ad;%( (p/( r_v*thetaa))*(wa/(wa+ep)) )/ phi_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))/((y(1)*c_d+m_v*c_v+m_l*c_l+m_s*c_s)/m2)-m2^2) );
       
    dy(3) =0.125*( rho_ab*y(2) )^0.5 * c_ab * thetaa - rho_ab*g*...
        ( ( (r_v*y(3) )/p) *(m_v+ep*y(1)) /(m2*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