% Seawat V4, Benchmark example, See Langevin e.a. 2008, p23ff

%  TO 091106

% Example problem (verbatim from Langeving e.a. 2008, p23ff)
% An example problem was developed to illustrate the new model features
% presented in this report (see above). The problem consists of a
% two-dimensional cross section of a confined coastal aquifer initially
% saturated with relatively cold seawater at a temperature of 5C. Warmer
% freshwater with a temperature of 25C is injected into the coastal aquifer
% along the left boundary to represent flow from inland areas (fig. 1). The
% warmer freshwater flows to the right, where it discharges into a vertical
% ocean boundary. The ocean boundary is represented with hydrostatic
% conditions base on a fluid density calculated from seawater salinities at
% 5C. No-flow conditions are assigned to the top and bottom boundaries. A
% complete list of the intput values used for the problem is gien below.
% This problem is a simplified representation of what might occur in a
% coastal carbonate aquifer.

basename='coastal_flow';

% The entire table of Langevin et al, has been stored in the excel workbook
[LanNams,LanVals]=getExcelData([basename, '.xls'],'TableLangevin','Vertical');

% Assing values from table for as far as needed for this model
top =LanVals(index('top' ,LanNams));
NCOL=LanVals(index('NCOL',LanNams));
NROW=LanVals(index('NROW',LanNams));
NLAY=LanVals(index('NLAY',LanNams));
Dx  =LanVals(index('Dx'  ,LanNams));
Dy  =LanVals(index('Dy'  ,LanNams));
Dz  =LanVals(index('Dz'  ,LanNams));
kh  =LanVals(index('kh'  ,LanNams));
kv  =LanVals(index('kv'  ,LanNams));
peff=LanVals(index('peff'  ,LanNams));
Ss  =LanVals(index('Ss'  ,LanNams));
h_ocean=LanVals(index('h_ocean',LanNams));
c_ocean=LanVals(index('c_ocean',LanNams));
T_ocean=LanVals(index('T_ocean',LanNams));
c_inflow=LanVals(index('c_inflow',LanNams));
T_inflow=LanVals(index('T_inflow',LanNams));
Q_inflow=LanVals(index('Qinflow',LanNams));

%% Mesh using table data

[xGr,yGr,zGr,xm,ym,zm,DX,DY,DZ,NX,NY,NZ]=...
    modelsize3([0 Dx*(1:NCOL)],[0 Dy*(1:NROW)],[top, -(1:NLAY)*Dz]);

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

%% Generate all other matrices using table data
M=ones(NROW,NCOL,NLAY);      % dummy array of ones, size of model

IBOUND   =M;  IBOUND(:,end,:)=-1;
ICBUND   =M;  % use CHD package
STRTHD   =M*h_ocean;
STCONC{1}=M*c_ocean;   % intially model fillee with ocean water
STCONC{2}=M*T_ocean; % * initialy model filled with ocean water

HK       =M*kh;
VK       =M*kv;
PEFF     =M*peff;
SS       =M*Ss;

%% Get stress periods from workbook to construct boundary conditions generically
[PERnams,PERvals]=getPeriods([basename, '.xls']);
NPER=size(PERvals,1);

%% Wels (left prescirbed flow boundary)
[XM,YM,ZM]=meshgrid(xm,ym,zm);

LRCleft=cellIndices(find(XM<xGr(2)),size(M),'LRC');  % cell cell indices L R C

q=Q_inflow.*squeeze(DZ)/sum(DZ);      % divide inflow over cells

u=ones(NZ,1);                       % we have only one stress period

WEL=[];
for iPer =1:NPER
    WEL=[WEL; [iPer*u LRCleft q(:)]];                   % construct well file
end
%% CHD right flow boundary

CHDDENSOPT=2;   % environmental head at ocean boundary, Langevin et al 2008, p22

% similar for CHD boundary conditions on right hand boundary of model
LRCright=cellIndices(find(XM>xGr(end-1)),size(M),'LRC');

CHD=[];
for iPer=1:NPER
    CHD=[CHD;[iPer*u LRCright u*[h_ocean h_ocean CHDDENSOPT]]];
end

%% PNTSRC for SSM source sink mixing module MT3DMS/Seawat

ITYPEleft =2;  % conc in prescribed flow at left boundary
ITYPEright=1;  % conc at prescribed head at right boundary     
dummy=0;
% Notice that we have two species, salinity and temperature

PNTSRC=[];
for iPer=1:NPER
    PNTSRC=[PNTSRC; [...
                        iPer*u LRCleft  u*[dummy ITYPEleft  c_inflow T_inflow];...
                        iPer*u LRCright u*[dummy ITYPEright c_ocean  T_ocean],...
                    ]...
           ];
end