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

% preamble
read_parameters_stratus;
thermo_constants
quickplot=true;
%clear t;
pathdata=way_raw_data_balloon;
rawimgpath=[pathdata(1:28) 'Raw_images\'];

% locate (new) raw sounding files
aa = dir([pathdata '*_LATLON.txt']);
pint=1030:-1:10; % pressure coordinate for interp. (hPa)
hz=0:10:25000;   % height coordinate (m)
hint=NaN+zeros(length(pint),length(aa));
tint=hint;
tdint=hint;
thetaint=hint;
qint=hint;
rhint=hint;
uint=hint;
vint=hint;
latint=hint;
lonint=hint;

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;

% construct title string, e.g. 'STRATUS 2007'
tit = cruise;
lt = length(cruise);
tit(lt+1) = ' ';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file

% LOOP HERE to analyze many soundings
fst=str2mat(aa.name);
start=0710181500; % the date string as a number for identifying times
endt =0711082559; % includ all of Stratus 2007
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=i1:i2;
    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

if ~exist([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2p.mat'],'file') || true % always do this until the end of the cruise
    for i=i1:i2;
        balfn=aa(i).name; % simple choice of file for now
        if exist([pathdata(1:end-4) 'mat\' balfn(1:end-3) 'mat'],'file') % read a .mat binary file
            disp(['File to read: ' balfn(1:end-3) 'mat'])
            load([pathdata(1:end-4) 'mat\' balfn(1:end-3) 'mat'])
        else % read with a not-so-fancy function
            disp(['File to read: ' balfn])
            [h,press,t,td,rh,w_spd,w_dir,lat,lon]=read_edt([pathdata balfn]);
            u=-w_spd.*sin(pi/180*w_dir);
            v=-w_spd.*cos(pi/180*w_dir);
            save([pathdata(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 reports 0 when wind speed is MISSING, and also when
        % wind speed is actually small.

        if quickplot      % quick plots of t,td, rh, u, v
            if (~exist([rawimgpath balfn(1:21) '.eps'],'file') || ...
                ~exist([rawimgpath 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',[rawimgpath balfn(1:21) '.eps']);
                print('-djpeg',[rawimgpath balfn(1:21) '.jpg']);
            end

            if 0
                subplot(1,2,1,'align')
                plot(u,h*1e-3,'k','linewidth',2)
                hold on
                set(gca,'ylim',[0 5])
                ylabel('height (km)'); xlabel('zonal velocity (m s^{-1})')
                subplot(1,2,2,'align')
                plot(v,h*1e-3,'k','linewidth',2)
                hold on
                set(gca,'ylim',[0 5])
                ylabel('height (km)'); xlabel('meridional velocity (m s^{-1})')
                orient landscape
            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);
        latz(:,i)=interp1(h(ii),lat(ii),hz);
        lonz(:,i)=interp1(h(ii),lon(ii),hz);
    end
    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');
    save([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2z.mat'], 'hz','pz','tz','tdz','qz','rhz','uz','vz','latz','lonz','thetaz');
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

% plot lat-height of interpolated values
la=nanmean(latint); lo=nanmean(lonint);
yday=datenum(str2num(year),month,day,hour,minute,zeros(size(minute)))-datenum('2007-0-0');
x=repmat(la,[3 1]); y=repmat([500;550;NaN],[1 length(la)]);

if 0 % plot latitude-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(yday); text(x(1,j),530,num2str(round(yday(j))),'color','w','horizontalalignment','right'); end
    title({'STRATUS 2007 (Oct 18-23)' 'Ecuador-Peru coastal 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(yday(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(yday(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])
    shading flat
    set(gca,'ydir','reverse')
    colorbar
    hold on
    plot(x,y,'k')
    for j=2:4:length(yday); text(x(1,j),530,num2str(round(yday(j))),'color','k','horizontalalignment','right'); end
    title('meridional velocity (m s^{-1})')
    xlabel('latitude (\circ)')
    orient tall
    %print('-depsc2',[rawimgpath 'latpres_thetaquv.eps'])
    %print('-djpeg',[rawimgpath 'latpres_thetaquv.jpg'])
end;

% plot time-height series
if 0
    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])
    shading flat
    set(gca,'ydir','reverse')
    colorbar
    hold on
    title('meridional velocity (m s^{-1})')
    xlabel('year day (UTC, Oct 18-31)')
    orient tall
    %print('-depsc2',[rawimgpath 'timepres_thetaquv.eps'])
    %print('-djpeg',[rawimgpath 'timepres_thetaquv.jpg'])
end

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

% inversion base pressure
mask=repmat(pint(1:end-1)'>800 & pint(1:end-1)'<990,[1 size(tint,2)]);
[dTdp,iip]=max(mask.*diff(tint));
for ni=1:size(tint,2)
    mask=(pint>=pint(iip(ni)))';
    [tinv(ni),iinvp(ni)]=min(1e12*~mask+tint(:,ni));
end
pinv=pint(iinvp);

if 0 % check the inversion heights
    plot(pinv); set(gca,'ydir','reverse')
    for ni=1:6:66
        figure; plot(tint(:,ni:ni+5),pint,'b');
        ylim([700 1010]); set(gca,'ydir','reverse')
        hold on
        plot(tinv(ni:ni+5),pinv(ni:ni+5),'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 0
    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,1)
    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')
    %print('-depsc2',[rawimgpath 'composite800wvp.eps'])
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(day(i1:i2)+hour(i1:i2)/24,pint,zp(:,i1:i2), -20:2:20);
    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(yday(j))),'color','k','horizontalalignment','right'); end
    title({'relative humidity (%)'})
end