% Compute the diurnal cycle of the sondes while on station at the WHOI buoy
% 20S, 85W for Stratus 2007.
%
% SPdeS 2007 November 1

if ~exist('way_raw_data_balloon','var')
    read_parameters_Stratus
end
pathdata=[way_raw_data_balloon(1:end-4) 'mat\'];
%rawimgpath=way_images_balloon;
rawimgpath='E:\STRATUS_2007\RHB\balloon\Raw_Images\';
aa = dir([pathdata '*_LATLON.mat']);
fst=str2mat(aa.name);
start=0710261200; % the date string as a number for identifying times
endt =0710312459;
i1=max([1;find(str2num(fst(:,1:10))>=start,1,'first')]);
i2=min([find(str2num(fst(:,1:10))<=endt,1,'last');length(aa)]);

for i=1:length(aa); % get times from filenames
    balfn=aa(i).name; % simple choice of file for now
    % read timestamp from filename
    yr(i) = str2double(balfn(1:2))+2000;
    month(i) = str2double(balfn(3:4));
    day (i)  = str2double(balfn(5:6));
    hour(i)  = str2double(balfn(7:8));
    minute(i) = str2double(balfn(9:10));
end
yday=datenum(str2num(year),month,day,hour,minute,zeros(size(minute)))-datenum('2007-0-0');
load([pathdata 'wtec_edt_interp2z.mat']);
la=nanmean(latz); lo=nanmean(lonz);

% inversion base height
mask=repmat(hz(1:end-1)'<1800 & hz(1:end-1)'>200,[1 size(tz,2)]);
[dTdz,iiz]=max(mask.*diff(tz));
for ni=1:size(tz,2)
    mask=(hz<=hz(iiz(ni)))';
    [tinv(ni),iinvz(ni)]=min(1e12*~mask+tz(:,ni));
end
zinv=hz(iinvz);

if 0 % test inversion heights
    plot(zinv);
    for ni=1:6:66
        figure; plot(tz(:,ni:ni+5),hz,'b');
        ylim([0 1600]);
        hold on
        plot(tinv(ni:ni+5),zinv(ni:ni+5),'rx')
    end
end

% composite zinv 4-hourly (6x a day) for the 20S,85W period
time=(0:4:20);
zinvcomp=NaN+zeros(6,1);
zi=zinv(i1:i2); % keep only desired period
for it=1:length(time)
    zinvcomp(it)=nanmean(zi(hour(i1:i2)==mod(time(it)-1,24) | hour(i1:i2)==time(it)));
    zinvcmpmat(hour==mod(time(it)-1,24) | hour==time(it))=zinvcomp(it);
end

% normalized inversion height moves with diurnal cycle
hnorm=NaN+zeros(size(thetaz));
thetacomp=NaN+zeros(size(thetaz));
rhcomp=thetacomp;
qcomp=thetacomp;
ucomp=thetacomp;
vcomp=thetacomp;
for it=i1:i2
    % stretch/shift z to match inversion base height (only good for i1:i2)
    hnorm(1:iinvz(it),it)=hz(1:iinvz(it))'*(zinvcmpmat(it)/zinv(it)); % normalize in BL
    hnorm(iinvz(it)+1:end,it)=hz(iinvz(it)+1:end)+zinvcmpmat(it)-zinv(it); % offset above
    % interp transformed variable to standard coordinate (only good for i1:i2)
    thetacomp(:,it)=interp1(hnorm(:,it),thetaz(:,it),hz);
    rhcomp(:,it)=interp1(hnorm(:,it),rhz(:,it),hz);
    qcomp(:,it)=interp1(hnorm(:,it),qz(:,it),hz);
    ucomp(:,it)=interp1(hnorm(:,it),uz(:,it),hz);
    vcomp(:,it)=interp1(hnorm(:,it),vz(:,it),hz);
end

% thetaz(hz) is the data interpolated to standard levels hz
% thetaz(hnorm(:,t)) is transformed s.t. the inversions line up for values of hnorm
% thetacomp(hz) is transformed so that thetacomp inversions line up with hz
%   as the independent variable.

% composite diurnal cycle 4-hourly (6x a day) for the 20S,85W period
thetadiel=NaN+zeros(length(thetaz),6);
qdiel=thetadiel;
rhdiel=qdiel;
udiel=thetadiel;
vdiel=thetadiel;
al=false(1,length(yday));
al(i1:i2)=true;
for it=1:length(time)
    thetadiel(:,it)=nanmean(thetacomp(:,(hour==mod(time(it)-1,24) | hour==time(it)) & al),2);
    rhdiel(:,it)=nanmean(rhcomp(:,(hour==mod(time(it)-1,24) | hour==time(it)) & al),2);
    qdiel(:,it)=nanmean(qcomp(:,(hour==mod(time(it)-1,24) | hour==time(it)) & al),2);
    udiel(:,it)=nanmean(ucomp(:,(hour==mod(time(it)-1,24) | hour==time(it)) & al),2);
    vdiel(:,it)=nanmean(vcomp(:,(hour==mod(time(it)-1,24) | hour==time(it)) & al),2);
end

subs=[6 1:6 1:6 1];
taxis=[-4 time 24+time 48+4];
    clf
    b2rcolormap(24);
    subplot(5,1,1,'align')
    [c,h]=contourf(taxis,hz/1e3,thetadiel(:,subs),286:1:310);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3],'xlim',[-2 46])
    set(gca,'xtick',[0:4:48],'xticklabel',[20 00 04 08 12 16 20 00 04 08 12 16])
    colorbar
    title({'STRATUS 2007 (Oct 26-31) 20\circ S, 85\circ W WHOI buoy location' 'potential temperature (K)'})
    subplot(5,1,2,'align')
    [c,h]=contourf(taxis,hz/1e3,qdiel(:,subs),0:0.5:8.5);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3],'xlim',[-2 46])
    set(gca,'xtick',[0:4:48],'xticklabel',[20 00 04 08 12 16 20 00 04 08 12 16])
    colorbar
    title('specific humidity (g/kg)');
    subplot(5,1,3,'align')
    [c,h]=contourf(taxis,hz/1e3,rhdiel(:,subs),0:5:100);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3],'xlim',[-2 46])
    set(gca,'xtick',[0:4:48],'xticklabel',[20 00 04 08 12 16 20 00 04 08 12 16])
    colorbar
    title('relative humidity (%)');
    ylabel('altitude (km)')
    subplot(5,1,4,'align')
    [c,h]=contourf(taxis,hz/1e3,udiel(:,subs),-8:0.5:2);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3],'xlim',[-2 46])
    set(gca,'xtick',[0:4:48],'xticklabel',[20 00 04 08 12 16 20 00 04 08 12 16])
    colorbar
    title('zonal velocity (m s ^{-1})')
    subplot(5,1,5,'align')
    [c,h]=contourf(taxis,hz/1e3,vdiel(:,subs),-2:0.5:8);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3],'xlim',[-2 46])
    set(gca,'xtick',[0:4:48],'xticklabel',[20 00 04 08 12 16 20 00 04 08 12 16])
    colorbar
    title('meridional velocity (m s^{-1})')
    xlabel('local hour')
    orient tall
    %print('-depsc2',[rawimgpath 'diurnal_sonde_20S85W.eps'])
    %print('-djpeg',[rawimgpath 'diurnal_sonde_20S85W.jpg'])

% bring in ceilo cloud base and profiler cloud top.