%% Saltwater Intrusion Package, example 4

basename='swiex3';

% Example 3,
% This example considers interface flow in a two-aquifer system. The plane
% of flow is 4000 m long and 2 m wide; both aquifers are 20 m thick and are
% separated by a leaky layers of 1 m thickness. The origin of the x,z
% coordinate system is chosen at the lower left-hand corner of the plane.
% The hydraulic conductivity of the top aquifer is 2 m/d and of the bottom
% aquifer 4 m/d; the effective porosity of both aquifers is 0.2. Flow is
% semi-confined at all times. The first 600 m of the top of aquifer 1
% represents the ocean floor; the freshwater head at the ocean bottom is
% 50 m. The leaky layer between aqfuiers 1 and 2 has a vertical resistance
% of flow of c = 100 days (Vcont=0.01). The vertical resistance to outflow
% in the ocean is c=50 days (Vcont=0.02). There is a compehensiver flow of
% 0.015 m2/d towards the ocean (eqaully distributed over both aqufiers at
% the buondary).
% Initially the interface between the frsh and sawater is straight and is
% at the top of each aquifer at x=500 m and at the bottom of each aquifer
% at x=900 m; i.e. the initial slope is 0.05. The dimensinless density
% nu=0.025.
% Each aquifer is discretized into 200 cells of 20 m long and the timestep
% is 1 year. Teh last cell in the top aquifer has a source of 0.01 m3/d,
% and the last cell in the bottom aquifer a source of 0.02 m3/d. The ocan
% bottom is modeled with the general head boundary (GHB) package. The head
% in the GHB cells is specified to teh quivalent frsh water head at the
% bottom of the ocean (50 m). The Cond value is computed as the Vont value
% (0.02) of the ocean bottom by the area of cell (40). The water in the GHB
% cells representing the ocean is salt. Hence, when the GHB cell acts as a
% source, the infiltrating water is zone 2. Whne the GHB acts as a sink,
% water that flows out of the cells is of the same type as the water at the
% top of the aquifer. Hence, the ISOURCE parameter is set to -2. The top
% and toe tracking parameters are a TOESLOPE of 0.02, a TIPSLOPE of 0.04, a
% DELZETA of 0.06 and a ZETAMIN of 0.006.

%% Specifics of this example
K     =[2 4];% hydraulic conductivity
C     =100;  % resistane between first and second aquifer
peff  =0.2;  % effective porosity

hOcean=50;   % ocean head specified
cOcean=50;   % resistance ocean bottom
%% Specify mesh

xGr=0:20:4000; yGr=[0 2]; zGr=permute([41 21 20 0],[1,3,2]);

%% Mesh housekeeping
[xGr,yGr,xm,ym,DELX,DELY,NCOL,NROW]=modelsize(xGr,yGr);

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

[LAYnams,LAYvals,LAYCBD,NLAY]=getLayers([basename,'.xls'],zGr);

DELZ=NaN(1,1,NLAY);
zm  =NaN(1,1,NLAY);

TRAN =NaN(NROW,NCOL,NLAY);
VCONT=NaN(NROW,NCOL,NLAY-1);
k=1;
for iLay=1:NLAY
    TRAN(:,:,iLay)=K(iLay)*(zGr(k)-zGr(k+1));
    DELZ(iLay)=zGr(k)-zGr(k+1);
    zm  (iLay)=0.5*(zGr(k)+zGr(k+1));
    if LAYCBD(iLay)
        VCONT(:,:,iLay)=1./C(iLay);
        k=k+1;
    end
    k=k+1;
end

[DX,DY,DZ]=meshgrid(DELX,DELY,DELZ);
[XM,YM,ZM]=meshgrid(xm  ,ym  ,zm  );

%% Matrices of computation grid
IBOUND=ones(NROW,NCOL,NLAY);

%% wel
Idx   =find(XM>xGr(end-1));
LRC=cellIndices(Idx,size(IBOUND),'LRC');

Qin=[0.01; 0.02];

u=ones(size(Idx));

iPer=1;
WEL    =[u*iPer LRC Qin];
%% GHB

[XM,YM,ZM]=meshgrid(xm,ym,zm);

Idx   =find(XM<600 & ZM==zm(1));
LRC=cellIndices(Idx,size(IBOUND),'LRC');
vcont=DY(Idx).*DX(Idx)/cOcean;
u=ones(size(vcont(:)));

iPer=1;
GHB    =[u*iPer LRC u*hOcean vcont(:)];

%% Head boundary conditions
STRTHD=ones(size(IBOUND))*hOcean;  % start using fixed head

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

%% ZETA planes

% ZETA(iSurf).values contains the elevations of surf 3 in all layers
% size(ZETA(iSurf).values is therefore (NROW,NCOL,NLAY)

zetax{1,1}=[0 500    900 xGr(end)];  % x points for interface
zetax{1,2}=[0 500    900 xGr(end)];  % x points for interface
zetaz{1,1}=squeeze(zGr([1 1 2 2]))';
zetaz{1,2}=squeeze(zGr([3 3 4 4]))';

NSURF=1;  % We have only one surface here
ZETA(NSURF).values=NaN(NROW,NCOL,NLAY);  % allocated struct

for iSurf=1:NSURF
    for iLay=1:NLAY
        ZETA(iSurf).values(:,:,iLay)=...
            interp1(zetax{iSurf,iLay},...
                    zetaz{iSurf,iLay},...
                    xm);  % z points for interface
    end
end

% Actual distribution over the layers is done by writeSWI