% plots and diagnostics I don't usually do
plot_sonde_timeheight
plot_sonde_timepres

% plot to check the inversion heights
plot(pinv); set(gca,'ydir','reverse')
for ni=4:5:nsonde
    figure; plot(tint(:,ni:ni+4),pint,'b');
    ylim([700 1010]); set(gca,'ydir','reverse')
    hold on
    plot(tinv(ni:ni+4),pinv(ni:ni+4),'rx')
end

% cumulative integral of water vapor
dp=100.0; % 1 hPa=100 Pa
rhow=1.0e3;
qv=qint*1.0e-3;
qv(isnan(qv))=0; % not to worry in a cumsum
wvpcumint=dp/rhow/g*cumsum(qv); % water vapor path in precipitable meters
normwvpcumint=wvpcumint./repmat(wvpcumint(end,:),[length(wvpcumint) 1]);

%plot time series of WVP
mailbox=load([way_raw_data_mailbox 'mailbox.mat']);
area(octday,1e2*wvpcumint(end,:),'linewidth',1.5,'facecolor',[.7 .7 1])
hold on
plot(mailbox.time+octday(1)-yday(1),mailbox.wvp,'.','markersize');
axis([20 32 0 7])
title('VOCALS 2008 radiosonde water vapor path')
ylabel('cm')
xlabel('2008 October day (UT)')
set(gca,'tickdir','out')
set(gca,'xtick',20:32,'xticklabel',num2str([20:1:31 1]'))
if prtit; print('-depsc2',[way_proc_images_balloon 'WVP.eps']); end

% composite on the fraction of cumulative WVP at 800 hPa
fi800=find(pint==800);
iwet=normwvpcumint(fi800,:)<0.64;
imed=normwvpcumint(fi800,:)>=0.64 & normwvpcumint(fi800,:)<=0.79;
idry=normwvpcumint(fi800,:)>0.79;

subplot(1,2,1,'align')
plot(wvpcumint(:,iwet),pint,'b')
hold on
plot(wvpcumint(:,imed),pint,'g')
plot(wvpcumint(:,idry),pint,'r')
set(gca,'ydir','reverse')
axis tight
ylabel('pressure')
xlabel('cumulative water vapor path (m liq. equiv.)')

subplot(1,2,2)
plot(normwvpcumint(:,iwet),pint,'b')
hold on
plot(normwvpcumint(:,imed),pint,'g')
plot(normwvpcumint(:,idry),pint,'r')
set(gca,'ydir','reverse','fontsize',14)
axis tight
ylabel('pressure')
xlabel('normalized cumulative water vapor path')
subplot(1,2,2)
plot(lo(iwet),la(iwet),'bo')
hold on
plot(lo(imed),la(imed),'g.')
plot(lo(idry),la(idry),'rx')
axis equal; % axis tight
set(gca,'fontsize',14)
xlabel('longitude')
ylabel('latitude')
if prtit;
    print('-depsc2',[way_raw_images_balloon 'composite800wvp.eps'])
end

% test the qrat approach
plot(qv(31,:)*1000,wvpcumint(end,:)*100,'.')
xlabel('q (g/kg)'); ylabel('WVP (cm liq. equiv.)');
% qrat not so good--better to use calculated qv for radiative transfer

% perturbation geopotential
[c,h]=contourf(yday(i1:i2)+hour(i1:i2)/24,pint,zp(:,i1:i2), -30:3:30);
set(h,'edgecolor','none')
set(gca,'ydir','reverse')
set(gca,'ylim',[200 1020])
colorbar
% need to adjust pressure obs for each sonde by ship barometer.

% relative humidity
[c,h]=contourf(la,pint,rhint);
set(h,'edgecolor','none')
set(gca,'ydir','reverse')
colorbar
set(gca,'ylim',[500 1020])
set(gca,'ydir','reverse')
colorbar
hold on
plot(x,y,'k')
for j=2:4:length(yday); text(x(1,j),530,num2str(round(octday(j))),'color','k','horizontalalignment','right'); end
title({'relative humidity (%)'})