function bouy_readinterp(filename)

% This script reads the data from the Spanish Wave-Buoy Net and 
% fills the date gaps with NaNs or interpolated data.
% The output are vector variables with a user selected extension 
% version 1.6 04/05/2010 (Thanks to clever Irene)
% Haris Plomaritis, University of Cadiz 

% detects null data code and headerline automatically 

% READ DATA FROM WAVE BUOY FILE 

% filename = input ('Enter data filename   ','s');

% disp(' ');
% disp('Enter the number of headline rows')
%%
headline_number= 84; % input('headline number=');


% disp('Enter the normal time step of the recorded data (in hours) ')
step_norm= 1; %input('time step=');

% READ DATA FROM THE FILE
 A = importdata(filename,' ',headline_number); % Delimeter is set as tab

PARAM = A.textdata(:,1); % make text variable with all the informations
data = A.data;
DATE = [data(:,1) data(:,2) data(:,3) data(:,4) 0*data(:,4) 0*data(:,4)];
var_names = A.colheaders;
[j,k] = size(data);

%% ******************************************
% FIND GAPS  MAKE DATES AND FILL DATA WITH NaNs
% ******************************************

dnum = datenum([DATE]);
step = diff(dnum);
step_d=double(vpa(step_norm/24,8));
indst = find(step>step_d);
gap_length = fix(step(indst)./step_d);


% disp('Enter the maximum data gap (in hours)')
% disp('that you want to be filled with interpolations.')
% disp('(1. Large data gaps is better not to be interpolated)')
% disp('(2. If you want no interpolation to be undertaken set value to 0)')

gap_inter=6; %input('Data gap to be interpolated (in hours) <=');
gap_inter=gap_inter./step_norm;
A = importdata(filename,' ',headline_number);

idx = strmatch(' DATO NULO  : Es representado por', PARAM); % find missing data code
[null] = textread(filename, 'DATO NULO  : Es representado por %f',1,'headerlines',idx-1); %read missing data code

data1 = data;
indnan = find(data1 == null);
data1(indnan)=NaN;
indst1= indst;
counter1 = 0;
for i=1:length(indst)
    gap_time(i) = {datevec([dnum(indst(i))+step_d:step_d:dnum(indst(i))+step(indst(i))])};
    if gap_length(i) > gap_inter

        gap_fill(i) = {[gap_time{1,i}(:,1:4) NaN*ones(gap_length(i),k-4)]};
    else
        warning('off')
        gap_fill(i) = {[gap_time{1,i}(:,1:4) interp1([dnum(indst(i):indst(i)+1)],...
                        data1(indst(i):indst(i)+1,5:k),...
                        [dnum(indst(i))+step_d:step_d:dnum(indst(i))+step(indst(i))])]};
        warning('on')
        counter1 = counter1+1;
    end
    
    data_full   = [data1(1:indst1(i),:); gap_fill{1,i}(:,:); data1(indst1(i)+1:end,:)];
    data1       = data_full;
    
    if i<length(indst)
        indst1(i+1) = indst1(i+1)+sum(gap_length(1:i)); % CALCULATE THE NEW INDICES 
    else
    end
end

% ** FINISH FILLING THE GAPS 
% ******************************************************

%% *****************************************************
%  ** CHECK FOR PRE-EXISTING NaN GAPS AND FILL THEM UP
nj = length(data1);                         % get new data file size


 for j=5:k                                   % operation for each vector-row separately 
    indnanj = find(isnan(data1(:,j)));       % find Nans
    if numel(indnanj)~=0;                %If the matrix is not empty (Irene)
    df = diff(indnanj);                      % calc differences
    indnanjj = find(df>1);                   % find start of nan
    if numel(indnanjj)~=0;               %If the matrix is not empty (Irene)
    dff = diff(indnanjj);                    % calculate length of NaNs blocks 
    last = flipud(indnanj);
    dflast= diff(last);
    indlast = find(dflast<-1);
    
% make a vector with gap lengths 
    if indnanj(end)==nj;
        new_gaplength = [indnanjj(1); dff];             % omit last gap if at the end of file
    else
        new_gaplength = [indnanjj(1); dff; indlast(1)];    
    end
    
    loc=0;
    newgap_fill = cell(k, length(new_gaplength));
    for jj=1:length(new_gaplength)
        jj
        if  new_gaplength(jj) <= gap_inter
            newgap_fill(j,jj)  = {interp1([indnanj(1+loc)-1 indnanj(new_gaplength(jj)+loc)+1],...
                                  data1([indnanj(1+loc)-1 indnanj(new_gaplength(jj)+loc)+1],j),...
                                  [indnanj(1+loc) :indnanj(new_gaplength(jj)+loc)])'};
                         
            data1([indnanj(1+loc) :indnanj(new_gaplength(jj)+loc)],j)   = newgap_fill{j,jj}(:,:); 
            counter1           =  counter1+1;
                            
            loc =sum(new_gaplength(1:jj));                 % CALCULATE THE NEW INDICES 
        else
            loc =sum(new_gaplength(1:jj));                 % CALCULATE THE NEW INDICES 
        end
    end
   end %Irene
   end %Irene

 end


% %%
% disp(' ');
% disp('Give an extension code for the data')
% ext = input('ext=','s')
% 
 eval(['DATE = [data1(:,1) data1(:,2) data1(:,3) data1(:,4) 0*data1(:,4) 0*data1(:,4)];']);
 
 for ii = 1:k
     eval([char(var_names(ii)) '= data1(:,' num2str(ii) ');']);  % make variables
     eval(['ind = find(' char(var_names(ii)) '== -99.9);']);     % find bad values 
     eval([char(var_names(ii)) '(ind) = NaN;']);                 % convert to NaN
 end

%% PRINT SCREEN REPORT 

disp(' ');
disp([ num2str(length(indst)) ' DATE GAPS FOUND AND FILLED'])

disp(' ');
disp([ num2str(counter1) ' DATA GAPS (NaNs) FOUND TO BE <' num2str(gap_inter) 'hours AND INTERPOLATED'])
% out =
% clear A
% clear -regexp  '^AA|^MM|^DD|^HH|^ata|^i|^j|^k|^d|^l|^counter|^step|^var_names|^gap|^file|^n|^ext|^head'


%% Save data
fichero_all=[filename(1:end-4) '_noGaps'];
savefile_all = [fichero_all,'.mat'];
save(savefile_all,'PARAM','DATE', var_names{5:end} );
