function  writeSSM(basename,ssm)
% writeSSM(basename,ssm) --- write sink-source mixing package file
% TO 0706030 081227


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

ssm.unit=100;  % incompatibility with UD2REL (MT3DMS p97)

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

%D0 HEADING 1+2 (<=80 chars) -- No header allowd in DSP file
fprintf(    '%s\n',['# MT3DMS writeSSM ' datestr(now)]);


%D1 FWEL FDRN FRCH FEVT FRIV FGHB FCHD FMNW (FNEW(n), n=1:4) 10I2
% option flags for WEL, DRN etc packages + FNEW (n=1:4) future options
if ssm.FWEL, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FDRN, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FRCH, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FEVT, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FRIV, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FGHB, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FCHD, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FMNW, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FNW2, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FNW3, fprintf(fid,' T'); else fprintf(fid,' F'); end
if ssm.FNW4, fprintf(fid,' T'); else fprintf(fid,' F'); end
fprintf(fid,'    FWEL FDRN FRCH FEVT FRIV FGHB FCHD FMNW (FNEW(n n=1:3)\n');

fprintf(fid,'%10d     MXSS max number of sink/source points in model\n',ssm.MXSS);
   
%% ======================== FOR EACH PERIOD ========================

ASPREVIOUS=-1;

for iPer=1:ssm.NPER
    
    if ssm.FRCH
        %D3 INCRCH (I10) read recharge flux for this stress period?
        fprintf(fid,'%10d     INCRCH(%d)\n',ssm.INCRCH(iPer),iPer);
        if ssm.INCRCH(iPer)>=0
            for iComp=1:ssm.NCOMP
                %D4 CRCH(NCOL,NROW) if FRCH=T and INCRCH>=0
                warray(fid,ssm.CRCH{iPer,iComp},ssm.unit,'(10E12.3)',sprintf('CRCH{%d}',iPer));
            end
        end
    end

    if ssm.FEVT
        %D5 INCEVT (I10) read evt flux for this stress period?
        fprintf(fid,'%10d     INCEVT(%d)\n',ssm.INCEVT(iPer),iPer);
        if ssm.INCEVT(iPer)>=0
            for iComp=1:ssm.NCOMP
                %D6 CEVT(NCOL,NROW) if FEVT=T and INCEVTRCH>=0
                warray(fid,ssm.CEVT{IP,iComp},ssm.unit,'(10E12.3)',sprintf('CEVT{%d}',iPer));
            end
        end
    end

    %D7 NNS (I10) max of point souces for which conc needs to be specified or
    %read for the current stress period. (MXSS is overall maximum of all sources and sinks)
    % By default, unspecified point sources are assumed zero concentration.
    % (The concentration of point sinks is always set equal
    %to the concentration of groundwater at the sink location). MT3DMS p121
    %
    % Not all sources and sinks need to be specified, only those for which
    % a concentration needs ot be specified. For constant head cells,
    % specifying the concentration means, specifying the concentration of
    % inflowing water, not that of the entire cell, which is the result of
    % mixing. This is the same as for wells and different for constand
    % concentration cells, for which the concentration of the entire cell
    % is fixed, irrespective of its size and flows in and out of it.
    %
    % So how do we specify elegantly only those point sinks for which
    % specifying the concentration of the incoming water is relevant?
    % We have to do this for each stress period.
    %
    % The specification here takes precedence over what has been specified
    % in the BTN file through ICBUND. So we may overrid ICBUND starting
    % concentrations here, if we so wish.
    
    %find(mnw.parvals(:,1)==iPer & mnw.parvals(:,6)==0) WEGHALEN
    
    if isempty(ssm.PNTSRC),   I=[];
    else                      I=find(abs(ssm.PNTSRC(:,1))==iPer);     end
    
    if ~isempty(I) && ssm.PNTSRC(I(1),1)<0
        fprintf(fid,'%10d\n',ASPREVIOUS);
    else
        NSS=length(I);
        fprintf(fid,['%10d     NSS(period=%d).', ...
            '  Following lines: L R C CSS ITYPE CSSMS(1..%d)\n'], NSS,iPer,ssm.NCOMP);

        fmt='%10d%10d%10d%10.3g%10d';
        for iComp=1:ssm.NCOMP,fmt=[fmt,' %12g']; end  % is free format acc. to manual

        % Even if NCOMP=1, use L R C CSS ITYP CSSMS(1), see MT3DMS manural p122    
        if ssm.NCOMP==1 && size(ssm.PNTSRC,2)==6
            ssm.PNTSRC=ssm.PNTSRC(:,[1:6, 5]);     % add CSS as 7th column, i.e. CSSMS(1)
        end

        % Make sure to specify the correct number of Species in PNTSRC
        if NSS>0
            if size(ssm.PNTSRC,2)~=6+ssm.NCOMP
                errstr=['The number of species in PNTSRC=%d, which does not match the number\n',...
                        'specified in the field NCOMP=%d the worksheet MT3D of your workbook!\n',...
                        'See MT3DMS manual version 5.2, page 103 and item D8 on page 122.'];
                    error(errstr,size(ssm.PNTSRC,2)-6,ssm.NCOMP);
            end

            fprintf(fid,[fmt,'\n'],ssm.PNTSRC(I,2:(6+ssm.NCOMP))');
        end
    end
end

fclose(fid);
