data=importdata('C:\Data\GOMECC\filter_total2.txt');
data2=importdata('C:\Data\GOMECC\GOMECC07_flux_10uf.txt');
% This file inclueds the GOMECC data filtered for reldir, rain, lagtime,
% std(O3>3), delta O3ppb>3, filterchange, calibration, count drop. There
% fore the column O3ppb is for each filter copied and filtered values set
% to NaN

jd=data(1101:3480,1);               %julianday
O3ppb=data(1101:3480,2);            %O3ppb
u=data(1101:3480,3);                %true wind speed m/s
reldir=data(1101:3480,4);           %NaN=wrong wind direcrion
rain=data(1101:3480,5);             %NaN=rain>5    
lag=data(1101:3480,6);              %NaN=lag>70,<40
std=data(1101:3480,7);              %NaN=std>3
delta=data(1101:3480,8);            %NaN=deltaO3ppb>3
ppb5=data(1101:3480,9);             %NaN=O3ppb<5,>80
filter=data(1101:3480,12);          %NaN=filterchange
cali=data(1101:3480,13);            %NaN=calibration
miss=data(1101:3480,14);            %NaN=missing data
countdrop=data(1101:3480,15);       %NaN=count drop
total=data(1101:3480,16);            %total filter
total_o=data(1101:3480,17);          %total filter without lag
slowO3=data2(1101:3480,56);

jj=find(isnan(lag));
lag(jj)=1;

ii=find(isnan(reldir));
reldir(ii)=2;

gg=find(isnan(delta));
delta(gg)=3;

hh=find(isnan(std));
std(hh)=4;

ll=find(isnan(filter));
filter(ll)=5;

ff=find(isnan(ppb5));
ppb5(ff)=6;

tt=find(isnan(miss));
miss(tt)=7;
 
ss=find(isnan(countdrop));
countdrop(ss)=8;

dd=find(isnan(cali));
cali(dd)=9;

ij=find(isnan(u));
u(ij)=10;

kk=find(isnan(rain));
rain(kk)=11;


figure;plot(jd(kk),rain(kk),'y.') ;hold on; plot(jd(ij),u(ij),'c.');plot(jd(dd),cali(dd),'g.');plot(jd(ss),countdrop(ss),'r.');
plot(jd(tt),miss(tt),'k.');plot(jd(ff),ppb5(ff),'m.');plot(jd(ll),filter(ll),'y.');plot(jd(hh),std(hh),'c.');
plot(jd(gg),delta(gg),'g.');plot(jd(ii),reldir(ii),'r.');plot(jd(jj),lag(jj),'m.');
xlabel('day of year');axis([190 220 0 12]);
legend('11:rain','10:wind speed','9:calibration', '8:count drop', '7:missing data', '6:O3<5ppb,>180ppb', '5:filter change','4:std O3', '3:delta O3','2:relative wind direction' , '1:lagtime');

l=find(isnan(total));
figure;plot(jd(l),O3ppb(l),'.');hold on;plot(jd,total,'r.');
xlabel('Day of Year'); ylabel('[O3]/ppb');title('blue: eliminated, red: remainig; including lag time filter');
axis([190 220 -1 90])

k=find(isnan(total_o));
figure;plot(jd(k),O3ppb(k),'.');hold on;plot(jd,total_o,'r.');
xlabel('Day of Year'); ylabel('[O3]/ppb');title('blue: eliminated, red: remainig; excluding lag time filter');
axis([190 220 -1 90])

j=find(isfinite(total_o));
i=find(isnan(total_o));
figure;plot(jd(i),slowO3(i),'.');hold on;plot(jd(j),slowO3(j),'r.');
xlabel('Day of Year'); ylabel('[O3]/ppb');title('slowO3;blue: eliminated, red: remainig; excluding lag time filter');
axis([190 220 -1 90])

figure;plot(jd,O3ppb,'.');hold on;plot(jd,slowO3,'.r')
xlabel('Day of Year'); ylabel('[O3]/ppb');
title('red:O3 ML8810, blue: O3 FRCI');
axis([190 220 -1 90])

