function lakedraw

global errmark
errmark = -999.99;

numexps = 1;

expname = cell(numexps);
exp_name = cell(numexps);
exp_name_short = cell(numexps);

%expname = 'essea1948-2011_hydrpres_inclinit(ref)/';
%expname = 'essea1948-2011_radius1mm/';
%expname = 'essea1948-2011_hydrpres/';
expname{1} = 'essea1948-2011_nospinup/';
%expname{2} = 'essea1948-2011/';
%expname = 'Valkea-Kotinen2006/';
%expname = 'Kivu2003-2011/';

%expname = 'essea1948-2011_prod0/';
%expname = 'essea1948-2011_ref/essea1948-2011/';
%expname = 'essea1948-2011_0/';
%expname = 'essea1948-2011_ts0/';
%expname = 'essea1948-2011_uv0/';

exp_name{1}='no_spinup'; %with hydrostatic pressure in soil';
%exp_name{2}='No nudging due to advection'; %with hydrostatic pressure in soil';
%exp_name='0lbc experiment'; %with hydrostatic pressure in soil';
%exp_name='ts0lbc experiment'; %with hydrostatic pressure in soil';
%exp_name='uv0lbc experiment'; %with hydrostatic pressure in soil';
%exp_name='meth-prod=0';

%exp_name='Experiment with salinity'; %with hydrostatic pressure in soil';
%exp_name='Bubble radius 1mm';
%exp_name='no-bubble-model'; %with hydrostatic pressure in soil';

exp_name_short{1}='no-spinup'; %with hydrostatic pressure in soil';
%exp_name_short{2}='No nudging'; %with hydrostatic pressure in soil';

%rootdir='/home/victor/Files/main/models/LAKE/2013/March/'; 
%'/home/victor/Files/main/models/LAKE/2013/February/';
rootdir='D:/models/LAKE/2013/April/lake/';
%rootdir='D:/models/LAKE/2012/April/lake/';
%myfiledir=strcat(rootdir,'tools/');

% Kivu
%dt_out = 0.5; % input data timestep, hours

% essea1948-2011
dt_out = 6; 

% Kivu2003-2011
%startyear = 2003; % The first year in input files
%outyear = 2010; % The year to display
%outmonth = 'dec'; % The first month of display
%ndays = 365; % number of days to display

% essea1948-2011
startyear = 1948; % The first year in input files
startmonth = 'jan'; % The first month in input files
outyear = 1965; % The year to display 1959
outmonth = 'jan'; % The first month of display
outdaymonth = 1; % The first day in a month to display
ndays = 40*365; %1035; % number of days to display 365
%ndays = 10; %30.42

nperiodaver = 12; % period of averaging, in stepout
numper = 40; % the number of averaging cycles

%stepout = 30.42; % Timestep of output, days
% 30.42 - mean month length for 2003-2011
%stepout = 91.25; % Timestep of output, days
% 91.25 - 3 months
stepout = 30.42; % Timestep of output, days

ndaysyear = 365;
outyearlast = outyear + floor(ndays/ndaysyear);

outday = daysact(strcat('1-',startmonth,'-',num2str(startyear)), ...
    strcat('1-',outmonth,'-',num2str(outyear))) + outdaymonth;

nstep1 = 24.*outday/dt_out+1;
nstepout = nearest(24.*stepout/dt_out); % output interval, model timesteps
nsteplast = nearest(24.*(outday+ndays)/dt_out)+1; % yeardays(outyear)

nleg = 12; % Number of legend items

legitem0 = cell(nleg);
%legitem0{1} = 'DJF';
%legitem0{2} = 'MAM';
%legitem0{3} = 'JJA';
%legitem0{4} = 'SON';

legitem0{1} = 'December';
legitem0{2} = 'January';
legitem0{3} = 'February';
legitem0{4} = 'March';
legitem0{5} = 'April';
legitem0{6} = 'May';
legitem0{7} = 'June';
legitem0{8} = 'July';
legitem0{9} = 'August'; %'August 20-30-th';
legitem0{10} = 'September';
legitem0{11} = 'October';
legitem0{12} = 'November';
%legitem0{13} = 'January + 1yr';
%legitem0{14} = 'February + 1yr';
%legitem0{15} = 'March + 1yr';


