function  writeLPF(basename,lpf)
% writeLPF(basename,lpf) --- write linear flow package file
% TO 0706030


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

NLAY=length(lpf.LAYCON);

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

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

%1.
fprintf(fid,'%10d%10.3g%10d     ILPFCB HDRY NPLPF\n',lpf.ILPFCB,lpf.HDRY,lpf.NPLPF);

%2 Layer convertibility (confined / unconfined)
warray(fid,lpf.LAYCON,lpf.unit,'(25I3)','LAYCON','noctrlrec');

%3 Interblock conductance computation method
warray(fid,lpf.LAYAVG,lpf.unit,'(25I3)','LAYAVG','noctrlrec');

%4 Horizontal layer anisotropy
warray(fid,lpf.CHANI,lpf.unit,'(10G10.3)','CHANI','noctrlrec');

%5 vertical anisotropy: preferably use 0 indicating we provide VK directly
warray(fid,lpf.LAYVKA,lpf.unit,'(10G10.3)','LAYVKA','noctrlrec');

%6 Layer wet-dry switch 
warray(fid,lpf.LAYWET,lpf.unit,'(25I3)','LAYWET','noctrlrec');

%7  Layer wetting options
if any(lpf.LAYWET)
    fprintf(fid,'%10.3g%10d%10d     WETFCT IWETIT IHDWET\n',lpf.WETFCT,lpf.IWETIT,lpf.IHDWET);
end

%8  we don't use parameters
%9  we don't use parameters

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

transient=any(lpf.isTran);

for iL=1:NLAY
    %10.
    warray(fid,lpf.HK(:,:,iL),lpf.unit,'(15G15.5)',sprintf(' HK{%d}',iL));  % hiorzontal conductivity
    %11 hor anisotropy through CHANI not through item 11
    if lpf.CHANI(iL)<=0
        warray(fid,lpf.HANI(:,:,iL),lpf.unit,'(15G15.5)',sprintf(' Hor anisotropy{%d}',iL));  % horizontal anistotropy
    end           
    %12 vertical conductivity
    warray(fid,lpf.VKA(:,:,iL),lpf.unit,'(15G15.5)',sprintf(' Vert Anisotropy{%d}',iL));  % vertical conductivity
    if transient
        %13 Ss
        if ~isfield(lpf,'SS')
           error('SS was not specified in mf_adapt for transient simulation');
        end
        warray(fid,lpf.SS(:,:,iL),lpf.unit,'(15G15.5)',sprintf(' SS{%d}',iL));  % Ss
        if lpf.LAYCON(iL)
            %14 Sy
            if ~isfield(lpf,'SY')
               error('SY was not specified in mf_adapt for transient simulation for converible layer %d',iL);
            end
            warray(fid,lpf.SY(:,:,iL),lpf.unit,'(15G15.5)',sprintf(' SY{%d}',iL));  %Sy
        end
    end
    
    %15 Resistance layer below model layers
    if lpf.LAYCBD(iL)
        warray(fid,lpf.VKCB(:,:,iL),lpf.unit,'(15G15.5)',sprintf(' VKCB{%d}',iL));  %VKCB
    end
    
    %16
    if lpf.LAYWET(iL)&& lpf.LAYCON(iL);
        if iscell(lpf.WETDRY)
            if length(lpf.WETDRY{iL}(:))==1
                warray(fid,lpf.WETDRY{iL}*PLANE,lpf.unit,'(15G12.5)',sprintf(' WETDRY{%d}',iL));  %WETDRY
            else
                warray(fid,lpf.WETDRY{iL}      ,lpf.unit,'(15G12.5)',sprintf(' WETDRY{%d}',iL));  %WETDRY
            end
        else
            warray(fid,    lpf.WETDRY(iL)*PLANE,lpf.unit,'(15G12.5)',sprintf(' WETDRY{%d}',iL));  %WETDRY
        end
    end
end    

fclose(fid);
