%     program parcel

%     TO CONVERT TO SI UNITS
      TICE = 273.15;
      pa_per_mb = 100.;

%     INITIAL PARCEL PROPERTIES

%     environment at initial height (NEW)
      z = 0;		% initial height of parcel (m) 
      [Te pe] = ?	% T (K) and p (Pa) in environment at initial height
      RHe1 = ?		% relative humidity in environment at all heights
      
%     parcel properties at initial height
      thC = ?     	% potential temp (C)
      qv = ?		% mixing ratio (kg/kg) 
      qc = 0	    % liquid water mixing ratio (kg/kg)
      p = pe      	% pressure (Pa)
      th = ?		% potential temperature (K)
      T = ?  		% temperature (K)

%     for entraining parcel  (NEW)
      dz = 100;
      ztop = 10.e3;
      
%     set up column arrays for storing results
      max_size = 1000;
      
      p_mb = zeros(max_size,1)*NaN ;
      th_K = zeros(max_size,1)*NaN;
      T_C = zeros(max_size,1)*NaN;
      qv_gkg = zeros(max_size,1)*NaN;
      qc_gkg = zeros(max_size,1)*NaN;      

%     (NEW)   
      qvs_gkg = zeros(max_size,1)*NaN; 
      Te_C = zeros(max_size,1)*NaN;
      RHe = zeros(max_size,1)*NaN;

%     store initial values (assuming no adjustment req'd)
      i = 1;
      p_mb(i) = p/pa_per_mb;    %  (NEW)
      th_K(i) = th;
      T_C(i) = T - TICE;
      qv_gkg(i) = qv * 1.e3;
      qc_gkg(i) = qc * 1.e3;

%     (NEW)      
      qvs_gkg(i) = ?  	% parcel saturation mixing ratio (g/kg)      
      Te_C(i) = ?		% temperature in environment (C)
      RHe(i) = ?		% relative humidity in environment

%     display initial values
      disp ( ['    p (mb)' '  theta(K)' '      T(C)'  ' q_v(g/kg)' '  q_c(g/kg)'])
      disp ( [p_mb(i) th_K(1) T_C(i) qv_gkg(i) qc_gkg(i)] )

%     DO ASCENT

%       while (p > ptop)
%                     
%       i = i + 1;  
%       p = p - dp;

%     for entraining parcel 

      while (z < ztop)	% (NEW)

      i = i + 1;  
      z = z + dz;		% (NEW)

%     environment T, p (US Standard Atmosphere)	% (NEW)

      [Te pe] = ?	% T (K) and p (Pa) in environment at height z
      
      p = pe;		% parcel pressure = environment pressure

%     environment theta, qv, qc		% (NEW)

      the = ?		% (K)
      qve = ?		% (kg/kg)
      qce = 0;		

%     FOR DRY ADIABATIC ASCENT TH, QV, QC DO NOT CHANGE

%     entrainment 

      th = ?
      qv = ?
      qc = ?         

%     ADJUST PERFORMS ISOBARIC SATURATION ADJUSTMENT

      [th qv qc qvs pi] = satadjust(th,qv,qc,p);

      p_mb(i) = p/pa_per_mb;
      th_K(1) = th;
      T_C(i) = th*pi-TICE;
      qv_gkg(i) = qv * 1.e3;
      qc_gkg(i) = qc * 1.e3;
      
%     (NEW)      
      qvs_gkg(i) = qvs * 1.e3;	% parcel saturation mixing ratio (g/kg) 
      
      Te_C(i) = ?				% temperature in environment (C)
      RHe(i) = ?				% relative humidity in environment
     
%     display  values       
      disp ( [p_mb(i) th_K(1) T_C(i) qv_gkg(i) qc_gkg(i) z/1.e3] )

      end
      
%     disp ( [p_mb(i) th_K(1) T_C(i) qv_gkg(i) qc_gkg(i)] )

%     (NEW)
      tskew_m(p_mb,T_C,qv_gkg./qvs_gkg)   	% parcel
      hold on
      tskew(p_mb,Te_C,RHe)	       		% environment
      axis([-20 20 250 1050])   		% zoom in
       
orient landscape
print -dpsc2 'Std_Atm_skewt.ps'
      