k = [4 9]; % switches for graph output [1 2 3 4 5 6 7 8 9]
h = zeros(max(k),1);

filename = cell(max(k),1);
filename_obs = cell(max(k),1); % names of files with observations
col_obs  = cell(max(k),1); % column numbers in files with observations
title = cell(max(k),1);
xtitle = cell(max(k),1);
for i = 1, size(col_obs)
    col_obs{i} = zeros(2,1); % depth and veriable columns will be specified
end

% Only two experiments are now allowed
switch numexps
    case 1
        linetype(1,1:2)='- ';
    case 2
        linetype(1,1:2)='- ';
        linetype(2,1:2)='-.';
end

legitem = cell(nleg,1);

for m = 1:numexps

for i=1:nleg
    legitem{i} = strcat(legitem0{i},'-',exp_name_short{m});
end
    
% Directory time_series
dirs = cell(max(k),1);
for i=k
   dirs{i}=strcat(rootdir,'results/',expname{m},'time_series/');
end

title_suffix = strcat(',',num2str(outyear),'-',num2str(outyearlast)); %exp_name{m},','
% Bubble fluxes
filename{1} = 'fbblflx_ch4  1  1.dat'; title{1} = strcat('Methane bubble flux',title_suffix); xtitle{1} = 'Flux, mol/(m^2*s)'; scale(1) = 1;
col(1,1:5) = 'green';
filename_obs{1} = ''; col_obs{1}(:) = errmark;

filename{2} = 'fbblflx_co2  1  1.dat'; title{2} = strcat('Carbon dioxide bubble flux',title_suffix); xtitle{2} = 'Flux, mol/(m^2*s)'; scale(2) = 1;
col(2,1:5) = 'green';
filename_obs{2} = ''; col_obs{2}(:) = errmark;

filename{3} = 'fbblflx_o2  1  1.dat'; title{3} = strcat('Oxygen bubble flux',title_suffix); xtitle{3} = 'Flux, mol/(m^2*s)'; scale(3) = 1;
col(3,1:5) = 'green';
filename_obs{3} = ''; col_obs{3}(:) = errmark;

% Gas concentration profiles
filename{4} = 'methane_water  1  1.dat'; title{4} = strcat('Methane concentration',title_suffix); xtitle{4} = 'Concentration, nmol/l'; scale(4) = 2;
col(4,1:5) = 'green';
filename_obs{4} = ''; col_obs{4}(:) = errmark;

filename{5} = 'co2_water  1  1.dat'; title{5} = strcat('Carbon dioxide concentration',title_suffix); xtitle{5} = 'Concentration, nmol/l'; scale(5) = 1;
col(5,1:5) = 'green';
filename_obs{5} = ''; col_obs{5}(:) = errmark;

filename{6} = 'oxygen_water  1  1.dat'; title{6} = strcat('Oxygen concentration',title_suffix); xtitle{6} = 'Concentration, nmol/l'; scale(6) = 1;
col(6,1:5) = 'green';
filename_obs{6} = ''; col_obs{6}(:) = errmark;

% Temperature profile
filename{7} = 'water_temp  1  1.dat'; title{7} = strcat('Water temperature',title_suffix); xtitle{7} = 'Temperature, Celsius'; scale(7) = 1;
col(7,1:3) = 'red';
filename_obs{7} = 'HX189_aug95.xls'; %'HX189_aug95.xls'; 
% Note xlsread interprets xls file differently in Linux and Windows
%col_obs{7}(:) = errmark; 
col_obs{7}(1) =1; %4; % depth 
col_obs{7}(2) =2; %5; % temperature

% Salinity profile
filename{8} = 'sal_water  1  1.dat'; title{8} = strcat('Water salinity',title_suffix); xtitle{8} = 'Salinity, n/d'; scale(8) = 1;
col(8,1:4) = 'blue';
filename_obs{8} = ''; col_obs{8}(:) = errmark;

