% MF_ADAPT --  mt3dms injeciton i n uniform flow field, example
% Zie Zheng & Wang (1999) p138ff
%
%  TO 091112

basename='2D-Uniform';

% The problem deals with an injection well in a uniform flow field whose
% direction is at 45 degrees in NE direction. The analytical solution was
% given by Wilson and Miller (1978). The analytical solution is applicable
% under the condtion that the aquifer is of infinite exent and relatively
% thin in vertical direction, so that instantaneous vertical mixing can be
% assumed and the injection rate is insignificant with respect to the
% ambient uniform flow.

% The numerical model consists of 31 rows and 46 columns and 1 layer. The
% direction of the uniform ambient flow is along the rows of the model
% (x-direction). Teh model parameters are

%% Model parameters specified by Zheng

delx=10;    % m
dely=10;    % m
delz=10;    % m
peff=0.30;  % effective porosity
aL=10;      % m set in worsheet LAY
            % aT/aL = 0.1; set in worksheet LAY
vx=1/3;     % m/d groundwater seepage velocity = qx/peff
qx=vx*peff; % m/d Darcy flux i.e. specific discharge in x-direction
Q=1;        % m3/d injection well flow
cInj=1000;  % ppm concentration in injection water
            %simulation time 365 days

kh= 10; kv=0;  % values are irrelevant as the flow field is determined entirely
% by the given uniform seepage velocity and the injection flow
% however we will use a kh to create the desired uniform flow
% and, by the way, modflow needs it

%% Mesh definition

NLAY=1; NROW=31; NCOL=46;
DELX=ones(1,NCOL)* delx; xGr=            [0 cumsum(DELX)];
DELY=ones(NROW,1)* dely; yGr=sum(DELY) - [0;cumsum(DELY(:))];
DELZ=ones(NLAY,1)* delz; z= - [0; cumsum(DELZ(:))];

[xGr,yGr,z,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,z);

Z=NaN(NROW,NCOL,NLAY+1); for iLay=1:size(Z,3), Z(:,:,iLay)=z(iLay); end

%% Generate all other matrices
IBOUND=ones(NROW,NCOL,NLAY);  IBOUND(:,end,:)= -1;  % fix right-hand head
ICBUND=ones(NROW,NCOL,NLAY);   % no fixed concentration points needed
STRTHD=zeros(NROW,NCOL,NLAY);  % starthead can be anything if we use a steady state run

STCONC{1}=ones(NROW,NCOL,NLAY)*0;  % Initial concentration of species 1
STCONC{2}=ones(NROW,NCOL,NLAY)*0;  % Initial concentration of species 2

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

%% East and west boundary cells

% To get the a list of Layer Row Col triples for the cells of the west
% boundary of the model, we first generate a matrix of NaNs of the size of
% our model
G=NaN(NROW,NCOL,NLAY);
G(:,1,:)=0;             % Then mark the left boundary
I=find(G==0);           % Then get the indices of these cells
WEST=cellIndices(I,size(G),'LRC'); % Then use CellIndices to get the list for you
  % type help cellIndices to see how it works. Quite useful when dealing
  % with finite difference grids.
  
%% Wells at tehe west boundaries

q=vx*DELZ.*DELY*peff;   % compute the inflow of the west boundary cells

% then put the stress period number in front of the WEST list and the q
% behind it to obtain the list of wells in the form that MODFLOW wants
% them.

%     iLRC Q
WEL=[ones(size(q)) WEST q];  % stress peirod 1

%Finally add the injection well at LRC=[1 16 16] to the well list
WEL=[WEL; [1 1 16 16 Q]];
% ---------^-- stress period, we have only one in this example


%%                            CSS ITYPE CSSM(1..4)

% The inflow form the west size has concentration zero. Hence we need not
% specify point source concentration for it. But we do have to specify th
% concentration of the inflow of the well at [1 16 16]

itype=2; % ITYPE is 2 for wells

%       per L R  C CSS  ITYPE CSS(dummy if NCOMP=1)
PNTSRC=[1  1 16 16 cInj itype cInj]; % well at 1 16 16, itype 2
