function s = readMKLtil2012(path)
%

% OPT.path2      = 'https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/kustlijnkaartenboek/';
OPT.path2        = path;
OPT.Filename2    = {'Amel'; 'Delfl'; 'Goeree'; 'NBev'; 'NH'; 'Rijnl'; 'Schier'; 'Schouw'; 'Tersch'; 'Texel'; 'Vlie'; 'Voorne'; 'Walch'; 'ZVla'};
OPT.areaname2    = {'Ameland'; 'Delfland'; 'Goeree'; 'Noord-Beveland'; 'Noord-Holland'; 'Rijnland'; 'Schiermonnikoog'; ...
    'Schouwen'; 'Terschelling'; 'Texel'; 'Vlieland'; 'Voorne'; 'Walcheren'; 'Zeeuws-Vlaanderen'};

OPT.ncpath      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
OPT.Filename    = 'transect.nc';
id       = ncread([OPT.ncpath OPT.Filename],'id');
cross    = ncread([OPT.ncpath OPT.Filename],'cross_shore');
along    = ncread([OPT.ncpath OPT.Filename],'alongshore');
areaname = ncread([OPT.ncpath OPT.Filename],'areaname');
x        = ncread([OPT.ncpath OPT.Filename],'x');
y        = ncread([OPT.ncpath OPT.Filename],'y');
time     = ncread([OPT.ncpath OPT.Filename],'time');
OPT.refdatenum   = datenum(1970,1,1);
jarkusyear = str2num(datestr(OPT.refdatenum+ncread([OPT.ncpath OPT.Filename],'time'),'yyyy')); %#ok<ST2NM>


for j = 1:length(OPT.Filename2)
fidMKL    = fopen(fullfile(OPT.path2, [char(OPT.Filename2(j)), ' - Kustlijn.txt']),'r');
disp(['reading MKL until 2012 - ' OPT.areaname2{j}])

% first lines
for i = 1:18; fgetl(fidMKL); end
i = 0;

% read all area information
while ~feof(fidMKL)
    info = fscanf(fidMKL, '%f', [1 6]);
    
    if ~(isempty(info))
        if ~(length(info)==4)     % ** values
            if ~(info(3) == 0 || info(4) == 0)     % issues w/ day = 0, month = 0
                
                i    = i + 1;
                
                area.raai(i)     = info(1);
                area.year(i)     = info(2);
                area.time_MKL(i) = datenum([num2str(info(4)) '-' num2str(info(3)) '-' ...
                    num2str(info(2))],'dd-mm-yyyy') - OPT.refdatenum;
                area.hoog(i)     = info(5);
                area.laag(i)     = info(6);

                area.XzeeMKL(i)  = str2double(fscanf(fidMKL, '%s', [1 1]));
                area.Xmkl(i)     = str2double(fscanf(fidMKL, '%s', [1 1]));
                
            else
                fscanf(fidMKL, '%s', [1 1]);
                fscanf(fidMKL, '%s', [1 1]);
            end
        else
            fscanf(fidMKL, '%s', [1 1]);
            fscanf(fidMKL, '%s', [1 1]);
            fscanf(fidMKL, '%s', [1 1]);
            fscanf(fidMKL, '%s', [1 1]);
        end
        
    end
end

area.years  = unique(area.year);
area.raaien = unique(area.raai);

