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

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

watlq3q4 = zeros(sum(numtijd(7:12))-5,1);
saliq3q4 = zeros(sum(numtijd(7:12))-5,numlaag);
tempq3q4 = zeros(sum(numtijd(7:12))-5,numlaag);

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

fid1 = fopen(['sal_q3q4_' runid '.txt'],'w');
fid2 = fopen(['temp_q3q4_' runid '.txt'],'w');

fprintf (fid1,['%s \n'],'Statistische grootheden (gemiddelde, standaarddeviatie, minimum en maximum) saliniteit (oppervlak, dieptegemiddeld en bodem) over tweede half jaar 2004-2005; droog = 1 betekent punt drooggevallen');
fprintf (fid2,['%s \n'],'Statistische grootheden (gemiddelde, standaarddeviatie, minimum en maximum) temperatuur (oppervlak, dieptegemiddeld en bodem) over tweede half jaar 2004-2005; droog = 1 betekent punt drooggevallen');
fprintf (fid1,['%s \n'],'mstat nstat xstat ystat namstat Sopp_gem Sopp_std Sopp_min Sopp_max Sdpg_gem Sdpg_std Sdpg_min Sdpg_max Sbod_gem Sbod_std Sbod_min Sbod_max');
fprintf (fid2,['%s \n'],'mstat nstat xstat ystat namstat Topp_gem Topp_std Topp_min Topp_max Tdpg_gem Tdpg_std Tdpg_min Tdpg_max Tbod_gem Tbod_std Tbod_min Tbod_max');

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

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

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

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

        watlq3q4(sum(numtijd(7:9))-2:sum(numtijd(7:10))-3) = nc_varget(['SDS-' runid10 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(10) 1]);
        saliq3q4(sum(numtijd(7:9))-2:sum(numtijd(7:10))-3,:) = nc_varget(['SDS-' runid10 '_hist.nc'],'GRO',[0 0 m-1 0],[numtijd(10) 1 1 numlaag]);
        tempq3q4(sum(numtijd(7:9))-2:sum(numtijd(7:10))-3,:) = nc_varget(['SDS-' runid10 '_hist.nc'],'GRO',[0 1 m-1 0],[numtijd(10) 1 1 numlaag]);

        watlq3q4(sum(numtijd(7:10))-3:sum(numtijd(7:11))-4) = nc_varget(['SDS-' runid11 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(11) 1]);
        saliq3q4(sum(numtijd(7:10))-3:sum(numtijd(7:11))-4,:) = nc_varget(['SDS-' runid11 '_hist.nc'],'GRO',[0 0 m-1 0],[numtijd(11) 1 1 numlaag]);
        tempq3q4(sum(numtijd(7:10))-3:sum(numtijd(7:11))-4,:) = nc_varget(['SDS-' runid11 '_hist.nc'],'GRO',[0 1 m-1 0],[numtijd(11) 1 1 numlaag]);

        watlq3q4(sum(numtijd(7:11))-4:sum(numtijd(7:12))-5) = nc_varget(['SDS-' runid12 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(12) 1]);
        saliq3q4(sum(numtijd(7:11))-4:sum(numtijd(7:12))-5,:) = nc_varget(['SDS-' runid12 '_hist.nc'],'GRO',[0 0 m-1 0],[numtijd(12) 1 1 numlaag]);
        tempq3q4(sum(numtijd(7:11))-4:sum(numtijd(7:12))-5,:) = 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(7:12))-5
            if (watlq3q4(t) == watlq3q4(t-1))
                ngelijk = ngelijk + 1;
            end
        end

        % analyse tijdseries en wegschrijven
        if (ngelijk < 10)
            % dieptegemiddeld
            saliq3q4dg = saliq3q4*lagen;
            tempq3q4dg = tempq3q4*lagen;

            % wegschrijven m,n x,y naam (gem,stdev,min,max)_opp (gem,stdev,min,max)_dpgem (gem,stdev,min,max)_bod
            fprintf (fid1,['%4d %4d %10.3f %10.3f %s ',...
                '%10.3f %10.3f %10.3f %10.3f ',...
                '%10.3f %10.3f %10.3f %10.3f ',...
                '%10.3f %10.3f %10.3f %10.3f \n'],...
                mstat(m), nstat(m), xstat(m), ystat(m), namstat(m,1:20),...
                mean(saliq3q4(:,1)), std(saliq3q4(:,1),1), min(saliq3q4(:,1)), max(saliq3q4(:,1)),...
                mean(saliq3q4dg), std(saliq3q4dg,1), min(saliq3q4dg), max(saliq3q4dg),...
                mean(saliq3q4(:,numlaag)), std(saliq3q4(:,numlaag),1), min(saliq3q4(:,numlaag)), max(saliq3q4(:,numlaag)));

            fprintf (fid2,['%4d %4d %10.3f %10.3f %s ',...
                '%10.3f %10.3f %10.3f %10.3f ',...
                '%10.3f %10.3f %10.3f %10.3f ',...
                '%10.3f %10.3f %10.3f %10.3f \n'],...
                mstat(m), nstat(m), xstat(m), ystat(m), namstat(m,1:20),...
                mean(tempq3q4(:,1)), std(tempq3q4(:,1),1), min(tempq3q4(:,1)), max(tempq3q4(:,1)),...
                mean(tempq3q4dg), std(tempq3q4dg,1), min(tempq3q4dg), max(tempq3q4dg),...
                mean(tempq3q4(:,numlaag)), std(tempq3q4(:,numlaag),1), min(tempq3q4(:,numlaag)), max(tempq3q4(:,numlaag)));

        end
    end
end
fclose(fid1);
fclose(fid2);