function series

rootdir='/home/victor/Files/main/models/LAKE/2014/06/lake/';
%expname='essea1948-2011_nospinup/';
expname='Seida1_2007-2008/';
dir=strcat(rootdir,'results/',expname,'time_series/');

monthdays0 = 30.42;

work = cell(2,1);

[work{1},m,n]=getlmipascii(dir,'methane_series  1  1.dat');
nsurfconc = 7;
nbotconc = 8;
nbotebul = 18;
nbotdiff = 20;
%nsurfebul = 
nsurfdiff = 21;
ntime(1) = 5;

[work{2},m,n]=getlmipascii(dir,'T_fluxes  1  1.dat');
nsignwaveheight = 20;
nsurftemp = 6;
nbottemp = 9;   
ntime(2) = 5;

hr2yr = 365*24;
hr2mnth = 30.42*24;
molm3tonM = 1E+6;
mcmol3tonM = 1E+3;

nmeas = 6;
% Observed values - water surface temperature (Seida1)
tsurfmeas = zeros(nmeas,4);
ch4conc    = zeros(nmeas,4);
ch4dflx     = zeros(nmeas,4);

% Lake 1
tsurfmeas(1,1) = 2007; ch4conc(1,1) = 2007; ch4dflx(1,1) = 2007; 
tsurfmeas(1,2) = 7; ch4conc(1,2) = 7; ch4dflx(1,2) = 7;
tsurfmeas(1,3) = 22; ch4conc(1,3) = 22; ch4dflx(1,3) = 22; 
tsurfmeas(1,4) = 19.70; ch4conc(1,4) = 0.22*mcmol3tonM; ch4dflx(1,4) = NaN; 

tsurfmeas(2,1) = 2007; ch4conc(2,1) = 2007; ch4dflx(2,1) = 2007; 
tsurfmeas(2,2) = 7; ch4conc(2,2) = 7; ch4dflx(2,2) = 7; 
tsurfmeas(2,3) = 28; ch4conc(2,3) = 28; ch4dflx(2,3) = 28; 
tsurfmeas(2,4) = 18.60; ch4conc(2,4) = 0.28*mcmol3tonM; ch4dflx(2,4) = NaN; 

tsurfmeas(3,1) = 2007; ch4conc(3,1) = 2007; ch4dflx(3,1) = 2007; 
tsurfmeas(3,2) = 8; ch4conc(3,2) = 8; ch4dflx(3,2) = 8; 
tsurfmeas(3,3) = 4; ch4conc(3,3) = 4; ch4dflx(3,3) = 4; 
tsurfmeas(3,4) = 16.80; ch4conc(3,4) = 0.30*mcmol3tonM; ch4dflx(3,4) = 4.79; 

tsurfmeas(4,1) = 2007; ch4conc(4,1) = 2007; ch4dflx(4,1) = 2007; 
tsurfmeas(4,2) = 8; ch4conc(4,2) = 8; ch4dflx(4,2) = 8; 
tsurfmeas(4,3) = 11; ch4conc(4,3) = 11; ch4dflx(4,3) = 11; 
tsurfmeas(4,4) = 12.40; ch4conc(4,4) = 0.42*mcmol3tonM; ch4dflx(4,4) = 9.47; 

tsurfmeas(5,1) = 2007; ch4conc(5,1) = 2007; ch4dflx(5,1) = 2007; 
tsurfmeas(5,2) = 8; ch4conc(5,2) = 8; ch4dflx(5,2) = 8; 
tsurfmeas(5,3) = 18; ch4conc(5,3) = 18; ch4dflx(5,3) = 18; 
tsurfmeas(5,4) = 15.70; ch4conc(5,4) = 0.44*mcmol3tonM; ch4dflx(5,4) = 9.75; 

tsurfmeas(6,1) = 2007; ch4conc(6,1) = 2007; ch4dflx(6,1) = 2007; 
tsurfmeas(6,2) = 8; ch4conc(6,2) = 8; ch4dflx(6,2) = 8; 
tsurfmeas(6,3) = 27; ch4conc(6,3) = 27; ch4dflx(6,3) = 27; 
tsurfmeas(6,4) = 11.90; ch4conc(6,4) = 0.08*mcmol3tonM; ch4dflx(6,4) = 0.93; 

