%% ===== mf_setup  --- mfLab SETUP SCRIPT ===================
% This generic script generates inpt files for MODFLOW, MT3D and SEAWAT
% and an increasing number of other packages for the MODFLOW suite.
% First the paths are set to the directories containing the mfiles needed
% by mflab.
% In a second step mf_adapt in run locally, which is a script that
% generates the actual model.
% The matrices generated in mf_adapt are stored in basename.mat for
% subsequent retrieval during the generation of the input files for the
% grounwater models and as a backup.
% The basename of the curent model is stored in basename.mat. Default is
% the name of the local directory.
% By default all geneated filename consist of basename plus an extension
% denoting the function of the file.
% mf_setup will use basename.xls, an excel file as a parameter container.
% It is used to read the packages )to be invoked and the parameters that
% need to be set in the simulations. These
% parameters are necessary for the input files.
% See the documentaion for more detailed information.
% TO 090101 091010 091204

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

clear variables

%% Set paths, as mf_adapt may need them

fprintf('======= Start construction of your model ============\n');

fprintf('Your current HOME directory is\n%s\n',pwd);

% make a shortcut in the shortcut menu to set the path to the
% directories where the mfiles of mflab reside
% by clicking that shortcut after loading modflow
% Matlab is aware of these folders and the functions in them.
%     
% if ismac
%     MODELS='/Users/Theo/GRWMODELS/';
% else
%     MODELS='Z:\tolsthoorn On My Mac\GRWMODELS\';
% end
%
% MFLAB=[MODELS,'mflab' filesep]; % mflab home directory
% addpath([MFLAB 'write']);
% addpath([MFLAB 'read' ]);
% addpath([MFLAB 'gridcoords']);
% addpath([MFLAB 'etc']);
% addpath([MFLAB 'fdm']);

%% =====RUN MMF_ADAPT ===================================================
%  MF_ADAPT is a local m-file that may be used to the model matrices
%  created or to make a model completely
%  also in mf_adapt you set the paths for your local model

fprintf('Running model construction script mf_adapt..\n');

mf_adapt    % local script to assemble or adapt your conceptual model

fprintf('... mf_adapt finished, model has been constructed.\n');

%% Save the model so we can generate input files

fprintf('Saving data to <<%s.mat>>\n',basename);

save('name.mat','basename');

S='save([basename,''.mat''],''IBOUND'',''STRTHD''';  % always

% Grid, expect you to use xGr, yGr
if exist('xGr'   ,'var'), S=[S,',''xGr'''  ]; end  % always use xGr, vector!
if exist('yGr'   ,'var'), S=[S,',''yGr'''  ]; end  % always use yGr, vector!
if exist('zGr'   ,'var'), S=[S,',''zGr'''  ]; end  % always use zGr, vector!
if exist('Z'     ,'var'), S=[S,',''Z'''    ]; end  % full 3D Z matrix!!

% Grid, preferably use DELX, DELY or DX, DY to prevent confusion
if exist('DELC'  ,'var'), DELC=DELC(1,:,1); S=[S,',''DELC'''  ]; end  % DX not DY!!
if exist('DELR'  ,'var'), DELR=DELR(:,1,1); S=[S,',''DELR'''  ]; end  % DY not DX!!
if exist('DELX'  ,'var'), DELX=DELX(1,:,1); S=[S,',''DELX'''  ]; end  % clear alternative
if exist('DELY'  ,'var'), DELY=DELY(:,1,1); S=[S,',''DELY'''  ]; end  % clear alternative
if exist('DX'    ,'var'), DX  =DX(  1,:,1); S=[S,',''DX'''    ]; end  % clear alternative
if exist('DY'    ,'var'), DY  =DY(  :,1,1); S=[S,',''DY'''    ]; end  % clear alternative
if exist('Dx'    ,'var'), Dx  =Dx(  1,:,1); S=[S,',''Dx'''    ]; end  % clear alternative
if exist('Dy'    ,'var'), Dy  =Dy(  :,1,1); S=[S,',''Dy'''    ]; end  % clear alternative
if exist('dx'    ,'var'), dx  =dx(  1,:,1); S=[S,',''dx'''    ]; end  % clear alternative
if exist('dy'    ,'var'), dy  =dy(  :,1,1); S=[S,',''dy'''    ]; end  % clear alternative

if exist('LAYCBD','var'), S=[S,',''LAYCBD''']; end  % confinging layer config

% Check DX and DY as the might be 2D after using meshgrid inaccidentally
if exist('DX','var'), DX=DX(1,:); end
if exist('DY','var'), DY=DY(:,1); end

% Hydraulic parameters for LPF package
if exist('HK'    ,'var'), S=[S,',''HK'''    ]; end  % LPF (alias for HY in BCF)
if exist('HANI'  ,'var'), S=[S,',''HANI'''  ]; end  % LPF ky/kx
if exist('VK'    ,'var'), S=[S,',''VK'''    ]; end  % LPF
if exist('VKA'   ,'var'), S=[S,',''VKA'''   ]; end  % LPF kv/kh if LAYVKA ~=0
if exist('SS'    ,'var'), S=[S,',''SS'''    ]; end  % for LPF
if exist('SY'    ,'var'), S=[S,',''SY'''    ]; end  % for LPF
if exist('VKCB'  ,'var'), S=[S,',''VKCB'''  ]; end  % for LPF
if exist('WETDRY','var'), S=[S,',''WETDRY''']; end  % for LPF

% Hydraulic parameters for BCF package
if exist('HY'    ,'var'), S=[S,',''HY'''    ]; end  % for BCF=KH
if exist('TRAN'  ,'var'), S=[S,',''TRAN'''  ]; end  % for BCF
if exist('VCONT' ,'var'), S=[S,',''VCONT''' ]; end  % for BCF
if exist('SF1'   ,'var'), S=[S,',''SF1'''   ]; end  % for BCF
if exist('SF2'   ,'var'), S=[S,',''SF2'''   ]; end  % for BCF

% Recharge and evaptranspiration packages
if exist('RECH'   ,'var'), S=[S,',''RECH''' ]; end
if exist('IRCH'   ,'var'), S=[S,',''IRCH''' ]; end
if exist('SURF'   ,'var'), S=[S,',''SURF''' ]; end
if exist('EVTR'   ,'var'), S=[S,',''EVTR''' ]; end
if exist('EXDP'   ,'var'), S=[S,',''EXDP''' ]; end
if exist('IEVT'   ,'var'), S=[S,',''IEVT''' ]; end

% Boundary conditions for MODFLOW
if exist('WEL'   ,'var'), S=[S,',''WEL'''   ]; end
if exist('MNW'   ,'var'), S=[S,',''MNW'''   ]; end
if exist('DRN'   ,'var'), S=[S,',''DRN'''   ]; end
if exist('RIV'   ,'var'), S=[S,',''RIV'''   ]; end
if exist('GHB'   ,'var'), S=[S,',''GHB'''   ]; end
if exist('CHD'   ,'var'), S=[S,',''CHD'''   ]; end

% Boundary conditions transport MT3DMS and Seawat
if exist('ICBUND','var'), S=[S,',''ICBUND''']; end  % MT3D
if exist('STCONC','var'), S=[S,',''STCONC''']; end  % MT3D BTN
if exist('PNTSRC','var'), S=[S,',''PNTSRC''']; end  % MT3D SSM
if exist('BTNOBS','var'), S=[S,',''BTNOBS''']; end  % MT3D BTN

% Recharge and EVT concentrations in SSM
if exist('CRCH','var'), S=[S,',''CRCH''']; end  % MT3D SSM
if exist('CEVT','var'), S=[S,',''CEVT''']; end  % MT3D SSM

% Necessary for transport, effective porosity
if exist('POROSITY','var') && ~exist('PEFF','var'), PEFF=POROSITY; end
if exist('PEFF'  ,'var'), S=[S,',''PEFF'''  ]; end  % MT3D BTN
if exist('DMCOEF','var'), S=[S,',''DMCOEF''']; end  % MT3D DSP

