%% Saltwater Intrusion Packge, example 1
% See 

% Rotating interface in a vertical cross section 40 m high and 250 m long
% and 2 m wide. The origin of the x,z coordinate systems is in the upper
% left corner. Flow in confined at all times and the top and bottom of the
% aquifer are impermeabile and hoirzontal. There is a constant flow Qxo
% of 1 m2/d from left to right, gradient of .0125. The hydraulic
% conductivity k= 2 m/d and the effective porosity is 0.2. Aquifer and
% water are considered incompressible. The interface is initially at a
% 45 degrees angle from (x,z)=(80,0) to (x,z)=(120,-40).
% The density of fresh water is 1000 kg/m3 and of saltwater 1025 kg/m3. 
% The domain is discretized into 50 cells of 5 m wide, and the timestep 
% is 2 d. A source of 2 m3/d is specified in cell 1. Modflow requires the
% specificaiton of 1 head (the rotation of the interface is independent
% of this choice); here, the head is fixed to 0.05 m at the top of the last
% cell.
% There is a freshwater zone and a saltwater zone and thus one surface
% between them (NPNL=1). Flow is treated as stratified (ISTRAT=1) and the
% elevations of the interface is writeen to a file called swieex1.zta on
% unit number ISWIZT=52 every NPRN=100 times. Thet total number of time
% steps is 200, and thus the position of the interface is printed at 200
% and 400 days. The stope of th einterface varies from 1 (initial position)
% to about 0.2 (when it reaches the corner of the domain). The maximum slop
% of the toe an dtop is chosen to be 0.2 (TOESLOPE=TOPSLOPE=02), and the
% deltaGamma and deltaGammamin parameters are computed arccoding to
% equation 16 (ZETAMIN=0.01, DELZETA=0.1). Since this is stratified flow
% with 2 zones, 2 v values must be specified (0 and 0.025). The source/sink
% terms are specified to be fresh water (zone 1) everywhere (ISOURCE=1);
% the only source in this model is in cell 1.

basename='swiex1';

%% Specifics of this example
k     =2;    % hydraulic conductivity
peff  =0.2;  % effective porosity
strthd=0.05; % specified head in cell at right
Qo    =2;    % source in cell 1

%% Specify mesh
xGr=0:5:250; yGr=[0 2]; zGr=[0 -40];

%% Mesh housekeeping
[xGr,yGr,zGr,xm,ym,zm,DELX,DELY,DELZ,NCOL,NROW,NLAY]=modelsize3(xGr,yGr,zGr);

Z=zGr(ones(NROW,1),ones(1,NCOL),1:NLAY+1);

%% Matrices of computation grid
IBOUND=ones(NROW,NCOL,NLAY);
HK    =ones(size(IBOUND))*k;
VK    =HK;

%% Head boundary conditions
IBOUND(:,end,1)=-1;                % fixed head at right hand side
STRTHD=ones(size(IBOUND))*strthd;  % start using fixed head

%% Get the wells (left hand side of model) -- here done in a general way

[XM,YM,ZM]=meshgrid(xm,ym,zm);
LRC=cellIndices(find(XM(:,xm<xGr(2),:)),size(XM),'LRC');
u=ones(size(LRC(:,1)));

NPER=1; WEL=[];
for iPer=1:NPER
    WEL=[WEL; [iPer*u LRC Qo*u]];
end
% You might also do this as simple as
%WEL=[1    1 1 1    Qo];  % iper L R C Qo

%% SWI specific matrices
SSZ=ones(size(IBOUND))*peff;     % effective porosity called SSZ in SWI
ISOURCE=zeros(size(IBOUND));     % source zone incase of sinks 

%% ZETA planes
NSURF=1;  % We have only one surface here, so NSURF=1 in this case

% A plane with elevations is necessary for all interfaces
% We just specify this single interface as simply as possible using
% coordinates in space (x,z)
x=[0 80 120 250];  % x points for interface
z=[0  0 -40 -40];  % z points for interface

% ZETA must be specified for every interface in all layers, but we have
% only one layer in this example. We use the above specified surface and
% interpolate its elevations to the cell centers of all (1) layers
ZETA(NSURF).values=interp1(x,z,xm); % zeta all cells

% Actual distribution over the layers is done by writeSWI