for yj = 1:length(area.years) % from 1965 to current year
    [true_year, loc_year] = ismember(area.years(yj),area.year);
    
    for r = 1:length(area.raaien)
        [true_raai, loc_raai] = ismember(area.raaien(r),area.raai);
        
        if sum(area.year==area.years(yj) & area.raai==area.raaien(r)) == 1 % issues w/ values for the same year
            area.tableXmkl(yj,r)     = area.Xmkl(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
            area.tableXzeeMKL(yj,r)  = area.XzeeMKL(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
            area.tablehoog(yj,r)     = area.hoog(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
            area.tablelaag(yj,r)     = area.laag(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
            area.tabletime_MKL(yj,r) = area.time_MKL(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
            area.tableraai(yj,r)     = area.raai(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
        else
            area.tableXmkl(yj,r)     = NaN;
            area.tableXzeeMKL(yj,r)  = NaN;
            area.tablehoog(yj,r)     = NaN;
            area.tablelaag(yj,r)     = NaN;
            area.tabletime_MKL(yj,r) = NaN;
            area.tableraai(yj,r)     = NaN;
        end
    end
    
end

tableXmkl(j)     = {area.tableXmkl};
tableXzeeMKL(j)  = {area.tableXzeeMKL};
tablehoog(j)     = {area.tablehoog};
tablelaag(j)     = {area.tablelaag};
tabletime_MKL(j) = {area.tabletime_MKL};
tableraai(j)     = {area.tableraai};

clear area
fclose(fidMKL);

end

%% create the matrix for each parameter, with a structure similar to the jarkus file
disp('creating MKL until 2012 matrix as in Jarkus file ')

% define the final variables
hoog    = NaN(length(time),length(id)); laag = NaN(length(time),length(id)); 
XzeeMKL = NaN(length(time),length(id)); Xmkl = NaN(length(time),length(id));
timeMKL = NaN(length(time),length(id));

areaname = areaname'; % due to an erroneous Jarkus version
for i = 1:length(id)
    dum_area = ismember(OPT.areaname2,cellstr(areaname(i,:)));
    if any(dum_area)
    % dum_id = ismember(along(i),tableraai{ismember(OPT.areaname2,cellstr(areaname(i,:)))});
    dum_id = ismember(tableraai{dum_area},along(i));
    
        for t = 1:length(time)
            dum_time = tabletime_MKL{dum_area};
            [dum_t, index] = ismember(jarkusyear(t),str2num(datestr((dum_time(:,any(dum_id))+OPT.refdatenum),'yyyy'))); %#ok<ST2NM>
    %         if ismember(cellstr(areaname(i,:)),OPT.areaname2) && ...
    %                  any(any(dum_id))  
            if dum_t
               Xmkl(t,i)    = tableXmkl{dum_area}(index,any(dum_id));
               XzeeMKL(t,i) = tableXzeeMKL{dum_area}(index,any(dum_id));
               hoog(t,i)    = tablehoog{dum_area}(index,any(dum_id));
               laag(t,i)    = tablelaag{dum_area}(index,any(dum_id));
               timeMKL(t,i) = tabletime_MKL{dum_area}(index,any(dum_id));

            end

        end
    end

end

%% convert Xmkl in RD coordinates and in WGS84
disp('creating MKL until 2012 x, y, lat, lon ')
% MKL in RD coordinates
MKL_x = NaN(length(time),length(id)); MKL_y = NaN(length(time),length(id));

for i = 1:length(id)
    for t = 1:length(time)
        if ~isnan(Xmkl(t,i))
            MKL_x(t,i) = interp1(cross', x(:,i), Xmkl(t,i));
            MKL_y(t,i) = interp1(cross', y(:,i), Xmkl(t,i));
        else
            MKL_x(t,i) = NaN;
            MKL_y(t,i) = NaN;
        end
    end
end

% MKL in WGS84 coordinates
[MKL_lon, MKL_lat] = convertCoordinates(MKL_x,MKL_y,'CS1.code',28992,'CS2.code',4326);

MKL_Vol = NaN(size(Xmkl)); % not defined in the files, but defined in next versions

s.id = id;
s.areaname = areaname;
% s.along    = along;
s.time     = time;

s.MKL_lat = MKL_lat;
s.MKL_lon = MKL_lon;
s.MKL_x   = MKL_x;
s.MKL_y   = MKL_y;

s.Xmkl    = Xmkl;
s.Volume  = MKL_Vol;

s.timeMKL = timeMKL;
s.hoog    = hoog;
s.laag    = laag;
s.XzeeMKL = XzeeMKL;




