%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%STRATUS 2006 SOUNDING READER %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%READS EDT %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
warning off
read_parameters  % DEW 10092008 setup for Dan's Mapped H: on laptop

eval(['cd ', way_raw_data_balloon]);
ls([way_raw_data_balloon '*LAT*.txt'])

%%%%%%%%%%%%%START READING AND LOADING PTU FILES %%%%%%%%%%%%%%%%%%%%%%
aa = dir('*LON.txt');

k = length(aa);
%% DEW 5/16/2007  Uncomment the next 3 lines if you want to process 1 specific file at a time
%% Next 3 line should be commented if you want to process all the soundings
dt=input('Input date/time YYMMDDhhmm (ex:0705060111_EDT_LATLON.txt) ', 's');
aa(1).name(1:10) = dt(1:10);
k = 1;

disp('Number of radiosondes : '); disp(k);

%%%%%%%%% starting file is st%%%%%%%%%%%%%%%%%%%%%%%%
st =1;
%%%%%%%%%%%%%%%%%% last file == k %%%%%%%%%%%%%%%

for i = st:k;
    close all
    i;
    aa(i).name
    dt(i,1:10) = aa(i).name(1:10);
    %   0705060111_EDT_LATLON.txt
    yr(i) = str2num(aa(i).name(1:2))+2000;
    month(i) = str2num(aa(i).name(3:4));
    day (i)  = str2num(aa(i).name(5:6));
    hour(i)  = str2num(aa(i).name(7:8));
    minute(i) = str2num(aa(i).name(9:10));

    fid1=fopen([aa(i).name],'r');
    balfn = aa(i).name;

    if fid1 == -1
        error(['Cannot read ' aa(i).name]);
    end %if end

    %%%%%%%%%%%START READING HEADER %%%%%%%%%%%%%%%%%%%%%%%
    %                                                  RS Serial #: B2815335
    %                                                  Station Name: WTEC
    %                                                  Launch Time: 0111Z
    %                                                  Launch Date: 06 MAY 07
    %
    %EDT LEVEL OUTPUT with GPS Lat, Long, Alt
    %
    %    Time     Height  Pressure      Temp    Dew P.    RH    W Spd    W Dir             Lat            Long      GPS Alt
    %     sec       mtrs       hPa        øC       øC    Pct      m/s     Az ø       Decimal ø       Decimal ø         mtrs
    %
    %%%%
    for j=1:17
        ij=fgetl(fid1);
    end %for end

    clear ij;

    %%%%%%%%%%%END READING HEADER %%%%%%%%%%%%%%%%%%%%%%%%

    k=0;
    %   -0.88        7.0    1013.0     25.65     20.9     75     10.0     80.0      11.9248377     -44.2862418        -21.3
    %    0.12        7.5    1012.9     25.65     20.8     75      9.7     79.6      11.9247976     -44.2862551        -14.5
    while ~feof(fid1)
        temp=fgetl(fid1);
        if ~isempty(temp)
            k=k+1;
            [data]=sscanf(temp,'%f') ;
            h(k,i)     = data(2); %Ht m
            press(k,i) = data(3);  %Press mb
            t(k,i)     = data(4);  % T C
            dp(k,i)    = data(5);  %dew Pt C
            rh(k,i)    = data(6);  % RH %
            w_spd(k,i)=data(7);  % wnd spd m/s
            w_dir(k,i)=data(8);  % wnd dir true
            lat(k,i)=data(9); % tenths of degrees N = + S = -
            lon(k,i)=data(10);  % tenths of degrees W = - E = +
        end % if end
    end %while end
    balt(:,1) = h(:,1);
    balt(:,2) = press(:,1);
    balt(:,3) = t(:,1);
    balt(:,4) = rh(:,1);
    fclose(fid1);
    
    %write out thermo data that will be used to calculate pwv
    if 1
    f=[way_proc_data_balloon cruise '_sonde_' aa(i).name(1:2) '_' aa(i).name(3:4) '_' aa(i).name(5:6) '_' aa(i).name(7:8) aa(i).name(9:10) '.txt'];
    flist=fopen(f,'w');
    save( f, 'balt', '-ascii','-tabs')
    fclose('all');
