function [W]=fillGAPSindata(W0,fields1,maximumgapsize)

%     % using arguments (e.g. 'size',5 OR 'time',0.3 days)
%     if ~isempty(findstr(lower(varargin{2}),'size')) || ~isempty(findstr(lower(varargin{2}),'gap'))
%         maximumgapsize = varargin{3}
%     elseif ~isempty(findstr(lower(varargin{2}),'time')) || ~isempty(findstr(lower(varargin{2}),'dt'))
%         
%     else
%         fprintf('Error input not specified correctly\n')
%         return
%     end
    fprintf('  - Fill small gaps in data with  interpolation [');

    NRstations = length(W0);
    NRfields   = length(fields1);
    W = struct;
    for jj=1:NRstations
        IDnans={};
        if isfield(W0,'loc')
            W(jj).loc  = W0(jj).loc;
        end
        
        for kk=1:NRfields
            if isfield(W0(jj),fields1{kk})
                if isfield(W0(jj),'time')
                    if isfield(W0(jj).time,'datenum')
                        timenum = W0(jj).time(kk).datenum(:);
                        data1   = W0(jj).(fields1{kk}).data(:);
                    else
                        timenum = W0(jj).time(:);
                        data1   = W0(jj).(fields1{kk})(:);
                    end
                else
                    timenum = W0(jj).(fields1{kk}).time(:);
                    data1   = W0(jj).(fields1{kk}).val(:);
                end
                IDnans=findmultipleNANs(data1,maximumgapsize);
                if strcmpi(fields1{kk},'th0') || strcmpi(fields1{kk},'udir') || strcmpi(fields1{kk},'xdir') || strcmpi(fields1{kk},'dir') 
                    W(jj).time(kk).datenum = interpNANs(timenum);
                    W(jj).(fields1{kk}).data = interpNANsDIR(data1);
                else
                    W(jj).time(kk).datenum = interpNANs(timenum);
                    W(jj).(fields1{kk}).data = interpNANs(data1);
                end
                W(jj).time(kk).datenum(IDnans)=nan;
                W(jj).(fields1{kk}).data(IDnans)=nan;
            else
                W(jj).time(kk).datenum = nan;
                W(jj).(fields1{kk}).data = nan;
            end
        end
        fprintf('#');
    end
    fprintf(']\n');
end

%-----------------------------------------------------------------------------------------
%% SUB-FUNCTION FIND MULTIPLE NANS
%-----------------------------------------------------------------------------------------
function [IDnans,nrNANs]=findmultipleNANs(X,numberofNANs)

    X=X(:);

    %% determine the number of consecutive NAN's
    % ------- replaced by a faster alternative ----------
    % if two nan's after each other the ID of both NAn's will be 2
    % tic
    % Xnan = isnan(X);
    % nrNANs=[];
    % cntr=0;
    % for ii=1:length(Xnan)
    %     if Xnan(ii)==1
    %         cntr=cntr+1;
    %         for ii1=ii:-1:ii-(cntr-1)
    %             nrNANs(ii1)=cntr;
    %         end
    %     else
    %         nrNANs(ii)=0;
    %         cntr=0;
    %     end
    % end
    % toc

    X1 = [0;X(:);0];
    nrNANs = zeros(size(X));
    IDXnotnan   = find(~isnan(X1));
    NANlength   = diff(IDXnotnan)-1;
    IDNANlength = find(NANlength~=0);
    for ii=1:length(IDNANlength)
        value = NANlength(IDNANlength(ii));
        IDrange = [IDXnotnan(IDNANlength(ii)):IDXnotnan(IDNANlength(ii)+1)-2]';
        nrNANs(IDrange)=value;
    end

    IDnans = find(nrNANs>numberofNANs);
end

%-----------------------------------------------------------------------------------------
%% SUB-FUNCTION INTERP-NANS
%-----------------------------------------------------------------------------------------
function [data2]=interpNANs(data,useextrapolation)
    % function [data2]=interpNANs(data)
    %
    % INPUT:
    %     data                Array [Nx1] with data
    %     useextrapolation    switch to use extrapolation (use 'nearest' or 'linear' for extrapolation)
    %                         this does not affect the interpolation which will be 'linear' 
    %
    % data = [10, 9, nan, 8, nan, 7.5, nan, nan]

    if nargin<2
        useextrapolation='nearest';  %<- use extrapolation if it is set to 1
    end
    rotate=0;
    if max(size(data,1),2)<size(data,2)
        rotate=1;
        data=data';
    end
    if size(data,2)==1
        x = [1:length(data)]';
        y = data;
    else
        x = data(:,1);
        y = data(:,2);
    end

    % make sure that x is increasing
    [x,ids]  = sort(x);
    y        = y(ids);
    eps      = 1e-6;
    for tt=1:length(x)-1
        if x(tt)==x(tt+1)
            x(tt+1)=x(tt+1)+eps;
        end
    end

    %% remove nans
    id = find(~isnan(y));
    x2 = x(id);
    y2 = y(id);
    if length(id)==length(x)
         ynew = y2;
    else
        if length(x2)>1
            if strcmpi(useextrapolation,'linear')
                ynew = interp1(x2,y2,x,'linear','extrap');
            else
                y3 = interp1(x2,y2,x,'linear');
                id = find(~isnan(y3));
                x2 = x(id);
                y2 = y3(id);
                ynew = interp1(x2,y2,x,'nearest','extrap');
            end
        elseif length(x2)==1
            ynew=repmat(y2,size(x));
        else
            %fprintf('Need more data points than 1\n')
            ynew=nan(size(x));
        end
    end

    if size(data,2)==2
        data2 = [x,ynew];
    else
        data2 = [ynew];
    end
    if rotate==1
        data2 = data2';
    end
end

%-----------------------------------------------------------------------------------------
%% SUB-FUNCTION INTERP-NANS DIR
%-----------------------------------------------------------------------------------------
function [data2]=interpNANsDIR(data)
    % function [data2]=interpNANsDIR(data)
    %
    % INPUT:
    %     data                Array [Nx1] with data
    %     useextrapolation    switch to use extrapolation (use 'nearest' or 'linear' for extrapolation)
    %                         this does not affect the interpolation which will be 'linear' 
    %
    % [data2]=interpNANsDIR([nan;nan;30;nan;330;300;nan;320;300;nan;nan;nan;20;nan;nan;350;30;nan;285]);

    rotate=0;
    if max(size(data,1),2)<size(data,2)
        rotate=1;
        data=data';
    end

    data1 = data;
    dir0 = data(find(~isnan(data),1));
    for ii=1:length(data)
        if ~isnan(data(ii));
            dir = data(ii);
            diroffset = dir0-mod(dir0,360);
            dir0B = mod(dir0+180,360);
            if dir<dir0B && abs(dir+diroffset-dir0)>180
                diroffset = diroffset+360;
            elseif dir>dir0B && abs(dir+diroffset-dir0)>180
                diroffset = diroffset-360;
            end
            data1(ii) = dir+diroffset;
            dir0 = data1(ii);
        end
    end

    data2 = mod(interpNANs(data1,'nearest'),360);

    if rotate==1
        data2 = data2';
    end
end

