%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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)  = 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;
spd = spd_int
dir = dir_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);
u_wnd     = spd_int.*-sin((3.14*dir)/180)
v_wnd     = spd_int.*-cos((3.14*dir)/180)


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

%[mx,peak] = max(2:201,w)+2;
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)));

pd = find((peak-dummy)>0);
max_pd = max(pd);
dum_max_pd = dummy(max_pd);
inv_b(w) = dum_max_pd;
if ~isempty(pd)
    min_pd = min(pd);
    dum_min_pd = dummy(min_pd);
    inv_t(w) = dum_min_pd;
else
    inv_t(w) = 0
end

%inv_t(w) = dummy(min(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

fname = aa(57).name; 

%%%% Graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%qw = length(fname);

%for u=1:qw;
   figure(1)
   orient tall
   set(gcf,'DefaultAxesFonts',10)
   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',10);
   ylabel('Height (km)','FontName','times','FontSize',10);
   titlestr = sprintf('%s/%s/%s %s UTC',...
                 fname(3:4),fname(5:6),fname(1:2),fname(7:8)); 
   title([ '\theta(red), \thetaes(green), \thetav(blue) ' titlestr],'FontName','times','FontSize',11);
   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',10);
   %ylabel('Height (km)','FontName','times','FontSize',16);
   titlestr = sprintf('RH(Dotted), Mixing ratio(Solid), %s/%s/%s %s UTC',...
                 fname(3:4),fname(5:6),fname(1:2),fname(7:8)); 
   title(titlestr,'FontName','times','FontSize',11);
   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',10);
   ylabel('Height (km)','FontName','times','FontSize',10);
   titlestr = sprintf('RICO 05, Wspd, %s/%s/%s %s UTC',...
                 fname(3:4),fname(5:6),fname(1:2),fname(7:8)); 
   title(titlestr,'FontName','times','FontSize',11);
   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',10);
   %ylabel('Height (m)','FontName','times','FontSize',16);
   titlestr = sprintf('RICO 05, Wdir, %s/%s/%s %s UTC',...
                 fname(3:4),fname(5:6),fname(1:2), fname(7:8)); 
   title(titlestr,'FontName','times','FontSize',11);
   %print_buffer = ['-djpeg90 F:\rico\sonde\results\images\RIC02005.sounding.' str2num(aa(w).name(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(day(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

% figure(98)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(v_wnd),[-20 20]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Meridional Wind','FontName','times','FontSize',16);

% figure(99)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(u_wnd),[-20 20]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Zonal Wind','FontName','times','FontSize',16);

% figure(100)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(spd),[0 22]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Wind Speed','FontName','times','FontSize',16);
% 
% figure(101)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(rh),[0 100]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Relative Humidity','FontName','times','FontSize',16);
% 
% figure(102)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(r),[0 22]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% set(gca,'ztick',[0:22]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Mixing Ratio','FontName','times','FontSize',16);
% 
% figure(103)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(t_int),[0 25]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% set(gca,'ztick',[0:22]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Temperature','FontName','times','FontSize',16);
% 
% figure(104)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(theta),[290 320]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% set(gca,'ztick',[0:22]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Theta (K)','FontName','times','FontSize',16);
% 
% figure(105)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(thetav),[295 315]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% set(gca,'ztick',[0:22]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Virtual Potential Temperature (K)','FontName','times','FontSize',16);
% 
% figure(106)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(thetae),[305 350]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% set(gca,'ztick',[0:22]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Equivalent Potential Temperature (K)','FontName','times','FontSize',16);

% figure(107)
% set(gcf,'DefaultAxesFonts',16);
% imagesc(real(thetaes),[325 365]);
% axis xy; colorbar;
% hold on
% set(gca,'xtick',[1:10:w]);
% set(gca,'xticklabel',[ floor(day(1:10:w))]);
% set(gca,'ytick',[51 101 151 201 251 301 351 401]);
% set(gca,'yticklabel',[ h_int(51:50:401)]);
% set(gca,'ztick',[0:22]);
% xlabel('Day','FontName','times','FontSize',16);
% ylabel('Height (m)','FontName','times','FontSize',16);
% title('RICO R/V Seward Johnson Equivalent Saturation Potential
% Temperature (K)','FontName','times','FontSize',16);

end