%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Project Name: RICO 2005 
%Purpose: Reads the ptu and wnd files and plots the sounding data
%
%Input Data: 
%
%Last Modified: Kristen Rasmussen 6/16/05
%
%%%%%%%%%%READS THE PTU AND WND FILES %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
fclose all;

addpath 'C:\rico'
addpath 'C:\rico\raw'

pathdata = 'C:\rico\raw';
outpath  = 'C:\rico\sonde\results';
eval(['cd ', pathdata]);

%%%%%%%%%%%%%START READING AND LOADING PTU FILES %%%%%%%%%%%%%%%%%%%%%%
aa = dir('*.ptu');
k = length(aa);
disp('Number of rawinsonde *.ptu files : '); disp(k); 
%%%%%%%%% starting file is st%%%%%%%%%%%%%%%%%%%%%%%%
st =1;
%%%%%%%%%%%%%%%%%% last file == k %%%%%%%%%%%%%%%

for w = st:k;
w

year(w) = 2000 + str2num(aa(w).name(1:2));
month(w) = str2num(aa(w).name(3:4));
day (w)  = str2num(aa(w).name(5:6));
hour(w)  = str2num(aa(w).name(7:8));

fid1=fopen(aa(w).name,'r');

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

%%%%%%%%%%%START READING HEADER %%%%%%%%%%%%%%%%%%%%%%%

for i=1:46
    ij=fgets(fid1);
    if i==6
        mmm(w)=str2num(ij(37:38));
    end %if end
    
    if i==9
        lat(w) = str2num(ij(17:21));
        lat_dir(w) = ij(23);
        lon(w) = str2num(ij(26:30));
        lon_dir(w) = ij(32);
    end %if end
   
end %for end

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

i=0;
while ~feof(fid1)
    temp=fgetl(fid1);num=length(findstr(temp,'/'));
    if num < 10
    i=i+1;    
        [data]=sscanf(temp,'%d %d %f %f %f %f') ;
       press(i,w) = data(3);
       h(i,w)     = data(4);
       t(i,w)     = data(5);
       rh(i,w)    = data(6);
    end %if end
end %while end
fclose(fid1);
%%%%%%%%%%%%%%%%%FINISH READING THE PTU FILES %%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%START READING THE WND FILES %%%%%%%%%%%%%%%%%%%%%%%%

fid2=fopen([aa(w).name(1:8) '.wnd']);

if fid2 == -1
error(['Cannot read ' ceil_filename]);
end %if end


for es=1:46
    sv=fgets(fid2);
end %for end

i=0;
while ~feof(fid2)
    temp=fgetl(fid2);num=findstr(temp,'/');temp(num)='9';
    i=i+1;    
        [data]=sscanf(temp,'%d %d %f %f %f %f %f') ;
       press(i,w) = data(3);
       h(i,w)  = data(4);
       spd(i,w)    = data(6);
       dir(i,w)    = data(7);
end %while end
 
end %for end

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



h_int = [0:10:4000];

for i=1:w;
    i
    if i==17
        g=find(h(1:550,i)>0 & h(1:550,i)<4000);
    elseif i==14
        g=find(h(3:530,i)>0 & h(3:530,i)<4000);  
    elseif i==31
        g=find(h(850:end,i)>0 & h(850:end,i)<4000);
    else
    g=find(h(:,i)>0 & h(:,i)<4000);
    end
    %First get height between 0 to 4000m
    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 4000
    clear g;
    g=c; %g has the positive row and column vectors between 0 to 4000m 
    clear c gk;
%    plot(h(g,i));pause(.01)
    
    gk=find(diff(h(g,i))>0); %same thing again 
    c = g(gk);
    clear g; 
    g=c;
    clear c gk;
%    plot(h(g,i));pause(.01)
     
    gk=find(diff(h(g,i))>0); %same thing again 
    c = g(gk);
    clear g; 
    g=c;
    clear c gk;
%    plot(h(g,i));pause(.01)
    
    gk=find(diff(h(g,i))>0); %same thing again 
    c = g(gk);
    clear g; 
    g=c;
    clear c gk;
%    plot(h(g,i));pause(.01)
    
    gk=find(diff(h(g,i))>0); %same thing again 
    c = g(gk);
    clear g; 
    g=c;
    clear c gk;
%    plot(h(g,i));pause(.01)
   % return
    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)<10000);
    if isempty(gg)==0
        press_int(:,i) =  interp1(h(g(gg),i),press(g(gg),i),h_int)';
    end
        clear gg; gg = find(spd(g,i)<10000);
    if isempty(gg)==0
        spd_int(:,i) =  interp1(h(g(gg),i),spd(g(gg),i),h_int)';
    end
        clear gg; gg = find(dir(g,i)<10000);
    if isempty(gg)==0
        dir_int(:,i) =  interp1(h(g(gg),i),dir(g(gg),i),h_int)';
    end
    clear g;
end

save 'C:\rico\sonde\results\RICO_sounding press_int t_int rh_int spd_int dir_int h_int year month day hour mmm lat lat_dir lon lon_dir';

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

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

load 'C:\rico\sonde\results\RICO_sounding press_int t_int rh_int spd_int dir_int h_int year month day hour mmm lat lat_dir lon lon_dir';
lat_deg = floor(lat);
lon_deg = floor(lon);
t = t_int;
rh = rh_int;
press = press_int;
h_int = [0:10:4000];
h = h_int;
k = length(lat);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 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-0.0001)/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

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-0.01/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



