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