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);
droog = zeros(numstat,1);

for m = [1865 4158 2012 4355 1995 5143 6977 1809 4310 1465 3189 1203 2901 1188 2260 1281 2580 1166 2723 4692 1147 1888]+49-1000 % numstat
    % naam station weergeven
    namstat(m,1:20)

    % coordinaten history stations bepalen
    xstat = xgrid(nstat(m),mstat(m));
    ystat = 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 tijdseries saliniteit en temperatuur (oppervlakte, dieptegemiddeld
    % en bodem), elk half 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, ystat);
    fprintf (fid1,['%s \n'], 'Tijdseries (half uur) S_opp, S_dpg, S_bod, T_opp, T_dpg, T_bod, tijd in minuten na 1-jun-2004');

%    fprintf (fid1,['%s %4d %4d \n'], 'locatie_mn', mstat(m), nstat(m));
    droog(m) = (ngelijk >= 10);
    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);

    ncfile = ['saltemp_' runid '_' strrep(strtrim(namstat(m,1:20)),' ','') '_hist.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','Tijdseries (half uur en uur) en tijdgemiddelden (over 2 getijperioden, 2 weken en een half jaar) S_opp, S_dpg, S_bod, T_opp, T_dpg, T_bod, tijd in minuten na 1-jun-2004');

    tijden_hfu = 30.*((1:sum(numtijd(1:12))-11)-1);
    [a nrhfu] = size(tijden_hfu);
    nc_adddim(ncfile,'series_hfu',nrhfu);

    tijden_uur = 30.*((1:2:sum(numtijd(1:12))-11)-1);
    [a nruur] = size(tijden_uur);
    nc_adddim(ncfile,'series_uur',nruur);

    nr2tg = 0;
    for t = 1:50:sum(numtijd(1:12))-11
        if (t+50<sum(numtijd(1:12))-11)
            nr2tg = nr2tg + 1;
            tijden_2tg(nr2tg) = 30*mean(t-1:t+49-1);
            sal_2tg(nr2tg,1) = mean(saliq1q4(t:t+49,1));
            sal_2tg(nr2tg,2) = mean(saliq1q4dg(t:t+49));
            sal_2tg(nr2tg,3) = mean(saliq1q4(t:t+49,numlaag));
            tmp_2tg(nr2tg,1) = mean(tempq1q4(t:t+49,1));
            tmp_2tg(nr2tg,2) = mean(tempq1q4dg(t:t+49));
            tmp_2tg(nr2tg,3) = mean(tempq1q4(t:t+49,numlaag));
        end
    end
    nc_adddim(ncfile,'series_2tg',nr2tg);

    nr2wk = 0;
    for t = 1:672:sum(numtijd(1:12))-11
        if (t+672<sum(numtijd(1:12))-11)
            nr2wk = nr2wk + 1;
            tijden_2wk(nr2wk) = 30*mean(t-1:t+671-1);
            sal_2wk(nr2wk,1) = mean(saliq1q4(t:t+671,1));
            sal_2wk(nr2wk,2) = mean(saliq1q4dg(t:t+671));
            sal_2wk(nr2wk,3) = mean(saliq1q4(t:t+671,numlaag));
            tmp_2wk(nr2wk,1) = mean(tempq1q4(t:t+671,1));
            tmp_2wk(nr2wk,2) = mean(tempq1q4dg(t:t+671));
            tmp_2wk(nr2wk,3) = mean(tempq1q4(t:t+671,numlaag));
        end
    end
    nc_adddim(ncfile,'series_2wk',nr2wk);

    nrhjr = 2;
    tijden_hjr(1) = 30*mean(0:sum(numtijd(1:6))-5-1);
    sal_hjr(1,1) = mean(saliq1q4(1:sum(numtijd(1:6))-5,1));
    sal_hjr(1,2) = mean(saliq1q4dg(1:sum(numtijd(1:6))-5));
    sal_hjr(1,3) = mean(saliq1q4(1:sum(numtijd(1:6))-5,numlaag));
    tmp_hjr(1,1) = mean(tempq1q4(1:sum(numtijd(1:6))-5,1));
    tmp_hjr(1,2) = mean(tempq1q4dg(1:sum(numtijd(1:6))-5));
    tmp_hjr(1,3) = mean(tempq1q4(1:sum(numtijd(1:6))-5,numlaag));

    tijden_hjr(2) = 30*mean(sum(numtijd(1:6))-5-1:sum(numtijd(1:12))-11-1);
    sal_hjr(2,1) = mean(saliq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,1));
    sal_hjr(2,2) = mean(saliq1q4dg(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11));
    sal_hjr(2,3) = mean(saliq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,numlaag));
    tmp_hjr(2,1) = mean(tempq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,1));
    tmp_hjr(2,2) = mean(tempq1q4dg(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11));
    tmp_hjr(2,3) = mean(tempq1q4(sum(numtijd(1:6))-5:sum(numtijd(1:12))-11,numlaag));

    nc_adddim(ncfile,'series_hjr',nrhjr);

    nc_adddim(ncfile,'een',1);
    nc_adddim(ncfile,'naam',20);

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

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

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

    varstruct.Name = 'droog';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'een'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'droog',droog(m));
    clear varstruct

    varstruct.Name = 'tijd_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'tijd_half_uur',tijden_hfu);
    clear varstruct

    varstruct.Name = 'S_opp_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_opp_half_uur',saliq1q4(:,1));
    clear varstruct

    varstruct.Name = 'S_dpg_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_dpg_half_uur',saliq1q4dg(:));
    clear varstruct

    varstruct.Name = 'S_bod_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_bod_half_uur',saliq1q4(:,numlaag));
    clear varstruct

    varstruct.Name = 'T_opp_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_opp_half_uur',tempq1q4(:,1));
    clear varstruct

    varstruct.Name = 'T_dpg_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_dpg_half_uur',tempq1q4dg(:));
    clear varstruct

    varstruct.Name = 'T_bod_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_bod_half_uur',tempq1q4(:,numlaag));
    clear varstruct

    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, ystat);
    fprintf (fid2,['%s \n'], 'Tijdseries (half uur) waterstand, tijd in minuten na 1-jun-2004');
