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

watlq1q2 = zeros(sum(numtijd(1:6))-5,1);
saliq1q2 = zeros(sum(numtijd(1:6))-5,numlaag);
tempq1q2 = zeros(sum(numtijd(1:6))-5,numlaag);

% coordinaten history stations dimensioneren
xstat = zeros(numstat,1);
ystat = zeros(numstat,1);
naam = namstat(1:numstat,1:20);
sal_opp = zeros(numstat,4);
sal_dpg = zeros(numstat,4);
sal_bod = zeros(numstat,4);
temp_opp = zeros(numstat,4);
temp_dpg = zeros(numstat,4);
temp_bod = zeros(numstat,4);
droog = zeros(numstat,1);

fid1 = fopen(['sal_q1q2_' runid '.txt'],'w');
fid2 = fopen(['temp_q1q2_' runid '.txt'],'w');

fprintf (fid1,['%s \n'],'Statistische grootheden (gemiddelde, standaarddeviatie, minimum en maximum) saliniteit (oppervlak, dieptegemiddeld en bodem) over eerste 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 eerste 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');

stat = 0;

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
        stat = stat + 1;
        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
        watlq1q2(1:numtijd(1)) = nc_varget(['SDS-' runid1 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(1) 1]);
        saliq1q2(1:numtijd(1),:) = nc_varget(['SDS-' runid1 '_hist.nc'],'GRO',[0 0 m-1 0],[numtijd(1) 1 1 numlaag]);
        tempq1q2(1:numtijd(1),:) = nc_varget(['SDS-' runid1 '_hist.nc'],'GRO',[0 1 m-1 0],[numtijd(1) 1 1 numlaag]);

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

        watlq1q2(sum(numtijd(1:2))-1:sum(numtijd(1:3))-2) = nc_varget(['SDS-' runid3 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(3) 1]);
        saliq1q2(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]);
        tempq1q2(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]);

        watlq1q2(sum(numtijd(1:3))-2:sum(numtijd(1:4))-3) = nc_varget(['SDS-' runid4 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(4) 1]);
        saliq1q2(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]);
        tempq1q2(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]);

        watlq1q2(sum(numtijd(1:4))-3:sum(numtijd(1:5))-4) = nc_varget(['SDS-' runid5 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(5) 1]);
        saliq1q2(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]);
        tempq1q2(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]);

        watlq1q2(sum(numtijd(1:5))-4:sum(numtijd(1:6))-5) = nc_varget(['SDS-' runid6 '_hist.nc'],'ZWL',[0 m-1],    [numtijd(6) 1]);
        saliq1q2(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]);
        tempq1q2(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]);

        % analyse droogval
        ngelijk = 0;
        for t = 2:sum(numtijd(1:6))-5
            if (watlq1q2(t) == watlq1q2(t-1))
                ngelijk = ngelijk + 1;
            end
        end

        % analyse tijdseries en wegschrijven

        droog(stat) = (ngelijk >= 10);
        %        if (ngelijk < 10)
        % dieptegemiddeld
        saliq1q2dg = saliq1q2*lagen;
        tempq1q2dg = tempq1q2*lagen;

        sal_opp(stat,1) = mean(saliq1q2(:,1));
        sal_opp(stat,2) = std(saliq1q2(:,1),1);
        sal_opp(stat,3) = min(saliq1q2(:,1));
        sal_opp(stat,4) = max(saliq1q2(:,1));
        sal_dpg(stat,1) = mean(saliq1q2dg);
        sal_dpg(stat,2) = std(saliq1q2dg,1);
        sal_dpg(stat,3) = min(saliq1q2dg);
        sal_dpg(stat,4) = max(saliq1q2dg);
        sal_bod(stat,1) = mean(saliq1q2(:,numlaag));
        sal_bod(stat,2) = std(saliq1q2(:,numlaag),1);
        sal_bod(stat,3) = min(saliq1q2(:,numlaag));
        sal_bod(stat,4) = max(saliq1q2(:,numlaag));

        temp_opp(stat,1) = mean(tempq1q2(:,1));
        temp_opp(stat,2) = std(tempq1q2(:,1),1);
        temp_opp(stat,3) = min(tempq1q2(:,1));
        temp_opp(stat,4) = max(tempq1q2(:,1));
        temp_dpg(stat,1) = mean(tempq1q2dg);
        temp_dpg(stat,2) = std(tempq1q2dg,1);
        temp_dpg(stat,3) = min(tempq1q2dg);
        temp_dpg(stat,4) = max(tempq1q2dg);
        temp_bod(stat,1) = mean(tempq1q2(:,numlaag));
        temp_bod(stat,2) = std(tempq1q2(:,numlaag),1);
        temp_bod(stat,3) = min(tempq1q2(:,numlaag));
        temp_bod(stat,4) = max(tempq1q2(:,numlaag));

        % 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 %3d \n'],...
            mstat(m), nstat(m), xstat(m), ystat(m), namstat(m,1:20),...
            sal_opp(stat,1), sal_opp(stat,2), sal_opp(stat,3), sal_opp(stat,4),...
            sal_dpg(stat,1), sal_dpg(stat,2), sal_dpg(stat,3), sal_dpg(stat,4),...
            sal_bod(stat,1), sal_bod(stat,2), sal_bod(stat,3), sal_bod(stat,4), droog(stat));

        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 %3d \n'],...
            mstat(m), nstat(m), xstat(m), ystat(m), namstat(m,1:20),...
            temp_opp(stat,1), temp_opp(stat,2), temp_opp(stat,3), temp_opp(stat,4),...
            temp_dpg(stat,1), temp_dpg(stat,2), temp_dpg(stat,3), temp_dpg(stat,4),...
            temp_bod(stat,1), temp_bod(stat,2), temp_bod(stat,3), temp_bod(stat,4), droog(stat));

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

ncfile = ['sal_q1q2_' 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','Statistische grootheden (gemiddelde, standaarddeviatie, minimum en maximum) saliniteit (oppervlak, dieptegemiddeld en bodem) over eerste half jaar 2004-2005; 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 = 'S_opp_gem';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'S_opp_gem',sal_opp(1:stat,1));
clear varstruct

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

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

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

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

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

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

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

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

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

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

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

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

ncfile = ['temp_q1q2_' 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','Statistische grootheden (gemiddelde, standaarddeviatie, minimum en maximum) temperatuur (oppervlak, dieptegemiddeld en bodem) over eerste half jaar 2004-2005; 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 = 'T_opp_gem';
varstruct.Nctype = 'double';
varstruct.Dimension = {'stations'};
nc_addvar(ncfile,varstruct);
nc_varput(ncfile,'T_opp_gem',temp_opp(1:stat,1));
clear varstruct

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

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

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

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

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

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

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

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

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

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

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

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