disp('DYNAMO_2011_hr_ave');

dta = 60; % averaging time in minutes
clear jda y1 yvar1 nn1 y2 yvar2 nn2 y3;

%% straight average of all data
[jda,y1,yvar1,nn1] = timave_1(jdy,dta,x);
[jda,y3,yvar3,nn1] = timave_1(jdy,dta,yxx);

%% straight average of 'best' bulk met data
best = [jdy,u_best,ts_best,qs_best,ta_best,rh_best,qa_best,rs_best,rl_best,P_best,rain_best];
%        1    2      3       4       5       6       7       8       9      10        11
[jda,y4,yvar4,nn1] = timave_1(jdy,dta,best);

%% selected turbulence averages and data quality parameters
% apply selection criteria
% kkk=find((reldir<60 & reldir>-60) & ushp_im<5 & sig_h<4 & sig_u<0.2 & sp2<0.8  & sqrt(ww)./ugw<.7+.0015*U.^2 & sqrt(vv)./ugu<1.4 & isfinite(usib));% & isledist>3.5
% kkk=find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<1.4 & rain_best<2 & sp2<1.2  & (sqrt(wwj)./ugw<2+.0015*U.^2) & (sqrt(vvj)./ugu<2) & isfinite(usib) & missing<300 & badSon<300);% & isledist>3.5
kkk = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<0.6 & sp2<0.8 & rain_best<5 & (sqrt(wwj)./ugw<2+.0015*u_best.^2) & (sqrt(vvj)./ugu<2) & isfinite(usib) & badSon<50);
zzz = setdiff(1:length(jdy), kkk); % indices for data NOT passing criteria

% ship plume index from wind variances
jplume = ones(length(jdy),1);
ii = find(sqrt(wwj)./ugw<2+.0015*U.^2 & sqrt(vvj)./ugu<2);
jplume(ii) = 0;
% ship maneuver index from stdev of ship speed
jmanuv = 3*ones(length(jdy),1);
ii = find(sig_h<5 & sig_u<0.6 & sp2<0.8);
jmanuv(ii) = sig_u(ii);

% average values which pass criteria
clear y2 yj jym xx
xx = x; xx(zzz,2:end) = NaN;  % remove data not passing criteria
jplume(zzz) = NaN; jmanuv(zzz) = NaN;
yxxx = yxx; yxxx(zzz,2:end) = NaN;
[jdax,y2,yvar2,nn2] = timave_1(jdy,dta,xx);
[jdax,yjx,yjvar,nn2] = timave_1(jdy,dta,jplume);
[jdax,yjmx,yjmvar,nn2] = timave_1(jdy,dta,jmanuv);
[jdax,ycr,ycrvar,nn2] = timave_1(jdy,dta,yxxx);

%% identify turbulence variables and insert filtered result into y1/y3
yj = zeros(1,length(jda));
yjm = zeros(1,length(jda));

y1(:,116:128) = y2(:,116:128);  % cov and var
y1(:,134:141) = y2(:,134:141);  % cx^2
y3(:,2:16) = ycr(:,2:16);       % .cr files

y1(:,44) = y2(:,44);     % tilt
yj = yjx(:);             % plume
yj(isnan(yj)) = 1;
yjm = yjmx(:);           % maneuver
yjm(isnan(yjm)) = 3;

%% humidity spectrum
clear ycq;
% average filtered structure function parameter cql from licor
[jdar,ycqfl,ycqvarl,nnql] = timave_1(jdy(kkk),dta,xx(kkk,140));

%% Other averages
clear yu yts yqs yq
[jdaa,yu,yvaru,nnu] = timave_1(jdy,dta,u10n);
[jdaa,yts,tsvar,nnts] = timave_1(jdy,dta,dt(:,4));
[jdaa,yqs,qsvar,nnts] = timave_1(jdy,dta,dt(:,6));  % average qs
yq = qair_p([y4(:,5),y4(:,6),y4(:,10)]);

%% vector average winds
clear v1 theta1 v2_scs thet1_scs v1rel thet1rel yushp yhead v1x thet1x jdam
[jda,v1,thet1] = timave_vec1(jdy,dta,x(:,16),x(:,41)); v1 = v1'; thet1 = thet1';                        % psd true
[jda,v1_scs,thet1_scs] = timave_vec1(jdy,dta,x(:,2),x(:,34)); v1_scs = v1_scs'; thet1_scs = thet1_scs'; % scs true
[jda,v1rel,thet1rel] = timave_vec1(jdy,dta,x(:,38),x(:,42)); v1rel = v1rel'; thet1rel = thet1rel';      % relative
[jdam,v1x,thet1x] = timave_vec1(jdy,dta,xx(:,16),xx(:,41)); v1x = v1x'; thet1x = thet1x';               % psd turb select
[jdax,yushp,yhead] = timave_vec1(jdy,dta,ushp,head); yushp = yushp'; yhead = yhead';                    % ship dopler log, RL gyro
[jdax,ysog,ycog] = timave_vec1(jdy,dta,sog_im,cog_im); ysog = ysog'; ycog = ycog';                      % scs sog and cog

%August 4, 2008 - minor corrections to direction averaging.
i = find(thet1<0);     thet1(i) = thet1(i)+360;
i = find(thet1_scs<0); thet1_scs(i) = thet1_scs(i)+360;
i = find(thet1rel<-180);  thet1rel(i) = thet1rel(i)+360;
i = find(thet1rel>180);  thet1rel(i) = thet1rel(i)-360;
i = find(ycog<0);      ycog(i) = ycog(i)+360;
i = find(yhead<0);     yhead(i) = yhead(i)+360;

v1f = zeros(1,length(jda))';
ii = find(jda==jdam(1));
v1f(ii:length(jda)) = v1x;
if ii>1
   v1f(1:ii-1) = v1(1:ii-1);
end;

%% ave 10-min bulk fluxes
[jda,hsbar,hsvar,nnhs] = timave_1(jdy,dta,hsbb);
[jda,hlbar,hlvar,nnhl] = timave_1(jdy,dta,hlbb);
[jda,taubar,tauvar,nntau] = timave_1(jdy,dta,taubb);
[jda,dtwrmbar,dtwrmvar,nndtwrm] = timave_1(jdy,dta,dt(:,34));
[jda,hl_webar,hl_webvar,nnhl_webb] = timave_1(jdy,dta,hl_webb);
[jda,rhbar,rhvar,nnhs] = timave_1(jdy,dta,rh);

%% licor
hlc_licbar = zeros(1,length(jda));
if isempty(hlc_lic)
   hlc_licbar(:) = NaN;
else
    hlc_licx = hlc_lic; hlc_licx(zzz) = NaN; % remove bad values
  	[jdax,hlx,hlc_licvar,nnlic] = timave_1(jdy,dta,hlc_licx);
	hlc_licbar= hlx; % cov and var
end;
