% Read and plot EDT file (STRATUS 2007, SPdeS)

% preamble
read_parameters;
thermo_constants
quickplot=true; % 4-panel T,Td, RH, U, V plots for raw images
plotit=true;
prtit=false;
%clear t;

% construct title string, e.g. 'STRATUS 2007'
tit = [cruise ' ' yyyy]; 
% locate raw sounding files
aa = dir([way_raw_data_balloon '*_LATLON.txt']);
nsonde=length(aa);
% select subset based on yymmddhhMM dates
fst=str2mat(aa.name);
start=0810180000; % the date string as a number for identifying times
endt =0811040000; % include all of VOCALS 2008 % sonde 35 is farthest south
i1=max([1;find(str2num(fst(:,1:10))>=start,1,'first')]);
i2=min([find(str2num(fst(:,1:10))<=endt,1,'last');length(aa)]);

% preallocate arrays for interpolation
pint=(1030:-1:10)'; % pressure coordinate for interp. (hPa)
hint=NaN+zeros(length(pint),nsonde);
tint=hint;
tdint=hint;
thetaint=hint;
qint=hint;
rhint=hint;
uint=hint;
vint=hint;
wsint=hint;
wdint=hint;
latint=hint;
lonint=hint;

hz=(0:10:25000)';   % height coordinate (m)
pz=NaN+zeros(length(hz),length(aa));
tz=pz; tdz=pz; qz=pz; rhz=pz; uz=pz; vz=pz; latz=pz; lonz=pz; thetaz=pz;
wsz=pz; wdz=pz;

ppop=NaN+zeros(nsonde,1);
%tinv=ppop;
tinvtop=ppop;
tinvbase=ppop;
%iinvp=ppop;
iinvptop=ppop;
iinvpbase=ppop;
yr=ppop;
month=ppop;
day=ppop;
hour=ppop;
minute=ppop;
% convert filename timestamps to time variables
for i=i1:i2;
    balfn=aa(i).name; % simple choice of file for now
    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

