function writeBCN(basename,bcn,bcnp)
% writeBCN(basename,bcn,bcnp) -- write modflow boundary condition file
% the bcn.name must be one of 'WEL', 'MNW', 'GHB', 'DRN', 'RIV', 'CHD'
% TO 070702 090807 100530 (write format change)


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

if isempty(bcn.parvals)  % check if any data are provided for this package
    fprintf(['No data are provided for package <<%s>> !\n',...
        'Provide data or switch off this package in the NAM worksheet.\n'],...
        bcn.type);
end

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

AsPrevious=-1; % flag indicating to reuse drains of previous stress period

%% Number of bcn per period
nDPer=zeros(bcn.NPER,1);    % the number of input lines per stress period or -1 (as previous) if not found
                            % to remove all boundaries in a given stress period use a single dummy one
for iper=1:bcn.NPER
    if any(bcn.parvals(:,1)==-iper)
        nDPer(iper)=AsPrevious;
        if iper==1,
            error('AsPrevious (negative stress period) cannot be used for first stress period !\n');
        end
    else
        nDPer(iper)=length(find(bcn.parvals(:,1)==iper)); % then count the number
    end
end
MXACTD=max(nDPer);  % max number of active drains at any one time

%% Generate the output

%% 0
fprintf(fid,'# MATLAB writeBCN for %s,  %s\n',bcn.type,datestr(now));
fprintf(    '# MATLAB writeBCN for %s,  %s\n',bcn.type,datestr(now));

%% 1
if bcnp.NPar~=0
    fprintf(fid,'PARAMETER%10d%10d\n',[bcnp.NP,bcnp.MXL]);
end

%% 2

% printing the simulation header allowing for auxiliary variables
% while taking the exception for CHD into account

% throw away empty AUX cells
for i=1:length(bcn.AUX)
    if isempty(bcn.AUX{i}),
        bcn.AUX(i:end)=[];
        break;
    end
end

% ===== first print the variables =================================
fprintf(fid,'%10d',MXACTD);  % always

if ~strcmp(bcn.type,'CHD'), fprintf(fid,'%10d',bcn.ICB   ); end % all but CHD

if  strcmp(bcn.type,'MNW'), fprintf(fid,'%10d',bcn.IWELPT); end % only MNW

for i=1:length(bcn.AUX)
    fprintf(fid,'  AUX %s',bcn.AUX{i});  % the variables
end

% ===== next print information on these variables =================
fprintf(fid,'    <===== MXACTD');
if ~strcmp(bcn.type,'CHD')
    fprintf(fid,'     I%sICB',bcn.type);
end
if  strcmp(bcn.type,'MNW') % only MNW
    fprintf(fid,'     IWELPT');
end

if ~isempty(bcn.AUX)
    fprintf(fid,' +  Auxiliary variables');
end
fprintf(fid,'\n');

if  strcmp(bcn.type,'MNW') % only MNW
    if bcn.LOSSTYPE==1
        fprintf(fid,'%10s\n','SKIN');
    end
    if bcn.LOSSTYPE==2
        fprintf(fid,'%10s\n','LINEAR');
    end
    if bcn.LOSSTYPE==3
        fprintf(fid,'%10s','NONLINEAR');
        fprintf(fid,'%10s\n',bcn.PLOSSMNW)
    end
end
    
%% 3

 fmt='%10.3g'; % works for numbers >-10e100, otherwise the size occupied will be 11 characters instead of 10 (=' %9g');
 
 for iPar=1:bcnp.NPar
    bcnp.NLST(iPar)=length(find(bcnp.parvals(:,1)==iPar));
    fprintf(fid,'%11s%11s% 10d%10d\n',char(bcnp.PARNAM(iPar)),char(bcnp.PARTYP(iPar)),bcnp.Parval(iPar),bcnp.NLST(iPar));
    L=size(bcnp.parnams,2)-4;
    fprintf(fid,['%10d%10d%10d',repmat(fmt,[1,L]),'\n'],bcnp.parvals(bcnp.parvals(:,1)==iPar,2:end));
 end

%% The non parameter (ordinary) values for this boundary condition
L=size(bcn.parvals,2)-4;

switch bcn.type
    case 'WEL', n=1; info='Q';
    case 'MNW', n=9; info='Qdes MN QWval RW Skin Hlim Href DD LWGRP';
    case 'DRN', n=2; info='Elev Cond';
    case 'RIV', n=3; info='Stage Conc RivBotElev';
    case 'GHB', n=2; info='Bhead Cond';
    case 'CHD', n=2; info='Shead Ehead <<denseopt in Seawat>>';
    otherwise,
        error('Illegal/unknown bcn.type=<<%s>>\n',bcn.type);
end
if L<n+length(bcn.AUX)
     error(['Writing %s file:\n',...
         'You use %d auxiliary parameters\n',...
         'The number of items per line must be >= %d'\n,...
         'See if you forgot to prepend the stress period number.'],...
         bcn.type,length(bcn.AUX),4+n+length(bcn.AUX));
end


for iPer=1:bcn.NPER   %length(bcn.ITMP)
   I2=find(bcn.parvals(:,1)==iPer);  % which drains on in this period?
   if iPer>1 &&  ...
           (nDPer(iPer)==AsPrevious || ...
           (nDPer(iPer)==nDPer(iPer-1) && all(all(bcn.parvals(I2,:)==bcn.parvals(I1,:)))))
                fprintf(fid,'%10d%10d\n',AsPrevious,bcnp.NPar);
   else
       %5
       fprintf(fid,'%10d%10d     ITMP(%d)  NP(%d). Following lines: L R C %s\n',...
           nDPer(iPer),bcnp.NPar,iPer,iPer,info); % ITMP NP  should be NP(iPer)
       %6
       if nDPer(iPer)~=AsPrevious,
           if ~strcmp(bcn.type,'MNW')
               fprintf(fid,['%10d%10d%10d',repmat(fmt,[1,L]),'\n'],bcn.parvals(I2,2:end)');
           else % only MNW, because of strings ('MN' and 'DD'
               for irec=I2'
                   fprintf(fid,'%10d%10d%10d',bcn.parvals(irec,2:4)');
                   fprintf(fid,fmt,bcn.parvals(irec,5)');
                   if bcn.parvals(irec,6)==1
                       fprintf(fid,'%10s','MN');
                   else
                       fprintf(fid,'%10s','');
                   end
                   fprintf(fid,repmat(fmt,[1,5]),bcn.parvals(irec,7:11)');
                   if bcn.parvals(irec,12)==1
                       fprintf(fid,'%10s','DD');
                   else
                       fprintf(fid,'%10s','');
                   end
                   fprintf(fid,[fmt,'\n'],bcn.parvals(irec,13)');
               end
           end
       else
           fprintf(fid,'%10d\n',AsPrevious);
       end
   end
   I1=I2;
end

fclose(fid);
