function writeDIS(basename,dis)
% writeDIS(basename,dis)  --- writing discretization file
% TO 070630


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later

if ~exist('ITMUNI','var'), ITMUNI = 4; end   % time units (4=days)     Used as default here
if ~exist('LENUNI','var'), LENUNI = 2; end  % length units (2=meters) Used as default here

fid=fopen([basename,'.',dis.ext],'wt');

%0.
fprintf(fid,'# MATLAB  writeDIS %s\n',datestr(now));
fprintf(    '# MODFLOW writeDIS %s\n',datestr(now));

%1.
fprintf(fid,'%10i%10i%10i%10i%10i%10i      %s\n',...
    dis.NLAY,dis.NROW,dis.NCOL,dis.NPER,dis.ITMUNI,dis.LENUNI,...
    'NLAY, NROW NCOL NPER ITMUNI LENUNI');

%2  Resistance at bottom of layer?? -- Quasi 3D flag for each layer
warray(fid,dis.LAYCBD,dis.unit,'(40I2)','LAYCBD','norec');

%3.
warray(fid,dis.DELC(:)',dis.unit,'(10E12.3)','DELX or DELC');

%4.
warray(fid,dis.DELR(:)',dis.unit,'(10E12.3)','DELY or DELR');

%5
warray(fid,dis.Z(:,:,1),dis.unit,'(10E15.6)','TOP of model');

%6.
k=1;
for i=1:dis.NLAY
    k=k+1;
    warray(fid,dis.Z(:,:,k),dis.unit,'(10E15.6)',sprintf('BOT{%d}',k-1));
    if dis.LAYCBD(i) && i<dis.NLAY
        k=k+1;
        warray(fid,dis.Z(:,:,k),dis.unit,'(10E15.6)',sprintf('BOT{%d}',k-1));
    end
end

%% FOR EACH STRESS PERIOD
%7. 
for iP=1:dis.NPER
    fprintf(fid,'%10s%10d%10g',numsqueeze(dis.PERLEN(iP),'%g',10),dis.NSTP(iP),dis.TSMULT(iP));
    if dis.isTran(iP),
        fprintf(fid,'     TR');
    else
        fprintf(fid,'     SS');
    end;
    fprintf(fid,'     PERLEN NSTP TSMULT SS/TR\n');
end

fclose(fid);
