clear;

load /sirius2/zhup/wrk/tower_data/Everglads/out_10m_dat_1.mat;
% starting time 2012-09-28, 12:00 LST
% ending time   2012-09-30, 04:30 LST
% 40.5 hours
time1=[0:(40.5*4+1)*18000-1]/3600/20/24+28.5;

%remove spikes
for n=1:163
var=squeeze(t(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
t(:,n)=var;
var=squeeze(q(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
q(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
u(:,n)=var;
var=squeeze(v(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
v(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
w(:,n)=var;
end

p1=reshape(p,1,[]);
t1=reshape(t,1,[]);
q1=reshape(q,1,[]);
u1=reshape(u,1,[]);
v1=reshape(v,1,[]);
w1=reshape(w,1,[]);

pp(1:12*3600*20,1)=nan;
tt(1:12*3600*20,1)=nan;
qq(1:12*3600*20,1)=nan;
uu(1:12*3600*20,1)=nan;
vv(1:12*3600*20,1)=nan;
ww(1:12*3600*20,1)=nan;

pp(12*3600*20+1:24*3600*20,1)=p1(1:12*3600*20);
tt(12*3600*20+1:24*3600*20,1)=t1(1:12*3600*20);
qq(12*3600*20+1:24*3600*20,1)=q1(1:12*3600*20);
uu(12*3600*20+1:24*3600*20,1)=u1(1:12*3600*20);
vv(12*3600*20+1:24*3600*20,1)=v1(1:12*3600*20);
ww(12*3600*20+1:24*3600*20,1)=w1(1:12*3600*20);

pp(1:24*3600*20,2)=p1(12*3600*20+1:36*3600*20);
qq(1:24*3600*20,2)=q1(12*3600*20+1:36*3600*20);
tt(1:24*3600*20,2)=t1(12*3600*20+1:36*3600*20);
uu(1:24*3600*20,2)=u1(12*3600*20+1:36*3600*20);
vv(1:24*3600*20,2)=v1(12*3600*20+1:36*3600*20);
ww(1:24*3600*20,2)=w1(12*3600*20+1:36*3600*20);

pp(1:24*3600*20,3)=nan;
tt(1:24*3600*20,3)=nan;
qq(1:24*3600*20,3)=nan;
uu(1:24*3600*20,3)=nan;
vv(1:24*3600*20,3)=nan;
ww(1:24*3600*20,3)=nan;

pp(1:4.75*3600*20,3)=p1(36*3600*20+1:end);
tt(1:4.75*3600*20,3)=t1(36*3600*20+1:end);
qq(1:4.75*3600*20,3)=q1(36*3600*20+1:end);
uu(1:4.75*3600*20,3)=u1(36*3600*20+1:end);
vv(1:4.75*3600*20,3)=v1(36*3600*20+1:end);
ww(1:4.75*3600*20,3)=w1(36*3600*20+1:end);

load /sirius2/zhup/wrk/tower_data/Everglads/out_10m_dat_2.mat;
% starting time 2012-09-30, 10:30 LST
% ending time   2012-09-30, 18:45 LST
% 8.25 hours
time2=[0:(8.25*4+1)*18000-1]/3600/20/24+30.4375;

%remove spikes
for n=1:34
var=squeeze(t(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
t(:,n)=var;
var=squeeze(q(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
q(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
u(:,n)=var;
var=squeeze(v(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
v(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
w(:,n)=var;
end

p2=reshape(p,1,[]);
t2=reshape(t,1,[]);
q2=reshape(q,1,[]);
u2=reshape(u,1,[]);
v2=reshape(v,1,[]);
w2=reshape(w,1,[]);

pp(1:24*3600*20,4)=nan;
tt(1:24*3600*20,4)=nan;
qq(1:24*3600*20,4)=nan;
uu(1:24*3600*20,4)=nan;
vv(1:24*3600*20,4)=nan;
ww(1:24*3600*20,4)=nan;

pp(10.5*3600*20+1:19*3600*20,4)=p2;
tt(10.5*3600*20+1:19*3600*20,4)=t2;
qq(10.5*3600*20+1:19*3600*20,4)=q2;
uu(10.5*3600*20+1:19*3600*20,4)=u2;
vv(10.5*3600*20+1:19*3600*20,4)=v2;
ww(10.5*3600*20+1:19*3600*20,4)=w2;

load /sirius2/zhup/wrk/tower_data/Everglads/out_10m_dat_3.mat;
% starting time 2012-10-01, 09:00 LST
% ending time   2012-10-01, 18:15 LST
% 9.25 hours
time3=[0:(9.25*4+1)*18000-1]/3600/20/24+31.375;

%remove spikes
for n=1:38
var=squeeze(t(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
t(:,n)=var;
var=squeeze(q(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
q(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
u(:,n)=var;
var=squeeze(v(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
v(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
w(:,n)=var;
end

p3=reshape(p,1,[]);
t3=reshape(t,1,[]);
q3=reshape(q,1,[]);
u3=reshape(u,1,[]);
v3=reshape(v,1,[]);
w3=reshape(w,1,[]);

pp(1:24*3600*20,5)=nan;
tt(1:24*3600*20,5)=nan;
qq(1:24*3600*20,5)=nan;
uu(1:24*3600*20,5)=nan;
vv(1:24*3600*20,5)=nan;
ww(1:24*3600*20,5)=nan;

pp(9*3600*20+1:18.5*3600*20,5)=p3;
tt(9*3600*20+1:18.5*3600*20,5)=t3;
qq(9*3600*20+1:18.5*3600*20,5)=q3;
uu(9*3600*20+1:18.5*3600*20,5)=u3;
vv(9*3600*20+1:18.5*3600*20,5)=v3;
ww(9*3600*20+1:18.5*3600*20,5)=w3;

load /sirius2/zhup/wrk/tower_data/Everglads/out_10m_dat_4.mat;
% starting time 2012-10-02, 09:00 LST
% ending time   2012-10-02, 11:45 LST
% 2.75 hours
time4=[0:(2.75*4+1)*18000-1]/3600/20/24+32.375;

%remove spikes
for n=1:12
var=squeeze(t(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
t(:,n)=var;
var=squeeze(q(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
q(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
u(:,n)=var;
var=squeeze(v(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
v(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
w(:,n)=var;
end

p4=reshape(p(:,1:12),1,[]);
t4=reshape(t(:,1:12),1,[]);
q4=reshape(q(:,1:12),1,[]);
u4=reshape(u(:,1:12),1,[]);
v4=reshape(v(:,1:12),1,[]);
w4=reshape(w(:,1:12),1,[]);

% starting time 2012-10-03, 09:00 LST
% ending time   2012-10-03, 20:30 LST
% 11.5 hours
time5=[0:(11.5*4+1)*18000-1]/3600/20/24+33.375;

%remove spikes
for n=13:59
var=squeeze(t(:,n));
if n==24 
var(abs(var-nanmean(var))>0.5*nanstd(var))=nan;
else
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
end
t(:,n)=var;
var=squeeze(q(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
q(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
u(:,n)=var;
var=squeeze(v(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
v(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
w(:,n)=var;
end

p5=reshape(p(:,13:end),1,[]);
t5=reshape(t(:,13:end),1,[]);
q5=reshape(q(:,13:end),1,[]);
u5=reshape(u(:,13:end),1,[]);
v5=reshape(v(:,13:end),1,[]);
w5=reshape(w(:,13:end),1,[]);

pp(1:24*3600*20,6)=nan;
tt(1:24*3600*20,6)=nan;
qq(1:24*3600*20,6)=nan;
uu(1:24*3600*20,6)=nan;
vv(1:24*3600*20,6)=nan;
ww(1:24*3600*20,6)=nan;

pp(9*3600*20+1:12*3600*20,6)=p4;
tt(9*3600*20+1:12*3600*20,6)=t4;
qq(9*3600*20+1:12*3600*20,6)=q4;
uu(9*3600*20+1:12*3600*20,6)=u4;
vv(9*3600*20+1:12*3600*20,6)=v4;
ww(9*3600*20+1:12*3600*20,6)=w4;

pp(1:24*3600*20,7)=nan;
tt(1:24*3600*20,7)=nan;
qq(1:24*3600*20,7)=nan;
uu(1:24*3600*20,7)=nan;
vv(1:24*3600*20,7)=nan;
ww(1:24*3600*20,7)=nan;

pp(9*3600*20+1:20.75*3600*20,7)=p5;
tt(9*3600*20+1:20.75*3600*20,7)=t5;
qq(9*3600*20+1:20.75*3600*20,7)=q5;
uu(9*3600*20+1:20.75*3600*20,7)=u5;
vv(9*3600*20+1:20.75*3600*20,7)=v5;
ww(9*3600*20+1:20.75*3600*20,7)=w5;

load /sirius2/zhup/wrk/tower_data/Everglads/out_10m_dat_5.mat;
% starting time 2012-10-05, 09:00 LST
% ending time   2012-10-05, 21:30 LST
% 12.5 hours
time6=[0:(12.5*4+1)*18000-1]/3600/20/24+35.375;

%remove spikes
for n=1:51
var=squeeze(t(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
t(:,n)=var;
var=squeeze(q(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
q(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
u(:,n)=var;
var=squeeze(v(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
v(:,n)=var;
var=squeeze(u(:,n));
var(abs(var-nanmean(var))>3.0*nanstd(var))=nan;
w(:,n)=var;
end

p6=reshape(p,1,[]);
t6=reshape(t,1,[]);
q6=reshape(q,1,[]);
u6=reshape(u,1,[]);
v6=reshape(v,1,[]);
w6=reshape(w,1,[]);

pp(1:24*3600*20,8)=nan;
tt(1:24*3600*20,8)=nan;
qq(1:24*3600*20,8)=nan;
uu(1:24*3600*20,8)=nan;
vv(1:24*3600*20,8)=nan;
ww(1:24*3600*20,8)=nan;

pp(9*3600*20+1:21.75*3600*20,8)=p6;
tt(9*3600*20+1:21.75*3600*20,8)=t6;
qq(9*3600*20+1:21.75*3600*20,8)=q6;
uu(9*3600*20+1:21.75*3600*20,8)=u6;
vv(9*3600*20+1:21.75*3600*20,8)=v6;
ww(9*3600*20+1:21.75*3600*20,8)=w6;

time=[1:24*3600*20]/3600/20;

N=5;

time_N=reshape(time,N*60*20,[]);
ppp=reshape(pp,N*60*20,[],8);
ttt=reshape(tt,N*60*20,[],8);
qqq=reshape(qq,N*60*20,[],8);
uuu=reshape(uu,N*60*20,[],8);
vvv=reshape(vv,N*60*20,[],8);
www=reshape(ww,N*60*20,[],8);

m_p=squeeze(nanmean(ppp,1));
m_t=squeeze(nanmean(ttt,1));
m_q=squeeze(nanmean(qqq,1));
m_u=squeeze(nanmean(uuu,1));
m_v=squeeze(nanmean(vvv,1));
m_w=squeeze(nanmean(www,1));

for k=1:8
for n=1:24*60/N
wt(n,k)=nanmean(squeeze((www(:,n,k)-m_w(n,k)).*(ttt(:,n,k)-m_t(n,k))));
wq(n,k)=nanmean(squeeze((www(:,n,k)-m_w(n,k)).*(qqq(:,n,k)-m_q(n,k))));
wu(n,k)=nanmean(squeeze((www(:,n,k)-m_w(n,k)).*(uuu(:,n,k)-m_u(n,k))));
wv(n,k)=nanmean(squeeze((www(:,n,k)-m_w(n,k)).*(vvv(:,n,k)-m_v(n,k))));
end
end

m_time=squeeze(mean(time_N,1));
m_wt=squeeze(nanmean(wt,2));
m_wq=squeeze(nanmean(wq,2));
m_wu=squeeze(nanmean(wu,2));
m_wv=squeeze(nanmean(wv,2));

figure(1)
subplot(2,2,1);
plot(m_time,m_wt);
subplot(2,2,2);
plot(m_time,m_wq);
subplot(2,2,3);
plot(m_time,m_wu);
subplot(2,2,4);
plot(m_time,m_wv);

save data/flux_10m_dat.mat m_time m_wt m_wq m_wu m_wv;

---can

%hold on;
%plot(squeeze(mean(time_N,1)),squeeze(mean(ppp(:,:,1),1)),'r');
%hold off;

[n1 n2]=size(tt);

m_time=squeeze(nanmean(TT,1));
m_t1=squeeze(nanmean(tt,1));
m_q1=squeeze(nanmean(qq,1));
m_w1=squeeze(nanmean(ww,1));
m_u1=squeeze(nanmean(uu,1));
m_v1=squeeze(nanmean(vv,1));

for n=1:96
wq(n,1)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(qq(:,n)-m_q1(n))));
wu(n,1)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(uu(:,n)-m_u1(n))));
wv(n,1)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(vv(:,n)-m_v1(n))));
end

for n=97:163
wt(n-96,2)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(tt(:,n)-m_t1(n))));
wq(n-96,2)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(qq(:,n)-m_q1(n))));
wu(n-96,2)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(uu(:,n)-m_u1(n))));
wv(n-96,2)=nanmean(squeeze((ww(:,n)-m_w1(n)).*(vv(:,n)-m_v1(n))));
end

wt(68:96,2)=nan;
wq(68:96,2)=nan;
wu(68:96,2)=nan;
wv(68:96,2)=nan;


m_time2=mean(reshape(time2,[],34),1);
m_p2=squeeze(nanmean(p,1));
m_t2=squeeze(nanmean(t,1));
m_u2=squeeze(nanmean(u,1));
m_v2=squeeze(nanmean(v,1));
m_w2=squeeze(nanmean(w,1));

m_time3=mean(reshape(time3,[],38),1);
m_p3=squeeze(nanmean(p,1));
m_t3=squeeze(nanmean(t,1));
m_u3=squeeze(nanmean(u,1));
m_v3=squeeze(nanmean(v,1));
m_w3=squeeze(nanmean(w,1));

m_time4=mean(reshape(time4,[],12),1);
m_p4=squeeze(nanmean(p(:,1:12),1));
m_t4=squeeze(nanmean(t(:,1:12),1));
m_u4=squeeze(nanmean(u(:,1:12),1));
m_v4=squeeze(nanmean(v(:,1:12),1));
m_w4=squeeze(nanmean(w(:,1:12),1));

m_time5=mean(reshape(time5,[],47),1);
m_p5=squeeze(nanmean(p(:,13:59),1));
m_t5=squeeze(nanmean(t(:,13:59),1));
m_u5=squeeze(nanmean(u(:,13:59),1));
m_v5=squeeze(nanmean(v(:,13:59),1));
m_w5=squeeze(nanmean(w(:,13:59),1));

m_time6=mean(reshape(time6,[],51),1);
m_p6=squeeze(nanmean(p,1));
m_t6=squeeze(nanmean(t,1));
m_u6=squeeze(nanmean(u,1));
m_v6=squeeze(nanmean(v,1));
m_w6=squeeze(nanmean(w,1));

plot(time1,w1);
hold on;
plot(time2,w2);
plot(time3,w3);
plot(time4,w4);
plot(time5,w5);
plot(time6,w6);
plot(m_time1,m_w1,'r');
plot(m_time2,m_w2,'r');
plot(m_time3,m_w3,'r');
plot(m_time4,m_w4,'r');
plot(m_time5,m_w5,'r');
plot(m_time6,m_w6,'r');

hold off;


hold off;

