%%%% READ SOUNDING DATA FROM ptu and wnd FORMATs...
%  New edt file format wasn't in place yet for MW15
%%%% COMPOSITE TIME-HEIGHT DISPLAY ....wn
%%% pk
%% DEW 01/27/2005
%% DEW 3/19/2008 ICEALOT
clear

warning off MATLAB:divideByZero  % turn off warning message

pathdata = 'c:\Documents and Settings\krasmussen\My Documents\RICO Project\RICO Raw Data';
outpath  = 'c:\Documents and Settings\krasmussen\My Documents\RICO Project\RICO Sonde\Results';
pathdata = 'c:\data\icealot\data\sonde\ptuwnd';
outpath = 'c:\data\icealot\data\sonde\';


eval(['cd ', pathdata]);

aa = dir('*.wnd');
[m n] = size(aa);
p = length(aa);
[num] = num2str(p);
disp('Number of rawinsonde *.wnd files : '); disp( [num]);

% set up arrays  large enough to hold largest launche
balttmp(1:3600,1:p) = -999;
balwtmp(1:3600,1:p) = -999;
httmp(1:3600,1:p) = -999;
hwtmp(1:3600,1:p) = -999;
Hw(1:3600,1:p) = -999;
spd(1:3600,1:p) = -999;
dr(1:3600,1:p) = -999;
Ht(1:3600,1:p) = -999;
T(1:3600,1:p) = -999;
RH(1:3600,1:p) = -999;
P(1:3600,1:p) = -999;

