% read_son_day % VOCALS 2008 :: 2008-10-16 :: Simon de Szoeke % Read sonic anemometer data stream. disp(['read_son_day: ', cruise, ' ' yyyy]); fclose('all'); clear st1;clear st2;clear st3; clear jd_son; jd=sprintf('%03i',ddd); jhold=0; jax=1; st1=[];st2=[];st3=[]; jaxd=[];jade=[];temp=[]; for jam=0:23, %cycle thru 24 hourly stats files hr=sprintf('%02i',jam); e=[way_raw_data_flux, 'day',jd,'\son0' yyyy(3:4) jd hr '_raw.txt']; flist=fopen(e,'r'); st1=[];temp=[]; if flist>0, %if the file exists, disp(['Reading Sonic file ',e]) k=0; %SPdeS keep count while ~feof(flist) % end-of-file not reached k=k+1; try temp=textscan(flist,'%2f%2f%3f %*4c %f %f %f %f %f %*3c %*[^\n]','delimiter', ',', 'headerlines', 1, 'emptyvalue', NaN,'treatAsEmpty', 'Sonic'); st1=[st1 cell2mat(temp)']; catch % if one line is cut short, blank it out totally [l8,c8]=size(temp{1,8}); if l8==0 || c8==0; continue; end % read next batch of lines for ll=1:8 if length(temp{1,ll})~=l8 temp{1,ll}(l8+1)=[]; end; end; if isempty(temp{1,1}); temp{1,9}=temp{1,9}'; end; st1=[st1 cell2mat(temp)']; end; %end try/catch end; [nr1,nl1]=size(st1); st2(:,jax:jax+nl1-1)=st1(:,1:nl1); ii=find(~isnan(st1(1,:))); uu=st1(5,ii); vv=st1(6,ii); ww=st1(7,ii); jade=ddd+(jam+(st1(1,:)+(st1(2,:)+st1(3,:)/1000)/60)/60)/24; % time in decimal yearday 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); if plotit; 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, ' ' yyyy]); legend('5/3','Spd','W' ); xlabel('Frequency in Hz'); ylabel('Spectrum * Frequency'); if prtit; print('-djpeg90 ',[way_images_flux 'son_spectrum' jd hr '.jpg']); end; end 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 if hour==0 or 12 or 18 jax=jax+nl1; sonfac=3600*fsonic/60; if ~isempty(st2); jk=min(60,floor(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)=jade(jl); end; stx=st2; jhold=jhold+jk; st2=[]; st1=[]; jaxd=[jaxd jax]; jax=1; end; %if st2 isempty fclose(flist); % SPdeS end; %end if flist end; %end for jam=0:23 % cycle through hour clear u v w Tsonic U V dir u=st3(5,:); %m/s v=st3(6,:); %m/s w=st3(7,:); %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(8,:); %Temperature in degC else speedson=st3(8,:)/50; Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15; end; if rotationsonic switch sonicmodel case 'R3' %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 case 'VOCALS' U= u*cos(90+30/180*pi)-v*sin(90+30/180*pi); V= u*sin(90+30/180*pi)+v*cos(90+30/180*pi); otherwise 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); np=length(st3); if saveit; prt_jas_son; end; if plotit; %plot T profiles j_st=jd_son(1); j_en=jd_son(np); figure; plot(jd_son(1:np),Tsonic(1:np));xlabel(['JD', yyyy]);ylabel('Temperature (C)');title(['Sonic Temperature. ', cruise, yyyy]); if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_temp_' num2str(ddd) '.jpg']); end; 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 ' yyyy]);ylabel('Sonic Winds (m/s)');title(['Relative Wind Components: U=blue, V=green, W=red, Speed=cyan. ', cruise, yyyy]); if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_winds_' num2str(ddd) '.jpg']); end; % 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 dir(dir>180)=dir(dir>180)-360; 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 ' yyyy]); ylabel('Sonic Wind Direction (deg)'); title(['Relative Wind Direction. ', cruise, yyyy]); axis([j_st j_en -180 180]); if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_rel_wnd_dir_' num2str(ddd) '.jpg']); end; figure; plot(jd_son(1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));xlabel(['JD ' yyyy]);ylabel('Wind Speed m/s');title(['Sonic Relative Wind Speed. ', cruise, yyyy]); if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_rel_wnd_speed_' num2str(ddd) '.jpg']); end; end; % if plotit