% Lake 2
% tsurfmeas(1,1) = 2007; ch4conc(1,1) = 2007; ch4dflx(1,1) = 2007; 
% tsurfmeas(1,2) = 7; ch4conc(1,2) = 7; ch4dflx(1,2) = 7;
% tsurfmeas(1,3) = 22; ch4conc(1,3) = 22; ch4dflx(1,3) = 22; 
% tsurfmeas(1,4) = 19.70; ch4conc(1,4) = 0.52*mcmol3tonM; ch4dflx(1,4) = NaN; 
% 
% tsurfmeas(2,1) = 2007; ch4conc(2,1) = 2007; ch4dflx(2,1) = 2007; 
% tsurfmeas(2,2) = 7; ch4conc(2,2) = 7; ch4dflx(2,2) = 7; 
% tsurfmeas(2,3) = 28; ch4conc(2,3) = 28; ch4dflx(2,3) = 28; 
% tsurfmeas(2,4) = 18.60; ch4conc(2,4) = 0.3*mcmol3tonM; ch4dflx(2,4) = NaN; 
% 
% tsurfmeas(3,1) = 2007; ch4conc(3,1) = 2007; ch4dflx(3,1) = 2007; 
% tsurfmeas(3,2) = 8; ch4conc(3,2) = 8; ch4dflx(3,2) = 8; 
% tsurfmeas(3,3) = 4; ch4conc(3,3) = 4; ch4dflx(3,3) = 4; 
% tsurfmeas(3,4) = 16.80; ch4conc(3,4) = 0.29*mcmol3tonM; ch4dflx(3,4) = 0.97; 
% 
% tsurfmeas(4,1) = 2007; ch4conc(4,1) = 2007; ch4dflx(4,1) = 2007; 
% tsurfmeas(4,2) = 8; ch4conc(4,2) = 8; ch4dflx(4,2) = 8; 
% tsurfmeas(4,3) = 11; ch4conc(4,3) = 11; ch4dflx(4,3) = 11; 
% tsurfmeas(4,4) = 12.40; ch4conc(4,4) = 0.42*mcmol3tonM; ch4dflx(4,4) = 2.68; 
% 
% tsurfmeas(5,1) = 2007; ch4conc(5,1) = 2007; ch4dflx(5,1) = 2007; 
% tsurfmeas(5,2) = 8; ch4conc(5,2) = 8; ch4dflx(5,2) = 8; 
% tsurfmeas(5,3) = 18; ch4conc(5,3) = 18; ch4dflx(5,3) = 18; 
% tsurfmeas(5,4) = 15.70; ch4conc(5,4) = 0.29*mcmol3tonM; ch4dflx(5,4) = 2.11; 
% 
% tsurfmeas(6,1) = 2007; ch4conc(6,1) = 2007; ch4dflx(6,1) = 2007; 
% tsurfmeas(6,2) = 8; ch4conc(6,2) = 8; ch4dflx(6,2) = 8; 
% tsurfmeas(6,3) = 27; ch4conc(6,3) = 27; ch4dflx(6,3) = 27; 
% tsurfmeas(6,4) = 11.90; ch4conc(6,4) = 0.37*mcmol3tonM; ch4dflx(6,4) = 1.46; 