%%%%%%%%%% starting file %%%%%%%%%%%%%%%%%%%%%%%%
st=1;
%%%%%%%%%%%%%%%%%% last file == k %%%%%%%%%%%%%%%
ws=1;
we = m;
% set number of files to read in for testing
% we = 3;
% Loop through all soundings in directory
for w = ws:we
    w;
    [numf] = num2str(w);
    %disp('# of bal files in process :'); disp ( [numf]);
    sounding_filename = aa(w).name;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% READ THE name 04 ptu FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    year(w) = 2000 + str2num(sounding_filename(1:2));
    month(w) = str2num(sounding_filename(3:4));
    day (w)  = str2num(sounding_filename(5:6));
    hour(w)  = str2num(sounding_filename(7:8));

    fid1=fopen(sounding_filename,'r');
    sounding_filename
    %SEWARD JOHNSON name 2005 MW15 radiosonde data
    %rb04081117a
    %
    %
    %Complete listing of Fine Structure file
    %File: C:\Metgraph\archive\050123.001\05012315.59e
    %
    %Started at      23 January 2005 15:59 UTC
    %
    %Station  :      313
    %Location :      17.90 N  61.67 W         2 m/MSL
    %
    %Sounding type : PTU GPS
    %RS-number     : 4510161
    %Sounding processor serial number: 0
    %
    %
    %Ground check omitted
    %
    %Clouds    : /////
    %Special 1 : /////
    %Special 2 : /////
    %
    %Telemetry noise level  : 100 %
    %
    %                 P      T      U
    %Accepted (%) :  99.7   93.5   99.7
    %Replaced (%) :   0.3    6.5    0.3
    %Rejected (%) :   0.0    0.0    0.0
    %
    %GPS data quality
    %    Valid raw wind ratio        :   97.0 %
    %    Valid raw wind on 5 ratio   :   97.0 %
    %    Invalid raw wind on 4 ratio :    1.0 %
    %    Poor PDOP ratio             :    1.0 %
    %    Unidentified signals ratio  :    0.0 %
    %    Unused                      :    0.0 %
    %    Mean satellite track count  :    7.8
    %    Mean track count valid      :    8.0
    %    Gaps over 1 min             :    0.0 min
    %
    %
    %Reason for termination : Increasing pressure
    %
    %   Time Pressure      Height   Speed  Direction   Significance flags
    %
    % min  s      hPa     gpm     m   m/s    deg  Automatic  Operator
    %
    %   0  0   1016.0       2 /////   6.0    134       dfv         dfv
    %   0  2   1015.7       4 /////   5.5    121       d v         d v
    %   0  4   1015.0      11 /////   5.4    117
    %   0  6   1014.0      19 /////   5.4    113
    %   0  8   1012.7      30 /////   5.4    110


    if fid1 == -1
        error(['Cannot read ' sounding_filename])
    end
    % Read in all header lines
    for n = 1:46
        hd = fgetl(fid1);
        % print out Launch time
        if n == 6
            baltm = hd;
        end
    end
    nhts = 9999;

    k=1;
    %while ~feof(fid1)
    wt = 0;
    while ~feof(fid1)
        %for i=1:1800
        data = fgetl(fid1);
        if feof(fid1) == 1
            break
        end
        if (data(1,14:16) ~= '///' ) &  (data(1,41:43) ~= '///' )    % Test for flagged data
            wt = wt + 1; % count good wind pts
            balw(wt,1) = str2num(data(1,1:4)); %min
            balw(wt,3) = str2num(data(1,8:16));%p
            balw(wt,2) = str2num(data(1,17:24))/1000.;%ht msl km
            balw(wt,6) = str2num(data(1,37:43));%WD
            balw(wt,7) = str2num(data(1,31:36));%WS
        elseif (data(1,14:16) == '///' ) %disregard this level as Ht is flagged
            data
            wt = wt -1;
        else
            wt = wt + 1; % count good ht pt, but flag winds
            balw(wt,1) = str2num(data(1,1:4)); %min
            balw(wt,3) = str2num(data(1,8:16));%p
            balw(wt,2) = str2num(data(1,17:24))/1000.;%ht msl km
            balw(wt,6) = -999;%WD
            balw(wt,7) = -999;%WS
        end

    end
    lb = size(balw);
    if 0
        if lb(1) > 10
            %De-Bob balloon sounding
            bct = 1;
            % save 1st point
            bal1(bct,1) = balw(1,1); %min
            bal1(bct,2) = balw(1,2); %ht
            bal1(bct,3) = balw(1,3); %P
            bal1(bct,6) = balw(1,6); %wd
            bal1(bct,7) = balw(1,7); %ws
            bal1(bct,8) = balw(1,8); %dz
            blast = balw(1,2);
            for bb = 2:(length(balw))
                if balw(bb,2) > blast
                    bct = bct + 1;
                    bal1(bct,1) = balw(bb,1); %min
                    bal1(bct,2) = balw(bb,2); %ht
                    bal1(bct,3) = balw(bb,3); %P
                    bal1(bct,6) = balw(bb,6); %wd
                    bal1(bct,7) = balw(bb,7); %ws
                    bal1(bct,8) = balw(bb,8); %dz
                    blast = balw(bb,2);  % save last good pt assuming lower pt is good
                end
            end
            % put de-bobbed data back into ba array for later comparisons
            %clear bal
            %clear balwtmp
          
            %balw = bal1; % copy de-bobbed file into original name
        end % test > 10 pts to de-bob
    end %skip de-bobbing
    %save data in array names used for interpolating and plotting
    [r c ] = size(balw);
    sz(w,1) = r;
    sz(w,2) = c;
    hwtmp(1:r,w) = balw(:,2);
    balwtmp(1:r,w) = balw(:,7);
    if 0  %w > 1
        [a,c] = size(Hw);
        [b,c1] = size(balw);
        if a > b
            a = b
        end
        if b > a
            b = a
        end
        Hw(1:a,w) = balw(1:b,2)*1000; %convert ht to m
        spd(1:a,w) = balw(1:b,7);
        dr(1:a,w) = balw(1:b,6);
    end
    if 1  %w == 1
        [a,c] = size(balw);
        Hw(1:a,w) = balw(1:a,2)*1000; %convert ht to m
        spd(1:a,w) = balw(1:a,7);
        dr(1:a,w) = balw(1:a,6);
    end
    clear balt
    clear balw
    % READ IN PTU FILE fo same date/time as wnd file
    %at = dir('*.ptu');
    %[m n] = size(at);
    %p = length(at);
    %[num] = num2str(p);
    %disp('Number of rawinsonde *.ptu files : '); disp( [num]);
    %sounding_filename = aa(w).name
    sounding_filename(10:12) = 'ptu';
    fid3=fopen(sounding_filename,'r');
    sounding_filename
    %%%%%%%%%% starting file %%%%%%%%%%%%%%%%%%%%%%%%
    st=1;

    %
    %
    %Complete listing of Fine Structure file
    %File: C:\Metgraph\archive\001007.001\00100716.06e
    %
    %Started at       7 October 2000 16:06 UTC
    %
    %Station  :      313
    %Location :      33.24 N 128.89 W        10 m/MSL
    %
    %Sounding type : PTU GPS
    %RS-number     : 805103110
    %Sounding processor serial number: 0
    %
    %
    %Ground check omitted
    %
    %Clouds    : /////
    %Special 1 : /////
    %Special 2 : /////
    %
    %Telemetry noise level  : 80 %
    %
    %                 P      T      U
    %Accepted (%) :  35.3   64.2   60.7
    %Replaced (%) :  52.1   23.2   26.7
    %Rejected (%) :  12.6   12.6   12.6
    %
    %GPS data quality
    %    Valid raw wind ratio        :   73.0 %
    %    Valid raw wind on 5 ratio   :   73.0 %
    %    Invalid raw wind on 4 ratio :    1.0 %
    %    Poor PDOP ratio             :    1.0 %
    %    Unidentified signals ratio  :    0.0 %
    %    Unused                      :    0.0 %
    %    Mean satellite track count  :    4.7
    %    Mean track count valid      :    6.4
    %    Gaps over 1 min             :    0.0 min
    %
    %
    %Reason for termination : Manual stop
    %
    %   Time Pressure   Height  Temperature  RH   Dewp   Significance flags
    %
    % min  s      hPa      gpm     deg C      %  deg C  Automatic  Operator
    %
    % Read in all header lines
    for n = 1:46
        hd = fgetl(fid3);
        % print out Launch time
        if n == 6
            hd;
        end
    end
    % Read in data lines
    %   Time Pressure   Height  Temperature  RH   Dewp   Significance flags
    % min  s      hPa      gpm     deg C      %  deg C  Automatic  Operator
    %   0  0   1017.0       10      18.3     58    9.9 TU          TU
    %   0  2   1015.5       23      18.3     55    9.1
    %   0  4   1014.7       29      18.2     54    8.8
    tt = 0;
    while ~feof(fid3)
        %for i=1:1800
        data = fgetl(fid3);
        if feof(fid3) == 1
            break
        end