end;

    %%%%%%%%%%%%%%%%%FINISH READING THE EDT FILES %%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%QUALITY CONTROL %%%%%%%%%%%%%%%%%%%

    h_int = [0:10:25000];

    for i=1:length(yr)

        g=find(h(:,i)>0 & h(:,i)<25000);

        %First get height between 0 to 25000m
        gk=find(diff(h(g,i))>0); %% then get values with positive difference
        c = g(gk); %now we just have ascending values between 0 to 25000
        clear g;
        g=c; %g has the positive row and column vectors between 0 to 10000m
        clear c gk;

        gk=find(diff(h(g,i))>0); %same thing again
        c = g(gk);
        clear g;
        g=c;
        clear c gk;

        gk=find(diff(h(g,i))>0); %same thing again
        c = g(gk);
        clear g;
        g=c;
        clear c gk;

        gk=find(diff(h(g,i))>0); %same thing again
        c = g(gk);
        clear g;
        g=c;
        clear c gk;


        gk=find(diff(h(g,i))>0); %same thing again
        c = g(gk);
        clear g;
        g=c;
        clear c gk;

        if isempty(g)==0
            t_int(1:length(h_int),i)     =  interp1(h(g,i),t(g,i),h_int)';
            rh_int(1:length(h_int),i)    =  interp1(h(g,i),rh(g,i),h_int)';
        end


        gg = find(press(g,i)<1020);
        if isempty(gg)==0
            press_int(:,i) =  interp1(h(g(gg),i),press(g(gg),i),h_int)';
        end
        clear gg;

        gg = find(rh(g,i)<110);
        if isempty(gg)==0
            rh_int(:,i) =  interp1(h(g(gg),i),rh(g(gg),i),h_int)';
        end
        clear gg;

        gg = find(w_spd(g,i)<1000);
        if isempty(gg)==0
            spd_int(:,i) =  interp1(h(g(gg),i),w_spd(g(gg),i),h_int)';
        end
        clear gg;


        gg = find(w_dir(g,i)<1000);
        if isempty(gg)==0
            dir_int(:,i) =  interp1(h(g(gg),i),w_dir(g(gg),i),h_int)';
        end
        clear gg;

        gg = find(lat(g,i)<1000);
        if isempty(gg)==0
            lat_int(:,i) =  interp1(h(g(gg),i),lat(g(gg),i),h_int)';
        end
        clear gg;

        gg = find(lon(g,i)<1000);
        if isempty(gg)==0
            lon_int(:,i) =  interp1(h(g(gg),i),lon(g(gg),i),h_int)';
        end
        clear gg;

    end

end

ss = [way_raw_data_balloon, 'VOCALS08_sounding'];
save   ss press_int t_int rh_int spd_int dir_int h_int lat_int lon_int yr month day hour month minute dt;

%%%%%%%%%%%%%%%END OF LOADING AND SAVING %%%%%%%%%%%%%%%%%%

%%%%%%START OF ANALYSIS %%%%%%%%%%%%%%%%%%%%%
clear all;
read_parameters_VOCALS
load ss;

%lat_deg=floor(lat);
%lon_deg=floor(lon);
t=t_int;
rh=rh_int;
press=press_int;
h=h_int;
k=length(yr);


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


for w = 1:k

    if w==10 | w==16 | w==17  | w==18 | w==35
        continue;
    end


    i_lev = 10;

    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 = 10;
    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);


    i_1005 = 10;

    for i = i_1005:201;
        bnta(i,w) = 0.608*theta(i,w)/(1+0.608*r(i,w)/1000);
        dtheta_dp(i,w) = (theta(i+1,w)-theta(i,w))/(press(i+1,w)-press(i,w));
        dr_dp(i,w)     = (10^(-3)*(r(i+1,w)-r(i,w)))/(press(i+1,w)-press(i,w));
        mu(i,w) = -(dtheta_dp(i,w) - bnta(i,w)*dr_dp(i,w));
    end


end

jul_day_sonde = datenum(yr,month,day,hour,minute,0)-datenum(2006,12,31);

ss = [way_raw_data_balloon, 'VOCALS08_sound_results']
save ss
clear all

%%%% graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load ss;
%graphs = -1


qw = length(jul_day_sonde);

