function [th1, qv1, qc1, qvs1, pibar] = satadjust_dummy(th, qv, qc, pbar);
%
%  Performs an isobaric moist adiabatic adjustment.
%  The final state is either subsaturated with no liquid water
%  present, or exactly saturated with liquid water present.
%
%  This version iterates to obtain the adjustment, so accurate
%  guesses for TH and QV are not needed if ITTMAX=10. 
%
%  Units:  SI (MKS)
% 
%  Input -- 
%  TH: potential temperature, theta^* (K)
%  QV: mixing ratio of water vapor, q_v^* (kg/kg)
%  QC: mixing ratio of liquid water, q_c^* (kg/kg)
%  PBAR: pressure, p (Pa)
%
%  Output -- 
%  TH1: adjusted potential temperature, theta^{n+1} (K)
%  QV1: adjusted mixing ratio of water vapor, q_v^{n+1} (kg/kg)
%  QC1: adjusted mixing ratio of liquid water, q_c^{n+1} (kg/kg)
%  QVS1: saturation mixing ratio, q_vs^{n+1} (kg/kg) for TH1, PBAR
%  PIBAR: Exner function, pi (non-dimensional pressure)
%
% Required functions:
%  esat_Pa(T): saturation vapor pressure (Pa)
%  desdT_Pa(T): derivative of esat(T) (Pa/K)

%  physical constants
hlf = 2.500E+06;
cp = 1004;
rgas = 287;
pzero = 100000;

% iteration parameters
ittmax = 10;        % maximum number of iterations allowed
dT_crit = 0.001;    % stop iterating once T adjustment is less than dT_crit 

%  PBAR is a (hydrostatic) reference pressure

itt = 1;
pibar = (pbar/pzero)^(rgas/cp);
gamma = hlf/(cp*pibar);

thstar = th;
qvstar = qv;
qcstar = qc;

do_it = 1;
while (do_it == 1)

% Tstar = temperature before adjustment
Tstar = thstar*pibar;
% saturation vapor pressure at T=Tstar
es1 = esat_Pa(Tstar); 
% saturation mixing ration at T=Tstar, p=pbar
qvsat = 0.622/(pbar-es1)*es1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assume saturated after adjustment
% solve for adjusted th, qv (water vapor), qc (cloud water):
% th1 = ?
% qv1 = ?
% qc1 = ?
% saturation mixing ratio after adjustment = adjusted qv:
% qvs1 = qv1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TEMPORARY: NO ADJUSTMENT *** REMOVE THIS CODE ***
th1 = thstar;
qv1 = qvstar;
qc1 = qcstar;
qvs1 = qvsat;
% TEMPORARY: NO ADJUSTMENT *** REMOVE THIS CODE ***

% check saturation assumption used above:
% if qc1 < 0, saturation assumption must be incorrect
if (qc1 < 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% not saturated: solve for adjusted qc, qv, th   
   qc1 = 0;     % no cloud water if not saturated
%  qv1 = ?
%  th1 = ? 
% saturation mixing ratio after adjustment > adjusted qv:  
% qvs1= ?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


% T adjustment for this iteration
dT = (th1-thstar)*pibar;

if (abs(dT) < dT_crit) | (itt == ittmax)
  do_it = 0;
end

thstar = th1;
qvstar = qv1;
qcstar = qc1;	% added Sep 2008

itt = itt+1;

end

