% MF_ADAPT -- Henry problem, example comes with Seawat
%  TO 090101

basename='Henry';

%% The Henry problem (See Guo and Langevin (2002) p70:
% Henry (1964) presented an analytial solution for groundwater flowing
% toward a seawater boundary. Beause an analytical solution was available
% for the Henry problem, many numerical codes have been evaluated and
% tested with the Henry solution ...
% The basic Henry problem concerns va vertical cross section of 1 m high
% and 2 m wide where fresh water enters uniformly at the left boundary with
% a rate of 5.702 m2/d and concentration zero. A constant hydrostatic head
% boundary is applied at the right hand (seawater) side. The upper and
% lower boundaries are no-flow.

% Data
peff=0.35;   % effective porosity
k   =864;    % m/s conductivity (uniform)
rho1=1025;   % kg/m3 (see SEAWAT worksheet in workbook
rho0=1000;   % kg/m3 (fresh water)
QL  =5.702;  % m2/d or m3/d per m perpendicular to the model
Cfresh = 0;  % kg/m3 initial and fresh water salinity
Csea   =35;  % kg/m3 Sea water salinity
aL  =0;      % longitunidal dispersivity
aT  =0;      % transverse dispersivity
Dm1 =1.62925;% Molecular diffusion case 1
Dm2 =0.57024;% Molecular diffusion case 1

strthd =1.0; % initial heads (point water heads)

%% Generate a Henry benchmark problem (Guo and Langevin 2002, p70ff)

yGr= 0.0: 1   :1.00;
xGr= [0.0: 0.1:2.00,  2.01];
zGr= 1.0 :-0.1:0.00;

[xGr,yGr,zGr,xm,ym,zm,DX,DY,DZ,NX,NY,NZ]=modelsize3(xGr,yGr,zGr);

dz=reshape(DZ,[1,1,NZ]);

[XM ,YM ,ZM]=meshgrid(xm ,ym, zm ); % full 3D coordinate arrays
[XGR,YGR,Z ]=meshgrid(xGr,yGr,zGr); % full 3D coordiante arrays

%% Generate all other matrices

IBOUND=ones(NY,NX,NZ);  % we use CHD package don't need IBOUND==-1 at all
ICBUND=ones(NY,NX,NZ);  % fixed conc boundaries almost never necessary
STRTHD=ones(NY,NX,NZ)   *strthd;
STCONC=ones(NY,NX,NZ)   *Csea;

%SF1   =Ss*dz(ones(1:NY),ones(1,NX),1:NZ);
% Contrary to Langevin, we use steady state (he used SF1=0)

TRAN  =k *dz(ones(1:NY),ones(1,NX),1:NZ);
VCONT =k./dz(ones(1:NY),ones(1,NX),1:NZ);
VCONT=0.5*(VCONT(:,:,1:end-1)+VCONT(:,:,2:end));
PEFF  =ones(NY,NX,NZ)   *peff;

%% Specify the prescribed left-side inflow using the well package

IL=find(XM<xGr(2));  % global indices of left boundary cells
LRCleft=cellIndices(IL,size(IBOUND),'LRC');  % list of L R C indices of this boundary

IR=find(XM>xGr(end-1));  % global indices of left boundary cells
LRCright=cellIndices(IR,size(IBOUND),'LRC');  % list of L R C indices of this boundary

u=ones(size(IL));  % shorthand, is column IL long of ones

WEL=[u LRCleft u*QL/NZ];      % all wells stress period 1

%% Specify the right-side constant head using CHD and seawat option
% for hydrostaticy.

hSea=1.0;  % sea level

% In Seawat V4 CHD boundary options can be used to specify more accurately
% the head boundary when density changes during the simulation.
% See Langevin e.a. (2008), p15 (Seawat V4 manual)
% CHDDENSOPT = 0 -- use density in cell to compute ref density head
%              1 -- use CHDDENS to compute ref dens head
%              2 -- env head, use average density above point
%              3 -- use ref head h0, no converstion done
% if any line has option 1, CHDDENS must be provided for all lines
% the way its done below is quite transparent

% using 0 is in fact the default
CHDDENSOPT=0;  % used if parameter exists in MATLAB
CHDDENS=1025;  % used if parameter exists in MATLAB

CHD=[u LRCright u*[hSea hSea CHDDENSOPT CHDDENS] ];  % all CHD period 1

% if you desire not to use any CHDDENSOPT, then you must out-comment
% both parameters CHDDENSOPT and CHDDENS and use the line below
% CHD=[u LRCright u*[hSea hSea] ];  % ehead and shead are zero  

%% Specify point sources at right hand boundary for SSM package

LTYPE=1; % constant head
PNTSRC=[u LRCright u*[Csea LTYPE Csea] ];  % all pointsources period 1

% note that we can also specify the WELLS and CHD in the workbook sheets
% WEL and CHD respectively. But the way shown here is more powerfull