% Soil ice profile
filename{9} = 'wi_soil  1  1.dat'; title{9} = strcat('Methane hydrate',title_suffix); xtitle{9} = 'Methane hydrate, kg/kg'; scale(9) = 1;
col(9,1:4) = 'blue';
filename_obs{9} = ''; col_obs{9}(:) = errmark;

ntcols = 1; % number of non-data columns
nformat = 1; % format of profile time series

aver = 1; % averaging over output interval
avercalc = 0; % calculating of seasonal averages
for i = k
    if m == 1
        h(i) = figure; % h - figure handle
    end
    plotgraphs(dirs{i},filename{i},filename_obs{i},col_obs{i},title{i},xtitle{i},...
        legitem,nstep1,nstepout,nsteplast,nformat,ntcols,aver,h(i),scale(i),...
        col(i,1:size(col,2)),linetype(m,1:2),m,avercalc,nperiodaver,numper)
end

end

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

function[]=plotgraphs(dir,filename,filename_obs,col_obs,title1,xtitle,legitem,nstep1,...
    nstepout,nsteplast,nformat,ntcols,aver,h,scale,col,linetype,nexp,avercalc,nperiodaver,numper)

global errmark

figure(h)
hold on % holds the plot when adding new objects

switch deblank(col)
    case('red')
        ncol(1) = 1; ncol(2) = 0; ncol(3) = 0;
    case('green')
        ncol(1) = 0; ncol(2) = 1; ncol(3) = 0;
    case('blue')
        ncol(1) = 0; ncol(2) = 0; ncol(3) = 1;
end

if nperiodaver > 0
    work1ss = cell(nperiodaver,1);
    work2ss = cell(nperiodaver,1);
end

small_value = 1.e-20;

chdir(dir)

nper = 0;
nnper = 1;
fileexist = exist(filename,'file');
if (fileexist == 2)
    [work,m,n]=getlmipascii(dir,filename);
    if (nformat == 1)
        
       title(title1,'FontSize',20);
       set(gca,'XAxisLocation','top','Fontsize',18)
       xlabel(xtitle,'Fontsize',18)
       ylabel('Depth, m','Fontsize',18)            
       set(gca,'YDir','reverse','FontSize',16)    
       
       %xmax = errmark;
       %xmin = errmark;
       %ymax = errmark;
       %ymin = errmark;
       
       work1s(1:(n-ntcols)/2) = 0;
       work2s(1:(n-ntcols)/2) = 0;          
       nline = 0;
       for i=nstep1:min(nsteplast,m)
          k = ntcols + 1;
          work1(1:(n-ntcols)/2) = errmark;
          work2(1:(n-ntcols)/2) = errmark;
          while (k<n)
             work1((k+1-ntcols)/2) = work.data(i,k);
             work2((k+1-ntcols)/2) = work.data(i,k+1);
             k = k + 2;
          end
          if (aver) % averaging over output interval
              work1s = work1s + work1;
              work2s = work2s + work2;
          else
              work1s = work1;
              work2s = work2;
          end
          if (mod(i-nstep1,nstepout) == 0 && ~(i == nstep1 && aver))
              
            if (aver)
               work1s = work1s/nstepout;
               work2s = work2s/nstepout;              
            end
            
            if nperiodaver > 0
                nper = nper + 1;
                if nper > nperiodaver
                    nper = 1;
                    nnper = nnper + 1;
                end
                if i == nstep1 + nstepout
                    for kk = 1:nperiodaver
                        work1ss{kk} = zeros(size(work1s,1),1);
                        work2ss{kk} = zeros(size(work2s,1),1);
                    end
                end
                work1ss{nper} = work1ss{nper} + work1s;
                work2ss{nper} = work2ss{nper} + work2s;
                work1s = 0;
                work2s = 0;
            end 
            
           % Obsolete code for axis limits calculation
           % if (xmax == errmark)
           %    xmax = max(work1s);
           %     xmin = min(work1s);
           %    ymax = max(work2s);
           %     ymin = min(work2s);                
           % else
           %     xmax = max([xmax,max(work1s)]);
           %     xmin = min([xmin,min(work1s)]);
           %     ymax = max([ymax,max(work2s)]);
           %     ymin = min([ymin,min(work2s)]);                
           % end
           % if (ymin == ymax) 
           %     ymin = ymin - small_value;
           %     ymax = ymax + small_value;
           % end
           % axis([ymin ymax xmin xmax])
           
           if nperiodaver > 0 && nnper == numper
               work2wr = work2ss{nper}/numper;
               work1wr = work1ss{nper}/numper;
           elseif nperiodaver > 0
               continue
           else
               work2wr = work2s;
               work1wr = work1s;
           end                       
            
            if nperiodaver > 0
                color = 1 - 1.*(nper-1)/(nperiodaver-1);
            else
                color = 1 - 1.*(i-nstep1)/(nsteplast-nstep1);
            end
            
            % output beginning with second level from the surface, since
            % surface scalars in LAKE are not correct
            xr = ncol(1)*color;
            xg = ncol(2)*color;
            xb = ncol(3)*color;
            hline = line(work2wr(1:size(work2wr,2)),work1wr(1:size(work1wr,2)),'Color',...
                [xr xg xb],'LineStyle',linetype) ;
            axis tight %auto
            nline = nline + 1;
            set(findobj('Type','line'),'LineWidth',2)
            if scale == 2 % logarithmic scale
                set(gca,'XScale','log')
            end
            if nexp > 1
                % Get object handles
                [~,~,OUTH,OUTM] = legend;
                % Add object with new handle and new legend string to legend
                legend([OUTH;hline],OUTM{:}, legitem{nline},'Location','Best')
            end
