%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%SOUNDING READER %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%READS EDT %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display('sounding_plot')
clear all;
close all;
warning off
read_parameters_rhb
pathdata = way_raw_data_balloon %'F:\AMMA_2007\RHB\balloon\Raw\';
%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 %%%%%%%%%%%%%%%%%%%
sy = size(yr);
h_int = [0:10:25000];
for i=1:sy(1,2)
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 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,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 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
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:lt+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 '_' cruise '_' year '_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');