% set initial variables
clear all
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 = 5; %           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

%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(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);
      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 ;
      erupra = r^2*u*3.14159*rho_b ;
      thermflux = erupra*(thetap-thetaa)*c_b;
  % set initial condition for diff. equations, for indexes of equations see
%subroutine glaze
      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
      
    m = m_d + m_v +m_l +m_s;
 
  for k=0:10:4000

options = odeset('Refine', 1,'RelTol',1e-2,'AbsTol',1e-2) %'MaxStep',0.1);
[zh,y97n] = ode45(@glaze_97,[z z+k+hstep],y97, options); 

    
        m_d = y97n(end,1) %y97n(end,1);
        m_v = y97n(end,2);
         m_l = y97n(end,3);
         v   = y97n(end,4);
         t   = y97n(end,5);

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

  
 h(1)=subplot(5,1,1);
 plot(zh,y97n(:,1))
 ylabel('m_d')
 xlabel('z')
 h(2)=subplot(5,1,2);
 plot(zh,y97n(:,2))
 ylabel('m_v')
 xlabel('z')
 h(3)=subplot(5,1,3);
 plot(zh,y97n(:,3))
 ylabel('m_l')
 xlabel('z')
 h(4)=subplot(5,1,4);
 plot(zh,y97n(:,4))
 ylabel('v')
 xlabel('z')
 h(5)=subplot(5,1,5);
 loglog(zh,y97n(:,5))
 ylabel('log T')
 xlabel('log z')