function s = readMKL2016(path)
%

% jarkus database for transect info
OPT.ncpath      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
OPT.Filename   = 'transect.nc';

% general settings
OPT.refdatenum   = datenum(1970,1,1);

% MKL
% OPT.ncpath2      = 'https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/kustlijnkaartenboek/raw/MKL2013/';
OPT.ncpath2      = path; % local
OPT.Filename2    = {'03_BKL-MKL2016-TKL2017-trend export 2016 Ameland'; '09_BKL-MKL2016-TKL2017-trend export 2016 Delfland'; '12_BKL-MKL2016-TKL2017-trend export 2016 Goeree'; '15_BKL-MKL2016-TKL2017-trend export 2016 Noord-Beveland'; '07_BKL-MKL2016-TKL2017-trend export 2016 Noord-Holland'; '08_BKL-MKL2016-TKL2017-trend export 2016 Rijnland'; '02_BKL-MKL2016-TKL2017-trend export 2016 Schiermonnikoog'; '13_BKL-MKL2016-TKL2017-trend export 2016 Schouwen'; '04_BKL-MKL2016-TKL2017-trend export 2016 Terschelling'; '06_BKL-MKL2016-TKL2017-trend export 2016 Texel'; '05_BKL-MKL2016-TKL2017-trend export 2016 Vlieland'; '11_BKL-MKL2016-TKL2017-trend export 2016 Voorne'; '16_BKL-MKL2016-TKL2017-trend export 2016 Walcheren'; '17_BKL-MKL2016-TKL2017-trend export 2016 Zeeuws-Vlaanderen'};
OPT.areaname2    = {'Ameland'; 'Delfland'; 'Goeree'; 'Noord-Beveland'; 'Noord-Holland'; 'Rijnland'; 'Schiermonnikoog'; ...
    'Schouwen'; 'Terschelling'; 'Texel'; 'Vlieland'; 'Voorne'; 'Walcheren'; 'Zeeuws-Vlaanderen'};

% general settings
OPT.refdatenum   = datenum(1970,1,1);

%% Read data form jarkus 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');
time     = ncread([OPT.ncpath OPT.Filename],'time');

jarkusyear = str2num(datestr(OPT.refdatenum+ncread([OPT.ncpath OPT.Filename],'time'),'yyyy')); %#ok<ST2NM>

x        = ncread([OPT.ncpath OPT.Filename],'x');
y        = ncread([OPT.ncpath OPT.Filename],'y');

%% Read MKL, time_MKL, [High, Low, x_sea (boundaries)]
% cycle to read each area
for j = 1:length(OPT.Filename2)
    disp(['reading MKL 2016 - ' OPT.areaname2{j}])
    % fidMKL    = fopen([OPT.ncpath2, char(OPT.Filename2(j)), '.csv'],'r');
    % first line
    % for i = 1:1; fgetl(fidMKL); end
    % i = 0;
    % read all area information
    % while ~feof(fidMKL)
    %     info = fscanf(fidMKL, '%f,%f,%s,%f,%f', [1 5]);
    xlsxfile = fullfile(OPT.ncpath2, [char(OPT.Filename2(j)), '.xlsx']);
    arraytest      = xlsread( xlsxfile, 'MKL2016' );
    [nRowst,~] = size( arraytest );
    M2 = xlsread(xlsxfile, 'MKL2016', ['A1:I', num2str(nRowst+1)],'basic');
    
    [~,DatesXLS] = xlsread(xlsxfile, 'MKL2016', ['C2:C', num2str(nRowst+1)]);
    
    for ddd = 1:length(M2)
        % Dates(ddd-1) = datenum(DatesXLS(ddd,3),'dd-mmm-yy');
        Dates(ddd) = datenum(DatesXLS(ddd),'dd-mmm-yy');
    end
    
    % Argh. Same here.
    % M1 = xlsread([OPT.ncpath2, char(OPT.Filename2(j)), '.csv'],'A:B');
    M1_1 = M2(:,1);
    M1_2 = M2(:,2);
    M1 = [M1_1, M1_2];
    
    %     if ~isempty(info)
    %         if ~(info(3) == 0 || info(4) == 0)     % issues w/ day = 0, month = 0
    
    %             i    = i + 1;
    
    % area.raai     = M1(:,1)*1e6 + M1(:,2);
    area.raai     = M1(:,2);
    area.year     = str2num(datestr(Dates,'yyyy'));
    area.time_MKL = Dates - OPT.refdatenum;
    area.hoog     = M2(:,6);
    area.laag     = M2(:,7);
    
    area.XMaxSea  = M2(:,4);
    area.XMinLand = M2(:,5);
    
    area.XMkl     = M2(:,8);
    area.Volume   = M2(:,9);
    
    % avoid repetitions
    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.tableXMaxSea(yj,r)  = area.XMaxSea(find(area.year==area.years(yj) & area.raai==area.raaien(r)));
                area.tableXMinLand(yj,r) = area.XMinLand(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.tableVolume(yj,r)   = area.Volume(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.tableXMaxSea(yj,r)  = NaN;
                area.tableXMinLand(yj,r) = NaN;
                area.tablehoog(yj,r)     = NaN;
                area.tablelaag(yj,r)     = NaN;
                area.tabletime_MKL(yj,r) = NaN;
                area.tableVolume(yj,r)   = NaN;
                area.tableraai(yj,r)     = NaN;
            end
        end
        
    end
    
    tableXMkl(j)     = {area.tableXMkl};
    tableXMaxSea(j)  = {area.tableXMaxSea};
    tableXMinLand(j) = {area.tableXMinLand};
    tablehoog(j)     = {area.tablehoog};
    tablelaag(j)     = {area.tablelaag};
    tabletime_MKL(j) = {area.tabletime_MKL};
    tableVolume(j)   = {area.tableVolume};
    tableraai(j)     = {area.tableraai};
    
    clear area
    clear Dates
    
end

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

% define the final variables
hoog    = NaN(length(time),length(id)); laag = NaN(length(time),length(id));
XMaxSea = NaN(length(time),length(id)); XMinLand = NaN(length(time),length(id));
XMkl = NaN(length(time),length(id)); Volume = 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(:,dum_id)+OPT.refdatenum),'yyyy'))); %#ok<ST2NM>

            if dum_t
                XMkl(t,i)     = tableXMkl{dum_area}(index,dum_id);
                XMaxSea(t,i)  = tableXMaxSea{dum_area}(index,dum_id);
                XMinLand(t,i) = tableXMinLand{dum_area}(index,dum_id);
                hoog(t,i)     = tablehoog{dum_area}(index,dum_id);
                laag(t,i)     = tablelaag{dum_area}(index,dum_id);
                Volume(t,i)   = tableVolume{dum_area}(index,dum_id);
                timeMKL(t,i)  = tabletime_MKL{dum_area}(index,dum_id);
                
            end
            
        end
    end
    
end

%% convert Xmkl in RD coordinates and in WGS84
disp('creating MKL 2016 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)
    disp(['interpolate ' num2str(i) '/' num2str(length(id))])
    for t = 1:length(time)
        MKL_x(t,i) = interp1(cross', x(:,i), XMkl(t,i));
        MKL_y(t,i) = interp1(cross', y(:,i), XMkl(t,i));
    end
end

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

s.id = id;
s.areaname = areaname;
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  = Volume;

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


