function  writeBCF(basename,bcf)
% writeBCF(basename,bcf) --- write basic flow package file (is BCF6)
% TO 0706030 081225 090713 091206

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

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

%0. No header allowed in BCF file
 fprintf(    '%s\n',['# MODFLOW writeBCF ' datestr(now)]);

%1.
fprintf(fid,'%10d%10g%10d%10g%10d%10d     %s\n',...
    bcf.IBCFCB,bcf.HDRY,bcf.IWDFLG,bcf.WETFCT,bcf.IWETIT,bcf.IHDWET,...
    '     IBCFCB HDRY IWDFLG WETFCT IWETIT IHDWET');
%2
bcf.LAYCON(bcf.LAYCON<0)=0; bcf.LAYCON(bcf.LAYCON>3)=0;
bcf.LAYAVG(bcf.LAYAVG<0)=0; bcf.LAYAVG(bcf.LAYAVG>3)=0;

 for iL=1:bcf.NLAY
     fprintf(fid,'%1d%1d',bcf.LAYAVG(iL),bcf.LAYCON(iL));
     if rem(iL,40)==0 && iL<bcf.NLAY
         fprintf(fid,'\n');
     end
 end
fprintf(fid,'     LTYPE(NLAY) (40I2) LAYAVG|LAYCON\n');
%3 Anisotropy factor
warray(fid,bcf.TPRY,bcf.unit,'(10G12.5)','TPRY (hor anisotropy factor)');

%% computing VCONT required by the bcf package

PLANE=ones(bcf.NROW,bcf.NCOL); % to change a layer par into a layer-wide cell par

k=0;
for iL=1:bcf.NLAY
    k=k+1;
    %4 SF1
    if any(bcf.isTran)
        warray(fid,bcf.SF1(:,:,iL),bcf.unit,'(10G12.5)',...
            sprintf('SF1(%d) LAYCON=%d',iL,bcf.LAYCON(iL)));
    end
    %5 TRAN where transmissiviy is fixed
    if any(bcf.LAYCON(iL)==[0 2])  % i.e. fixed transmissivity
        warray(fid,bcf.TRAN(:,:,iL),bcf.unit,'(10E12.3)',...
            sprintf('TRAN(%d) LAYCON=%d',iL,bcf.LAYCON(iL)));
    %6 HY where layer is convertable
    else % Horizontal conductiviy
        if ~isfield(bcf,'HY'),
            error('Need HY for kh in BCF layer %d!',iL);
        end
        warray(fid,bcf.HY(:,:,iL),bcf.unit,'(10E12.3)',...
            sprintf('HY(%d) LAYCON=%d',iL,bcf.LAYCON(iL)));
    end
    %7 VCONT for all but the last layer
    if iL<bcf.NLAY  % write the vertical conductance between the two layers
        warray(fid,bcf.VCONT(:,:,iL),bcf.unit,'(10E12.3)',...
            sprintf('VCONT(%d) LAYCON=%d',iL,bcf.LAYCON(iL)));
    end
    %8 SF2 secondary storage coefficient = SY if LAYCON is 2 or 3
    if any(bcf.isTran) && any(bcf.LAYCON(iL)==[2 3])
        if ~isfield(bcf,'SF2'),
            error('Need SF2 for SY in BCF layer %d',iL);
        end 
        warray(fid,bcf.SF2(:,:,iL),bcf.unit,'(10G12.5)',...
            sprintf('SF2(%d) LAYCON=%d',iL,bcf.LAYCON(iL)));
    end
    %9 WETDRY rewettabiliy only if LAYCON is 1 or 3
    if bcf.IWDFLG~=0 && any(bcf.LAYCON(iL)==[1 3])
        if iscell(bcf.WETDRY)
            if length(bcf.WETDRY{iL}(:))==1
                warray(fid,bcf.WETDRY{iL}*PLANE,bcf.unit,'(10G12.5)',sprintf(' WETDRY{%d} LAYCON=%d',iL,bcf.LAYCON(iL)));  %WETDRY
            else
                warray(fid,bcf.WETDRY{iL}      ,bcf.unit,'(10G12.5)',sprintf(' WETDRY{%d} LAYCON=%d',iL,bcf.LAYCON(iL)));  %WETDRY
            end
        else
            warray(fid,bcf.WETDRY(iL)*PLANE,bcf.unit,'(10G12.5)', sprintf(' Wetdry(%d) LAYCON=%d',iL,bcf.LAYCON(iL)));  %WETDRY
        end
    end

end
    
fclose(fid);
