clear all
close all
clc

runid = 'kz04';

% for m=1:12
%    eval(['runid' int2str(m) '= ''' runid 'm' sprintf('%0.2d',m) ''''])
% end

runid1 = [runid 'm01']; % kz04m01
runid2 = [runid 'm02']; % kz04m02
runid3 = [runid 'm03']; % kz04m03
runid4 = [runid 'm04']; % kz04m04
runid5 = [runid 'm05']; % kz04m05
runid6 = [runid 'm06']; % kz04m06
runid7 = [runid 'm07']; % kz04m07
runid8 = [runid 'm08']; % kz04m08
runid9 = [runid 'm09']; % kz04m09
runid10 = [runid 'm10']; % kz04m10
runid11 = [runid 'm11']; % kz04m11
runid12 = [runid 'm12']; % kz04m12

dt = 30; % tijdsinterval in history

% inlezen parameters stations
namstat = nc_varget(['SDS-' runid1 '_hist.nc'],'NAMPOL');
mstat   = nc_varget(['SDS-' runid1 '_hist.nc'],'MPOL');
nstat   = nc_varget(['SDS-' runid1 '_hist.nc'],'NPOL');
% inlezen parameters grid
xgrid   = nc_varget(['SDS-' runid1 '_hist.nc'],'XZETA');
ygrid   = nc_varget(['SDS-' runid1 '_hist.nc'],'YZETA');
lagen   = nc_varget(['SDS-' runid1 '_hist.nc'],'HLAY');

numtijd = zeros(12,1);
% uitlezen aantal stations en tijdstippen in q1+q2 en q3+q4
dimgro = nc_varsize(['SDS-' runid1 '_hist.nc'],'GRO');
numtijd(1) = dimgro(1);
numconc = dimgro(2);
numstat = dimgro(3);
numlaag = dimgro(4);

dimgro = nc_varsize(['SDS-' runid2 '_hist.nc'],'GRO');
numtijd(2) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid3 '_hist.nc'],'GRO');
numtijd(3) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid4 '_hist.nc'],'GRO');
numtijd(4) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid5 '_hist.nc'],'GRO');
numtijd(5) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid6 '_hist.nc'],'GRO');
numtijd(6) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid7 '_hist.nc'],'GRO');
numtijd(7) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid8 '_hist.nc'],'GRO');
numtijd(8) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid9 '_hist.nc'],'GRO');
numtijd(9) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid10 '_hist.nc'],'GRO');
numtijd(10) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid11 '_hist.nc'],'GRO');
numtijd(11) = dimgro(1);
dimgro = nc_varsize(['SDS-' runid12 '_hist.nc'],'GRO');
numtijd(12) = dimgro(1);

watlduur = zeros(sum(numtijd(1:12))-11,1);
saliduur = zeros(sum(numtijd(1:12))-11,1);

% coordinaten history stations dimensioneren
xstat = zeros(numstat,1);
ystat = zeros(numstat,1);
naam = namstat(1:numstat,1:20);
duur0 = zeros(numstat,8);
droog = zeros(numstat,1);

fid1 = fopen(['sal_duur_' runid '.txt'],'w');
fprintf (fid1,['%s \n'],'Maximale onderschrijdingsduur bodemsaliniteit voor gegeven saliniteitswaarde, in dagen; droog = 1 betekent punt drooggevallen');
fprintf (fid1,['%s \n'],'mstat nstat xstat ystat namstat duur_S<=5 duur_S<=10 duur_S<=15 duur_S<=20 duur_S<=25 duur_S<=30 duur_S<=35 duur_S<=40 droog');
sali0 = [5 10 15 20 25 30 35 40];
stat = 0;

for m = 1:numstat
    % naam station weergeven
    disp(namstat(m,1:20));
    statin = strcmp(namstat(m,1:5),'point');
    if (statin==1)
        stat = stat + 1;
        % coordinaten history stations bepalen
        xstat(stat) = xgrid(nstat(m),mstat(m));
        ystat(stat) = ygrid(nstat(m),mstat(m));
        naam(stat,1:20)=namstat(m,1:20);

        % inlezen tijdseries waterstand, saliniteit en temperatuur eerste
        % half jaar per maand
        watlduur(1:numtijd(1)) = nc_varget(['SDS-' runid1 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(1) 1]);
        saliduur(1:numtijd(1)) = nc_varget(['SDS-' runid1 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(1) 1 1 1]);

        watlduur(numtijd(1):sum(numtijd(1:2))-1) = nc_varget(['SDS-' runid2 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(2) 1]);
        saliduur(numtijd(1):sum(numtijd(1:2))-1) = nc_varget(['SDS-' runid2 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(2) 1 1 1]);

        watlduur(sum(numtijd(1:2))-1:sum(numtijd(1:3))-2) = nc_varget(['SDS-' runid3 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(3) 1]);
        saliduur(sum(numtijd(1:2))-1:sum(numtijd(1:3))-2) = nc_varget(['SDS-' runid3 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(3) 1 1 1]);

        watlduur(sum(numtijd(1:3))-2:sum(numtijd(1:4))-3) = nc_varget(['SDS-' runid4 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(4) 1]);
        saliduur(sum(numtijd(1:3))-2:sum(numtijd(1:4))-3) = nc_varget(['SDS-' runid4 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(4) 1 1 1]);

        watlduur(sum(numtijd(1:4))-3:sum(numtijd(1:5))-4) = nc_varget(['SDS-' runid5 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(5) 1]);
        saliduur(sum(numtijd(1:4))-3:sum(numtijd(1:5))-4) = nc_varget(['SDS-' runid5 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(5) 1 1 1]);

        watlduur(sum(numtijd(1:5))-4:sum(numtijd(1:6))-5) = nc_varget(['SDS-' runid6 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(6) 1]);
        saliduur(sum(numtijd(1:5))-4:sum(numtijd(1:6))-5) = nc_varget(['SDS-' runid6 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(6) 1 1 1]);

        watlduur(sum(numtijd(1:6))-5:sum(numtijd(1:7))-6) = nc_varget(['SDS-' runid7 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(7) 1]);
        saliduur(sum(numtijd(1:6))-5:sum(numtijd(1:7))-6) = nc_varget(['SDS-' runid7 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(7) 1 1 1]);

        watlduur(sum(numtijd(1:7))-6:sum(numtijd(1:8))-7) = nc_varget(['SDS-' runid8 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(8) 1]);
        saliduur(sum(numtijd(1:7))-6:sum(numtijd(1:8))-7) = nc_varget(['SDS-' runid8 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(8) 1 1 1]);

        watlduur(sum(numtijd(1:8))-7:sum(numtijd(1:9))-8) = nc_varget(['SDS-' runid9 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(9) 1]);
        saliduur(sum(numtijd(1:8))-7:sum(numtijd(1:9))-8) = nc_varget(['SDS-' runid9 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(9) 1 1 1]);

        watlduur(sum(numtijd(1:9))-8:sum(numtijd(1:10))-9) = nc_varget(['SDS-' runid10 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(10) 1]);
        saliduur(sum(numtijd(1:9))-8:sum(numtijd(1:10))-9) = nc_varget(['SDS-' runid10 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(10) 1 1 1]);

        watlduur(sum(numtijd(1:10))-9:sum(numtijd(1:11))-10) = nc_varget(['SDS-' runid11 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(11) 1]);
        saliduur(sum(numtijd(1:10))-9:sum(numtijd(1:11))-10) = nc_varget(['SDS-' runid11 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(11) 1 1 1]);

        watlduur(sum(numtijd(1:11))-10:sum(numtijd(1:12))-11) = nc_varget(['SDS-' runid12 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(12) 1]);
        saliduur(sum(numtijd(1:11))-10:sum(numtijd(1:12))-11) = nc_varget(['SDS-' runid12 '_hist.nc'],'GRO',[0 0 m-1 numlaag-1],[numtijd(12) 1 1 1]);

        % analyse droogval
        ngelijk = 0;
        duur(1:8) = 0;
        for t = 2:sum(numtijd(1:12))-11
            if (watlduur(t) == watlduur(t-1))
                ngelijk = ngelijk + 1;
            end
            for k = 1:8
                if (saliduur(t-1) < sali0(k) && saliduur(t) <= sali0(k))
                    duur(k) = duur(k) + dt;
                elseif (saliduur(t-1) < sali0(k) && saliduur(t) > sali0(k))
                    duur(k) = duur(k) + dt*(sali0(k)-saliduur(t-1))/(saliduur(t)-saliduur(t-1));
                    duur0(stat,k) = max(duur(k),duur0(stat,k));
                    duur(k) = 0;
                elseif (saliduur(t-1) == sali0(k) && saliduur(t) > sali0(k))
                    duur0(stat,k) = max(duur(k),duur0(stat,k));
                    duur(k) = 0;
                elseif (saliduur(t-1) == sali0(k) && saliduur(t) <= sali0(k))
                    duur(k) = duur(k) + dt;
                elseif (saliduur(t-1) > sali0(k) && saliduur(t) < sali0(k))
                    duur(k) = dt*(sali0(k)-saliduur(t))/(saliduur(t-1)-saliduur(t));
                    % anders niets
                end
            end
        end
        for k = 1:8
            duur0(stat,k) = max(duur(k),duur0(stat,k));
            duur0(stat,k) = duur0(stat,k)./1440.0;
        end
        %         duur0 = 100.*duur0./(sum(numtijd(1:12))-11-1)./dt;
        % analyse tijdseries en wegschrijven
        droog(stat) = (ngelijk >= 10);

        %        if (droog(m) == 1)
        % wegschrijven m,n x,y naam duur onderschrijving in dagen
        fprintf (fid1,['%4d %4d %10.3f %10.3f %s %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %2d\n'],...
            mstat(m), nstat(m), xstat(stat), ystat(stat), naam(stat,1:20), duur0(stat,1), duur0(stat,2), duur0(stat,3),...
            duur0(stat,4), duur0(stat,5), duur0(stat,6), duur0(stat,7), duur0(stat,8), droog(stat) );
        %        end

    end
end
fclose(fid1);

ncfile = ['sal_duur_' runid '.nc'];
nc_create_empty(ncfile);

attval = sprintf('created on %s', datestr(now));
nc_attput(ncfile,nc_global,'history',attval);
nc_attput(ncfile,nc_global,'author','J. Adema');
nc_attput(ncfile,nc_global,'description','Maximale onderschrijdingsduur bodemsaliniteit voor gegeven saliniteitswaarde, in dagen; droog = 1 betekent punt drooggevallen');

nc_adddim(ncfile,'stations',stat);
nc_adddim(ncfile,'naam',20);

varstruct.Name = 'x_station';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'x_station',xstat(1:stat));
clear varstruct

varstruct.Name = 'y_station';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'y_station',ystat(1:stat));
clear varstruct

varstruct.Name = 'naam_station';
varstruct.Nctype = 'char';
varstruct.Dimension = {'stations' 'naam'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'naam_station',naam(1:stat,1:20));
clear varstruct

varstruct.Name = 'duur_S_5';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_5',duur0(1:stat,1));
clear varstruct

varstruct.Name = 'duur_S_10';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_10',duur0(1:stat,2));
clear varstruct

varstruct.Name = 'duur_S_15';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_15',duur0(1:stat,3));
clear varstruct

varstruct.Name = 'duur_S_20';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_20',duur0(1:stat,4));
clear varstruct

varstruct.Name = 'duur_S_25';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_25',duur0(1:stat,5));
clear varstruct

varstruct.Name = 'duur_S_30';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_30',duur0(1:stat,6));
clear varstruct

varstruct.Name = 'duur_S_35';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_35',duur0(1:stat,7));
clear varstruct

varstruct.Name = 'duur_S_40';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'duur_S_40',duur0(1:stat,8));
clear varstruct

varstruct.Name = 'droog';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'droog',droog(1:stat));
clear varstruct