%                  T                            HT                          RH                            
        if (data(1,33:35) ~= '///') &  (data(1,14:16) ~= '///' )  &  (data(1,40:42) ~= '///' )  % Test for good data
            tt = tt + 1; %count good data pts
            balt(tt,1) = str2num(data(1,1:4)); %min
            balt(tt,2) = str2num(data(1,17:25))/1000.;%ht msl km
            balt(tt,3) = str2num(data(1,8:16));%p
            balt(tt,4) = str2num(data(1,26:35));%Temp
            balt(tt,5) = str2num(data(1,36:42));%RH
        elseif (data(1,33:35) ~= '///') &  (data(1,40:42) == '///' )  % Test for good T and bad RH data
            tt = tt + 1; %count good data pts
            balt(tt,1) = str2num(data(1,1:4)); %min
            balt(tt,2) = str2num(data(1,17:25))/1000.;%ht msl km
            balt(tt,3) = str2num(data(1,8:16));%p
            balt(tt,4) = str2num(data(1,26:35));%Temp
            balt(tt,5) = -999;  %RH
        elseif (data(1,14:16) == '///' ) %disregard this level as Ht is flagged
            data
            tt = tt -1;
        else
            tt = tt + 1; % count good ht pt, but flag T RH
            balt(tt,1) = str2num(data(1,1:4)); %min
            balt(tt,2) = str2num(data(1,17:25))/1000.;%ht msl km
            balt(tt,3) = str2num(data(1,8:16));%p
            balt(tt,4) = -999;  %T
            balt(tt,5) = -999;  %RH
        end
    end
    %De-Bob balloon sounding
    clear bal1
    % Test for > 10 pts to de-bob
    if 0
        lb = size(balw);
        if lb(1) > 10
            bct = 1;
            %save 1st point
            bal1(bct,1) = balt(bb,1); %min
            bal1(bct,2) = balt(1,2); %ht
            bal1(bct,3) = balt(1,3); %P
            bal1(bct,4) = balt(1,4); %T
            bal1(bct,5) = balt(1,5); %RH
            blast = balt(1,2);
            for bb = 2:(length(balt))
                if balt(bb,2) > blast
                    bct = bct + 1;
                    bal1(bct,1) = balt(bb,1); %min
                    bal1(bct,2) = balt(bb,2); %ht
                    bal1(bct,3) = balt(bb,3); %P
                    bal1(bct,4) = balt(bb,4); %T
                    bal1(bct,5) = balt(bb,5); %RH
                    blast = balt(bb,2);  % save last good pt assuming lower pt is good
                end
            end
            % put de-bobbed data back into ba array for later comparisons
            %clear bal

            baltmpt = balt; % save original data for testing
            clear balt
            balt = bal1; % copy de-bobbed file into original name
        end % > 10 pts to de-bob
    end % skip de-bobbing
    [r c ] = size(balt);
    sz(w,3) = r;
    sz(w,4) = c;
    httmp(1:r,w) = balt(:,2);
    balttmp(1:r,w) = balt(:,4);

    %save data into arrays for interpolating and plotting
    if 0  %w > 1
        [a,c] = size(Ht);
        [b,c1] = size(balt);
        if a > b
            a = b;
        end
        if b > a
            b = a;
        end
        Ht(1:a,w) = balt(1:b,2)*1000; %convert ht to m
        P(1:a,w) = balt(1:b,3);
        T(1:a,w) = balt(1:b,4);
        RH(1:a,w) = balt(1:b,5);
    end

    if 1  %w == 1
        [a,c] = size(balt);
        Ht(1:a,w) = balt(1:a,2)*1000; %convert ht to m
        P(1:a,w) = balt(1:a,3);
        T(1:a,w) = balt(1:a,4);
        RH(1:a,w) = balt(1:a,5);
    end
    clear balt
