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

% preamble
read_parameters_stratus;
thermo_constants
quickplot=false;
%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;
rhint=hint;
uint=hint;
vint=hint;
latint=hint;
lonint=hint;

% 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 =0710230600; % 
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('wtec_edt_interp2p.mat','file')
    for i=i1:i2;
        balfn=aa(i).name; % simple choice of file for now
        if exist([pathdata balfn(1:end-3) 'mat'],'file') % read a .mat binary file
            disp(['File to read: ' balfn(1:end-3) 'mat'])
            load([pathdata 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 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.^(-Rd/Cp);
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);
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)]);

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'])

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

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