%%%% READ SOUNDING DATA FROM LATEST PSD LATLON format FORMAT...
%%%% COMPOSITE TIME-HEIGHT DISPLAY ....
%%% pk
display('sounding_composite')
%% DEW 9/20/2007
clear

read_parameters_Stratus
prt = 1;  % print figures to jpg files
warning off MATLAB:divideByZero  % turn off warning message

pathdata = way_raw_data_balloon;
outpath  = way_images_balloon;
fname = '*Lon.txt';
s1=[way_raw_data_balloon fname];
aa = dir(s1);
[m n] = size(aa);
p = length(aa);
[num] = num2str(p);
disp('Number of rawinsonde *LON.txt files : '); disp( [num]);

%%%%%%%%%% starting file %%%%%%%%%%%%%%%%%%%%%%%%
st=1;
%%%%%%%%%%%%%%%%%% last file == k %%%%%%%%%%%%%%%
ws=1;
we=m;
for w = ws:we  %m;
   w;
   [numf] = num2str(w);

   disp('# of WTEC*.txt file to process :'); disp ( [numf]);
   sounding_filename = aa(w).name
   %WTEC_200705060111_edt.txt
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %% READ THE AMMA 2007 EDT FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   %   0705060111_EDT_LATLON.txt
   yr(w) = str2num(aa(w).name(1:2))+2000;
   month(w) = str2num(aa(w).name(3:4));
   day (w)  = str2num(aa(w).name(5:6));
   hour(w)  = str2num(aa(w).name(7:8));
   minutes(w) = str2num(aa(w).name(9:10));


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

   if fid1 == -1
       error(['Cannot read ' sounding_filename])
   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)
       %for k = 1:30  %nhts
       %dout=fgetl(fid1); % read in extra space between each data line
       spout=fgetl(fid1);
       %if length(spout) > 1
       [data] =sscanf(spout,'%f');
       %end
       ss = size(data);
       if ss(1) == 11
           k = k + 1;
           H(k,w)     = data(2); %Ht m
           P(k,w) = data(3);  %Press mb
           T(k,w)     = data(4);  % T C
           Td(k,w)    = data(5);  %dew Pt C
           RH(k,w)    = data(6);  % RH %
           spd(k,w)=data(7);  % wnd spd m/s
           dir(k,w)=data(8);  % wnd dir true
           lat(k,w)=data(9); % tenths of degrees N = + S = -
           lon(k,w)=data(10);  % tenths of degrees W = - E = +
       end
       %%%%%%%%
   end
   k;
end  % Lopp through all the files loading the arrays
fclose(fid1);


%%%% read the sounding, now QC and interpolation %%%%%%
%%%% missing value is -32768
%%%%% Interpolate the sounding....from surface to 6 km %%%%%%%%%%%%%%%%%%%
%%% FIRST QUALITY CONTROL....

h_int = [0:25:25000];

for i=ws:we  % loop through all the profiles
   g = find (H(:,i)>0 & H(:,i)<=25000);
   gk = find(diff(H(g,i))>0);  % Check for increasing hts
   c = g(gk);clear g; g=c;clear c gk;
   f=find(T(g,i)>-273);


   if isempty(f)==0
       t_int(1:length(h_int),i)     =  interp1(H(g(f),i),T(g(f),i),h_int)'; clear f;
   end
   f=find(RH(g,i)>0);
   if isempty(f)==0
       rh_int(1:length(h_int),i)     =  interp1(H(g(f),i),RH(g(f),i),h_int)'; clear f;
   end
   f=find(P(g,i)>0);
   if isempty(f)==0
       press_int(1:length(h_int),i)     =  interp1(H(g(f),i),P(g(f),i),h_int)'; clear f;
   end
   f=find(dir(g,i)>0);
   if isempty(f)==0
       dir_int(1:length(h_int),i)     =  interp1(H(g(f),i),dir(g(f),i),h_int)'; clear f;
   end
   f=find(spd(g,i)>0);
   if isempty(f)==0
       spd_int(1:length(h_int),i)     =  interp1(H(g(f),i),spd(g(f),i),h_int)'; clear f;
   end
   % Put surface values in 0 ht interpolated arrays
   t_int(1,i) = T(1,i);
   rh_int(1,i) = RH(1,i);
   press_int(1,i) = P(1,i);
   spd_int(1,i) = spd(1,i);
   dir_int(1,i) = dir(1,i);


end
% Initialize RH profiles from 2004/12/16 1100 and 1152 flights
% move 1100 sounding into the bkup folder not used
%rh_int(:,48) = nan;
%rh_int(:,49) = nan;

save sounding_all press_int t_int rh_int spd_int dir_int h_int yr month day hour minutes T;


%%%%%% SOUNDINGS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all

read_parameters_Stratus
prt = 1; % print jpg of soundings

load sounding_all
t=t_int;%-273.15;
rh=rh_int;
press=press_int;
h_int = [0:25:25000];
h=h_int;
k=length(day);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 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);


%%%%%%%% SOUNDING QC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% LCL Calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



for w = 1:k;

   i_lev = 8;

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

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%% COMMON TIME COORDINATE %%%%%

for q=1:w
   jul_day_sonde(q) = datenum(yr(q),month(q),day(q),hour(q),minutes(q),00)-datenum(2006,12,31,00,00,00);
end;

save sound_resultsall

graphs = -1;
% pause
ymax = 25;  % used to scale images km
%%% Year Day 188 = JULY 6!!
figure(1)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(theta),[280,340]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+6)  = 'theta';
title(tit,'FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_theta.jpg'];
   print('-djpeg90 ', print_buffer);