peak  = find(mu(3:201,w)==max(mu(3:201,w)))+2;
dummy = find(mu(1:201,w)<0.3);
inv_b(w) = dummy(max(find(peak-dummy>0)));
inv_top_height(w) = h(inv_t(w));
inv_top_theta(w) = theta(inv_t(w),w);
inv_top_r(w)     = r(inv_t(w),w);
inv_base_height(w) = h(inv_b(w));
inv_base_theta(w) = theta(inv_b(w),w);
inv_base_r(w)     = r(inv_b(w),w);
end




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


%%%% COMMON TIME COORDINATE %%%%%
jul_day_sonde = get_julian_day(year,month,day,hour)-get_julian_day(2004,12,31);

save 'C:\rico\sonde\results\rico_sound_results';
clear all
load  'C:\rico\sonde\results\rico_sound_results';
graphs = -1
pause

qw = length(jul_day_sonde);

for u=1:qw;
    figure(u)
    orient tall
    set(gcf,'DefaultAxesFonts',14)
    subplot(2,2,1)
    h=plot(theta(:,u),h_int/1000,'r'); grid on
    set(h,'linewidth',2);xlim([290 360]);
    hold on
    h=plot(thetaes(:,u),h_int/1000,'g'); grid on
    set(h,'linewidth',2);xlim([290 360]);
    h=plot(thetav(:,u),h_int/1000,'b'); grid on
    set(h,'linewidth',2);xlim([290 360]);
    xlabel('\theta , \thetaes , \thetav (^{o}K)','FontName','times','FontSize',14);
    ylabel('Height (km)','FontName','times','FontSize',14);
    title([ '\theta(Red), \thetaes(Green), \thetav(Blue), Day: ' num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    subplot(2,2,2)
    h=plot(r(:,u),h_int/1000,'k'); grid on
    set(h,'linewidth',2);xlim([0 20]);
    hold on
    h=plot(rh(:,u)/10,h_int/1000,'k--'); grid on
    set(h,'linewidth',2);
    xlabel('0.1*RH(%), r (gr/kr)','FontName','times','FontSize',14);
    %ylabel('Height (km)','FontName','times','FontSize',16);
    title(['RH(Dotted), Mixing ratio(Solid), Day: ' num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    subplot(2,2,3)
    h=plot(spd_int(:,u),h_int/1000,'k'); grid on
    set(h,'linewidth',2);xlim([0 30]);
    xlabel('Wind Speed (ms^{-1})','FontName','times','FontSize',14);
    ylabel('Height (km)','FontName','times','FontSize',14);
    title([ 'RICO 05, Wspd, Day: ' num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    subplot(2,2,4)
    h=plot(dir_int(:,u),h_int/1000,'k.'); grid on
    set(h,'linewidth',2);xlim([0 360]);
    xlabel('Wind Direction (^{o})','FontName','times','FontSize',14);
    %ylabel('Height (m)','FontName','times','FontSize',16);
    title([ 'RICO 05, Wdir, Day: ' num2str(jul_day_sonde(u))],'FontName','times','FontSize',14);
    print_buffer = ['-djpeg90 F:\rico\sonde\results\images\RIC02005.sounding.' num2str(jul_day_sonde(u)) '.jpg'];
    eval(['print ',print_buffer]);
    %pause
       
end 

k = find(inv_top_r-inv_base_r<0);

figure(95)
set(gcf,'DefaultAxesFonts',16)
plot(inv_top_theta(k)-inv_base_theta(k),inv_top_r(k)-inv_base_r(k),'ks');
xlabel('\Delta \theta_{inversion} (^{o}K)','FontName','times','FontSize',16);
ylabel('\Delta r_{inversion} (g/kg)','FontName','times','FontSize',16);
title('RICO 2005, Soundings','FontName','times','FontSize',16);

print -djpeg90 F:\rico\sonde\041231.001\results\rico_dtdr.jpg



figure(96)
set(gcf,'DefaultAxesFonts',16)
plot((inv_top_theta(k)-inv_base_theta(k))./(inv_top_height(k)-inv_base_height(k))*1000,(inv_top_r(k)-inv_base_r(k))./(inv_top_height(k)-inv_base_height(k))*1000,'ks');
xlabel('(\Delta \theta)/(\Delta z) _{inversion} (^{o}K/Km)','FontName','times','FontSize',16);
ylabel('(\Delta r)/(\Delta z) _{inversion} (g/kg/Km)','FontName','times','FontSize',16);
title('RICO 2005, Soundings','FontName','times','FontSize',16);
print -djpeg90 F:\rico\sonde\results\images\rico_gradtgradr.jpg



figure(97)
set(gcf,'DefaultAxesFonts',16);
imagesc(real(mu(1:201,:)),[-.5 2.5]);
axis xy; colorbar;
hold on
plot(inv_b(k),'k.');
plot(inv_t(k),'k.');
set(gca,'xtick',[1:5:w]);
set(gca,'xticklabel',[ floor(jul_day_sonde(1:5:w))]);
set(gca,'ytick',[51 101 151 201]);
set(gca,'yticklabel',[ h_int(51:50:201)]);
xlabel('Day','FontName','times','FontSize',16);
ylabel('Height (m)','FontName','times','FontSize',16);
title('RICO 2005, \mu = \Delta\theta/\DeltaP - \beta * \Deltar/\DeltaP, \beta = 0.608\theta/(1+0.608r)','FontName','times','FontSize',16);
print -djpeg90 F:\rico\sonde\results\images\rico_mu.jpg

close('all');