%    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);

    varstruct.Name = 'waterstand_half_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hfu'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'waterstand_half_uur',watlq1q4);
    clear varstruct

    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, ystat);
    fprintf (fid3,['%s \n'], 'Tijdseries (uur) S_opp, S_dpg, S_bod, T_opp, T_dpg, T_bod, tijd in minuten na 1-jun-2004');
%    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);

    varstruct.Name = 'tijd_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'tijd_uur',tijden_uur);
    clear varstruct

    varstruct.Name = 'S_opp_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_opp_uur',saliq1q4(1:2:sum(numtijd(1:12))-11,1));
    clear varstruct

    varstruct.Name = 'S_dpg_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_dpg_uur',saliq1q4dg(1:2:sum(numtijd(1:12))-11));
    clear varstruct

    varstruct.Name = 'S_bod_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_bod_uur',saliq1q4(1:2:sum(numtijd(1:12))-11,numlaag));
    clear varstruct

    varstruct.Name = 'T_opp_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_opp_uur',tempq1q4(1:2:sum(numtijd(1:12))-11,1));
    clear varstruct

    varstruct.Name = 'T_dpg_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_dpg_uur',tempq1q4dg(1:2:sum(numtijd(1:12))-11));
    clear varstruct

    varstruct.Name = 'T_bod_uur';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_uur'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_bod_uur',tempq1q4(1:2:sum(numtijd(1:12))-11,numlaag));
    clear varstruct

    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, ystat);
    fprintf (fid4,['%s \n'], 'Tijdgemiddelden (over 2 getijperioden) S_opp, S_dpg, S_bod, T_opp, T_dpg, T_bod, tijd in minuten na 1-jun-2004');
