disp(['read_son_day_', cruise, '_', year]); fclose('all'); clear st1;clear st2;clear st3;clear jd_son; jd=num2str(ddd); if ddd<100 jd=['0' num2str(ddd)]; end; if ddd<10 jd=['00' num2str(ddd)]; end; jhold=0; jax=1; st1=[];st2=[];st3=[]; jaxd=[];jade=[]; for jam=0:23, %cycle thru 24 hourly stats files if jam<10, hr=['0' num2str(jam)]; else hr=num2str(jam); end; e=[way_raw_data_flux, 'day',jd,'\son0' year(3:4) jd hr '_raw.txt']; flist=fopen(e,'r'); if flist>0, %if the file exists, disp(['Reading Sonic file ',e]) % st1 =textscan(flist,'%f %f %f %f %f %*[^\n]','delimiter', ',', 'headerlines', 1, 'emptyvalue', NaN,'treatAsEmpty', 'Sonic'); st1 =textscan(flist,'%2f%2f%2f%3f %*3c %f %f %f %f %f %*3c %*[^\n]','delimiter', ',', 'headerlines', 1, 'emptyvalue', NaN,'treatAsEmpty', 'Sonic'); st1=cell2mat(st1)'; st1(5:9,any(isnan(st1),1))=NaN; %put NaN into any rows containing NaNs (to remove bad lines, gap...), except JD % % despiking used during GOMECC07... % sonraw = st1; % load sonic array into temp array % sonfix % run despiking routine % clear st1 % st1 = sonraw_dep; %copy despiked data into working arrays [nr1,nl1]=size(st1); st2(:,jax:jax+nl1-1)=st1(:,1:nl1); uu=st1(6,:); vv=st1(7,:); ww=st1(8,:); hr=floor(t1/1e7); mn=floor((t1-hr*1e7)/1e5); sc=(t1-hr*1e7-mn*1e5)/1e3; t2=hr*3600+mn*60+sc; jade=ddd + (st1(1,:)+(st1(2,:)+(st1(3,:)+st1(4,:)/1000)/60)/60)/24; jax=jax+nl1; if jam==0 | jam==12 | jam==18, disp(['Calculating sonic spectra for ',hr]) SS=sqrt(uu.^2+vv.^2); %SS=sqrt(UU.^2+VV.^2); [Su,F]=psd2(detrend(SS),length(SS),fsonic);[Sus,Fsu]=specsmoo(Su,fsonic); [Sw,F]=psd2(detrend(ww),length(ww),fsonic);[Sws,Fsw]=specsmoo(Sw,fsonic); figure; loglog(Fsu(length(Fsu)-31:length(Fsu)),((Fsu(length(Fsu)-31:length(Fsu))).^(-2/3))./(Fsu(length(Fsu)-18).^(-2/3)).*Sus(length(Fsu)-18).*Fsu(length(Fsu)-18),'r'); hold on loglog(Fsu,Sus.*Fsu,'b',Fsw,Sws.*Fsw,'g'); title(['Sonic Anemometer Spectra: Day ',jd,' Hour ',hr, ', ',cruise, ' ', year]); legend('5/3','Spd','W' ); xlabel('Frequency in Hz'); ylabel('Spectrum * Frequency'); print('-djpeg90 ',[way_images_flux 'son_spectrum' jd hr '.jpg']); if jam==12, % Print out spectra file for hr 12 specfile=[way_spec_data_flux, 'w_sp',jd,num2str(hr),'.txt']; fspec=fopen(specfile,'w'); specs=[Fsu;Sus;Fsw;Sws]; fprintf(fspec,'%11.5f %11.5f %11.5f %11.5f \r\n',specs); fclose(fspec); clear fspec specs specfile end;%end hhh=12 clear SS ww Su F Sus Fsu Sw Sws Fsw end;%end hhh=0 or 12 or 18 sonfac=3600*fsonic/60; if ~isempty(st2); jk=min(60,ceil(jax/sonfac)); for j=1:jk, jl=sonfac*(j-1)+1;jj=min(jax-1,sonfac*j); st3(:,j+jhold)=median(st2(:,jl:jj)')'; % jd_son(j+jhold)=st2(1,jl); jd_son(j+jhold)=jade(jl); end; %end for stx=st2; jhold=jhold+jk; st2=[]; st1=[]; jaxd=[jaxd jax]; jax=1; end; %if isempty end; %end if flist end; %end for clear u v w Tsonic U V dir u=st3(6,:); %m/s v=st3(7,:); %m/s w=st3(8,:); %m/s % % % Time in flux data files is the number of 100-nano second intervals since % % % Jan 1, 1601 divided by 10000. It was chosen to maintain compatibility % % % with an HP-UNIX version of the code form long ago % % % 0 'dd-mmm-yyyy HH:MM:SS' % % %01-Mar-2000 15:45:17 % % dates = datestr(datenum(1601,1,1) + datenum(jd_son(1,:)*10000*100e-9/86400)); % % hh = str2num(dates(:,13:14)); % % mm = str2num(dates(:,16:17)); % % ss = str2num(dates(:,19:20)); % % clear jd_son % % jd_son= ddd + (hh+((mm+(ss/60))/60))/24; % % % test to see if first point is from the previous day % % if str2num(dates(1,1:2)) ~= str2num(dates(2,1:2)) % % jd_son(1) = jd_son(1)-1; % % end if outputsonic==1 Tsonic=st3(9,:)/100; %Temperature in degC else speedson=st3(9,:)/50; Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15; end; if rotationsonic==1 C=strcmp(sonicmodel,'R3'); if C==1 %means sonicmodel=R3 U=u*cos(30/180*pi)-v*sin(30/180*pi); %to North V=u*sin(30/180*pi)+v*cos(30/180*pi); %to West else U=-u*cos(30/180*pi)-v*sin(30/180*pi); %to North V=-u*sin(30/180*pi)+v*cos(30/180*pi); %to West end; else U=u; V=v; end; S=sqrt(u.*u+v.*v+w.*w); dir=atan2(-V,U)*180/pi+180; %jeff coordinates, +boward, +portward, 180 means wind from dir=mod(dir+720,360); if prtit==1; np=length(st3); prt_jas_son_Stratus_07; end; if plotit==1; %plot T profiles j_st=jd_son(1); j_en=jd_son(np); figure;plot(jd_son(1:np),Tsonic(1:np));xlabel(['JD', year]);ylabel('Temperature (C)');title(['Sonic Temperature. ', cruise, year]); print('-djpeg90 ',[way_images_flux '\Sonic_temp_' num2str(ddd) '.jpg']); figure;plot(jd_son(1:np),U(1:np),jd_son(1:np),V(1:np),jd_son(1:np),w(1:np),jd_son(1:np),S(1:np));xlabel(['JD', year]);ylabel('Sonic Winds (m/s)');title(['Relative Wind Components: U=blue, V=green, W=red, Speed=cyan. ', cruise, year]); print('-djpeg90 ',[way_images_flux '\Sonic_winds_' num2str(ddd) '.jpg']); for i=1:np % Jeff Otten's mod to display direction -180 to +180. if dir(i) > 180 dir(i) = dir(i) - 360; end end figure;plot(jd_son(1:np),dir(1:np),'b', ... %data points [ddd+4/24 ddd+4/24],[-180 180],'k', ... %4 hour watch dividers [ddd+8/24 ddd+8/24],[-180 180],'k', ... [ddd+12/24 ddd+12/24],[-180 180],'k', ... [ddd+16/24 ddd+16/24],[-180 180],'k', ... [ddd+20/24 ddd+20/24],[-180 180],'k', ... [ddd ddd+1],[45 45],'g', ... %45 degree good zone [ddd ddd+1],[-45 -45],'g', ... [ddd ddd+1],[60 60],'r', ... %60 degree questionable zone [ddd ddd+1],[-60 -60],'r'); ... xlabel(['JD', year]); ... ylabel('Sonic Wind Direction (deg)'); ... title(['Relative Wind Direction. ', cruise, year]); ... axis([j_st j_en -180 180]); print('-djpeg90 ',[way_images_flux '\Sonic_rel_wnd_dir_' num2str(ddd) '.jpg']); figure;plot(jd_son(1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));xlabel(['JD', year]);ylabel('Wind Speed m/s');title(['Sonic Relative Wind Speed. ', cruise, year]); print('-djpeg90 ',[way_images_flux '\Sonic_rel_wnd_speed_' num2str(ddd) '.jpg']); end;