%           legitem{(i-nstep1)/nstepout+1} = num2str((i-nstep1)/nstepout+1);
            work1s = 0;
            work2s = 0;
          end                    
       end       
       
       % if LEGEND is not created in the figure, create default legend:
       if nexp == 1
           legend(legitem{:},'Location','Best');
       end       
       
       disp(['Graph plotted with data from_',dir,'_',filename]) 
    end
end

% Reading and plotting data from file with measurements
fileexist = exist(filename_obs,'file');
if fileexist == 2
   [num,txt,raw] = xlsread(filename_obs);    
   hscat = scatter(num(:,col_obs(2)),num(:,col_obs(1)),'filled'); %,'Color',...
                % [ncol(1) ncol(2) ncol(3)]) ; % Scatter plot with circle markers
   axis tight          
   % Adding new legend item
   [~,~,OUTH,OUTM] = legend;   
   legend([OUTH;hscat],OUTM{:}, 'Observations','Location','Best')
end

if avercalc
% Calculating averages
ncolmonth = 2;
ncol = 45; % 7
aver(1:2) = 0;
j = 0; k = 0;
dt_out = 6; %model output timestep, hours
istart = 3*365*24/dt_out; %start immediately after spinup period
for i = istart:m-1 % the last record is not complete!
    if work.data(i,ncolmonth) >= 6 && work.data(i,ncolmonth) <= 9 % summer
        aver(1) = aver(1) + work.data(i,ncol);
        j = j + 1;
    else %if work.data(i,ncolmonth) >= 3 && work.data(i,ncolmonth) <= 5 % late winter
        aver(2) = aver(2) + work.data(i,ncol);
%        if (isnan(aver(2)))
%            disp(strcat('NaN!',num2str(work.data(i,ncol)),num2str(i)))
%        end
        k = k + 1;        
    end
end
aver(1) = aver(1)/j;
aver(2) = aver(2)/k;
disp(strcat('Averages of_',title1))
disp('Average, summer')
disp(num2str(aver(1)))
disp('Average, winter')
disp(num2str(aver(2)))
end

% Printing figure to png file
print(h,'-dpng',strcat(title1,'.png'))

%close % Closing figure

end