%    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:nr2tg
        fprintf (fid4,frmt2,tijden_2tg(t),...
              sal_2tg(t,1),sal_2tg(t,2),sal_2tg(t,3),...
              tmp_2tg(t,1),tmp_2tg(t,2),tmp_2tg(t,3));
    end
    fclose(fid4);

    varstruct.Name = 'tijd_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'tijd_2tg',tijden_2tg);
    clear varstruct

    varstruct.Name = 'S_opp_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_opp_2tg',sal_2tg(:,1));
    clear varstruct

    varstruct.Name = 'S_dpg_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_dpg_2tg',sal_2tg(:,2));
    clear varstruct

    varstruct.Name = 'S_bod_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_bod_2tg',sal_2tg(:,3));
    clear varstruct

    varstruct.Name = 'T_opp_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_opp_2tg',tmp_2tg(:,1));
    clear varstruct

    varstruct.Name = 'T_dpg_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_dpg_2tg',tmp_2tg(:,2));
    clear varstruct

    varstruct.Name = 'T_bod_2tg';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2tg'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_bod_2tg',tmp_2tg(:,3));
    clear varstruct

    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, ystat);
    fprintf (fid5,['%s \n'], 'Tijdgemiddelden (over 2 weken) S_opp, S_dpg, S_bod, T_opp, T_dpg, T_bod, tijd in minuten na 1-jun-2004');
%    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:nr2wk
        fprintf (fid5,frmt2,tijden_2wk(t),...
              sal_2wk(t,1),sal_2wk(t,2),sal_2wk(t,3),...
              tmp_2wk(t,1),tmp_2wk(t,2),tmp_2wk(t,3));
    end
    fclose(fid5);

    varstruct.Name = 'tijd_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'tijd_2wk',tijden_2wk);
    clear varstruct

    varstruct.Name = 'S_opp_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_opp_2wk',sal_2wk(:,1));
    clear varstruct

    varstruct.Name = 'S_dpg_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_dpg_2wk',sal_2wk(:,2));
    clear varstruct

    varstruct.Name = 'S_bod_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_bod_2wk',sal_2wk(:,3));
    clear varstruct

    varstruct.Name = 'T_opp_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_opp_2wk',tmp_2wk(:,1));
    clear varstruct

    varstruct.Name = 'T_dpg_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_dpg_2wk',tmp_2wk(:,2));
    clear varstruct

    varstruct.Name = 'T_bod_2wk';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_2wk'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_bod_2wk',tmp_2wk(:,3));
    clear varstruct

    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, ystat);
    fprintf (fid6,['%s \n'], 'Tijdgemiddelden (over half jaar) S_opp, S_dpg, S_bod, T_opp, T_dpg, T_bod, tijd in minuten na 1-jun-2004');
%    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

    for t = 1:nrhjr
        fprintf (fid6,frmt2,tijden_hjr(t),...
          sal_hjr(t,1),sal_hjr(t,2),sal_hjr(t,3),...
          tmp_hjr(t,1),tmp_hjr(t,2),tmp_hjr(t,3));
    end
    fclose(fid6);

    varstruct.Name = 'tijd_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'tijd_hjr',tijden_hjr);
    clear varstruct

    varstruct.Name = 'S_opp_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_opp_hjr',sal_hjr(:,1));
    clear varstruct

    varstruct.Name = 'S_dpg_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_dpg_hjr',sal_hjr(:,2));
    clear varstruct

    varstruct.Name = 'S_bod_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'S_bod_hjr',sal_hjr(:,3));
    clear varstruct

    varstruct.Name = 'T_opp_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_opp_hjr',tmp_hjr(:,1));
    clear varstruct

    varstruct.Name = 'T_dpg_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_dpg_hjr',tmp_hjr(:,2));
    clear varstruct

    varstruct.Name = 'T_bod_hjr';
    varstruct.Nctype = 'double';
    varstruct.Dimension = {'series_hjr'};
    nc_addvar(ncfile,varstruct);
    nc_varput(ncfile,'T_bod_hjr',tmp_hjr(:,3));
    clear varstruct
end