% SEAWAT (direct reading in of DENSE and VISC probably never going to be used
if exist('DENSE' ,'var'), S=[S,',''DENSE''' ]; end  % SEAWAT
if exist('VISC'  ,'var'), S=[S,',''VISC'''  ]; end  % SEAWAT

% Additional options for CHD in relation to Seawat, see Langevin e.a. 2008,
% p12-15 & 22
if exist('CHDDENSOPT','var'), S=[S,',''CHDDENSOPT''']; end % CHD-SEAWAT
if exist('CHDDENS','var')   , S=[S,',''CHDDENS'''   ]; end % CHD-SEAWAT

% Necessary for SWI, the Saltwater Instrusion Package
% p12-15 & 22
if exist('ISOURCE','var'), S=[S,',''ISOURCE''']; end % SWI
if exist('ZETA'   ,'var'), S=[S,',''ZETA'''   ]; end % SWI
if exist('SSZ'    ,'var'), S=[S,',''SSZ'''    ]; end % SWI

% To allow execuring specific mfiles after mf_setup
if exist('AFTERMFSETUP','var'), S=[S,',''AFTERMFSETUP''']; end % PEST

S=[S,');'];

eval(S); fprintf('%s\n',S);

fprintf('.... Saved !\n');

%%  =====START GENERATING MODFLOW etc INPUT FILES ============================
% Cleaning up the workspace removes cluttering variables and is important
% to make sure later on that you are working with the correct variables
% saved in <<basename.mat>> and that you saved all variables necessary
% to run the model

fprintf('Cleanup workspace ....');

clear variables                % clean workspace, but leave debug settings
%close all                      % close all figures
load('name'); load(basename);  % load data computed in mf_adapt and just saved

fprintf('... workspace cleaned and <<%s.mat>> reloaded.\n',basename);

%% ======= IF NECESSRY ADAPT TO YOUR OWN SITUATION ====================

tic;   % tick is to measure elapsed time from this point toc to present it

% Path to the executables. I put them all in MODELS/bin
% because I have several differntly compiled versions. It's your
% choise, anyway set the parameters MODFLOW MT3D etc down below to
% their actual lcoations on your hard drive.

% NOTICE the " " in the paths below to manage spaces in file names
% in system command used to launch the external executable later on

fprintf('Defining paths to your excecutables\n');

setExecutables;
        
%% Start generation of input files

fprintf('\nStarting generation of input files ....\n');

%% NAM Namefile INFO
fprintf('Preparing name file struct.\n');

nam.MF=MODFLOW; % path to modflow model
nam.MT=MT3D;    % path to mt3d model
nam.SWT=SEAWAT;  % path to seawat model
nam.SWI=SWI;    % path to swi model

%% Legal packages, the list may be extented
% Specific for MODFLOW
nam.LegalMF={...
    'GLOBAL','LIST',...
    'BAS6','DIS',...
    'BCF6','LPF',...
    'RCH','EVT',...
    'WEL','MNW1','RIV','GHB','DRN','CHD',...
    'PES','SEN','OBS','HOB','LMT6','UMT',...
    'CFP','COC','CRCH',...
    'PCG','DE4','SIP','SOR', 'GMG',...
    'OC','DATA ','DATA(BINARY)'...
    };

% Specific legal packages for MT3DMS (explicitly defined as not all MF
% packages are  legal MT3D packages
transportPkgs={'BTN','ADV','DSP','SSM','GCG','RCT'};

nam.LegalMT ={'LIST','FTL',transportPkgs{1:end}};

% Specific for Seawat adds SWT specific packages to legal MF packages
nam.LegalSWT={nam.LegalMF{:},'VDF','VSC',transportPkgs{1:end}};

% Specific for Salt Water Intrusion (SWI), add 'SWI' package to legal MF packages
nam.LegalSWI={nam.LegalMF{:},'SWI'};

%% STARTING WITH THE ACTUAL DATA FOR THIS MODEL =============

% ===== The EXCEL workbook with all parameter settings ===============
fprintf('Basename current model is ''%s''\n',basename);

if     exist([basename,'.xls' ],'file'), XLSF=[basename,'.xls'];
elseif exist([basename,'.xlsx'],'file'), XLSF=[basename,'.xlsx'];
else
    error('Can''t find accompanying Excel workbook <<%s>> nor <<%s>>',...
    [basename,'.xls'],[basename,'.xlsx']);
end
fprintf('Getting nam file data from %s\n',XLSF);

%% xlsread mac is different from xlsread windows !!!

if ismac
    warning('off','MATLAB:xlsread:Mode')
    try
        [Nnum,Ntxt]=xlsread(XLSF,'NAM','','basic');
    catch ME
        if ~isempty(findstr(ME.identifier,'WorksheetNotFound'))
                error(['%s ''NAM''\n',...
                    'Save the worksheet as an Excel 5.0/95 Workbook (.xls)',...
                    ' so that Matlab can read it on non-Windows systems\n'],...
                    ME.message);
        end
    end
    warning('on','MATLAB:xlsread:Mode')
else
    [Nnum,Ntxt]=xlsread(XLSF,'NAM');    
end
% Get the data for the *.nam file from the NAM worksheet in the XSLF workbook 

%% Another incompatibility between xlsread mac and windows
if ismac   %  size(Nnum,2)==size(Ntxt,2)
    Nnum=Nnum(2:end,:);  % cut off header line
    Ntxt=Ntxt(2:end,:);  % cut off header line
else
    % on pc Nnum needs not be cut contrary to Mac
    Ntxt=Ntxt(2:end,:);  % cut off header line
end

%%
nam.PCKG=Ntxt(:,1);     % List of packages
nam.EXT =Ntxt(:,3);     % corresponding file extenstions
nam.UNIT=Nnum(:,1);     % corresponding file units
nam.SCEN=Nnum(:,3);
		% =0 package not inlcuded
		% >0 value is scenario number. This corresponds with data column
	    %     in the worksheets MFLOW, MT3D and SWT next to the package/parnam column
	    %     this column must be present when package is processed!!
        % <0  for VDF run seawat but do not run package
        % <0  for SWI run swi    but do not run package

% ====UNIT numbers for the output files (hds, ddn, bgt,zta) ==============
% These units are associate with the given extension "hds" "ddn" "bgt" "zta"
% and not with the specified package names

% Internally the model is given the unit to write values onto and the
% filename associated with the unit. There is no package connected to it.
% Therefore we must be careful to match the unit of values to write with
% the one in the name file. In MatlabModflow, i.e. this script, we use the
% given extension given in the nam sheet as an indication for the type of output
% to write. Therefore the extensions used in the name file for MODFLOW
% output are obligitory. Therefore the hds are always written to the file
% <basename>.hds, the interfaces produced by SWI are always output to the
% file <basename>.zta, the budgets are always written to <basename.bgt>.
% You may use other extensions, but then you must add extra lines below to
% link those exensions with the proper MODFLOW (MT3DMS etc) output. See for
% instance the dummy line with extension xxx to write budget to xxx if this
% xxx is included in the nam sheet.

%% REMOVE all PACKAGES THAT ARE OFF

% scenario switch (3rd col in NAM sheet) if <0 then run program without
% package switched on
if nam.SCEN(strmatchi('VDF',nam.PCKG,'noerr'))<0, nam.runSWT=1; else nam.runSWT=0; end

% remove all package that are off
nam.PCKG=nam.PCKG(abs(nam.SCEN)>0);
nam.EXT =nam.EXT (abs(nam.SCEN)>0);
nam.UNIT=nam.UNIT(abs(nam.SCEN)>0);
nam.SCEN=nam.SCEN(abs(nam.SCEN)>0);

%% Possibly use another MODFLOW program instead of default MF2K
%  This is done on the NAM worksheet using pseudo packages MF5 and ASP
i=strmatchi('MF5',nam.PCKG,'noerr');
if i>0 && nam.SCEN(i)>0,
    fprintf('NAM worksheet specifies use of MF2005 instead of MF2k\n');
    MODFLOW=MF2005;
    nam.PCKG(i)=[]; nam.EXT(i)=[]; nam.UNIT(i)=[]; nam.SCEN(i)=[];
end

i=strmatchi('ASP',nam.PCKG,'noerr');
if i>0 && nam.SCEN(i)>0
    if ismac
        error('MF2KASP specified in nam sheet, but John Doherty''s MF2KASP is only available under windows !');
    end
    fprintf('NAM worksheet specifies use of MF2KASP instead of MF2K\n');
    MODFLOW=MF2KASP;
    nam.PCKG(i)=[]; nam.EXT(i)=[]; nam.UNIT(i)=[]; nam.SCEN(i)=[];    
end


%% Check if unit numbers are legal

UNITS=NaN(length(nam.UNIT),2);
UNITS(:               ,1)=sort(nam.UNIT);  uniqueunits=unique(nam.UNIT);
UNITS(1:length(uniqueunits),2)=uniqueunits;
I=find(UNITS(:,1)~=UNITS(:,2));
if ~isempty(I)
    error('Unit number %d in NAM file is not unique !',UNITS(I(1),1));
end

%% Check uniqueness of file extensions

for i=1:length(nam.EXT)-1
    for j=i+1:length(nam.EXT)
        if strmatchi(nam.EXT{i},nam.EXT(j:end),'noerr'),
            error('File extension <<%s>> in nam file is not unique !',nam.EXT{i});
        end
    end
end

writeNAM(basename,nam)  % GENERATE THE NAM+BAT FILES FOR MODFLOW, MT3D AND SEAWAT

%% SIMULATION INFO MODFLOW (always needed)
fprintf('Getting simulation parameters from worksheet %-8s in %s\n','MFLOW',XLSF);
[MFLOWparnams,MFLOWparvals]=getExcelData(XLSF,'MFLOW','Vertical');

%% SIMULATION INFO MT3D (if  MT3D is to be run)
if strmatch('BTN',nam.PCKG)
    fprintf('Getting simulation parameters from worksheet %-8s in %s\n','MT3D',XLSF);
    [MT3Dparnams,MT3Dparvals]=getExcelData(XLSF,'MT3D','Vertical');
end

%% SIMULATION INFO SEAWAT (if  SEAWAT is to be run)
if strmatch('VDF',nam.PCKG)
    fprintf('Getting simulation parameters from worksheet %-8s in %s\n','SEAWAT',XLSF);
    [SEAWATparnams,SEAWATparvals]=getExcelData(XLSF,'SEAWAT','Vertical');
end

%% STRESS PERIOD INFO (always needed)
fprintf('Getting stress period data from    worksheet %-8s in %s\n','PER',XLSF);
[PERparnams,PERparvals]=getPeriods(XLSF);

%% LAYER INFO (always needed)
fprintf('Getting layer info from            worksheet %-8s in %s\n','LAY',XLSF);
if ~exist('Z','var')
    error(['You must define Z(NROW,NCOL,NLAY+sum(LAYBCD)+1) in mf_adapt!\n',...
     'Sum(LAYBCD) is the number of confined layers as defined in worksheet LAY,\n',...
     'which is needed if BCF or BCF6 is used']);
end
[LAYparnams,LAYparvals]=getLayers(XLSF,Z);

%% Check if the necessary DATA(BINARY) and DATA files are defined for hds
% and ddn (heads and drawdown)

ioc=strmatchi('OC',nam.PCKG,'noerr');

I___UN={'HED' ,'DDN'};   % see I*UN in MFLOW worksheet
SV    ={'HDSV','DDSV'};  % see PER worksheet
for i=1:length(I___UN)
    % do we want output for HEAD or DDN?
    if nam.UNIT(ioc)>0 && any(PERparvals(:,strmatchi(SV{i},PERparnams)))
        
        % Then look for the corresponding output file in the NAM fle/worksheet
        % Is any DATA or DATA(BINARY) output file defined in NAM ???
        Idata=strmatchi('DATA',nam.PCKG,'noerr');
        if ~Idata   % if none, throw error message and stop:
            error(['To write output as desired by the variable <<%s>> in the PER worksheet\n',...
                'you must specfy an output file for it in the NAM worksheet\n',...
                '(i.e. one line starting with DATA or DATA(BINARY))\n',...
                'which must be in accordance with the unit number for I%sUN in the\n',...
                'MFLOW worksheet. You may omit all output by switching off\n',...
                'the output package (OC) in the NAM file/worksheet.'],...
                SV{i},I___UN{i});
        end

        % What is the desired output unit specified in MOFLOW worksheet?
        iunit=MFLOWparvals(strmatchi(['I',I___UN{i},'UN'],MFLOWparnams,1));
        
        % is this the unit in one of the active data files in the NAM file/worksheet?
        if ~any(iunit==nam.UNIT(Idata))
            error(['Check the NAM sheet: Missing DATA output unit number <<%d>>\n',...
                'to ouput for <<I%sUN>> specified in worksheet MFLOW !'],iunit,I___UN{i});
        end
    end
end

%% Check if the necessary DATA(BINARY) and DATA files are defined

% These are the modules that may output cell-by-cell flow terms
% The concerned parameter names are obtained by preceding them with
% capital 'I' and appending 'CB'.
ICB={'LPF','BCF','RIV','DRN','GHB','CHD','WEL','MNW','EVT','RCH'};

% We will throw this error string in case an error is found

ioc=strmatchi('OC',nam.PCKG,'noerr');

% if any output of cell-by-cell flow terms is desired (see ICBCFL in PER worksheet) 
if nam.SCEN(ioc) && any(PERparvals(:,strmatchi('ICBCFL',PERparnams)))
    % Then we must find the desired unit number in one of the NAM file lines
    % specifying output with DATA or DATA(BINARY) in the first column.

    Idata=strmatchi('DATA',nam.PCKG,'noerr');  % find these DATA, DATA(BINARY) lines
    if ~Idata
        error(['To write cell-by-cell flow terms, you must specfy output files',...
            '(these are lines starting with DATA or DATA(BINARY))\n',...
            'in the NAM file/NAM worksheet in accordance with the unit numbers',...
            'for the packages I*CB in the MFLOW worksheet\n',...
            'If you do not want any cell by cell output, set ICBCFL=0 in PER worksheet.']);
    end
    
    % check for the individual packages that may write cell by cell output   
    for i=1:length(ICB)  % all possible such packages
        % get the desired output unit number from MFLOW worksheet
        idx=strmatchi(['I',ICB{i},'CB'],MFLOWparnams,1);
        if idx==0, icb_unit=0; else icb_unit=MFLOWparvals(idx); end
        if icb_unit>0
            if ~any(icb_unit==nam.UNIT(Idata))
                error(['Check the NAM data: Missing DATA output unit number <<%d>>\n',...
                    'to write cell-by-cell flow terms for <<%s>> as specified in\n',...
                    'worksheet MFLOW !'],icb_unit,ICB{i});
            end
        end
    end
end

%% ===== THE BAS FILE ====================
bas.SCEN=nam.SCEN(strmatch('BAS6',nam.PCKG));

if bas.SCEN
    fprintf('Generating basic struct\n');
    bas.unit=nam.UNIT(strmatch('BAS6',nam.PCKG));
    bas.ext =nam.EXT {strmatch('BAS6',nam.PCKG)};   
    bas.NLAY=size(IBOUND,3); 
    bas.HNOFLO =MFLOWparvals(strmatchi('HNOFLO',MFLOWparnams),bas.SCEN);
    bas.IBOUND=IBOUND;
    bas.STRTHD=STRTHD;        % initial heads zero

    writeBAS6(basename,bas);
end


%% ===== THE DIS-file struct ==============
% STRESS PERIODS
dis.SCEN=nam.SCEN(strmatch('DIS',nam.PCKG));

if dis.SCEN
    fprintf('Generating discretization struct\n');
    dis.unit=nam.UNIT(strmatch('DIS',nam.PCKG));
    dis.ext =nam.EXT {strmatch('DIS',nam.PCKG)};
    dis.LAYCBD = LAYparvals(:,strmatchi('LAYCBD',LAYparnams));  % confining bed below this layer
    
    % Allow user preference of using dx,Dx,DX or DELX instead of DELC
    % and                            dy,Dy,DY or DELR instead of DELR
    if exist('dx'  ,'var'), DELC=dx;   fprintf('Setting DELC=dx\n');   end
    if exist('dy'  ,'var'), DELR=dy;   fprintf('Setting DELR=dy\n');   end
    if exist('Dx'  ,'var'), DELC=Dx;   fprintf('Setting DELC=Dx\n');   end % overwrites if present
    if exist('Dy'  ,'var'), DELR=Dy;   fprintf('Setting DELR=Dy\n');   end
    if exist('DX'  ,'var'), DELC=DX;   fprintf('Setting DELC=DX\n');   end % overwrites if present
    if exist('DY'  ,'var'), DELR=DY;   fprintf('Setting DELR=DY\n');   end
    if exist('DELX','var'), DELC=DELX; fprintf('Setting DELC=DELX\n'); end % overwrites if present
    if exist('DELY','var'), DELR=DELY; fprintf('Setting DELR=DELY\n'); end

%% Get size of model from mf_adapt.m
    if length(size(DELR))==3,  % if DELR is given as full 3D array
        dis.DELR=squeeze(DELR(:,1,1)); % if given as vector
    else
        dis.DELR=DELR;
    end % DELR from mf_adapt
    
    if length(size(DELC))==3, % if DELC is given as full 3D array
        dis.DELC=squeeze(DELC(1,:,1)); % else, if given as vector
    else
        dis.DELC=DELC;
    end
    
    dis.NROW=length(dis.DELR);  % default matlab variable
    dis.NCOL=length(dis.DELC);  % default matlab variable
    dis.NLAY=size(Z,3)-1-sum(dis.LAYCBD(1:end-1));
    dis.NPER=size(PERparvals,1);
        
    dis.ITMUNI =MFLOWparvals(strmatchi('ITMUNI',MFLOWparnams),dis.SCEN);
    dis.LENUNI =MFLOWparvals(strmatchi('LENUNI',MFLOWparnams),dis.SCEN);  % confining bed below this layer

    dis.LAYCBD = LAYparvals(:,strmatchi('LAYCBD',LAYparnams));  % confining bed below this layer

    % Layer elevations
    NQ3D=sum(LAYparvals(:,strmatchi('LAYCBD',LAYparnams)));

    dis.Z=Z;  % Use Z from mf_adapt

    % Generating stress period info for discretizaton file: 
    dis.PERLEN = PERparvals(:,strmatchi('PERLEN',   PERparnams));
    dis.NSTP   = PERparvals(:,strmatchi('NSTP',     PERparnams));
    dis.TSMULT = PERparvals(:,strmatchi('TSMULT',   PERparnams));
    dis.isTran = PERparvals(:,strmatchi('Transient',PERparnams));
    
    writeDIS(basename,dis);
end


%% ===== THE BCF-file =====================
bcf.SCEN=nam.SCEN(strmatch('BCF6',nam.PCKG));

if bcf.SCEN    
    fprintf('Generating BCF struct and file\n');
    bcf.unit  =nam.UNIT(strmatch('BCF6',nam.PCKG));
    bcf.ext   =nam.EXT {strmatch('BCF6',nam.PCKG)};
    bcf.NROW  =dis.NROW;
    bcf.NCOL  =dis.NCOL;
    bcf.NLAY  =dis.NLAY;
    bcf.DELC  =dis.DELC;
    bcf.DELR  =dis.DELR;
    bcf.LAYCBD=dis.LAYCBD;
    bcf.isTran=dis.isTran;
    bcf.Z      =Z; % all top and bottom elevations of the model grid
    
    bcf.IBCFCB =MFLOWparvals(strmatchi('IBCFCB',MFLOWparnams),bcf.SCEN);
    bcf.HDRY   =MFLOWparvals(strmatchi('HDRY'  ,MFLOWparnams),bcf.SCEN);
    bcf.IWDFLG =MFLOWparvals(strmatchi('IWDFLG',MFLOWparnams),bcf.SCEN);
    bcf.WETFCT =MFLOWparvals(strmatchi('WETFCT',MFLOWparnams),bcf.SCEN);
    bcf.IWETIT =MFLOWparvals(strmatchi('IWETIT',MFLOWparnams),bcf.SCEN);
    bcf.IHDWET =MFLOWparvals(strmatchi('IHDWET',MFLOWparnams),bcf.SCEN);

    bcf.LAYCON =LAYparvals(:,strmatchi('LAYCON' ,LAYparnams));
    bcf.LAYAVG =LAYparvals(:,strmatchi('LAYAVG' ,LAYparnams));
    bcf.TPRY   =LAYparvals(:,strmatchi('TPRY'   ,LAYparnams));
    bcf.WETDRY =LAYparvals(:,strmatchi('WETDRY' ,LAYparnams));
    
    % allow specifiation of WETDRY as cell array also in mf_adapt
    if exist('WETDRY','var'),     % if not, we keep bcf.WETDRY as is
        if ~iscell(WETDRY),
            error('WETDRY specified in mf_adapt must be a cell array with one model layer WETDRY values per array cell');
        end
        wetdry=bcf.WETDRY;        % store WETDRY from spreadsheet LAY
        bcf.WETDRY=cell(size(WETDRY));
        for iL=1:size(bcf.WETDRY,1)  % may be smaller than number of layers
            if wetdry(iL)==0,        % this triggers use of mf_adapt values
                bcf.WETDRY{iL}=WETDRY{iL};  % default
            else
                bcf.WETDRY{iL}=wetdry(iL);
            end
        end
    end

    
    if exist('HY'  ,'var'), bcf.HY   =HY;  end  % if laycon is 1 or 3 for those layers HY   must exist
    if exist('TRAN','var'), bcf.TRAN=TRAN; end  % if laycon is 0 or 2 for those layers TRAN must exist

    if bcf.NLAY>1, bcf.VCONT=VCONT; end % must always exist
    if any(bcf.isTran)
        bcf.SF1=SF1;          % for confined and (the top) unconfined layer
        if any(bcf.LAYCON==2 | bcf.LAYCON==3)
            bcf.SF2=SF2;      % only for rewettable layers
        end
    end

    writeBCF(basename,bcf) 
end

%% ===== THE lpf-file ==================
lpf.SCEN=nam.SCEN(strmatch('LPF',nam.PCKG));

if lpf.SCEN
    fprintf('Generating LPF struct and file\n');
    lpf.unit=nam.UNIT(strmatch('LPF',nam.PCKG));
    lpf.ext =nam.EXT {strmatch('LPF',nam.PCKG)};
    lpf.NLAY=dis.NLAY;
    lpf.NROW=dis.NROW;
    lpf.NCOL=dis.NCOL;

    lpf.ILPFCB =MFLOWparvals(strmatchi('ILPFCB',MFLOWparnams),lpf.SCEN); % file unit number for saving cell by cell flow terms   

    lpf.HDRY   =MFLOWparvals(strmatchi('HDRY'  ,MFLOWparnams),lpf.SCEN); % recharge option, in what layer
    lpf.NPLPF  =MFLOWparvals(strmatchi('NPLPF' ,MFLOWparnams),lpf.SCEN); % recharge option, in what layer
    lpf.WETFCT =MFLOWparvals(strmatchi('WETFCT',MFLOWparnams),lpf.SCEN);
    lpf.IWETIT =MFLOWparvals(strmatchi('IWETIT',MFLOWparnams),lpf.SCEN);
    lpf.IHDWET =MFLOWparvals(strmatchi('IHDWET',MFLOWparnams),lpf.SCEN);

    
    lpf.LAYCON =LAYparvals(:,strmatchi('LAYCON' ,LAYparnams)); % Layer type 1=convertible
    lpf.LAYAVG =LAYparvals(:,strmatchi('LAYAVG' ,LAYparnams)); % Hor cond comp method (0=harmonic)
    lpf.LAYWET =LAYparvals(:,strmatchi('LAYWET' ,LAYparnams)); % Layer wettability flag
    
    lpf.LAYCBD =LAYparvals(:,strmatchi('LAYCBD' ,LAYparnams)); % Confining unit present flag
    if any(lpf.LAYCBD>0)
        try
            lpf.VKCB   =VKCB; % vertical conductivity of confining units must be specified in matlab
        catch  ME
            fprintf('If conf. unit present (any LAYCBD>0), then you Must specify VKCB in workspace');
            rethrow(ME);
        end
    end
    
    lpf.WETDRY =LAYparvals(:,strmatchi('WETDRY' ,LAYparnams)); % Layer wetting method flag
    
    % allow specifiation of WETDRY as cell array also in mf_adapt
    if exist('WETDRY','var'),     % if not, we keep bcf.WETDRY as is
        if ~iscell(WETDRY),
            error('WETDRY specified in mf_adapt must be a cell array with one model layer WETDRY values per array cell');
        end
        wetdry=lpf.WETDRY;        % store WETDRY from spreadsheet LAY
        lpf.WETDRY=cell(size(WETDRY));
        for iL=1:size(lpf.WETDRY,1)  % may be smaller than number of layers
            if wetdry(iL)==0,        % this triggers use of mf_adapt values
                lpf.WETDRY{iL}=WETDRY{iL};  % default
            else
                lpf.WETDRY{iL}=wetdry(iL);
            end
        end
    end            

    lpf.CHANI  =LAYparvals(:,strmatchi('CHANI'  ,LAYparnams)); % Layer wetting method flag
    lpf.LAYVKA =LAYparvals(:,strmatchi('LAYVKA' ,LAYparnams)); % 0=use kV
    
    % Vertical conductivity and anisotropy
    if exist('VK','var')
        if ~exist('VKA','var')
            VKA=VK;
        else  % if both, overwrite VKA with VK according to LAYVKA
            for iz=1:NZ
                if lpf.LAYVKA(iz)==0
                    lpf.VKA(:,:,iz)=HK(:,:,iz);
                end
            end
        end
        clear VK;
    else
        if ~exist('VKA','var')
            error('You must specify VKA or VK in mf_adapt');
        end
    end
    lpf.VKA =VKA;

    % horizontal conductivity in x-direction
    if exist('HK'  ,'var'),   lpf.HK  =HK;
    else error('You must dfined HK in mf_adapt if you use the LPF package');
    end
    
    if exist('HANI','var'), lpf.HANI=HANI;
    elseif any(lpf.CHANI<=0)
        error('You must define HANI if any CHANI<0 in the lay worksheet');
    end
    
    lpf.isTran = dis.isTran;  % whether or not transient flow
    
    if any(lpf.isTran)
        if ~exist('SS','var'), error('Transient model needs SS'); end
        lpf.SS=SS;             % always if transient
        if any(lpf.LAYCON),
            if ~exist('SY','var'), error('Transient model needs SY'); end
            lpf.SY=SY;         % only if convertable (LAYCON=[1,2,3])
        end
    end
   

    lpf.LAYWET(lpf.LAYCON==0)=0; % LAYWET must be 0 if LAYCON is 0

    writeLPF(basename,lpf) 
end

%% ===== the VDF-file for SEAWAT ================================
vdf.SCEN=nam.SCEN(strmatch('VDF',nam.PCKG));

if vdf.SCEN>0
    fprintf('Generating VDF struct and file\n');
    vdf.unit=nam.UNIT(strmatch('VDF',nam.PCKG));
    vdf.ext =nam.EXT {strmatch('VDF',nam.PCKG)};
    vdf.NLAY=dis.NLAY;
    vdf.NPER=dis.NPER;
    
    % for each simulation
    %1  Note alternative parameter names depending on version
    if strmatch('MT3DRHOFLG',SEAWATparnams)
        vdf.MT3DRHOFLG=SEAWATparvals(strmatchi('MT3DRHOFLG',SEAWATparnams),vdf.SCEN); % recharge option, in what layer
    else
        vdf.MT3DRHOFLG=SEAWATparvals(strmatchi('MTDNCONC'   ,SEAWATparnams),vdf.SCEN); % recharge option, in what layer
    end
    vdf.MFNADVFD=SEAWATparvals(strmatchi('MFNADVFD',SEAWATparnams),vdf.SCEN); 
    vdf.NSWTCPL =SEAWATparvals(strmatchi('NSWTCPL' ,SEAWATparnams),vdf.SCEN); 
    vdf.IWTABLE =SEAWATparvals(strmatchi('IWTABLE' ,SEAWATparnams),vdf.SCEN); 

    %2
    vdf.DENSEMIN =SEAWATparvals(strmatchi('DENSEMIN',SEAWATparnams),vdf.SCEN); 
    vdf.DENSEMAX =SEAWATparvals(strmatchi('DENSEMAX',SEAWATparnams),vdf.SCEN); 
 
    %3
    vdf.DNSCRIT =SEAWATparvals(strmatchi('DNSCRIT' ,SEAWATparnams),vdf.SCEN);
    
    if vdf.MT3DRHOFLG>=0
        %4
        vdf.DENSEREF=SEAWATparvals(strmatchi('DENSEREF',SEAWATparnams),vdf.SCEN); 
        vdf.DRHODC  =SEAWATparvals(strmatchi('DRHODC(1)',SEAWATparnams),vdf.SCEN);
        if vdf.MT3DRHOFLG==0
            %6
            vdf.INDENSE =PERparvals(:,strmatchi('INDENSE',PERparnams)); % if >0 read dense for each period
            if any(vdf.INDENSE)>0
                %7
                vdf.DENSE=DENSE;  % get DENS{NPER}(NROW,NCOL) from mf_adapt
            end
        end

    elseif vdf.MT3DRHOFLG==-1
    %4a
        vdf.DENSEREF  =SEAWATparvals(strmatchi('DENSEREF' ,SEAWATparnams),vdf.SCEN);
        vdf.DRHODPRHD =SEAWATparvals(strmatchi('DRHODPRHD',SEAWATparnams),vdf.SCEN); 
        vdf.PRHDREF   =SEAWATparvals(strmatchi('PRHDREF'  ,SEAWATparnams),vdf.SCEN); 
        %4b
        vdf.NSRHOEOS  =SEAWATparvals(strmatchi('NSRHOEOS' ,SEAWATparnams),vdf.SCEN); 
        %4c
        vdf.MTRHOSPEC  =SEAWATparvals(strmatchi('MTRHOSPEC',SEAWATparnams),1:vdf.NSRHOEOS);
        vdf.DRHODC     =SEAWATparvals(strmatchi('DRHODC'   ,SEAWATparnams,'exact'),1:vdf.NSRHOEOS); 
        vdf.CRHOREF    =SEAWATparvals(strmatchi('CRHOREF'  ,SEAWATparnams),1:vdf.NSRHOEOS); 
    end
    %5
    vdf.FIRSTDT =SEAWATparvals(strmatchi('FIRSTDT' ,SEAWATparnams),vdf.SCEN); 
     
    writeVDF(basename,vdf) 
end

%% ===== the VSC-file for SEAWAT ================================
vsc.SCEN=nam.SCEN(strmatch('VSC',nam.PCKG));

if vsc.SCEN
    fprintf('Generating VSC struct and file\n');
    vsc.unit=nam.UNIT(strmatch('VSC',nam.PCKG));
    vsc.ext =nam.EXT {strmatch('VSC',nam.PCKG)};
    vsc.NLAY=dis.NLAY;
    vsc.NPER=dis.NPER;
    
    % for each simulation
    %1
    vsc.MT3DMUFLG=SEAWATparvals(strmatchi('MT3DMUFLG',SEAWATparnams),vsc.SCEN); % recharge option, in what layer

    %2
    vsc.VISCMIN =SEAWATparvals(strmatchi('VISCMIN',SEAWATparnams),vsc.SCEN); 
    vsc.VISCMAX =SEAWATparvals(strmatchi('VISCMAX',SEAWATparnams),vsc.SCEN); 
    
    if vsc.MT3DMUFLG>=0
    %3
        vsc.VISCREF =SEAWATparvals(strmatchi('VISCREF',SEAWATparnams),vsc.SCEN); 
        vsc.DMUDC   =SEAWATparvals(strmatchi('DMUDC'  ,SEAWATparnams),vsc.SCEN); 
        vsc.CMUREF  =SEAWATparvals(strmatchi('CMUREF' ,SEAWATparnams),vsc.SCEN);
        
        if vsc.MT3DMUFLG==0
            %4
                vsc.INVISC = PERparvals(:,strmatchi('INVISC',   PERparnams));
                if any(vsc.INVISC)>0
                    %5
                    vsc.VISC=VISC; % get VISC{NPER}(NROW,NCOL) from mf_adapt
                end
        end

    elseif vsc.MT3DMUFLG==-1
        %3a
        vsc.VISCREF   =SEAWATparvals(strmatchi('VISCREF',SEAWATparnams),vsc.SCEN); 
        vsc.NSMUEOS   =SEAWATparvals(strmatchi('NSMUEOS'  ,SEAWATparnams),vsc.SCEN); 
        vsc.MUTEMPOPT =SEAWATparvals(strmatchi('MUTEMPOPT',SEAWATparnams),vsc.SCEN); 
        %3c
        vsc.MTMUSPEC=SEAWATparvals(strmatchi('MTMUSPEC' ,SEAWATparnams),1:vsc.NSMUEOS);
        vsc.DMUDC   =SEAWATparvals(strmatchi('DMUDC'    ,SEAWATparnams),1:vsc.NSMUEOS);
        vsc.CMUREF  =SEAWATparvals(strmatchi('CMUREF'   ,SEAWATparnams),1:vsc.NSMUEOS);
        if vsc.MUTEMPOPT>0
            %3d
            vsc.MTMUTEMPSPEC=SEAWATparvals(strmatchi('MTMUTEMPSPEC',SEAWATparnams),vsc.SCEN);
            vsc.AMUCOEF     =SEAWATparvals(strmatchi('AMUCOEF'     ,SEAWATparnams),vsc.MUTEMPOPT);
            vsc.AMUCOEF(isnan(vsc.AMUCOEF))=[];
        end
    end

    writeVSC(basename,vsc) 
end

%% ===== THE rch-file ============================
rch.SCEN=nam.SCEN(strmatch('RCH',nam.PCKG));

if rch.SCEN
    fprintf('Generating RCH struct and file\n');
    rch.unit=nam.UNIT(strmatch('RCH',nam.PCKG));
    rch.ext =nam.EXT {strmatch('RCH',nam.PCKG)};
    rch.NLAY=dis.NLAY;
    rch.NROW=dis.NROW;
    rch.NCOL=dis.NCOL;
    rch.NPER=size(PERparvals,1);
    % for each simulation
    %1
    rch.NPRCH=0;         % we will use no parameters
    %2    
    rch.NRCHOP =MFLOWparvals(strmatchi('NRCHOP',MFLOWparnams),rch.SCEN); % recharge option, in what layer
    rch.IRCHCB =MFLOWparvals(strmatchi('IRCHCB',MFLOWparnams),rch.SCEN); % file unit number for saving cell by cell flow terms
    
    %3
    %4
    
    %% for each stress period
    %5  
    
    rch.INRECH =PERparvals(:,strmatchi('INRECH',PERparnams)); % if >0 read recharge layer
    rch.INIRCH =PERparvals(:,strmatchi('INIRCH',PERparnams)); % skipped unless recharge in specific layers

    rch.rech   =PERparvals(:,strmatchi('RECH'  ,PERparnams));
    rch.irch   =PERparvals(:,strmatchi('IRCH'  ,PERparnams));

    %6
    if any(rch.INRECH==0) && ~exist('RECH','var'); % RECH must by defined in mf_adapt
        error('If any IN%s in the PER worksheet = 0 then you must define %s in mf_adapt !','RECH','RECH');
    end
    for iPer=1:rch.NPER
        if rch.INRECH(iPer)==0
            if iscell(RECH)
                rch.RECH{iPer}=RECH{iPer};
            else
                rch.RECH{iPer}=RECH(:,:,iPer);
            end
        else
            rch.RECH{iPer}=rch.rech(iPer);
        end
    end

    %7
    %8
    if any(rch.INIRCH==0) && ~exist('IRCH','var')
        error('If any IN%s in the PER worksheet = 0 then you must define %s in mf_adapt !',...
                'IRCH','IRCH');
    end
    for iPer=1:rch.NPER
        if rch.INIRCH(iPer)==0
            if iscell(rch.IRCH)
                rch.IRCH{iPer}=IRCH{iPer};
            else
                rch.IRCH{iPer}=IRCH(:,:,iPer);
            end
        else
            rch.IRCH{iPer}=rch.irch(iPer);
        end
    end
    
    writeRCH(basename,rch) 
end

%% s===== THE evt-file (the Evaporation package) ==============
evt.SCEN=nam.SCEN(strmatch('EVT',nam.PCKG));

if evt.SCEN
    fprintf('Generating EVT struct and file\n');
    evt.unit=nam.UNIT(strmatch('EVT',nam.PCKG));
    evt.ext =nam.EXT {strmatch('EVT',nam.PCKG)};
    evt.NLAY=dis.NLAY;
    evt.NROW=dis.NROW;
    evt.NCOL=dis.NCOL;
    evt.NPER=dis.NPER;
    evt.NPER=size(PERparvals,1);
    
    evtp.NPEVT=0;
    %1
    %2
    evt.NEVTOP =MFLOWparvals(strmatchi('NEVTOP',MFLOWparnams),evt.SCEN);
    evt.IEVTCB =MFLOWparvals(strmatchi('IEVTCB',MFLOWparnams),evt.SCEN); % file unit number for saving cell by cell flow terms
    %3
    %4
    %5
    evt.INSURF =PERparvals(:,strmatchi('INSURF',PERparnams));
    evt.INEVTR =PERparvals(:,strmatchi('INEVTR',PERparnams));
    evt.INEXDP =PERparvals(:,strmatchi('INEXDP',PERparnams));
    evt.INIEVT =PERparvals(:,strmatchi('INIEVT',PERparnams));
    
    evt.surf   =PERparvals(:,strmatchi('SURF',PERparnams));
    evt.evtr   =PERparvals(:,strmatchi('EVTR'  ,PERparnams));
    evt.exdp   =PERparvals(:,strmatchi('EXDP'  ,PERparnams));
    evt.ievt   =PERparvals(:,strmatchi('IEVT'  ,PERparnams));
    
    evt.SURF   =cell(evt.NPER,1);
    evt.EVTR   =cell(evt.NPER,1);
    evt.EXDP   =cell(evt.NPER,1);
    evt.IEVT   =cell(evt.NPER,1);
    
    NPER=dis.NPER;
    
    %6
    if any(evt.INSURF==0) && ~exist('SURF','var')
                   error('If any IN%s in the PER worksheet = 0 then you must define EVT%s in mf_adapt !','SURF');
    end
    if evt.INSURF(1)==-2; evt.SURF{1}=evt.surf(1)+Z(:,:,1); evt.INSURF(1)=2; end
    if evt.INSURF(1)==-1; evt.SURF{1}=evt.surf(1);          evt.INSURF(1)=1; end
    for iPer=1:NPER
        if evt.INSURF(iPer)==0
            if iscell(SURF)
                evt.SURF{iPer}=SURF{iPer};
            else
                evt.SURF{iPer}=SURF(:,:,iPer);
            end
        else
            evt.SURF{iPer}=evt.surf(iPer);
        end
    end
    %7
    if any(evt.INEVTR==0) && ~exist('EVTR','var')
        error('If any IN%s in the PER worksheet = 0 then you must define %s in mf_adapt !','EVTR');
    end
    for iPer=1:NPER
        if evt.INEVTR(iPer)==0
            if iscell(EVTR)
                evt.EVTR{iPer} =EVTR{iPer};
            else
                evt.EVTR{iPer} =EVTR(:,:,iPer);
            end
        else
            evt.EVTR{iPer}=evt.evtr(iPer);
        end
    end
    %8 parameters not yet implemented
    %9
    if any(evt.INEXDP==0) && ~exist('EXDP','var'); % EVTR must by defined in mf_adapt
        error('If any IN%s in the PER worksheet = 0 then you must define %s in mf_adapt !','EXDP');
    end
    if evt.INEXDP(1)==-1; evt.EXDP{1}=evt.exdp(1); evt.INEXDP(1)=1; end
    for iPer=1:NPER
        if evt.INEXDP(iPer)==0
            if iscell(EXDP)
                evt.EXDP{iPer}=EXDP{iPer};
            else
                evt.EXDP{iPer}=EXDP(:,:,iPer);
            end
        else
            evt.EXDP{iPer}=evt.exdp(iPer);
        end
    end
    
    %10
    if evt.NEVTOP==2
        if any(evt.INIEVT==0) && ~exist('IEVT','var'); % EVTR must by defined in mf_adapt
            error('If any IN%s in the PER worksheet = 0 then you must define %s in mf_adapt !','IEVT');
        end
        for iPer=1:evt.NPER
            if evt.INIEVT(iPer)==0
                if iscell(IEVT)
                    evt.IEVT{iPer}=IEVT(iPer);
                else
                    evt.IEVT{iPer}=IEVT(:,:,iPer);
                end
            else
                evt.IEVT{iPer}=evt.ievt(iPer);
            end
        end
    end

    writeEVT(basename,evt,evtp) 
end


%% ===== THE BOUNDADRY CONDTION files GHB, DRN, RIV, CHD ====================================
bctypes={'WEL','MNW','GHB','DRN','RIV','CHD'};
for ibc=1:length(bctypes)
    BCtype=bctypes{ibc};
    bcn.type=BCtype;
    bcn.SCEN=nam.SCEN(strmatch(BCtype,nam.PCKG));  % check if this boundary package is on

    if bcn.SCEN  % if this boundary package is "on" then ...
        fprintf('Generating %s struct and file\n',BCtype);
        BCN=[];
        
        % if AUX variables are used, the number of columns of the
        % respective bcn  variable must match them
        bcn.AUX=cell(1,5);
        
        % UNIT numbers for cell-by-cell flow term output files, specified in
        % worksheet MFLOW
        switch BCtype
            case 'WEL', if exist('WEL','var'), BCN=WEL; end
                bcn.ICB =MFLOWparvals(strmatchi('IWELCB',MFLOWparnams),bcn.SCEN);
            case 'MNW', if exist('MNW','var'), BCN=MNW; end
                bcn.ICB      =MFLOWparvals(strmatchi('IMNWCB',  MFLOWparnams),bcn.SCEN);
                bcn.IWELPT   =MFLOWparvals(strmatchi('IWELPT',  MFLOWparnams),bcn.SCEN);
                bcn.LOSSTYPE =MFLOWparvals(strmatchi('LOSSTYPE',MFLOWparnams),bcn.SCEN);
                bcn.PLOSSMNW =MFLOWparvals(strmatchi('PLOSSMNW',MFLOWparnams),bcn.SCEN);
            case 'GHB', if exist('GHB','var'), BCN=GHB; end
                bcn.ICB =MFLOWparvals(strmatchi('IGHBCB',MFLOWparnams),bcn.SCEN);
            case 'DRN', if exist('DRN','var'), BCN=DRN; end
                bcn.ICB =MFLOWparvals(strmatchi('IDRNCB',MFLOWparnams),bcn.SCEN);                
            case 'RIV', if exist('RIV','var'), BCN=RIV; end
                bcn.ICB =MFLOWparvals(strmatchi('IRIVCB',MFLOWparnams),bcn.SCEN);            
            case 'CHD', if exist('CHD','var'), BCN=CHD; end
                bcn.ICB =MFLOWparvals(strmatchi('ICHDCB',MFLOWparnams),bcn.SCEN);
                % possible seawat options
                if exist('CHDDENSOPT','var'), bcn.AUX{1}='CHDDENSOPT'; end
                if exist('CHDDENS'   ,'var'), bcn.AUX{2}='CHDDENS'   ; end
        end
        
        idx=strmatchi(BCtype,nam.PCKG,'noerr');
        bcn.unit=nam.UNIT(idx); % unit number of input file of this package
        bcn.ext =nam.EXT {idx}; % extention of input file of this package

        [bcn.parnams,bcnfromxls]=getExcelData(XLSF,BCtype,'Horizontal');  % sheet must exist in xls file
        if exist(BCtype,'var')        % BCtype data in spreadsheet and in mf_adapt
            ncolxls  =length(bcn.parnams);
            ncoladapt=size(BCN,2);
            bcnfromadapt=BCN;
            bcn.parvals=[bcnfromxls;bcnfromadapt];
        else                          % BCtype data only in spreadsheet (may be empty)
            bcn.parvals=bcnfromxls;
        end
        
        % Inform user about where he specified the boundaries, in XLSF, in
        % mf_adapt or both or nowhere, while the package in NAM is on.
        if size(bcnfromxls,1)>0
            fprintf('You specified %s in %s file',BCtype,XLSF);
            if exist(BCtype,'var')
                fprintf(' AND you specified %s also in mf_adapt\n',BCtype);
            else
                fprintf('\n');
            end
        elseif exist(BCtype,'var')
          fprintf('You specified %s only in mf_adapt.\n',BCtype);
        else
            fprintf('You neither specified %s in %s nor in mf_adapt!\n',BCtype,XLSF);
            fprintf('However, the %s package is on in the NAM sheet\n',BCtype);
        end
         
        bcnp.NPar=0; % drain parameters not yet implemented
        bcn.NPER =dis.NPER;

        writeBCN(basename,bcn,bcnp)
    end
    % replace bcn by a varible reflecting its role by name
    switch BCtype
        case 'WEL', wel=bcn;
        case 'MNW', mnw=bcn;
        case 'GHB', ghb=bcn;
        case 'DRN', drn=bcn;
        case 'RIV', riv=bcn;
        case 'CHD', chd=bcn;
    end
    clear('bcn'); % and do away with bcn
end



%% ==== THE OC-file (output control package) =======================
oc.SCEN=nam.SCEN(strmatch('OC',nam.PCKG));

if oc.SCEN
    % for each simulation
    fprintf('Generating OC struct and file\n');
    oc.unit=nam.UNIT(strmatch('OC',nam.PCKG));
    oc.ext =nam.EXT {strmatch('OC',nam.PCKG)};
    oc.IHEDFM=MFLOWparvals(strmatchi('IHEDFM',MFLOWparnams),oc.SCEN);
    oc.IDDNFM=MFLOWparvals(strmatchi('IDDNFM',MFLOWparnams),oc.SCEN);
    oc.IHEDUN=MFLOWparvals(strmatchi('IHEDUN',MFLOWparnams),1); % file unit number for saving heads
    oc.IDDNUN=MFLOWparvals(strmatchi('IDDNUN',MFLOWparnams),1); % file unit number for saving drawdowns

    % for each stres period
    oc.INCODE =PERparvals(:,strmatchi('INCODE',PERparnams));
    oc.IHDDFL =PERparvals(:,strmatchi('IHDDFL',PERparnams));
    oc.IBUDFL =PERparvals(:,strmatchi('IBUDFL',PERparnams));
    oc.ICBCFL =PERparvals(:,strmatchi('ICBCFL',PERparnams));
    oc.NPER=size(oc.INCODE,1);
    oc.NSTP=dis.NSTP;
    
    oc.Hdpr =PERparvals(:,strmatchi('Hdpr',PERparnams));
    oc.Ddpr =PERparvals(:,strmatchi('Ddpr',PERparnams));
    oc.Hdsv =PERparvals(:,strmatchi('Hdsv',PERparnams));
    oc.Ddsv =PERparvals(:,strmatchi('Ddsv',PERparnams));

    writeOC(basename,oc);
end


%% ==== THE CFP-file (Conduit FLow Process package) =======================
% under development
% cfp.SCEN=nam.SCEN(strmatch('CFP',nam.PCKG));
% 
% if cfp.SCEN
%     % for each simulation
% fprintf('Generating CFP struct and file\n');
% cfp.unit=nam.UNIT(strmatch('CFP',nam.PCKG));
% cfp.ext =nam.EXT {strmatch('CFP',nam.PCKG)};
% cfp.MODE       =MFLOWparvals(strmatchi('IHEDFM',MFLOWparnams),cfp.SCEN);
% cfp.NNODES     =MFLOWparvals(strmatchi('IDDNFM',MFLOWparnams),oc.SCEN);
% cfp.NPIPES     =
% cfp.NLAYERS    =
% cfp.TEMPERATURE=
% cfp.Connections=
% cfp.GEOHEIGHT=
% cfp.SA_EXCHANGE=
% cfp.EPSILON=
% cfp.NITER=
% cfp.RELAX=
% cfp.P_NR=
% cfp.PIPES= ...
% cfp.NO_N N_HEAD
% cfp.NO_N K_EXCHANGE
% cfp.NCL
% cfp.CL
% cfp.LTEMP
% cfp.VOID LCRITREY_L TCRITREY_L
% cfp.NPER=dis.NPER;
% 
% writeCFP(basename,cfp);
% end

%% ===== THE coc-file (Conduit Flow Process) =====
% coc.SCEN=nam.SCEN(strmatch('COC',nam.PCKG));
% coc.unit=nam.UNIT(strmatch('COC',nam.PCKG));
% coc.ext =nam.EXT {strmatch('COC',nam.PCKG)};

% in mf_adapt prepare this struct

% coc.NODE_NUMBERS (node numbers for output)
% coc.N_NTS = node output time interval (heads and flows)
% coc.PIPE_NUMBERS (pipe numbers for output)
% coc.T_NTS = tube output time interval 

% writeCOC(basename,coc);
% 
% %% ===== THE crch-file (Conduit Flow Process) =====
% crch.SCEN=nam.SCEN(strmatch('CRCH',nam.PCKG));
% crch.unit=nam.UNIT(strmatch('CRCH',nam.PCKG));
% crch.ext =nam.EXT {strmatch('CRCH',nam.PCKG)};

% in mf_adapt prepare this struct
% it has data for each stress period

% crch.IFLAG_CRCH  (1:NPER)); if >0 get crch.PCRCH.CRCH for this stress period

% crch.P_CRCH(1:NPER).CRCH(1:NNODES,1:2)
% the first value is the node number and the second the fraction of diffuse
% recharge routed directy to the nodes
% NNODES is total number of CFP nodes

% writeCRCH(basename,crch);

%% ===== THE pcg-file (P Conjugate solver Package) =====
pcg.SCEN=nam.SCEN(strmatch('PCG',nam.PCKG));

if pcg.SCEN
    fprintf('Generating PCG struct and file\n');
    pcg.unit=nam.UNIT(strmatch('PCG',nam.PCKG));
    pcg.ext =nam.EXT {strmatch('PCG',nam.PCKG)};
    pcg.MXITER=MFLOWparvals(strmatchi('MXITER',MFLOWparnams),pcg.SCEN);
    pcg.ITER1 =MFLOWparvals(strmatchi('ITER1' ,MFLOWparnams),pcg.SCEN);
    pcg.NPCOND=MFLOWparvals(strmatchi('NPCOND',MFLOWparnams),pcg.SCEN);
    pcg.HCLOSE=MFLOWparvals(strmatchi('HCLOSE',MFLOWparnams),pcg.SCEN);
    pcg.RCLOSE=MFLOWparvals(strmatchi('RCLOSE',MFLOWparnams),pcg.SCEN);
    pcg.RELAX =MFLOWparvals(strmatchi('RELAX' ,MFLOWparnams),pcg.SCEN);
    pcg.NBPOL =MFLOWparvals(strmatchi('NBPOL' ,MFLOWparnams),pcg.SCEN);
    pcg.IPRPCG=MFLOWparvals(strmatchi('IPRPCG',MFLOWparnams),pcg.SCEN);
    pcg.MUTPCG=MFLOWparvals(strmatchi('MUTPCG',MFLOWparnams),pcg.SCEN);
    pcg.DAMP  =MFLOWparvals(strmatchi('DAMP'  ,MFLOWparnams),pcg.SCEN);

    writePCG(basename,pcg);
end

%% ===== THE de4-file (P Conjugate solver Package) =====
de4.SCEN=nam.SCEN(strmatch('DE4',nam.PCKG));

if de4.SCEN
    fprintf('Generating DE4 struct and file\n');
    de4.unit=nam.UNIT(strmatch('DE4',nam.PCKG));
    de4.ext =nam.EXT {strmatch('DE4',nam.PCKG)};
    de4.ITMX  =MFLOWparvals(strmatchi('ITMX'  ,MFLOWparnams),de4.SCEN);
    de4.MXUP  =MFLOWparvals(strmatchi('MXUP'  ,MFLOWparnams),de4.SCEN);
    de4.MXLOW =MFLOWparvals(strmatchi('MXLOW' ,MFLOWparnams),de4.SCEN);
    de4.MXBW  =MFLOWparvals(strmatchi('MXBW'  ,MFLOWparnams),de4.SCEN);
    de4.IFREQ =MFLOWparvals(strmatchi('IFREQ' ,MFLOWparnams),de4.SCEN);
    de4.MUTD4 =MFLOWparvals(strmatchi('MUTD4' ,MFLOWparnams),de4.SCEN);
    de4.ACCL  =MFLOWparvals(strmatchi('ACCL'  ,MFLOWparnams),de4.SCEN);
    de4.HCLOSE=MFLOWparvals(strmatchi('HCLOSE',MFLOWparnams),de4.SCEN);
    de4.IPRD4 =MFLOWparvals(strmatchi('IPRD4' ,MFLOWparnams),de4.SCEN);
    
    writeDE4(basename,de4);
end

%% ===== THE SIP-file (Strongly Implicit Procedure Package) =====
sip.SCEN=nam.SCEN(strmatch('SIP',nam.PCKG));

if sip.SCEN
    fprintf('Generating SIP struct and file\n');
    sip.unit=nam.UNIT(strmatch('SIP',nam.PCKG));
    sip.ext =nam.EXT {strmatch('SIP',nam.PCKG)};
    sip.MXITER =MFLOWparvals(strmatchi('MXITER',MFLOWparnams),sip.SCEN);
    sip.NPARM  =MFLOWparvals(strmatchi('NPARM' ,MFLOWparnams),sip.SCEN);
    sip.ACCL  =MFLOWparvals(strmatchi('ACCL'   ,MFLOWparnams),sip.SCEN);
    sip.HCLOSE=MFLOWparvals(strmatchi('HCLOSE' ,MFLOWparnams),sip.SCEN);
    sip.IPCALC=MFLOWparvals(strmatchi('IPCALC' ,MFLOWparnams),sip.SCEN);
    sip.WSEED =MFLOWparvals(strmatchi('WSEED'  ,MFLOWparnams),sip.SCEN);
    sip.IPRSIP=MFLOWparvals(strmatchi('IPRSIP' ,MFLOWparnams),sip.SCEN);
    
    writeSIP(basename,sip);
end

%% ===== THE SOR-file (Sliced Successive Over Relaxation Package) =====
sor.SCEN=nam.SCEN(strmatch('SOR',nam.PCKG));

if sor.SCEN
    fprintf('Generating SOR struct and file\n');
    sor.unit=nam.UNIT(strmatch('SOR',nam.PCKG));
    sor.ext =nam.EXT {strmatch('SOR',nam.PCKG)};
    sor.MXITER =MFLOWparvals(strmatchi('MXITER',MFLOWparnams),sor.SCEN);
    sor.ACCL  =MFLOWparvals(strmatchi('ACCL'   ,MFLOWparnams),sor.SCEN);
    sor.HCLOSE=MFLOWparvals(strmatchi('HCLOSE' ,MFLOWparnams),sor.SCEN);
    sor.IPRSOR=MFLOWparvals(strmatchi('IPRSOR' ,MFLOWparnams),sor.SCEN);
    
    writeSOR(basename,sor);
end


%% ===== THE GMG-file (Geometric Multigrid solverPackage) =====
gmg.SCEN=nam.SCEN(strmatch('GMG',nam.PCKG));

if gmg.SCEN
    fprintf('Generating GMG struct and file\n');
    gmg.unit=nam.UNIT(strmatch('GMG',nam.PCKG));
    gmg.ext =nam.EXT {strmatch('GMG',nam.PCKG)};
    gmg.RCLOSE =MFLOWparvals(strmatchi('RCLOSE',MFLOWparnams),gmg.SCEN);
    gmg.IITER  =MFLOWparvals(strmatchi('IITER'   ,MFLOWparnams),gmg.SCEN);
    gmg.HCLOSE=MFLOWparvals(strmatchi('HCLOSE' ,MFLOWparnams),gmg.SCEN);
    gmg.MXITER=MFLOWparvals(strmatchi('MXITER' ,MFLOWparnams),gmg.SCEN);
    gmg.DAMP=MFLOWparvals(strmatchi('DAMP' ,MFLOWparnams),gmg.SCEN);
    gmg.IADAMP=MFLOWparvals(strmatchi('IADAMP' ,MFLOWparnams),gmg.SCEN);
    gmg.IOUTGMG=MFLOWparvals(strmatchi('IOUTGMG' ,MFLOWparnams),gmg.SCEN);
    gmg.ISM=MFLOWparvals(strmatchi('ISM' ,MFLOWparnams),gmg.SCEN);
    gmg.ISC=MFLOWparvals(strmatchi('ISC' ,MFLOWparnams),gmg.SCEN);
    gmg.RELAX=MFLOWparvals(strmatchi('RELAX' ,MFLOWparnams),gmg.SCEN);
    
    writeGMG(basename,gmg);
end


%% ===== THE LMT6 file (Link file generation LMT6) =====
lmt.SCEN=nam.SCEN(strmatch('LMT',nam.PCKG));

if lmt.SCEN
    lmt.unit=nam.UNIT(strmatch('LMT',nam.PCKG));  % must not be used !!
    lmt.ext =nam.EXT {strmatch('LMT',nam.PCKG)};
    lmt.OUTPUT_FILE_NAME='';
    lmt.OUTPUT_FILE_HEADER='EXTENDED';
    lmt.OUTPUT_FILE_FORMAT='UNFORMATTED';

    writeLMT(basename,lmt);
end

%% ===== END OF INPUT FOR THE FLOW PROCESs ==================

%% ==== MT3DSM input ==================================
% partly from the MT3D sheet in the spreadsheet


%% THE BTN-file (Basic Transprot Process for MT3D)
btn.SCEN=nam.SCEN(strmatch('BTN',nam.PCKG));

if btn.SCEN    
    btn.NCOMP=MT3Dparvals(strmatchi('NCOMP',MT3Dparnams),btn.SCEN);
    btn.MCOMP=MT3Dparvals(strmatchi('MCOMP',MT3Dparnams),btn.SCEN);
    
    fprintf('Generating Basic Transport Process struct\n');
    btn.unit=nam.UNIT(strmatch('BTN',nam.PCKG));
    btn.ext =nam.EXT {strmatch('BTN',nam.PCKG)};
    %3
    btn.NLAY=dis.NLAY;   % # layers
    btn.NROW=dis.NROW;   % # rows
    btn.NCOL=dis.NCOL;   % # columns
    btn.NPER=size(PERparvals,1);   % # stress periods
    %4
    btn.TUNIT='DAY';
    btn.LUNIT='M';
    btn.MUNIT='KG';
    
    %5
    %Packge-use flags: (ADV DSP SSM RCT GCG XXX XXX XXX XXX XXX)';
    i=strmatch('ADV',nam.PCKG); if ~isempty(i), adv.SCEN=nam.SCEN(i); else adv.SCEN=0; end
    i=strmatch('DSP',nam.PCKG); if ~isempty(i), dsp.SCEN=nam.SCEN(i); else dsp.SCEN=0; end
    i=strmatch('SSM',nam.PCKG); if ~isempty(i), ssm.SCEN=nam.SCEN(i); else ssm.SCEN=0; end
    i=strmatch('RCT',nam.PCKG); if ~isempty(i), rct.SCEN=nam.SCEN(i); else rct.SCEN=0; end
    i=strmatch('GCG',nam.PCKG); if ~isempty(i), gcg.SCEN=nam.SCEN(i); else gcg.SCEN=0; end

    % active flags then become
    btn.TRNOP=[adv.SCEN,dsp.SCEN,ssm.SCEN,rct.SCEN,gcg.SCEN,0,0,0,0,0];
    
    %6
    btn.LAYCON=LAYparvals(:,strmatchi('LAYCON',LAYparnams));
    %7
    btn.DELR=dis.DELR;
    %8
    btn.DELC=dis.DELC;
    %9 %10
    btn.Z=Z;
    %11
    btn.PRSITY=PEFF;
    %12
    btn.ICBUND=ICBUND;  % 0 inactive, 1 active, -1 fixed concentration 
    %13
    btn.STCONC=STCONC;
    %14
    btn.CINACT=MT3Dparvals(strmatchi('CINACT',MT3Dparnams),btn.SCEN);
    btn.THKMIN=MT3Dparvals(strmatchi('THKMIN',MT3Dparnams),btn.SCEN);
    %15
    btn.IFMTCN=MT3Dparvals(strmatchi('IFMTCN',MT3Dparnams),btn.SCEN);
    btn.IFMTNP=MT3Dparvals(strmatchi('IFMTNP',MT3Dparnams),btn.SCEN);
    btn.IFMTRF=MT3Dparvals(strmatchi('IFMTRF',MT3Dparnams),btn.SCEN);
    btn.IFMTDP=MT3Dparvals(strmatchi('IFMTDP',MT3Dparnams),btn.SCEN);
    btn.SAVUCN=MT3Dparvals(strmatchi('SAVUCN',MT3Dparnams),btn.SCEN);
    %16 %17
    % Note if TIMPRS(1)==0 it will be truncated in writeBTN
    btn.NPRS  =MT3Dparvals(strmatchi('NPRS',MT3Dparnams),btn.SCEN);
    btn.TIMPRS=MT3Dparvals(strmatchi('TIMPRS',MT3Dparnams),:); % record with times to print output
    %18
    btn.NPROBS=MT3Dparvals(strmatchi('NPROBS',MT3Dparnams),btn.SCEN);
     %19  skipped, are the locations of the observation points
    %% OBSERVATION POINTS FOR CONCENTRATION
    %if btn.NOBS>0
    fprintf('Getting observation points\n');
    [btn.OBSparnams,btn.OBS]=getExcelData(XLSF,'BTNOBS','Horizontal');
    if size(btn.OBS,1)>0
        fprintf('%d conc observation points read from worksheet BTNOBS\n',size(btn.OBS,1));
    end
    if exist('BTNOBS','var')
        fprintf('%d conc observation points obtained from mf_adapt\n',size(BTNOBS,1));
        if size(btn.OBS,1)>0  % get extra BNTOBS dirctly from mf_adapt
            btn.OBS=[btn.OBS;BTNOBS];
        else
            btn.OBS=BTNOBS; % read directly from mf_adapt
        end
    end
    %end
    %20
    btn.CHKMAS=MT3Dparvals(strmatchi('CHKMAS',MT3Dparnams),btn.SCEN);
    btn.NPRMAS=MT3Dparvals(strmatchi('NPRMAS',MT3Dparnams),btn.SCEN);
    %21
    btn.PERLEN = dis.PERLEN;
    btn.NSTP   = dis.NSTP;
    btn.TSMULT = dis.TSMULT;
    %22 if TSMULT<0 skip
    %23
    btn.DT0    =PERparvals(:,strmatchi('DT0'    ,PERparnams)); % i.e. let the model choose
    btn.MXSTRN =PERparvals(:,strmatchi('MXSTRN' ,PERparnams)); % maximum number of transport steps within one flow time step
    btn.TTSMULT=PERparvals(:,strmatchi('TTSMULT',PERparnams)); % succuessive transport time step multiplier
    btn.TTSMAX =PERparvals(:,strmatchi('TTSMAX' ,PERparnams)); %  use default, set no maximum
    
    writeBTN(basename,btn);
end


%% ===== THE ADV-file (Advection process) ==========
adv.SCEN=nam.SCEN(strmatch('ADV',nam.PCKG));

if adv.SCEN   
    fprintf('Generating Advection Process struct\n');
    adv.unit=nam.UNIT(strmatch('ADV',nam.PCKG));
    adv.ext =nam.EXT {strmatch('ADV',nam.PCKG)};
    %3
    adv.MIXELM=MT3Dparvals(strmatchi('MIXELM' ,MT3Dparnams),adv.SCEN);   % standard finite difference method
    adv.PERCEL=MT3Dparvals(strmatchi('PERCEL' ,MT3Dparnams),adv.SCEN);
    adv.MXPART=MT3Dparvals(strmatchi('MXPART' ,MT3Dparnams),adv.SCEN);
    adv.NADVFD=MT3Dparvals(strmatchi('NADVFD' ,MT3Dparnams),adv.SCEN);
    %2
    adv.ITRACK=MT3Dparvals(strmatchi('ITRACK' ,MT3Dparnams),adv.SCEN);   % particle tracking algorith flag 3=good compromise
    adv.WD    =MT3Dparvals(strmatchi('WD'     ,MT3Dparnams),adv.SCEN);   % concentration weighting factor, 05 should be adequat
    %3
    adv.DCEPS =MT3Dparvals(strmatchi('DCEPS'  ,MT3Dparnams),adv.SCEN);
    adv.NPLANE=MT3Dparvals(strmatchi('NPLANE' ,MT3Dparnams),adv.SCEN);
    adv.NPL   =MT3Dparvals(strmatchi('NPL'    ,MT3Dparnams,'exact'),adv.SCEN);
    adv.NPH   =MT3Dparvals(strmatchi('NPH'    ,MT3Dparnams),adv.SCEN);
    adv.NPMIN =MT3Dparvals(strmatchi('NPMIN'  ,MT3Dparnams),adv.SCEN);
    adv.NPMAX =MT3Dparvals(strmatchi('NPMAX'  ,MT3Dparnams),adv.SCEN);
    %4
    adv.INTERP =MT3Dparvals(strmatchi('INTERP',MT3Dparnams),adv.SCEN);
    adv.NLSINK =MT3Dparvals(strmatchi('NLSINK',MT3Dparnams),adv.SCEN);
    adv.NPSINK =MT3Dparvals(strmatchi('NPSINK',MT3Dparnams),adv.SCEN);
    
    adv.DCHMOC =MT3Dparvals(strmatchi('DCHMOC',MT3Dparnams),adv.SCEN);
    writeADV(basename,adv);
end


%% ===== THE DSP-file (Dispersion process package) =============
dsp.SCEN=nam.SCEN(strmatch('DSP',nam.PCKG));

if dsp.SCEN
   fprintf('Generating Dispersion Process struct\n');
    dsp.unit=nam.UNIT(strmatch('DSP',nam.PCKG));
    dsp.ext =nam.EXT {strmatch('DSP',nam.PCKG)};
    dsp.NCOMP=btn.NCOMP;
    %3
    dsp.NLAY=dis.NLAY;
    dsp.NROW=dis.NROW;
    dsp.NCOL=dis.NCOL;
    
    dsp.MultiDiffusion=MT3Dparvals(strmatchi('MULTIDIFFUSION',MT3Dparnams),dsp.SCEN);
    
    % if AL exists in mf_adapt, use it. It must be a cell array iLay long
    % each cell has the data for one layer (NROW,NCOL) or a single value
    if exist('AL','var'), dsp.AL=AL; else dsp.AL=cell(dsp.NLAY,1); end

    dsp.al    =LAYparvals(:,strmatchi('AL'    ,LAYparnams));
    for iLay=1:dsp.NLAY
        if dsp.al(iLay)>=0,  % use spreadsheet value if value >=0
            dsp.AL{iLay}=dsp.al(iLay);
        end
    end
    
    dsp.TRPT  =LAYparvals(:,strmatchi('TRPT'  ,LAYparnams));  % aTH/aL
    dsp.TRPV  =LAYparvals(:,strmatchi('TRPV'  ,LAYparnams));  % aTV/aL
        dmcoef=LAYparvals(:,strmatchi('DMCOEF',LAYparnams));  % effctief molecular diffusion coefficient
    if size(dmcoef,2)>dsp.NCOMP
        dmcoef=dmcoef(:,1:dsp.NCOMP);
    end
    if exist('DMCOEF','var'), dsp.DMCOEF=DMCOEF; else dsp.DMCOEF=cell(dsp.NLAY,dsp.NCOMP); end
    for iComp=1:dsp.NCOMP
        for iLay=1:dsp.NLAY
            if dmcoef(iLay,iComp)>=0,
                dsp.DMCOEF{iLay,iComp}=dmcoef(iLay,iComp);
            end
        end
    end
    writeDSP(basename,dsp);
end


%% ===== THE SSM-file (Source-sink mixing process package) =====
ssm.SCEN=nam.SCEN(strmatch('SSM',nam.PCKG));

if ssm.SCEN
    fprintf('Generating Source-Sink Mixing Process struct\n');
    ssm.unit=nam.UNIT(strmatchi('SSM',nam.PCKG));
    ssm.ext =nam.EXT {strmatchi('SSM',nam.PCKG)};
    
    ssm.NPER=size(PERparvals,1);
    ssm.NROW=dis.NROW;
    ssm.NCOL=dis.NCOL;
    ssm.NLAY=dis.NLAY;
    ssm.NCOMP=btn.NCOMP;
    ssm.MCOMP=btn.MCOMP;  % don't know if needed here ## check 

    %D1 FWEL FDRN FRCH FEVT FRIV FGHB FCHD FMNW (FNEW(n), n=1:3) 10I2
    % option flags (T F) for WEL, DRN etc packages
    % FNEW (n=1:4) future option flags
    %3 FLAGS for packages used by FLOW process
    %conversion from  0 1 to 'F' 'T' is done in writeSSM
    ssm.FWEL  =~isempty(strmatch('WEL',nam.PCKG));
    ssm.FDRN  =~isempty(strmatch('DRN',nam.PCKG));
    ssm.FRCH  =~isempty(strmatch('RCH',nam.PCKG));
    ssm.FEVT  =~isempty(strmatch('EVT',nam.PCKG));
    ssm.FRIV  =~isempty(strmatch('RIV',nam.PCKG));
    ssm.FGHB  =~isempty(strmatch('GHB',nam.PCKG));
    ssm.FCHD  =~isempty(strmatch('CHD',nam.PCKG));  % added by Theo
    ssm.FMNW  =~isempty(strmatch('MNW1',nam.PCKG)); % added by Ruben
    ssm.FNW2  =0;
    ssm.FNW3  =0;
    ssm.FNW4  =0;   
     
    if ssm.FRCH % recharge concentration has to be defined
        ssm.INCRCH =PERparvals(:,strmatchi('INCRCH',PERparnams));
        if any(ssm.INCRCH==0), % Then we need to define CRCH in mf_adapt
            try
                ssm.CRCH=CRCH;  % has CRCH been defined in mf_adapt?
            catch ME
               error('If any IN%s in the PER worksheet = 0 then you must define CRCH{NLAY,NCOMP} in mf_adapt','CRCH');
            end
            for iComp=1:ssm.NCOMP
                crch=PERparvals(:,strmatchi(sprintf('CRCH_%d',iComp),PERparnams));
                for iPer=1:ssm.NPER
                  if ssm.INCRCH(iPer)>0, ssm.CRCH{iPer,iComp}=crch(iPer); end
                end
            end
        else
            ssm.CRCH=cell(ssm.NPER,ssm.NCOMP);
            for iComp=1:ssm.NCOMP
                crch=PERparvals(:,strmatchi(sprintf('CRCH_%d',iComp),PERparnams));
                for iPer=1:ssm.NPER
                  ssm.CRCH{iPer,iComp}=crch(iPer);
                end
            end
        end
    end
    
    if ssm.FEVT
        ssm.INCEVT =PERparvals(:,strmatchi('INCEVT',PERparnams));
        if any(ssm.INCEVT)==0
            try
                ssm.CEVT=CEVT;  % CRCH must be defined in mf_adapt 
            catch ME
               error('If any IN%s in the PER worksheet = 0 then you must define %s in mf_adapt !','CEVT');
            end
            for iComp=1:ssm.NCOMP
              cevt=PERparvals(:,strmatchi(strcmp('CEVT_%d',iComp),PERparnams));
              for iPer=1:ssm.NPER
                  if ssm.INCEVT(iPer)>0, ssm.CEVT{iPer,iComp}=cevt(iPer); end
              end
            end
        else
            ssm.CEVT=cell(ssm.NPER,ssm.NCOMP);
            for iComp=1:ssm.NCOMP
                cevt=PERparvals(:,strmatchi(sprintf('CEVT_%d',iComp),PERparnams));
                for iPER=1:ssm.NPER
                    ssm.CEVT{iPer,iComp}=cevt(iPer);
                end
            end
        end
    end

    %if ssm.FMNW                       WEGHALEN
    %    ssm.mnw.parvals=mnw.parvals;  WEGHALEN
    %end                       WEGHALEN
        
%% point sources, are specified in spreadsheet or in mf_adapt or both
    fprintf('Getting point sourses of different types, see MT3MS manual, D8, p122\n');
    [dummy,ssm.PNTSRC]=getExcelData(XLSF,'PNTSRC','Horizontal');
    
    if exist('PNTSRC','var') && ~isempty(PNTSRC)  % have PNTSRC been additionally defined in mf_adapt?
        ssm.PNTSRC=[ssm.PNTSRC; PNTSRC];  % if so, add them to those from spreadsheet
    end

    % count total number of sources and sinks in any one period include in
    % the flow model
    ssm.MXSS=sum(IBOUND(:)<0 | ICBUND(:)<0);  % maximum number of fixed head or conc cells
    mxwel=0; mxdrn=0; mxriv=0; mxghb=0; mxchd=0; mxmnw=0;
    for iPer=1:ssm.NPER
        if ssm.FWEL, mxwel=max(mxwel,sum(wel.parvals(:,1)==iPer)); end
        if ssm.FDRN, mxdrn=max(mxdrn,sum(drn.parvals(:,1)==iPer)); end
        if ssm.FRIV, mxriv=max(mxriv,sum(riv.parvals(:,1)==iPer)); end
        if ssm.FGHB, mxghb=max(mxghb,sum(ghb.parvals(:,1)==iPer)); end
        if ssm.FCHD, mxchd=max(mxchd,sum(chd.parvals(:,1)==iPer)); end
        if ssm.FMNW, mxmnw=max(mxmnw,sum(mnw.parvals(:,1)==iPer)); end
    end
    mxpntsrc=0;
    if isfield(ssm,'PNTSRC') && ~isempty(ssm.PNTSRC)
        for iPer=1:ssm.NPER
            mxpntsrc=max(mxpntsrc,sum(ssm.PNTSRC(:,1)==iPer));
        end
    end
    ssm.MXSS=ssm.MXSS+mxwel+mxdrn+mxriv+mxghb+mxchd+mxmnw+mxpntsrc;
    
    writeSSM(basename,ssm);
end


%% ===== THE GCG-file (Generalize conjugate gradient solver package)
gcg.SCEN=nam.SCEN(strmatch('GCG',nam.PCKG));

if gcg.SCEN
    fprintf('Generating Generalized Conjugate Gradient Solver Process struct\n');
    gcg.unit=nam.UNIT(strmatch('GCG',nam.PCKG));
    gcg.ext =nam.EXT {strmatch('GCG',nam.PCKG)};
    %3
    gcg.MXITER=MT3Dparvals(strmatchi('MXITER',MT3Dparnams),gcg.SCEN);
    gcg.ITER1 =MT3Dparvals(strmatchi('ITER1' ,MT3Dparnams),gcg.SCEN);
    gcg.ISOLVE=MT3Dparvals(strmatchi('ISOLVE',MT3Dparnams),gcg.SCEN);
    gcg.NCRS  =MT3Dparvals(strmatchi('NCRS'  ,MT3Dparnams),gcg.SCEN);
    gcg.ACCL  =MT3Dparvals(strmatchi('ACCL'  ,MT3Dparnams),gcg.SCEN);
    gcg.CCLOSE=MT3Dparvals(strmatchi('CCLOSE',MT3Dparnams),gcg.SCEN);
    gcg.IPRGCG=MT3Dparvals(strmatchi('IPRGCG',MT3Dparnams),gcg.SCEN);

    writeGCG(basename,gcg);
end


%% ===== THE rct-file (chemical reaction package) =======
rct.SCEN=nam.SCEN(strmatch('RCT',nam.PCKG));

if rct.SCEN
    rct.NCOMP=btn.NCOMP;
    rct.MCOMP=btn.MCOMP;  % Don't know if this is needed in rct ## check
    rct.NCOL=dis.NCOL;
    rct.NROW=dis.NROW;
    rct.NLAY=dis.NLAY;
        
    fprintf('Generating Chemical Reaction Process struct\n');
    rct.unit=nam.UNIT(strmatch('RCT',nam.PCKG));
    rct.ext =nam.EXT {strmatch('RCT',nam.PCKG)};
    %E1 ISOTHM IREACT IRCTOP IGETSC
    %  ISOTHM sorption type
    %    0=no sorption
    %    1=linear sorption
    %    2=freundlich isotherm
    %    3=langmuir
    %    4=first order kinetic sorption (non-equilibrium)
    %    5=dual domain mass transfer (without sorption)
    %    6=dual domain mass transfer (with sorption)
    rct.ISOTHM=MT3Dparvals(strmatchi('ISOTHM',MT3Dparnams),rct.SCEN);

    % IREACT type of kinetic rate reaction flag
    %    0=no kinetic rate reaction
    %    1=first-order kinetic rate reaction
    %      chemical reaction simulations need add-on package
    rct.IREACT=MT3Dparvals(strmatchi('IREACT',MT3Dparnams),rct.SCEN);
  
    % IRCTOP reaction variable entry method flag
    %  >=2 all variables are entered as 3D array (using RARRAY)
    %  < 2 all variables are entered as 1D with one value per layer
    rct.IRCTOP=MT3Dparvals(strmatchi('IRCTOP',MT3Dparnams),rct.SCEN);
    if rct.IRCTOP<2
        error('IRCTOP = %d, this is deprecated and will produce no output! Set IRCTOP=2 in MT3D worksheet!',...
            rct.IRCTOP);
    end
    
    % IGETSC initial conc for the adsorbed phase reading flag (ISOTHM=4,5,6)
     rct.IGETSC=MT3Dparvals(strmatchi('IGETSC',MT3Dparnams),rct.SCEN);
     
     rct.NLAY=size(LAYparvals,1);
     
     % generate cell arrays to store data for layers and 
     rct.RHOB   =cell(rct.NLAY,1);
     rct.PRSITY2=cell(rct.NLAY,1);
     rct.SRCONC =cell(rct.NLAY,rct.NCOMP);
     rct.SP1    =cell(rct.NLAY,rct.NCOMP);
     rct.SP2    =cell(rct.NLAY,rct.NCOMP);
     rct.RC1    =cell(rct.NLAY,rct.NCOMP);
     rct.RC2    =cell(rct.NLAY,rct.NCOMP);
     
     % in which columns of worksheet, note that these may be multiple
     % columns because we use more than one species
     J_RHO=strmatchi('RHOB'   ,LAYparnams);
     J_PRS=strmatchi('PRSITY2',LAYparnams);
     J_SRC=strmatchi('SRCONC' ,LAYparnams);
     J_SP1=strmatchi('SP1'    ,LAYparnams);
     J_SP2=strmatchi('SP2'    ,LAYparnams);
     J_RC1=strmatchi('RC1'    ,LAYparnams);
     J_RC2=strmatchi('RC2'    ,LAYparnams);
     
     % get those multiple columns, but only those for our species
     rhob   =LAYparvals(:,J_RHO);
     prsity2=LAYparvals(:,J_PRS);
     srconc =LAYparvals(:,J_SRC(1:rct.NCOMP));
     sp1    =LAYparvals(:,J_SP1(1:rct.NCOMP));
     sp2    =LAYparvals(:,J_SP2(1:rct.NCOMP));
     rc1    =LAYparvals(:,J_RC1(1:rct.NCOMP));
     rc2    =LAYparvals(:,J_RC2(1:rct.NCOMP));
          
     % if any of these parameters in worksheet < 0 a cell array
     % must be defined in mf_adapt, so replace cell aray with one from mf_adapt
     if any(rhob   <0), rct.RHOB   =RHOB;    end  % cell array {NLAY,1}
     if any(prsity2<0), rct.PRSITY2=PRSITY2; end  % cell array {NLAY,1}
     if any(srconc <0), rct.SRCONC =SRCONC;  end  % cell array {NLAY,NCOMP}
     if any(sp1    <0), rct.SP1    =SP1;     end  % cell array {NLAY,NCOMP}
     if any(sp2    <0), rct.SP2    =SP2;     end  % cell array {NLAY,NCOMP}
     if any(rc1    <0), rct.RC1    =RC1;     end  % cell array {NLAY,NCOMP}
     if any(rc2    <0), rct.RC2    =RC2;     end  % cell array {NLAY,NCOMP}
     
     
     % if value in spreadsheet>0 use that instead of the one defined in
     % mf_adapt
     for iLay=1:rct.NLAY,
         if rhob(   iLay)>=0, rct.RHOB{   iLay}=rhob(   iLay); end
         if prsity2(iLay)>=0, rct.PRSITY2{iLay}=prsity2(iLay); end
         for iComp=1:rct.NCOMP
             if srconc( iLay,iComp)>=0, rct.SRCONC{iLay,iComp}=srconc(iLay,iComp); end
             if sp1(    iLay,iComp)>=0, rct.SP1{   iLay,iComp}=sp1(   iLay,iComp); end
             if sp2(    iLay,iComp)>=0, rct.SP2{   iLay,iComp}=sp2(   iLay,iComp); end
             if rc1(    iLay,iComp)>=0, rct.RC1{   iLay,iComp}=rc1(   iLay,iComp); end
             if rc2(    iLay,iComp)>=0, rct.RC2{   iLay,iComp}=rc2(   iLay,iComp); end
         end
     end
        
    writeRCT(basename,rct);
end

% ===== END of MT3DSM input ==============================

%% ======================== SWI input ============================
swi.SCEN=nam.SCEN(strmatch('SWI',nam.PCKG));

if swi.SCEN
    fprintf('Generating Salt Water Intrusion struct\n');
    swi.unit=nam.UNIT(strmatch('SWI',nam.PCKG));
    swi.ext =nam.EXT {strmatch('SWI',nam.PCKG)};
    
    swi.NCOL     =dis.NCOL;
    swi.NROW     =dis.NROW;
    swi.NLAY     =dis.NLAY;
    swi.LAYCBD   = dis.LAYCBD;
    
    swi.ISWIZT   = MFLOWparvals(strmatchi('ISWIZT'  , MFLOWparnams),swi.SCEN);
    swi.ISTRAT   = MFLOWparvals(strmatchi('ISTRAT'  , MFLOWparnams),swi.SCEN);
    swi.NPRN     = MFLOWparvals(strmatchi('NPRN'    , MFLOWparnams),swi.SCEN);
    swi.TOESLOPE = MFLOWparvals(strmatchi('TOESLOPE', MFLOWparnams),swi.SCEN);
    swi.TIPSLOPE = MFLOWparvals(strmatchi('TIPSLOPE', MFLOWparnams),swi.SCEN);
    swi.ZETAMIN  = MFLOWparvals(strmatchi('ZETAMIN' , MFLOWparnams),swi.SCEN);
    swi.DELZETA  = MFLOWparvals(strmatchi('DELZETA' , MFLOWparnams),swi.SCEN);
    swi.NU       = MFLOWparvals(strmatchi('NU'      , MFLOWparnams),       :);

    swi.Z        =Z;
    swi.ZETA     =ZETA;
    swi.SSZ      =SSZ;     % SWI uses SSZ in place of PEFF
    swi.ISOURCE  =ISOURCE;
    swi.NPLN     =length(ZETA);
    swi.NZONES   =swi.NPLN+1;
    try
        if swi.ISTRAT
            swi.NU     =swi.NU(1:swi.NZONES);
        else
            swi.NU     =swi.NU(1:swi.NZONES+1);
        end
    catch ME
        error('length swi.NU (%d) incompatible with number of planes in ZETA (%d), see manual',...
            length(NU),swi.NZONES);
    end
    
    writeSWI(basename,swi);
end

%% ============ End of input ============

fclose('all');  % just to make sure all files are closed
fprintf('.... finished generation of input files, all files closed! Time= %s!\n\n',datestr(now));

toc;

%% ========= RUN model(s): Seawat MT3DMS or mf2k =====================

if ismac, quote=''; else quote='"'; end


if    (nam.runSWT || ~isempty(vdf.SCEN)), system([quote,SEAWAT,quote,' swt_v4.nam']);
elseif           ~isempty(swi.SCEN), system([quote,SWI,quote,' < SWInamfilename.txt']);
elseif           ~isempty(btn.SCEN), system(['start /BELOWNORMAL ', quote,MODFLOW,quote,' mf2k.nam']);
%   fprintf('Pauzing... press enter to coninue!\n'); pause;
                                     system([quote,MT3D,quote    ,' mt3dms5.nam']);
else                                 system([quote,MODFLOW,quote ,' mf2k.nam']);
end

%% now run mf_analyze

%% Set AFTERMFSETUP to run this file at the end of mf_setup, used for calibration

if exist('AFTERMFSETUP','var')
    eval(AFTERMFSETUP);
end
