clear all
close all
clc

runid = 'kz04';

frmt1 = ['%10.1f' repmat(' %10.3f',1,2) ' \n'];
frmt2 = ['%10.1f' repmat(' %10.3f',1,6) ' \n'];

% 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

% 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);

watlq1q4 = zeros(sum(numtijd(1:12))-11,1);
saliq1q4 = zeros(sum(numtijd(1:12))-11,numlaag);
tempq1q4 = zeros(sum(numtijd(1:12))-11,numlaag);

% coordinaten history stations dimensioneren
xstat = zeros(numstat,1);
ystat = zeros(numstat,1);


for m = 1:numstat
    % naam station weergeven
    namstat(m,1:20)

    % coordinaten history stations bepalen
    xstat(m) = xgrid(nstat(m),mstat(m));
    ystat(m) = ygrid(nstat(m),mstat(m));


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

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

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

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

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

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

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

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

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

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

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

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

    % analyse droogval
    ngelijk = 0;
    for t = 2:sum(numtijd(1:12))-11
        if (watlq1q4(t) == watlq1q4(t-1))
            ngelijk = ngelijk + 1;
        end
    end

    % analyse tijdseries en wegschrijven
    % dieptegemiddeld
    saliq1q4dg = saliq1q4*lagen;
    tempq1q4dg = tempq1q4*lagen;

    % wegschrijven tijdserie saliniteit oppervlakte, dieptegemiddeld
    % en bodem elk uur
    fid1 = fopen(['saltemp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_hist.txt'],'w');
    fprintf (fid1,['%s %s \n'], 'naam station', namstat(m,1:20));
    fprintf (fid1,['%s %10.3f %10.3f \n'], 'locatie_xy', xstat(m), ystat(m));
    fprintf (fid1,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    if (ngelijk >= 10)
        fprintf (fid1,['%s \n'], 'let op: station is drooggevallen');
    end
    for t = 1:sum(numtijd(1:12))-11
        fprintf (fid1,frmt2, 30*(t-1), ...
            saliq1q4(t,1), saliq1q4dg(t), saliq1q4(t,numlaag), ...
            tempq1q4(t,1), tempq1q4dg(t), tempq1q4(t,numlaag));
    end
    fclose(fid1);

    fid2 = fopen(['wl_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_hist.txt'],'w');
    fprintf (fid2,['%s %s \n'], 'naam station', namstat(m,1:20));
    fprintf (fid2,['%s %10.3f %10.3f \n'], 'locatie_xy', xstat(m), ystat(m));
    fprintf (fid2,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    if (ngelijk >= 10)
        fprintf (fid2,['%s \n'], 'let op: station is drooggevallen');
    end
    for t = 1:sum(numtijd(1:12))-11
        fprintf (fid2,['%10.1f %10.3f \n'], 30*(t-1), watlq1q4(t,1));
    end
    fclose(fid2);

    fid3 = fopen(['saltemp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_uur.txt'],'w');
    fprintf (fid3,['%s %s \n'], 'naam station', namstat(m,1:20));
    fprintf (fid3,['%s %10.3f %10.3f \n'], 'locatie_xy', xstat(m), ystat(m));
    fprintf (fid3,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    if (ngelijk >= 10)
        fprintf (fid3,['%s \n'], 'let op: station is drooggevallen');
    end
    for t = 1:2:sum(numtijd(1:12))-11
        fprintf (fid3,frmt2, 30*(t-1), ...
            saliq1q4(t,1), saliq1q4dg(t), saliq1q4(t,numlaag), ...
            tempq1q4(t,1), tempq1q4dg(t), tempq1q4(t,numlaag));
    end
    fclose(fid3);

    fid4 = fopen(['saltemp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_2tg.txt'],'w');
    fprintf (fid4,['%s %s \n'], 'naam station', namstat(m,1:20));
    fprintf (fid4,['%s %10.3f %10.3f \n'], 'locatie_xy', xstat(m), ystat(m));
    fprintf (fid4,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    if (ngelijk >= 10)
        fprintf (fid4,['%s \n'], 'let op: station is drooggevallen');
    end
    for t = 1:50:sum(numtijd(1:12))-11
        if (t+50<sum(numtijd(1:12))-11)
            fprintf (fid4,frmt2, 30*mean(t-1:t+49-1), ...
                mean(saliq1q4(t:t+49,1)), mean(saliq1q4dg(t:t+49)), mean(saliq1q4(t:t+49,numlaag)), ...
                mean(tempq1q4(t:t+49,1)), mean(tempq1q4dg(t:t+49)), mean(tempq1q4(t:t+49,numlaag)));
        end
    end
    fclose(fid4);

    fid5 = fopen(['saltemp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_2wk.txt'],'w');
    fprintf (fid5,['%s %s \n'], 'naam station', namstat(m,1:20));
    fprintf (fid5,['%s %10.3f %10.3f \n'], 'locatie_xy', xstat(m), ystat(m));
    fprintf (fid5,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    if (ngelijk >= 10)
        fprintf (fid5,['%s \n'], 'let op: station is drooggevallen');
    end
    for t = 1:672:sum(numtijd(1:12))-11
        if (t+672<sum(numtijd(1:12))-11)
            fprintf (fid5,frmt2, 30*mean(t-1:t+671-1), ...
                mean(saliq1q4(t:t+671,1)), mean(saliq1q4dg(t:t+671)), mean(saliq1q4(t:t+671,numlaag)),...
                mean(tempq1q4(t:t+671,1)), mean(tempq1q4dg(t:t+671)), mean(tempq1q4(t:t+671,numlaag)));
        end
    end
    fclose(fid5);

    fid6 = fopen(['saltemp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_hjr.txt'],'w');
    fprintf (fid6,['%s %s \n'], 'naam station', namstat(m,1:20));
    fprintf (fid6,['%s %10.3f %10.3f \n'], 'locatie_xy', xstat(m), ystat(m));
    fprintf (fid6,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    if (ngelijk >= 10)
        fprintf (fid6,['%s \n'], 'let op: station is drooggevallen');
    end
    fprintf (fid6,frmt2, ...
        30*mean(0:sum(numtijd(1:6))-5-1), ...
        mean(saliq1q4(1:sum(numtijd(1:6))-5,1)), ...
        mean(saliq1q4dg(1:sum(numtijd(1:6))-5)), ...
        mean(saliq1q4(1:sum(numtijd(1:6))-5,numlaag)), ...
        mean(tempq1q4(1:sum(numtijd(1:6))-5,1)), ...
        mean(tempq1q4dg(1:sum(numtijd(1:6))-5)), ...
        mean(tempq1q4(1:sum(numtijd(1:6))-5,numlaag)));
    fprintf (fid6,frmt2, ...
        30*mean(sum(numtijd(1:6))-5-1:sum(numtijd(1:12))-11-1), ...
        mean(saliq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,1)), ...
        mean(saliq1q4dg(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11)), ...
        mean(saliq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,numlaag)), ...
        mean(tempq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,1)), ...
        mean(tempq1q4dg(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11)), ...
        mean(tempq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,numlaag)));
    fclose(fid6);
    %    fid3 = fopen(['temp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_uur.txt'],'w');
    %    fclose(fid3);
end