end % end of looping through ptu and wnd files
%%%%%%%%%%%%%%%%%
fnamet = sounding_filename(1:9);
% plotting data
%temp plot
if 0
    figure(1)
    %subplot(1,1,1),plot(balt(:,4),balt(:,2),'b')
    subplot(1,2,1),plot(balt(:,4),balt(:,2),'b');
    axis([-25,30,0,7])
    xlabel('Temp C')
    ylabel('Ht km-MSL')
    %hold

    subplot(1,2,2),plot(balt(:,5),balt(:,2),'b');
    title(fnamet)
    axis([0,100,0,7])
    xlabel('RH %')
    ylabel('Ht km-MSL')
    %hold
    fclose(fid3); % ptu file
end % ptu file found

%%%% read the sounding, now QC and interpolation %%%%%%
%%%% missing value is -32768
%%%%% Interpolate the sounding....from surface to 6 km %%%%%%%%%%%%%%%%%%%
%%% FIRST QUALITY CONTROL....
h_int = [0:15:50000];

for i=ws:we;
    aa(i).name
    ll = length(Ht(:,i));
    g = find (Ht(:,i)>0 & Ht(:,i)<=50000);
    gk = find(diff(Ht(g,i))>0);  % checks for balloon going down
    c = g(gk);clear g; g=c;clear c gk;
    
    f=find(T(g,i)>-273); %find good temperatures
    if isempty(f)==0
        t_int(1:length(h_int),i)     =  interp1(Ht(g(f),i),T(g(f),i),h_int)'; clear f;
    end
    f=find(RH(g,i)>=0);
    if isempty(f)==0
        rh_int(1:length(h_int),i)     =  interp1(Ht(g(f),i),RH(g(f),i),h_int)'; clear f;
    end
    f=find(P(g,i)>0);
    if isempty(f)==0
        press_int(1:length(h_int),i)     =  interp1(Ht(g(f),i),P(g(f),i),h_int)'; clear f;
    end
    f=find(dr(g,i)>0);
    if isempty(f)==0
        dr_int(1:length(h_int),i)     =  interp1(Hw(g(f),i),dr(g(f),i),h_int)'; clear f;
    end
    f=find(spd(g,i)>0);
    if isempty(f)==0
        spd_int(1:length(h_int),i)     =  interp1(Hw(g(f),i),spd(g(f),i),h_int)'; clear f;
    end
    % Put surface values in 0 ht interpolated arrays
    t_int(1,i) = T(1,i);
    rh_int(1,i) = RH(1,i);
    press_int(1,i) = P(1,i);
    spd_int(1,i) = spd(1,i);
    dr_int(1,i) = dr(1,i);

end

%save c:\data\name\balloon\results\name_2005_sounding press_int t_int rh_int spd_int dr_int h_int year month day hour T;


%%%%%% SOUNDINGS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%clear all

ymax = 20;  % used to scale images km

%load c:\data\rico\balloon\results\ICEALOT_2008_sounding
t=t_int;%-273.15;
rh=rh_int;
press=press_int;
h_int = [0:15:50000];
h=h_int;
k=length(day);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

es        = 6.112.*exp((17.67.*t)./(t+243.5));
e         = (rh.*es)/100;
r         = 622*(e./(press-e));
rs        = 622*(es./(press-es));
tk        = t+273.15;
a         = 1./(tk-55);
b         = log(rh./100)/2480;
tl        = 1./(a-b)+55;
theta     = tk.*(1000./press).^0.2854.*(1-0.28*10^(-3).*r);
thetav    = theta.*(1+0.608.*r/1000);
thetae    = theta.*exp(((3.376./tl)-0.00254).*r.*(1+0.81*10^(-3).*r));
thetaes   = theta.*exp(2.675.*rs./tk);


