%% Saltwater Intrusion Package, example 4

basename='swiex5';

% Exampe 5, shape of brakish zone below a sqaure island
% Consider a square island of 525 by 525 metrers. Teh origin of he
% coordinate system is chosen at the center of the island. Below the island
% is a 50 m thick aquifer with a hydraulic conductivity of 10 m/d. There is
% ann infiltration on the island of 2 mm/s (specified with the recharge
% package). In addition, there is an areas of 100 by 100 meter wehere there
% is a combined net discharge of 80 m3/d (the dark gray area in Figure 9a).
% The island is surrounded by an ocean with an equivalent freshwater level
% of 0.05 m; the vertical resistance to ooutflow in the ocean is c=2 days
% (Vcont=0.5).
% The domain is discretized into 21 by 21 cells of 50 by 50 m. A timestep
% of 50 days is used. The tip and toe tracking parameters are not too
% important (there is only a top where the brackish zone intersects the top
% of the aqufier); top and toeslope are 0.05; delzeta and zetamin are
% computed with equation 16. The ocean is represented by a strop of 150 m
% (6 cells) wiede around the island (the light gray area in Figure 9a).
% Groundwater below the island consists of frshwater (nu=0), brdakish water
% (nu=0.0125), and seawater (nu-0.025). Initially the brackish zone below
% the island is 10 m thick and extends from 25 to 35 m below the surface;
% below the strip respresentind the ocean, the surfaces bounding the
% brackish zone are at the top of the aquifer. Example 5 is summarized in
% Appendix E.

%% Specifics of this example
K     =10;    % hydraulic conductivity
peff  =0.2;   % effective porosity
hOcean=0.05;  % ocean head specified
cOcean=2;     % resistance ocean bottom
IOCEAN    =1; % zone number ocean floor
IDISCHARGE=2; % zone number discharging square

%% Specify mesh

xGr=-525:50:525; yGr=-525:50:525; zGr=permute([0 -50],[1,3,2]);

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

zone=NaN(NCOL,NROW);
zone(:,:)            = IOCEAN;
zone(4:end-3,4:end-3)=~IOCEAN;
zone(8:9    ,8:9)    = IDISCHARGE;

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

% Data arrays
TRAN =NaN(NROW,NCOL); TRAN(:,:)=K*(zGr(1)-zGr(end));
IBOUND=ones(NROW,NCOL);

%% WEL

Qdisch =-80;

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

Idx   =find(zone==IDISCHARGE);
LRC=cellIndices(Idx,size(IBOUND),'LRC');
Q=Qdisch*DY(Idx).*DX(Idx)/sum(DY(Idx).*DX(Idx));
u=ones(size(Q(:)));

iPer=1;
WEL    =[u*iPer LRC Q(:)];
%% GHB

Idx   =find(zone==IOCEAN);
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));     % source zone incase of sinks 

%% ZETA planes

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

ZETA(1).values(zone~=IOCEAN)=-25;
ZETA(2).values(zone~=IOCEAN)=-35;

% Finisch