%Fluxes and deposition velocity
jday=data2(1101:3480,1);
lagtime=data2(1101:3480,43);            % lagtime between w and O3 calculated by crosscovriance, 10 min lagtime
O3w_o=data2(1101:3480,44);              % w'O3' calculated with crosscovariance, 10 min lagtime, no webb correction
O3w_h=data2(1101:3480,45);              % w'O3' calculated with crosscovariance, flux at hourly mean lagtime, no webb correction
O3w_d=data2(1101:3480,46);              % w'O3' calculated with crosscovariance, flux at dayly mean lagtime, no webb correction
O3mean_o=data2(1101:3480,47);           % mean O3 concentration cross (ppb)
O3dev_o=data2(1101:3480,48);            % stdev ozone cross (ppb)
O3w=data2(1101:3480,49);                % O3'w' with constant lagtime of 5.1 s, and function cov, no webb correction
O3mean=data2(1101:3480,50);             % mean O3 concentration constant lagtime (ppb)
O3dev=data2(1101:3480,51);              % stdev ozone constant lagtime (ppb)
Vd_o=data2(1101:3480,52);               % deposition velocity 10 min lagtime (cm/s)
Vd_h=data2(1101:3480,53);               % deposition velocity hourly lagtime (cm/s)
Vd_d=data2(1101:3480,54);               % deposition velocity dayly lagtime (cm/s)
Vd=data2(1101:3480,55);                 % deposition velocity 5.1s lagtime (cm/s)
slowO3=data2(1101:3480,56);
lat=data2(1101:3480,23);
lon=data2(1101:3480,24);


%including lag time filter
p=find(isfinite(total));

Vdmean=mean(Vd(p))
Vdmean_o=mean(Vd_o(p))
Vdmedian=median(Vd(p))
Vdmedian_o=median(Vd_o(p))

ij=find(Vd(p)>0);
Vdmean2=mean(Vd(p(ij)))
Vdmedian2=median(Vd(p(ij)))

ij=find(Vd_o(p)>0);
Vdmean_o2=mean(Vd_o(p(ij)))
Vdmedian_o2=median(Vd_o(p(ij)))

figure; plot(jday(p),Vd_o(p),'.');
xlabel('Day of Year');ylabel('Vd / cm*s-1');title('Vd cross correl including lag time filter');
figure; plot(jday(p),Vd(p),'.');title('Vd constant lag time icluding lag time filter');
xlabel('Day of Year');ylabel('Vd / cm*s-1');

%Histogram
x=-5:0.005:5;
figure;hist(Vd_o(p),x);
xlabel('Vd / cm*s-1'); ylabel('Number of Data Points');title('Vd cross correl including lag time filter');
figure;hist(Vd(p),x);title('Vd constant lag time including lag time filter');
xlabel('Vd / cm*s-1'); ylabel('Number of Data Points');

%no lagtime as filter
g=find(isfinite(total_o));

Vdmean3=mean(Vd(g))
Vdmean_o3=mean(Vd_o(g))
Vdmedian3=median(Vd(g))
Vdmedian_o3=median(Vd_o(g))

ij=find(Vd(g)>0);
Vdmean4=mean(Vd(g(ij)))
Vdmedian4=median(Vd(g(ij)))

ij=find(Vd_o(g)>0);
Vdmean_o4=mean(Vd_o(g(ij)))
Vdmedian_o4=median(Vd_o(g(ij)))

figure; plot(jday(g),Vd_o(g),'.');
xlabel('Day of Year');ylabel('Vd / cm*s-1');title('Vd cross correl excluding lag time filter');

figure; plot(jday(g),Vd(g),'.');
xlabel('Day of Year');ylabel('Vd / cm*s-1');title('Vd const lag time excluding lag time filter');

%Histogramm data filtered
x=-5:0.005:5;
figure;hist(Vd_o(g),x);
xlabel('Vd / cm*s-1'); ylabel('Number of Data Points');title('Vd cross correl excluding lag time filter');
figure;hist(Vd(g),x);
xlabel('Vd / cm*s-1'); ylabel('Number of Data Points');title('Vd constant lag time excluding lag time filter');