for u=1:qw;
    close all
    %u = 1;
    figure(1)
    %orient tall
    set(gcf,'DefaultAxesFonts',14)
    subplot(1,2,1)
    h=plot(theta(:,u),h_int/1000,'r'); grid on
    set(h,'linewidth',2);xlim([280 360]);ylim([0 20])
    hold on
    % h=plot(thetaes(:,u),h_int/1000,'g'); grid on
    % set(h,'linewidth',2);xlim([280 380]);
    h=plot(thetav(:,u),h_int/1000,'b'); grid on
    set(h,'linewidth',2);xlim([280 400]);
    % h=plot(thetae(:,u),h_int/1000,'k'); grid on
    % set(h,'linewidth',2);xlim([280 380]);
    xlabel('T (K)','FontName','times','FontSize',14);
    ylabel('Ht (km)','FontName','times','FontSize',14);
    title([ 'VOCALS 2008 Day: ' dt(u,1:10)],'FontName','times','FontSize',14);
    legend('theta','thetav')

    subplot(1,2,2)
    set(gcf,'DefaultAxesFonts',14)
    %h=plot(r(:,u),h_int/1000,'k'); grid on
    set(h,'linewidth',2);
    ylim([0 20])
    %hold on
    h=plot(rh(:,u),h_int/1000,'b'); grid on
    set(h,'linewidth',2);
    xlabel('RH(%)','FontName','times','FontSize',14);
    xlim([0 100]);
    %ylabel('Height (km)','FontName','times','FontSize',14);
    %title(['VOCALS 2008 Day: ' num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    print_buffer = [way_proc_images_balloon, cruise, '.sounding_trh_' dt(u,1:10), '.png'];
    print('-dpng ', print_buffer);


    figure(2)
    subplot(1,2,1)
    h=plot(spd_int(:,u),h_int/1000,'k'); grid on
    set(h,'linewidth',2);xlim([0 20]);ylim([0 20])
    xlabel('Wind Spd (ms^{-1})','FontName','times','FontSize',14);
    ylabel('Ht (km)','FontName','times','FontSize',14);
    title([ 'VOCALS 2008 Day: ' dt(u,1:10)],'FontName','times','FontSize',14);

    subplot(1,2,2)
    h=plot(dir_int(:,u),h_int/1000,'k.'); grid on
    set(h,'linewidth',2);xlim([0 360]);ylim([0 20])
    xlabel('Wind Dir (True)','FontName','times','FontSize',14);
    %ylabel('Height (m)','FontName','times','FontSize',14);
    %title([ 'VOCALS 2008, Wdir, Day: '
    %num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    print_buffer = [way_proc_images_balloon, cruise, '.sounding_wnd_' dt(u,1:10), '.png'];
    print('-dpng ', print_buffer);


    figure(3)
    %orient tall
    set(gcf,'DefaultAxesFonts',14)
    subplot(1,2,1)
    h=plot(t(:,u),h_int/1000,'r'); grid on
    set(h,'linewidth',2);xlim([-10 35]);ylim([0 5])
    xlabel('T (C)','FontName','times','FontSize',14);
    ylabel('Ht (km)','FontName','times','FontSize',14);
    title([ 'VOCALS 2008 Day: ' dt(u,1:10)],'FontName','times','FontSize',14);

    subplot(1,2,2)
    set(gcf,'DefaultAxesFonts',14)
    h=plot(rh(:,u),h_int/1000,'b'); grid on
    set(h,'linewidth',2);;xlim([0 100]);ylim([0 5])
    xlabel('RH(%)','FontName','times','FontSize',14);
    print_buffer = [way_proc_images_balloon, cruise, '.sounding_trha_' dt(u,1:10), '.png'];
    print('-dpng ', print_buffer);

    figure(4)
    subplot(1,2,1)
    h=plot(spd_int(:,u),h_int/1000,'k'); grid on
    set(h,'linewidth',2);xlim([0 20]);ylim([0 5])
    xlabel('Wind Spd (ms^{-1})','FontName','times','FontSize',14);
    ylabel('Ht (km)','FontName','times','FontSize',14);
    title([ 'VOCALS 2008 Day: ' dt(u,1:10)],'FontName','times','FontSize',14);

    subplot(1,2,2)
    h=plot(dir_int(:,u),h_int/1000,'k.'); grid on
    set(h,'linewidth',2);xlim([0 360]);ylim([0 5])
    xlabel('Wind Dir (True)','FontName','times','FontSize',14);
    %ylabel('Height (m)','FontName','times','FontSize',14);
    %title([ 'VOCALS 2008, Wdir, Day: '
    %num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    print_buffer = [way_proc_images_balloon, cruise, '.sounding_wnda_' dt(u,1:10), '.png'];
    print('-dpng ', print_buffer);
end %plotting loop
