display(['pwv_sonde_07'])
% DEW 5/12/2007

clear all;
close all;


load f:/AMMA_2007/RHB/balloon/raw/AMMA07_sounding_all.mat

clear lat* lon* spd_int dir_int

jd=datenum(year,month,day,hour,minutes,0)-datenum(2006,12,31);

clear year month day hour minute

h=h_int;clear h_int
press=press_int;clear press_int;
t=t_int;clear t_int
rh=rh_int;clear rh_int



e=6.112.*exp(17.67.*t./(t+243.5)).*rh./100; %vapor pressure pa
w=e.*622./(press-0.378*e); %mixing ratio g/kg
es=6.112.*exp(17.67.*t./(t+243.5));%pa


%now, adjust for vapor pressure over ice to determine where get ice clouds.
ice=find(t<=0);
es(ice)=(1.0003+4.18e-6*press(ice)).*6.1115.*exp(22.452.*t(ice)./(t(ice)+272.55));


wm=0.5*(w(1:end-1,:)+w(2:end,:));%take each layer and avg mixing ratio

wm=0.001*wm;%convert from g/kg to g/g to be dimensionless

ql=wm./(wm+1);% change to specific humidity

dpp=1000.*(press(1:end-1,:)-press(2:end,:));
%find pressure difference for each layer in dynes/cm^2, starting from mb

g=980.616;%accln due to gravity cm/s2

pw=(ql./g).*dpp; %% Precipitable water cm

pw=nansum(pw); %cm

RHcorr=rh;
RHcorr(ice)=e(ice)./es(ice).*100;

figure()
set(gcf,'defaultaxesfonts',18)
set(gca,'fontname','time','fontsize',18);
plot(jd,pw,'r')
hold
plot(jd,pw,'r*')
ylim([0 8])
xlabel('Julian Day (UTC)')
ylabel('Sonde PWV (cm)')
title('AMMA 2007')

print_buffer=['-djpeg90 f:/AMMA_2007/RHB/balloon/Processed_images/AMMA07_PWV_sonde.jpg']
eval(['print ' print_buffer])

save f:/AMMA_2007/RHB/balloon/raw/AMMA_PWV_sonde jd pw t rh h