%%%%%%%% SOUNDING QC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% LCL Calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



for w = 1:k;

    i_lev = 8;

    init_tk = tk(i_lev,w);
    init_p  = press(i_lev,w);
    init_r  = r(i_lev,w);
    init_rh = rh(i_lev,w);
    init_thetae = thetae(i_lev,w);
    init_theta = theta(i_lev,w);
    i_1005 = 8;
    init_a  = 1/(init_tk-55);
    init_b  = log(init_rh/100)/2480;
    t_lcl(w)   = 1/(init_a-init_b)+55;
    p_lcl(w) = 1000*(t_lcl(w)/init_theta)^(1/0.2854);

    kke = find(press(:,w)>=400);
    kke = max(kke);
    zz = find(abs(p_lcl(w)-press(1:kke,w))==min(abs(p_lcl(w)-press(1:kke,w))));

    if isempty(zz)==1
        zp = i_1005;
    else
        zp = zz;
    end


    i_lcl(w) = zp;
    h_lcl(w) = h(zp);
    p_lcl(w) = press(zp);




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%% SOUNDING CLASSIFICATION %%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear buoyancy;
    clear intersect;

    buoyancy = init_thetae-thetaes(i_lcl(w):kke,w);

    rec_buoyancy(i_lcl(w):kke,w)=init_thetae-thetaes(i_lcl(w):kke,w);

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%% COMMON TIME COORDINATE %%%%%

for q=1:w
    jul_day_sonde(q) = datenum(year(q),month(q),day(q),hour(q),00,00)-datenum(year(q)-1,12,31,00,00,00);
end;



%save c:\data\rico\balloon\results\ICEALOT_2008_sound_results


graphs = -1;
% pause
if 1  % skip graphs
s = 'c:\data\icealot\data\sonde\Raw_Images\';
figure(1)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(theta(1:1667,:)),[270,390]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, \theta (k)','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_theta.png'];
print('-dpng ', print_buffer);

figure(2)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(thetav(1:1667,:)),[270 390]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, \theta_v (k)','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_thetav.png'];
print('-dpng ', print_buffer);

figure(3)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(rh_int(1:1667,:)),[0 100]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, Relative Humidity (%)','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_RH.png'];
print('-dpng ', print_buffer);


figure(4)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(r(1:1667,:)),[0 10]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, Mixing Ratio (g/kg)','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_MR.png'];
print('-dpng ', print_buffer);


figure(5)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(dr_int(1:1667,:)),[0 360]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, Direction','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_DIR.png'];
print('-dpng ', print_buffer);


figure(6)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(spd_int(1:1667,:)),[0 75]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, Wind Speed (m/s)','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_SPD.png'];
print('-dpng ', print_buffer);

figure(7)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int(1,1:1667)/1000,real(t_int(1:1667,:)),[-90 25]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Year Day (UTC)','FontName','times','FontSize',18);
ylabel('Ht  (km)','FontName','times','FontSize',18);
title('ICEALOT 2008, Temperature (C)','FontName','times','FontSize',18);
print_buffer = [s 'ICEALOT_2008_T.png'];
print('-dpng ', print_buffer);


% find mean profile for all soundings

for z = 1:(length(h_int))
    j = find(~isnan(t_int(z,:)));
    p_intbar(z,1) = mean(t_int(z,j));
    j = find(~isnan(theta(z,:)));
    p_intbar(z,2) = mean(theta(z,j));
    j = find(~isnan(rh(z,:)));
    p_intbar(z,3) = mean(rh(z,j));
    j = find(~isnan(t_int(z,:)));
    p_intbar(z,4) = mean(spd_int(z,j));
end
figure(8)
ylim('manual')
subplot(2,2,1),plot(p_intbar(:,1),h_int/1000)
ylim([0,ymax])
title('Mean Temp')
ylabel('Ht (km)')
subplot(2,2,2),plot(p_intbar(:,2),h_int/1000)
ylim([0,ymax])
title('Mean Theta')
ylabel('Ht (km)')
subplot(2,2,3),plot(p_intbar(:,3),h_int/1000)
ylim([0,ymax])
title('Mean RH')
ylabel('Ht (km)')
subplot(2,2,4),plot(p_intbar(:,4),h_int/1000)
ylim([0,ymax])
title('Mean Speed')
ylabel('Ht (km)')
%print_buffer = ['c:\data\rico\balloon\results\ICEALOT_2008_Mean_profiles.png'];
%print('-dpng ', print_buffer);
end  % skip graphs