end
figure(2)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(thetav),[280 340]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+8)  = 'theta-v';
title(tit,'FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_thetav.jpg'];
   print('-djpeg90 ', print_buffer);
end
figure(3)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(rh_int),[0 100]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+3)  = 'RH';
title(tit,'FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_RH.jpg'];
   print('-djpeg90 ', print_buffer);
end

figure(4)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(r),[0 15]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+13)  = 'Mixing-Ratio';
title(tit,'FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_MR.jpg'];
   print('-djpeg90 ', print_buffer);
end

figure(5)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(dir_int),[0 360]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+10)  = 'Direction';
title(tit,'FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_DIR.jpg'];
   print('-djpeg90 ', print_buffer);
end

figure(6)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(spd_int),[0 30]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+6)  = 'Speed';
title(tit,'FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_SPD.jpg'];
   print('-djpeg90 ', print_buffer);
end
figure(7)
set(gcf,'DefaultAxesFonts',16);
imagesc(jul_day_sonde,h_int/1000,real(t_int),[-90 25]);
axis xy; ylim([0 ymax]);
hc = colorbar;
xlabel('Day','FontName','times','FontSize',18);
ylabel('Ht (km)','FontName','times','FontSize',18);
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+12)  = 'Temperature';
title('AMMA 2007, Temperature (C)','FontName','times','FontSize',18);
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_T.jpg'];
   print('-djpeg90 ', print_buffer);
end

% find mean profile for all soundings

for z = 1:(length(h_int))
   j = find(~isnan(t_int(z,:)));
   p_intbar(z,1) = mean(t_int(z,j));
   j = find(~isnan(theta(z,:)));
   p_intbar(z,2) = mean(theta(z,j));
   j = find(~isnan(rh(z,:)));
   p_intbar(z,3) = mean(rh(z,j));
   j = find(~isnan(t_int(z,:)));
   p_intbar(z,4) = mean(spd_int(z,j));
end
figure(8)
subplot(2,2,1),plot(p_intbar(:,1),h_int/1000)
title('Mean Temp')
ylabel('Ht (km)')
subplot(2,2,2),plot(p_intbar(:,2),h_int/1000)
title('Mean Theta')
ylabel('Ht (km)')
subplot(2,2,3),plot(p_intbar(:,3),h_int/1000)
title('Mean RH')
ylabel('Ht (km)')
subplot(2,2,4),plot(p_intbar(:,4),h_int/1000)
title('Mean Speed')
ylabel('Ht (km)')
if prt
   print_buffer = [way_images_balloon '_' cruise '_' year '_Mean_profiles.jpg'];
   print('-djpeg90 ', print_buffer);
end
save meanp p_intbar


% PWV calculations
%%%%%%%%%%%%%%%%%%%%%%

display(['pwv_sonde_07'])
% DEW 5/12/2007

clear all;
read_parameters_Stratus
load sounding_all.mat

%clear lat* lon* spd_int dir_int

jd=datenum(yr,month,day,hour,minutes,0)-datenum(2006,12,31);

clear yr month day hour minute

h=h_int;clear h_int
press=press_int;clear press_int;
t=t_int;clear t_int
rh=rh_int;clear rh_int



e=6.112.*exp(17.67.*t./(t+243.5)).*rh./100; %vapor pressure pa
w=e.*622./(press-0.378*e); %mixing ratio g/kg
es=6.112.*exp(17.67.*t./(t+243.5));%pa


%now, adjust for vapor pressure over ice to determine where get ice clouds.
ice=find(t<=0);
es(ice)=(1.0003+4.18e-6*press(ice)).*6.1115.*exp(22.452.*t(ice)./(t(ice)+272.55));


wm=0.5*(w(1:end-1,:)+w(2:end,:));%take each layer and avg mixing ratio

wm=0.001*wm;%convert from g/kg to g/g to be dimensionless

ql=wm./(wm+1);% change to specific humidity

dpp=1000.*(press(1:end-1,:)-press(2:end,:));
%find pressure difference for each layer in dynes/cm^2, starting from mb

g=980.616;%accln due to gravity cm/s2

pw=(ql./g).*dpp; %% Precipitable water cm

pw=nansum(pw); %cm

RHcorr=rh;
RHcorr(ice)=e(ice)./es(ice).*100;

figure()
set(gcf,'defaultaxesfonts',18)
set(gca,'fontname','time','fontsize',18);
plot(jd,pw,'r')
hold
plot(jd,pw,'r*')
ylim([0 8])
xlabel('Year Day (UTC)')
ylabel('Sonde PWV (cm)')
clear tit
tit = cruise;
lt = length(cruise);
tit(lt+1) = '-';
tit((lt+2):lt+5) = year; % This year is read in from the parameter file
ltt = length(tit);
tit(ltt+1) = '-';
tit(ltt+2:ltt+10)  = 'PWV-Sonde';
title(tit)
print_buffer = [way_images_balloon '_' cruise '_' year '_PWV_sonde.jpg'];
print('-djpeg90 ', print_buffer);

% save data to be read in by AMMA07_allcruise_mb with a combined pwv plot
pwv(:,1) = jd;
pwv(:,2) = pw;
s1=[way_proc_data_balloon];
   flist1=fopen([s1 cruise '_' year '_sonde_pwv.txt'],'w');
   fprintf(flist1,'%12.6f %10.5f \r',pwv');
   fclose(flist1);