% activate this to add sounding to data.mat file
if true % ~exist([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2p.mat'],'file')
    % LOOP over sounding files
    for i=i1:i2;
        balfn=aa(i).name;
        if exist([way_raw_data_balloon(1:end-4) 'mat\' balfn(1:end-3) 'mat'],'file') % read a .mat binary file if possible
            disp(['File to read: ' balfn(1:end-3) 'mat'])
            load([way_raw_data_balloon(1:end-4) 'mat\' balfn(1:end-3) 'mat'])
        else % read .txt with a function
            disp(['File to read: ' balfn])
            [h,press,t,td,rh,w_spd,w_dir,lat,lon]=read_edt([way_raw_data_balloon balfn]);
            u=-w_spd.*sin(pi/180*w_dir);
            v=-w_spd.*cos(pi/180*w_dir);
            save([way_raw_data_balloon(1:end-4) 'mat\' balfn(1:end-3) 'mat'],'h','press','t','td','rh','w_spd','w_dir','lat','lon','u','v')
        end
        % problem: DigiCORA reported 0 when wind speed is MISSING, and also when
        % wind speed is actually small. --perhaps resolved by Dan Wolfe
        % after Stratus 2007.

        if quickplot      % quick plots of t,td, rh, u, v
            if (~exist([way_raw_images_balloon balfn(1:21) '.eps'],'file') || ...
                ~exist([way_raw_images_balloon balfn(1:21) '.jpg'],'file') );
                figure;%(1); clf;
                subplot(2,2,1,'align')
                plot(t,h*1e-3,'k','linewidth',2)
                hold on
                plot(td,h*1e-3,'k','linewidth',2)
                set(gca,'ylim',[0 5])
                ylabel('height (km)'); xlabel('temperature, dew point (C)')
                ht=title([cruise ' ' datestr(datenum(yr(i),month(i),day(i)+hour(i)/24+minute(i)/(24*60)),'yyyy-mmm-dd HH:MM UTC')]);
                subplot(2,2,2,'align')
                plot(rh,h*1e-3,'k','linewidth',2)
                set(gca,'ylim',[0 5])
                ylabel('height (km)'); xlabel('relative humidity (%)')
                title(sprintf('latitude: %5.2f, longitude: %5.2f', [lat(find(isfinite(lat),1)) lon(find(isfinite(lon),1))]))
                subplot(2,2,3,'align')
                plot(u,h*1e-3,'k','linewidth',2)
                set(gca,'ylim',[0 5])
                ylabel('height (km)'); xlabel('zonal velocity (m s^{-1})')
                subplot(2,2,4,'align')
                plot(v,h*1e-3,'k','linewidth',2)
                set(gca,'ylim',[0 5])
                ylabel('height (km)'); xlabel('meridional velocity (m s^{-1})')
                orient tall
                print('-deps2',[way_raw_images_balloon balfn(1:21) '.eps']);
                print('-djpeg',[way_raw_images_balloon balfn(1:21) '.jpg']);
            end
        end

        ppop(i)=min(press);
        q=qair3([press t rh]);

        %interpolate to constant & uniform vertical pressure levels
        % to save time, only do to append new data
        %if i>size(hint,2)
        %pint=1030:-1:10;
        [tmp,irev]=unique(press);
        ii=irev(end:-1:1);
        hint(:,i)=interp1(press(ii),h(ii),pint);
        tint(:,i)=interp1(press(ii),t(ii),pint);
        tdint(:,i)=interp1(press(ii),td(ii),pint);
        rhint(:,i)=interp1(press(ii),rh(ii),pint);
        qint(:,i)=interp1(press(ii),q(ii),pint);
        uint(:,i)=interp1(press(ii),u(ii),pint);
        vint(:,i)=interp1(press(ii),v(ii),pint);
        latint(:,i)=interp1(press(ii),lat(ii),pint);
        lonint(:,i)=interp1(press(ii),lon(ii),pint);
        %end
        %interpolate to constant/uniform vertical height levels
        % hz=0:10:25000;
        [tmp,ii]=unique(h);
        pz(:,i)=interp1(h(ii),press(ii),hz);
        tz(:,i)=interp1(h(ii),t(ii),hz);
        tdz(:,i)=interp1(h(ii),td(ii),hz);
        rhz(:,i)=interp1(h(ii),rh(ii),hz);
        qz(:,i)=interp1(h(ii),q(ii),hz);
        uz(:,i)=interp1(h(ii),u(ii),hz);
        vz(:,i)=interp1(h(ii),v(ii),hz);
        wsz(:,i)=interp1(h(ii),w_spd(ii),hz);
        latz(:,i)=interp1(h(ii),lat(ii),hz);
        lonz(:,i)=interp1(h(ii),lon(ii),hz);
    end
    wdz=angle(vz+i*uz)*180/pi+180;
    wdint=angle(vint+i*uint)*180/pi+180;
    thetaint=(273.15+tint).*repmat(pint/1000,[1 length(aa)]).^(-Rd/Cp);
    thetaz=(273.15+tz).*(pz/1000).^(-Rd/Cp);

    save([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2p.mat'], 'pint','hint','tint','tdint','rhint','uint','vint','latint','lonint','qint','thetaint','wdint','wsint','yr','month','day','hour','minute');
    save([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2z.mat'], 'hz','pz','tz','tdz','qz','rhz','uz','vz','latz','lonz','thetaz','wdz','wsz','yr','month','day','hour','minute');
else
    load([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2p.mat'])
    load([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2z.mat'])
end

%perturbation geopotential
zp=hint-repmat(nanmean(hint,2),[1 length(aa)]);

% inversion base and top temperature, pressure, and height
mask=repmat(pint(1:end-1)>800 & pint(1:end-1)<990,[1 size(tint,2)]);
[dTdp,iip]=max(mask.*diff(tint)); % strongest gradient
for ni=1:size(tint,2)
    mask=(pint>=pint(iip(ni))); % true below
    [tinvbase(ni),iinvpbase(ni)]=min(1e12*~mask+tint(:,ni));
    hinvbase(ni)=hint(iinvpbase(ni),ni); % cloud top (inversion base) height
    [tinvtop(ni) ,iinvptop(ni)]=max(-1e12*mask+tint(:,ni));
    hinvtop(ni)=hint(iinvptop(ni),ni); % inversion top height
end
pinvbase=pint(iinvpbase); % cloud top (inversion base) pressure
pinvtop=pint(iinvptop);

la=nanmean(latint); lo=nanmean(lonint);
yday=datenum(yr,month,day,hour,minute,zeros(size(minute)))-datenum(yr,0,0,0,0,0);
octday=datenum(yr,month,day,hour,minute,zeros(size(minute)))-datenum(yr,10,0,0,0,0);
x=repmat(la,[3 1]); y=repmat([500;550;NaN],[1 length(la)]); % for time ticks of spatial plots

% map sonde locations
if plotit
    map_sondes(lo,la)
    m_text(-83,3,'20','verticalalignment','middle','horizontalalignment','center');
    m_text(-87,-20,'25','verticalalignment','middle','horizontalalignment','center')
    m_text(-85,-22,'28','verticalalignment','middle','horizontalalignment','center');
    m_text(-75,-22,'30','verticalalignment','middle','horizontalalignment','center')
    if prtit; print('-depsc2',[way_proc_images_balloon 'map_sonde_VOCALS.eps']); end
end

% load ceilo cloud height to overplot
load([way_proc_data_ceilo 'Sc_base.mat']);
% interpolate lat,lon from sonde to ceilometer by time
lon_cb=interp1(yday(isfinite(yday)),lo(isfinite(yday)),time_cb);
lat_cb=interp1(yday(isfinite(yday)),la(isfinite(yday)),time_cb);

if plotit % plot latitude-pressure (vertical) section
    clf
    b2rcolormap(24);
    subplot(4,1,1,'align')
    [c,h]=contourf(la(i1:i2),pint,thetaint(:,i1:i2),286:2:320);
    set(h,'edgecolor','none')
    set(gca,'ylim',[500 1020])
    set(gca,'ydir','reverse')
    colorbar
    hold on
    plot(x,y,'w')
    for j=2:4:length(octday); text(x(1,j),530,num2str(round(octday(j))),'color','w','horizontalalignment','right'); end
    title({[cruise ' ' yyyy ' southward section'] 'potential temperature (K)'})
    subplot(4,1,2,'align')
    [c,h]=contourf(la(i1:i2),pint,qint(:,i1:i2),0:1:13);
    set(h,'edgecolor','none')
    set(gca,'ylim',[500 1020])
    set(gca,'ydir','reverse')
    colorbar
    hold on
    plot(x,y,'w')
    for j=2:4:length(yday); text(x(1,j),530,num2str(round(octday(j))),'color','w','horizontalalignment','right'); end
    title('specific humidity (g/kg)');
    subplot(4,1,3,'align')
    [c,h]=contourf(la(i1:i2),pint,uint(:,i1:i2),-10:10);
    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('zonal velocity (m s ^{-1})')
    ylabel('pressure (hPa)')
    subplot(4,1,4,'align')
    [c,h]=contourf(la(i1:i2),pint,vint(:,i1:i2),-10:10);
    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(octday); text(x(1,j),530,num2str(round(octday(j))),'color','k','horizontalalignment','right'); end
    title('meridional velocity (m s^{-1})')
    xlabel('latitude (\circ)')
    orient tall
    if prtit
        print('-depsc2',[way_proc_images_balloon 'latpres_thetaquv.eps'])
        print('-djpeg',[way_proc_images_balloon 'latpres_thetaquv.jpg'])
    end
end

% plot 85W latitude-height section
start85=0810180000;
endt85 =0810250900;
i1=max([1;find(str2num(fst(:,1:10))>=start85,1,'first')]);
i2=min([find(str2num(fst(:,1:10))<=endt85,1,'last');length(aa)]);

x=repmat(la(i1:i2),[3 1]); y=repmat([4.8;5.0;NaN],[1 i2-i1+1]);
y(1,hour(i1:i2)==23 | hour(i1:i2)==0)=4.6;
j=(hour(i1:i2)==11 | hour(i1:i2)==12); % midday for oct day labels
iz=hz<=5e3;
if plotit
    clf
    b2rcolormap(24);
    subplot(4,1,1,'align')
    [c,h]=contourf(la(i1:i2),hz(iz)*1e-3,thetaz(iz,i1:i2),286:2:320);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 5])
    colorbar
    hold on
    plot(lat_cb(1:end-1),cb2_85/1e3,'y.','markersize',2.5)
    plot(la(i1:i2),hinvbase(i1:i2)/1e3,'m.')
    plot(la(i1:i2),hinvtop(i1:i2)/1e3,'m.')
    plot(x,y,'w')
    text(x(1,j),4.8+zeros(1,sum(j)),num2str(floor(octday(j))),'color','w','horizontalalignment','right');
    title({[cruise ' ' yyyy ' southward section'] 'potential temperature (K)'})
    subplot(4,1,2,'align')
    [c,h]=contourf(la(i1:i2),hz(iz)*1e-3,qz(iz,i1:i2),0:1:15);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 5])
    colorbar
    hold on
    plot(x,y,'w')
    text(x(1,j),4.8+zeros(1,sum(j)),num2str(floor(octday(j))),'color','w','horizontalalignment','right');
    title('specific humidity (g/kg)');
    subplot(4,1,3,'align')
    [c,h]=contourf(la(i1:i2),hz(iz)*1e-3,uz(iz,i1:i2),-12:12);
    set(h,'edgecolor','none')
    colorbar
    set(gca,'ylim',[0 5]);
    colorbar
    hold on
    plot(x,y,'k')
    text(x(1,j),4.8+zeros(1,sum(j)),num2str(floor(octday(j))),'color','k','horizontalalignment','right');
    title('zonal velocity (m s ^{-1})')
    ylabel('height (km)')
    subplot(4,1,4,'align')
    [c,h]=contourf(la(i1:i2),hz(iz)*1e-3,vz(iz,i1:i2),-12:12);
    set(h,'edgecolor','none')
    colorbar
    set(gca,'ylim',[0 5])
    colorbar
    hold on
    plot(x,y,'k')
    text(x(1,j),4.8+zeros(1,sum(j)),num2str(floor(octday(j))),'color','k','horizontalalignment','right');
    title('meridional velocity (m s^{-1})')
    xlabel('latitude (\circ)')
    orient tall
    if prtit
        print('-depsc2',[way_proc_images_balloon 'latheight_85W_thetaquv.eps'])
        print('-djpeg',[way_proc_images_balloon 'latheight_85W_thetaquv.jpg'])
    end
end

% plot time-height series for the whole cruise
if false
    clf
    b2rcolormap(24);
    subplot(4,1,1,'align')
    [c,h]=contourf(yday(i1:i2),hz/1e3,thetaz(:,i1:i2),286:2:320);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3])
    colorbar
    hold on
    plot(time_cb(1:end-1),cb2_85/1e3,'y.','markersize',2) % cloud base from ceilo
    title({[cruise ' ' yyyy ' time-height section'] 'potential temperature (K)'})
    subplot(4,1,2,'align')
    [c,h]=contourf(yday(i1:i2),hz/1e3,qz(:,i1:i2),0:1:13);
    set(h,'edgecolor','none')
    set(gca,'ylim',[0 3])
    colorbar
    hold on
    title('specific humidity (g/kg)');
    subplot(4,1,3,'align')
    [c,h]=contourf(yday(i1:i2),hz/1e3,uz(:,i1:i2),-10:10);
    set(h,'edgecolor','none')
    colorbar
    set(gca,'ylim',[0 3]);
    colorbar
    hold on
    title('zonal velocity (m s ^{-1})')
    ylabel('height (km)')
    subplot(4,1,4,'align')
    [c,h]=contourf(yday(i1:i2),hz/1e3,vz(:,i1:i2),-10:10);
    set(h,'edgecolor','none')
    colorbar
    set(gca,'ylim',[0 3])
    colorbar
    hold on
    title('meridional velocity (m s^{-1})')
    xlabel(sprintf('year day (UTC, Oct %i-%i)',day([i1 i2])))
    orient tall
    %print('-depsc2',[way_raw_images_balloon 'timez_thetaquv.eps'])
    %print('-djpeg',[way_raw_images_balloon 'timez_thetaquv.jpg'])
end

% plot time-pressure series for the whole cruise
if false
    clf
    b2rcolormap(24);
    subplot(4,1,1,'align')
    [c,h]=contourf(yday(i1:i2),pint,thetaint(:,i1:i2),286:2:320);
    set(h,'edgecolor','none')
    set(gca,'ylim',[500 1020])
    set(gca,'ydir','reverse')
    colorbar
    hold on
    title({'STRATUS 2007 time-pressure section' 'potential temperature (K)'})
    subplot(4,1,2,'align')
    [c,h]=contourf(yday(i1:i2),pint,qint(:,i1:i2),0:1:13);
    set(h,'edgecolor','none')
    set(gca,'ylim',[500 1020])
    set(gca,'ydir','reverse')
    colorbar
    hold on
    title('specific humidity (g/kg)');
    subplot(4,1,3,'align')
    [c,h]=contourf(yday(i1:i2),pint,uint(:,i1:i2),-10:10);
    set(h,'edgecolor','none')
    set(gca,'ydir','reverse')
    colorbar
    set(gca,'ylim',[500 1020]);
    set(gca,'ydir','reverse')
    colorbar
    hold on
    title('zonal velocity (m s ^{-1})')
    ylabel('pressure (hPa)')
    subplot(4,1,4,'align')
    [c,h]=contourf(yday(i1:i2),pint,vint(:,i1:i2),-10:10);
    set(h,'edgecolor','none')
    set(gca,'ydir','reverse')
    colorbar
    set(gca,'ylim',[500 1020])
    set(gca,'ydir','reverse')
    colorbar
    hold on
    title('meridional velocity (m s^{-1})')
    xlabel('year day (UTC, Oct 18-Nov 3)')
    orient tall
    % not to end of cruise on scarecrow or LT
    %print('-depsc2',[way_raw_images_balloon 'timepres_thetaquv.eps'])
    %print('-djpeg',[way_raw_images_balloon 'timepres_thetaquv.jpg'])
end

%%% end of commonly-used processing, other diagnostics and plots follow

if false % 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
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]);

% 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;

if false
    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
end

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

if 0
    % 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 (%)'})
end