disp('HiWinGS_2013_hr_ave');

dta = 60; % averaging time in minutes
clear jda y1 yvar1 nn1 y2 yvar2 nn2 y3 yj jym xx yu yts yqs yq;
clear v1 theta1 v2_scs thet1_scs v1rel thet1rel yushp yhead v1x thet1x jdam;

%% 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_best<rdir_hi & reldir_best>rdir_lo) & sig_h<sigh_lim & sig_u<sigu_lim & ...
%             sp2_best<sp2_lim & rain_best<rain_lim & (sqrt(ww_best)./ugw<2+.0015*u_best.^2) & ...
%             (sqrt(vv_best)./ugu<2) & isfinite(usib) & badSon<badSon_lim);
kkk = find((reldir_best<rdir_hi & reldir_best>rdir_lo) & sig_h<sigh_lim & sig_u<sigu_lim & ...
            sp2_best<sp2_lim);
zzz = setdiff(1:length(jdy), kkk); % indices for data NOT passing criteria

% ship plume index from wind variances, 0 == good;
jplume = ones(length(jdy),1);
ii = find(sqrt(ww_best)./ugw<2+.0015*u_best.^2 & sqrt(vv_best)./ugu<2);
jplume(ii) = 0;
% ship maneuver index from stdev of ship speed, <1 is good...
jmanuv = 3*ones(length(jdy),1);
ii = find(sig_h<sigh_lim & sig_u<sigu_lim & sp2<sp2_lim);
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
y1(:,116:128) = y2(:,116:128);  % cov and var, PSD
y1(:,134:149) = y2(:,134:149);  % cx^2 & usib, tsib, qsib, lid, PSD
y1(:,185:191) = y2(:,185:191);  % cov and var, UH
y1(:,195:206) = y2(:,195:206);  % cx^2 & usib, tsib, qsib, lid, UH
y3(:,2:16) = ycr(:,2:16);       % .cr files

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

%% humidity spectrum
clear ycq;
% choose filtered structure function parameter cql from licor
switch IDflag
    case 'a'
%         [jdar,ycqfl,ycqvarl,nnql] = timave_1(jdy(kkk),dta,xx(kkk,140));
        ycqfl = y2(:,140);
    case 'b'
%         [jdar,ycqfl,ycqvarl,nnql] = timave_1(jdy(kkk),dta,xx(kkk,141));
        ycqfl = y2(:,141);
end

ii = ycqfl<1e-4 | ycqfl>1;   ycqfl(ii) = NaN;

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

% ustar and Cd10N
[~,yusr_psd,~,~] = timave_1(jdy,dta,usr_psd);
[~,yusr_uh,~,~] = timave_1(jdy,dta,usr_uh);
[~,yusr_id_psd,~,~] = timave_1(jdy,dta,usr_id_psd);
[~,yusr_id_uh,~,~] = timave_1(jdy,dta,usr_id_uh);
[~,ycd10n_id_psd,~,~] = timave_1(jdy,dta,cd10n_id_psd);
[~,ycd10n_id_uh,~,~] = timave_1(jdy,dta,cd10n_id_uh);

%% vector average winds, uship, heading, sog, cog
clear v1 theta1 v2_scs thet1_scs v1rel thet1rel yushp yhead v1x thet1x jdam
[jda,v1,thet1] = timave_vec1(jdy,dta,u_best,dir_best); v1 = v1'; thet1 = thet1';                        % best 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,urel_best,reldir_best); 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 from latest C35 run
[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_best);

%% licor hourly mean covariance hl
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,hlc_licbar,hlc_licvar,nnlic] = timave_1(jdy,dta,hlc_licx);
end
