%% Saltwater Intrusion Package, example 4

basename='swiex4';

% Example 4, upconing of an interface below a pumpoing well in a
% two-aquifer system

% Consider three-dimensional flow in the same two-aquifer system as in teh
% previous example. The area of interest extends 1300 m in the x-direction,
% and 1000 m in the y-direction. The upper aquifer is bounded on top by the
% ocean floor in the west (Figure 8z); a 300 m wide section of ocean floor
% is modeled. The flow entering the domain inthe east is 0.12 m2/d, so that
% the total flow entering the system is 120 m3/d. A well of discharge Q=70
% m3/d is situated in teh top aquifer at (x,y)=(975,525).
% Initially, the interface betwen the fresh and seawater is straight and is
% at the top of each aquifer at x=250 m and at the bottom of each aquifer
% at x=450 m; i.e. the initial slope is 0.1. The dimensionless density of
% the seawater is nu=0.025.
% Each aquifer is discretized in 20*26 square cells of 50 by 50 m; the well
% s at the center of cell (10,20). The timestep is one year. The top and
% toe tracking parameters are a TOESLOPE and TOPSLOPE of 0.05, a DELZETA of
% 0.25 and a ZETAMIN of 0.025. Example 4 is summarized in table 4; all
% input files are in Appendix D (of the SWI manual, part 1).
% Teh position of the interface after 50 years is to be shown as wel as the
% position after 400 years.

%% 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:50:1300; yGr=0:50:1000; zGr=permute([40 20.01 19.99 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

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

%% WEL
Qtot =120;
Qwell=-70;

[DX,DY,DZ]=meshgrid(DELX,DELY,DELZ);

% inflow from the right implemented as wells
kDB=TRAN.*DY; q=Qtot/sum(sum(sum(kDB(:,end,:))));
Qw=kDB.*q; Qw(:,1:end-1,:)=NaN;

Idx=find(~isnan(Qw));
LRC=cellIndices(Idx,size(IBOUND),'LRC');
u=ones(size(LRC(:,1)));

iPer=1;

WEL=[u*iPer LRC Qw(Idx)];

% the actual well, as a 50x50 (4 cell area) see fig manual
ix=find(xm>970 & xm<980);
iy=find(ym>520 & ym<530);
Q=NaN(size(IBOUND)); Q(iy,ix,1)=1;

Idx=find(~isnan(Q)); Q(Idx)=Qwell/length(Idx);

LRC=cellIndices(Idx,size(IBOUND),'LRC');
u=ones(size(LRC(:,1)));
WEL=[WEL; [u LRC Q(Idx)]];

%% GHB

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

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

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 250    450 xGr(end)];  % x points for interface
zetax{1,2}=[0 250    450 xGr(end)];  % x points for interface
zetay{1,1}= yGr([1 end]);
zetay{1,2}= yGr([1 end]);
zetaz{1,1}=ones(size(zetay{1,1}))*squeeze(zGr([1 1 2 2]))';
zetaz{1,2}=ones(size(zetay{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)=...
            interp2(zetax{iSurf,iLay},...
                    zetay{iSurf,iLay},...
                    zetaz{iSurf,iLay},...
                    xm,ym);  % z points for interface
    end
end

% Actual distribution over the layers is done by writeSWI