% MF_ADAPT -- Benchmark Problems and Application Example
% One-dimensional transport in a uniform flow field
% Zheng & Wang (1999) -- p130
%  TO 091106 091201


basename='1D-Uniform';

%% Problem

% The problem considered here is transport due to uniform through a linear
% homogeneous model subject to dispersion, linear sorption and decay.
% Model has cells 10 m wide and is 101 cells long with fixed head at its
% right hand boundary and prescribed flow at the left end of the model.
% As the orginal problem only states the seepage velocity, prescribed flow
% at the left hand side of the model is the natural boundary.
% For this problem, conductivity is immaterial, any value is ok.
% We have the following cases
%
% aL= 0 m, R=0, lambda (decay) = 0
% aL=10 m, R=0, lambda (decay) = 0
% aL=10 m, R=5, lambda (decay) = 0
% aL=10 m, R=5, lambda (decay) = 0
%
% We can simulate the four case in a single model run by means of 4 layers
% that are separated by making sure that vertical dispersion is zero for
% each layer.
% Linear sorption is switched on with by setting ISOTHM=1 in worksheet
% MT3D. Decay is switched on by setting IREACT=1 in the same worksheet.
% The parameters to be set in the LAY worksheet are
%
%DSP    DSP    DSP    DSP      RCT   RCT    RCT    RCT    RCT
% AL    TRPT   TRPV   DMCOEF_1 RHOB  SP1_1  SP2_1  RC1_1  RC2_1
%  0     0.1      0      0        1      0      0      0      0
% 10     0.1      0      0        1      0      0      0      0
% 10     0.1      0      0        1      1      0      0      0
% 10     0.1      0      0        1      1      0  0.002      0
%
%
% The other parameters are defined below

%% Conductivities and effective porosity

kh= 10; kv=0;   % m/d
peff=0.25;        % []


%% Mesh

%% Get the actual number of layers from the workbook
[LAYparnams,LAYparvals]=getExcelData([basename '.xls',],'LAY','Horizontal');

%% Use this number to further assemble the model grid

NLAY=size(LAYparvals,1); NROW=1; NCOL=101;
DELX=ones(1,NCOL) * 10;
DELY=1;
DELZ=ones(NLAY,1) * 1;
z=[0; -cumsum(DELZ(:))];

xGr=[-DELX(1)/2 cumsum(DELX)];   % center of first cell at xGr=0
yGr=[0 cumsum(DELY)];
Z=NaN(NROW,NCOL,NLAY+1);

for i=1:size(Z,3), Z(:,:,i)=z(i); end

%% and to generate all other matrices

IBOUND=ones(NROW,NCOL,NLAY);  IBOUND(:,end,:)= -1;
ICBUND=ones(NROW,NCOL,NLAY);  %ICBUND(:, 1 ,:)= -1; % use CHD package
STRTHD=zeros(NROW,NCOL,NLAY);

STCONC{1}=ones(NROW,NCOL,NLAY)*0;

HK    =ones(NROW,NCOL,NLAY)   *kh;
VK    =ones(NROW,NCOL,NLAY)   *kv;
PEFF  =ones(NROW,NCOL,NLAY)   *peff;

%% East and west boundary cell index triples + period number (only 1 period)

%       IPER         LAY        ROW            COL
WEST=[ones(NLAY,1) (1:NLAY)' ones(NLAY,1) ones(NLAY,1)     ];  % iPer Lay Row Col
EAST=[ones(NLAY,1) (1:NLAY)' ones(NLAY,1) ones(NLAY,1)*NCOL];  % iPer Lay Row Col

%% Wells at tehe west boundaries

vx=0.24;                 % m/d,  % given seepage velocity = vx/peff=

q=vx.*DELZ.*DELY*peff;  % Dary flux or specific discharge

%     iLRC Q
WEL=[ WEST q];

%% PNTSRC (point sources for the transport model (SSM module)

%                            CSS ITYPE CSSM(1..4)
%PNTSRC=[WEST ones(NLAY,1)*[0     2     1 1 1 1] ];
PNTSRC=[WEST  ones(NLAY,1)*[1     2     1] ];

% That's all, now run mf_setup and then mf_analyze