%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%SOUNDING READER %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%READS EDT %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display('sounding_plot')
clear all;
close all;
warning off
read_parameters_rhb

pathdata = way_raw_data_balloon 

tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file

%eval(['cd ', pathdata]);

%%%%%%%%%%%%%START READING AND LOADING PTU FILES %%%%%%%%%%%%%%%%%%%%%%
fname = '*Lon.txt';
s1=[way_raw_data_balloon fname];

aa = dir(s1);

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
    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([way_raw_data_balloon 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

    fclose(fid1);
    %%%%%%%%%%%%%%%%%FINISH READING THE EDT FILE %%%%%%%%%%%%%%%%%%%%%%

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

    h_int = [0:10:25000];

    for i=1: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

save sounding press_int t_int rh_int spd_int dir_int h_int lat_int lon_int yr month day hour month minute dt tit aa;


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

%%%%%%START OF ANALYSIS %%%%%%%%%%%%%%%%%%%%%
clear all;

load sounding;

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


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

save sound_results
clear all



%%%% graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
read_parameters_rhb
load sound_results
%graphs = -1


qw = length(jul_day_sonde);


for u=1:qw;
    close all
       % Create title for figures
    tit = cruise;
    lt = length(cruise);
    tit(lt+1) = '-';
    tit((lt+2):lt+5) = year; % This year is read in from the parameter file
    tit(lt+6) = '-';
    tit(lt+7:lt+8) = aa(u).name(3:4);
    tit(lt+9:lt+10) = aa(u).name(5:6);
    tit(lt+11:lt+12) = aa(u).name(7:8);
    tit(lt+13+14= aa(u).name(9:10);
    %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( tit, '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(tit,'FontName','times','FontSize',14);

        print_buffer = [way_images_balloon '_' tit '_sounding_trh.jpg'];
        print('-djpeg90 ', 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(tit,'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([ 'AMMA 2007, Wdir, Day: '
        %num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
        print_buffer = [way_images_balloon '_' cruise '_' year '_sounding_wnd.jpg'];
        print('-djpeg90 ', 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(tit,'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_images_balloon '_' cruise '_' year '_sounding_trha.jpg'];
        print('-djpeg90 ', 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(tit,'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([ 'AMMA 2007, Wdir, Day: '
        %num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
        print_buffer = [way_images_balloon '_' cruise '_' year '_sounding_wnda.jpg'];
        print('-djpeg90 ', print_buffer);

    end %plotting loop

%Save sounding data to use in calculating PW
save sondetmp t rh press theta h_int jul_day_sonde tit


%close('all');