display(['pwv_sonde'])
% DEW 5/12/2007
% DEW 3/22/2008

clear all;
close all;

%%%%%%%%%%%  Read parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
read_parameters_CalNEX_2010;     %reads the parameters for the cruise.

% % cruise='CALNEX';             % Cruise name  
% % year='2010';               % year (deprecated--make year into a double) 
% % yyyy='2010';               % year (string)
% % ship='Atlantis';               % Research vessel
% % 
% % %  balloon files
% % way_raw_images_balloon='D:\CALNEX_2010\Atlantis\balloon\Raw_Images\';  %daily saved images
% % way_proc_images_balloon='D:\CALNEX_2010\Atlantis\balloon\Processed_Images\';  %processed images
% % way_raw_data_balloon='D:\CALNEX_2010\Atlantis\balloon\Raw\edt\';  % Raw Data 
% % way_proc_data_balloon='D:\CALNEX_2010\Atlantis\balloon\Processed\'; % Processed Data

%ls([way_raw_data_balloon '*LAT*.txt'])

ls([way_proc_data_balloon '*.txt'])

fname=['*.txt'];

s1=[way_proc_data_balloon fname];
s='m:\';
%s1=['c:\data\icealot\data\Mailbox\' fname];
dr=dir(s1);
[n m]=size(dr);
mx=n;

for i=1:mx  %loop through temp profile text files
    clear balt
   flnm=dr(i).name
%c:\data\icealot\data\sonde\VOCALS_sonde_08_03_19_1831.txt
f=[way_proc_data_balloon flnm];
ss = [f];
balt = load(ss);

yr = str2num(flnm(19:20)) + 2000;
month = str2num(flnm(22:23));
day = str2num(flnm(25:26));
hour = str2num(flnm(28:29));
minutes = str2num(flnm(30:31));
jd=datenum(yr,month,day,hour,minutes,0)-datenum(yr-1,12,31);

clear yr month day hour minute
grh = find(balt(:,4) >= 0 & balt(:,4) <= 100); %find good rh values

h=balt(grh,1);
press=balt(grh,2);
t=balt(grh,3);
rh=balt(grh,4);


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

tpw(i,1) = jd;
tpw(i,2) = nansum(pw); %cm

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

end

figure(1)
clf
set(gcf,'defaultaxesfonts',18)
set(gca,'fontname','time','fontsize',18);
plot(tpw(:,1),tpw(:,2),'r')
hold
plot(tpw(:,1),tpw(:,2),'r*')
ylim([0 8])
xmin = fix(tpw(1,1));
xmax = fix(tpw(end,1)+1);
xlim([xmin xmax])
xlabel('Year Day (UTC)')
ylabel('Sonde PWV (cm)')
title('CALNEX 2010')

if 1
 print_buffer = [way_proc_images_balloon, cruise, '_sonde_pwv.png'];
 print('-dpng ', print_buffer);
end

eval(['cd ', way_proc_data_balloon]);
save CALNEX_PWV_sonde tpw