% Lake 3
% tsurfmeas(1,1) = 2007; ch4conc(1,1) = 2007; ch4dflx(1,1) = 2007; 
% tsurfmeas(1,2) = 7; ch4conc(1,2) = 7; ch4dflx(1,2) = 7;
% tsurfmeas(1,3) = 22; ch4conc(1,3) = 22; ch4dflx(1,3) = 22; 
% tsurfmeas(1,4) = 19.70; ch4conc(1,4) = 0.49*mcmol3tonM; ch4dflx(1,4) = 6.19; 
% 
% tsurfmeas(2,1) = 2007; ch4conc(2,1) = 2007; ch4dflx(2,1) = 2007; 
% tsurfmeas(2,2) = 7; ch4conc(2,2) = 7; ch4dflx(2,2) = 7; 
% tsurfmeas(2,3) = 28; ch4conc(2,3) = 28; ch4dflx(2,3) = 28; 
% tsurfmeas(2,4) = 18.60; ch4conc(2,4) = 0.92*mcmol3tonM; ch4dflx(2,4) = 22.54; 
% 
% tsurfmeas(3,1) = 2007; ch4conc(3,1) = 2007; ch4dflx(3,1) = 2007; 
% tsurfmeas(3,2) = 8; ch4conc(3,2) = 8; ch4dflx(3,2) = 8; 
% tsurfmeas(3,3) = 4; ch4conc(3,3) = 4; ch4dflx(3,3) = 4; 
% tsurfmeas(3,4) = 16.80; ch4conc(3,4) = 0.17*mcmol3tonM; ch4dflx(3,4) = 2.38; 
% 
% tsurfmeas(4,1) = 2007; ch4conc(4,1) = 2007; ch4dflx(4,1) = 2007; 
% tsurfmeas(4,2) = 8; ch4conc(4,2) = 8; ch4dflx(4,2) = 8; 
% tsurfmeas(4,3) = 11; ch4conc(4,3) = 11; ch4dflx(4,3) = 11; 
% tsurfmeas(4,4) = 12.40; ch4conc(4,4) = NaN; ch4dflx(4,4) = NaN; 
% 
% tsurfmeas(5,1) = 2007; ch4conc(5,1) = 2007; ch4dflx(5,1) = 2007; 
% tsurfmeas(5,2) = 8; ch4conc(5,2) = 8; ch4dflx(5,2) = 8; 
% tsurfmeas(5,3) = 18; ch4conc(5,3) = 18; ch4dflx(5,3) = 18; 
% tsurfmeas(5,4) = 15.70; ch4conc(5,4) = 0.07*mcmol3tonM; ch4dflx(5,4) = 2.5; 
% 
% tsurfmeas(6,1) = 2007; ch4conc(6,1) = 2007; ch4dflx(6,1) = 2007; 
% tsurfmeas(6,2) = 8; ch4conc(6,2) = 8; ch4dflx(6,2) = 8; 
% tsurfmeas(6,3) = 27; ch4conc(6,3) = 27; ch4dflx(6,3) = 27; 
% tsurfmeas(6,4) = 11.90; ch4conc(6,4) = 0.44*mcmol3tonM; ch4dflx(6,4) = 18.55; 

timemeas = zeros(nmeas,1);
for i=1:nmeas
    timemeas(i) = (tsurfmeas(i,1) - 2007)*12 + ...
        tsurfmeas(i,2)-1 + tsurfmeas(i,3)/monthdays0;        %
end

% Methane time series
h1 = figure;
line(work{1}.data(:,ntime(1))/hr2mnth,work{1}.data(:,nsurfconc)*molm3tonM,'Color','b')
line(work{1}.data(:,ntime(1))/hr2mnth,work{1}.data(:,nbotconc)*molm3tonM,'Color','r')
line(timemeas,ch4conc(:,4),'Color','k','Marker','o')
legend('Surface concentration','Bottom concentration','Surface concentration measured')
xlabel('Time, months'); ylabel('Concentration, nM')

h2 = figure;
line(work{1}.data(:,ntime(1))/hr2mnth,work{1}.data(:,nbotebul),'Color','g')
line(work{1}.data(:,ntime(1))/hr2mnth,work{1}.data(:,nbotdiff),'Color','r')
line(work{1}.data(:,ntime(1))/hr2mnth,work{1}.data(:,nsurfdiff),'Color','b')
line(timemeas,ch4dflx(:,4),'Color','k','Marker','o')
legend('Bottom ebullition','Bottom diffusion','Surface diffusion','Surface diffusion measured')
xlabel('Time, months'); ylabel('Flux, mg/(m^2 day)')

% Temperature time series
h3 = figure;
line(work{2}.data(:,ntime(2))/hr2mnth,work{2}.data(:,nsurftemp),'Color','r')
line(work{2}.data(:,ntime(2))/hr2mnth,work{2}.data(:,nbottemp),'Color','b')
line(timemeas,tsurfmeas(:,4),'Color','k','Marker','o')
% scatter(timemeas,tsurfmeas(:,4),'Marker','o')
legend('Surface temperature','Bottom temperature','Surface temperature measured')
xlabel('Time, months'); ylabel('Temperature, Celsius')

%line(4+work{2}.data(:,ntime(2))/hr2mnth,work{2}.data(:,nsignwaveheight),'Color','b') % 4 - starting month

% Averaging
firstyearaver = 0; %15; % The first year of averaging period
ndata = size(work{1}.data,1);
for i = 1:ndata
    if work{1}.data(i,ntime(1))/hr2yr > firstyearaver
        nstart = i;
        break
    end
end

disp(['Mean bottom diffusion=' num2str(mean(work{1}.data(nstart:ndata,nbotdiff)))])
disp(['Mean surface diffusion=' num2str(mean(work{1}.data(nstart:ndata,nsurfdiff)))])

end

function[work,m,n]=getlmipascii(dir,filename)
% n - number of columns in a file
chdir(dir)
work = importdata(strcat('./',filename));  %import raw-column data from ASCII
m = size(work.data,1);
n = size(work.data,2);
end
