function OUT = readOutput(folderName,fileName,includeInit)


if nargin==0
    folderName = pwd;
end

% determine the number of nodes per transverse line:
fid = fopen(fullfile(folderName,'SWMS_3D.IN','GRID.IN'),'rt');
tmp = textscan(fid,'%n',1,'headerlines',2);
nNodes = tmp{1};
fclose(fid);
clear tmp

% determine the number of nodes:
fid = fopen(fullfile(folderName,'SWMS_3D.IN','GRID.IN'),'rt');
tmp = textscan(fid,'%*n%*n%n',1,'headerlines',2);
nIJ = tmp{1};
fclose(fid);
clear tmp

% determine the number of print times:
fid = fopen(fullfile(folderName,'SWMS_3D.IN','SELECTOR.IN'),'rt');
fid = goToBlockStart(fid,'C');
fid = skipLines(fid,2);
tmp = textscan(fid,'%*f%*f%*f%*f%*f%n');
if includeInit
    nPrintTimes = tmp{1}+1;
else
    nPrintTimes = tmp{1};
end
fclose(fid);

% determine the print times:
fid = fopen(fullfile(folderName,'SWMS_3D.IN','SELECTOR.IN'),'rt');
fid = goToBlockStart(fid,'C');
fid = skipLines(fid,4);
tmp = textscan(fid,'%f',nPrintTimes);
if includeInit
    printTimes = union(tmp{1},0);
else
    printTimes = tmp{1};
end
fclose(fid);



formatStr = ['%n%f%f%f',repmat('%f',[1,nIJ])];
var = repmat(NaN,[ceil(nNodes/nIJ)*nIJ,nPrintTimes]);
iPrintTime = 1;
fid = fopen(fullfile(folderName,'SWMS_3D.OUT',fileName),'rt');
while true
    [fid,status] = goToTimeBlockStart(fid);
    if status == -1
        break
    end
    tmp = textscan(fid,'Time  ***      %f ***',1);
    
    if ismember(tmp{1},printTimes)
        varTimes(1,iPrintTime) = tmp{1};

        fid = skipLines(fid,3);
        tmp = textscan(fid,formatStr,ceil(nNodes/nIJ));

        if ~isempty(tmp{5})
            for iIJ = 1:nIJ
                var(iIJ:nIJ:ceil(nNodes/nIJ)*nIJ,iPrintTime) = tmp{:,4+iIJ};
            end

            if iPrintTime==nPrintTimes
                break
            else
                iPrintTime = iPrintTime+1;
            end
        else
            disp(['Encountered nans during t=',num2str(varTimes(1,iPrintTime)),' ',swmsTimeUnits,'.'])
        end
    else
        disp(['Skipping printed variable for t=',num2str(tmp{1})])
    end
    
end
fclose(fid);

OUT.values = var(1:nNodes,:);
OUT.times = varTimes;

switch fileName
    case 'Q.out'
        OUT.units = ['[',swmsLengthUnits,'^{3}.',swmsTimeUnits,'^{-1}]'];
    case 'h.out'
        OUT.units = ['[',swmsLengthUnits,']'];
    otherwise
        disp('no units assigned')
end
