%     program parcel
%     9-9-08 -- changed th_K(1) to th_K(i)

%     INITIAL PARCEL PROPERTIES
      pMB = 900. % pressure (mb)
      TC = 10.    % temperature (deg C)
      qv = 6.e-3 % water vapor mixing ratio
      qc = 0.	  % liquid water mixing ratio

%     ASCENT PARAMETERS   
      dpMB = 10.     % pressure interval (mb)
      ptopMB = 700.  % ending pressure (mb)

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

      p = pMB * pa_per_mb;
      T = TC + TICE;
      dp = dpMB  * pa_per_mb;
      ptop = ptopMB * pa_per_mb;

%     INPUT FOR ADJUST: TH, QV, QC
      th = theta ( T, p );  % theta(T,p) is a function

%     set up column arrays for storing results

      max_size = 1000;
      
      p_mb = zeros(max_size,1);
      th_K = zeros(max_size,1);
      T_C = zeros(max_size,1);
      qv_gkg = zeros(max_size,1);
      qc_gkg = zeros(max_size,1);

%     store initial values (assuming no adjustment req'd)
      i = 1;
      p_mb(i) = pMB;
      th_K(i) = th;
      T_C(i) = TC;
      qv_gkg(i) = qv * 1.e3;
      qc_gkg(i) = qc * 1.e3;
      
%     display initial values
      disp ( ['    p (mb)' '  theta(K)' '      T(C)'  ' q_v(g/kg)' '  q_c(g/kg)'])
      disp ( [p_mb(i) th_K(i) T_C(i) qv_gkg(i) qc_gkg(i)] )

%     DO ASCENT

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

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

%     ADJUST PERFORMS ISOBARIC SATURATION ADJUSTMENT

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

      p_mb(i) = p/pa_per_mb;
      th_K(i) = th;
      T_C(i) = th*pi-TICE;
      qv_gkg(i) = qv * 1.e3;
      qc_gkg(i) = qc * 1.e3;
%     display  values       
      disp ( [p_mb(i) th_K(i) T_C(i) qv_gkg(i) qc_gkg(i)] )

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

