% sounding_plot_VOCALS2008
% VOCALS 2008 :: 2008-10-31 :: Simon de Szoeke
% Read and plot EDT LATLON files. Interpolate to regular pressure and height
% coordinates. Save data.

% preamble
%read_parameters;
%thermo_constants
Rd=
quickplot=true; % 4-panel T,Td, RH, U, V plots for raw images
plotit=true;
prtit=true; %Keeps from saving the file image
warning('off','MATLAB:interp1:NaNinY')
%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;
hinvtop=ppop;
hinvbase=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);
novday=datenum(yr,month,day,hour,minute,zeros(size(minute)))-datenum(yr,11,0,0,0,0);
x=repmat(la,[3 1]); y=repmat([500;550;NaN],[1 length(la)]); % for time ticks of spatial plots

% 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
    % map sonde locations
    map_sondes(lo,la)
    if prtit; print('-depsc2',[way_proc_images_balloon 'map_sonde_VOCALS.eps']); 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)]);
    plot_sonde_latheight
    if prtit
        print('-depsc2',[way_proc_images_balloon 'latheight_85W_thetaquv.eps'])
        print('-djpeg',[way_proc_images_balloon 'latheight_85W_thetaquv.jpg'])
    end
 
    % plot 20S longitude-height section
    start20=0810271900;
    endt20 =0811040000;
    i1=max([1;find(str2num(fst(:,1:10))>=start20,1,'first')]);
    i2=min([find(str2num(fst(:,1:10))<=endt20,1,'last');length(aa)]);
    ii=(1:size(x,2))'>=i1 & (1:size(x,2))'<=i2;
    plot_sonde_lonheight
    if prtit
        print('-depsc2',[way_proc_images_balloon 'latheight_20S_thetaquv.eps'])
        print('-djpeg',[way_proc_images_balloon 'latheight_20S_thetaquv.jpg'])
    end

    % plot latitude-pressure (vertical) section
    plot_sonde_latpres
    if prtit
        print('-depsc2',[way_proc_images_balloon 'latpres_thetaquv.eps'])
        print('-djpeg',[way_proc_images_balloon 'latpres_thetaquv.jpg'])
    end
end

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