%% recharge on a rectangular area with fixed head boundary at left and right
%  you don't have to do this, but it is genrally good to clear the
%  enrvironment before starting a new task

clear variables; close all;

mf_settings;

%% Read in GIF geological map

% [MAP, cmap] = imread( 'GoblinValley1_rotated_flowmodel_new.gif' );
% MAP = double(MAP);



% make model grid coordinates
xGr= 0 : cellsize : (size(MAP,2) * cellsize); 
yGr= 0 : cellsize : (size(MAP,1) * cellsize);

%% Some Housekeeping, see function modelsize in mflab/mfiles/fdm
% the grid coordinates are used to compute cell centers, cell size, grid
% size. xGr and yGr will be sorted aligned and freed of double coordinates
% Note that DELC are the column width (alterantive DX) and DELR are the row
% widths (alternative DY). Also notice that the MODFLOW manuals define DELR
% as width along rows, which are actually column widts and DELC as width
% along columns, which are actually row widths. I fear this terminology
% causes schizophrenia in the minds of modelers, so I just use DELC for
% column widths and DELR for row widths. But you can equally well use
% DX and DY as alternatives, which are recognized as such by mflab.

[xGr,yGr,xm,ym,DELC,DELR,NCOL,NROW]=modelsize(xGr,yGr);

% modelsize sorts the coordinates and eliminates repeated coordinates (so
% you may put any grids and subgrids arbitrarily together, modelsize will
% apply the funciton unique to sort them and eliminate doubles. 

%% lowerase z is used here for the elevation of the top and bottom of all layers
% mflab expects a 3D uppercase Z variable to contain the elevations of the
% all model planes at the (xm,ym)-cell centers. This is the genearal from
% of the 3D fdm-grid, where columns and rows each have a single width, but
% where adjacent cells in the horizontal plane may have different
% elevations. Hence Z(:,:,1) the top plane of the Z matrix corresponds to
% the top elevation of the model and Z(:,:,end) with the bottom elevation.
% If there are aquitards between some of the aquifers (model layers) then
% these planes must be included in Z. The number of planes in Z is
% therefore (top + number of aquifers + number of aquitards between
% aquifers). Aquitards may be defined in BCF package with the parameters LAYCBD, where
% LAYCBD is 1 for each aquifer having an aquitard underneath.
%
% Notice that all 3D arrays used in Matlab obey the natural orientation of
% Matlab, that is, (NROW,NCOL,NLAY). The first dimension in the direction
% of the rows, the second along te columns and the third along the layers. 
% As we will see, it's trivial to reorientate arrays (usign the permute
% function) where we need them, for instance to show cross sections in any
% direction.

% simple flat model, 1 high
z=[0 , -1 ];  % Plane elevation vector
z=reshape(z,[1,1,numel(z)]);       % Orientate into z direction
Z=z(ones(NROW,1),ones(1,NCOL),:);  % Construct 3D Z matrix directly


%% IBOUND array
NLAY=1;

% This example uses the BCF6 package which, like the BCF package allows
% aquitards (LAYCBD) between aquifers (layers). In this case we have 3
% aqufiers and 2 aquitards between, hence NLAY=3 (aquifers) according to
% the BCF package. Notice that the alternative LPF package defines each
% layer as  regular model layer, aquitards as well as aquifers. If we choose
% using LPF instead of BCF, then NLAY would be 5 in this example.

IBOUND=ones(NROW,NCOL,NLAY);  % all ones, i.e. ordinary cells
STRTHD=zeros(size(IBOUND));  % all zeros, matrix same size as IBOUND



% Set up head gradient:
if strcmpi( FlowDirection, 'horizontal')
    % E-W (Horizontal)
    IBOUND(:, 1  , 1) = -1;             % then change left column layers 1 & 2 into fixed heads
    IBOUND(:, end, 1) = -1;
    STRTHD(:, 1  , :) = 1;
    STRTHD(:, end, :) = 0;
elseif strcmpi( FlowDirection, 'vertical')
    % N-S (Vertical flow)              
    IBOUND(1  , :, :) = -1;             % then change left column layers 1 & 2 into fixed heads
    IBOUND(end, inflowrange, :) = -1;
    STRTHD(1  , :, :) = 1;
    STRTHD(end, inflowrange, :) = 0;
else
    error('No flow direction defined')
end

%% TRAN and VCONT
% 
% Saltwash : 263- 813 mD. Fine to medium grained sst, sometimes
% conglomeratic.
% Brushy Basin bentonitic shale and siltstones, Average perm below 10 mD.
% Cutler: Fluvial sst, silt and shale. sst : 2-900 mD, 

% Assign permeabilities (mD) to the integer classes from the GIF image:
PERM = nan( size(MAP));

for idx = 1 : length( imagevaluesused ),
    PERM( MAP == imagevaluesused( idx) ) = permeabilities( idx );
end

% PERM( MAP == 0) = 102570; % SLIP SURFACES 10e5 mD for 0.1 cm grid blocks, 25000mD upscaled for 0.5cm grid blocks .1cm aperture 86m/day Matthai et al., 1998 
% PERM( MAP == 1) = 5e-9;   % Shale smear 5 picodarcy
% PERM( MAP == 2) = 50;  % Fine SST
% PERM( MAP == 4) = 0.07;   % Cemented SST



if sum(isnan(PERM(:)))>0,  error('not all grid cells have been assigned permeability'), end

% save( 'permeability.mat', 'PERM');

HK = PERM .* 2.75 .* 0.001 .* (1/3.28);
%      Darcy to ft/day, mD to D, ft/day to m/day   
% Hydrology Handbook, American Society of Civil Engineers, p264
VK = zeros(size(MAP));

%% RCH (RECH)

% No recharge here:
RECH= zeros( NROW, NCOL, 1) ; %{3.0e-8*ones(NROW,NCOL,1)};  % recharge in first stress period as cell array


if exist( 'inflowrange' ) == 1,
    save('adapt_settings.mat','FlowDirection', 'inflowrange', 'PERM');
else
    save('adapt_settings.mat','FlowDirection', 'PERM');
end

%% Effective porosity
% not used now, same as above, it if provided mflab only uses it when it
% needs it for transport modelling (MT3DMS, MODPATH, SEWAT, SWI etc)
% PEFF=0.30 *ones(NROW,NCOL,NLAY);

%% Boundaries conditions, if necessary, merged with values from excel file
% MODFLOW may make use of the following boundary condtions that are given
% on a cell-by-cell basis.
% They are "wells (WEL)", "general head boundaries(GHB)", "drains (DRN)", "rivers (RIV)" and
% "time varying constant head boundaries (CHD)".
% For each of these boundary types there is a worksheet in the workbook in
% which the list of of values can be provided in the general form
% PeriodNR Layer Row Col values ...
% (mflab requires the period number to make each line uniqued and to allow
% a uniform input where missing periods simply mean all such boundaries are
% off).
% While specifying such boundaries in whe workbook may be convenient when
% the number of cells is small, it can be an enormous burden to define
% hundreds of thousands of boundary cells through the worksheet. Therefore,
% you may defined these boundary condtions both in the worksheet and
% directly in mflab or both whatever is most convenient. mflab simply
% combines both inputs into a sinble list for modflow.
% The same method is applied for the definiction of point sources in MT3D
% thorugh the spreadsheet and or variable PNTSRC in mf_adapt.
%
% The DRN worksheet for this problem looks like this 
%PERIOD  LAYER  ROW COLUMN  ELEVATION  COND
%     1      1    8      2          0     1
%     1      1    8      3          0     1
%     1      1    8      4         10     1
%     1      1    8      5         20     1
%     1      1    8      6         30     1
%     1      1    8      7         50     1
%     1      1    8      8         70     1
%     1      1    8      9         90     1
%     1      1    8     10        100     1
%
% It is staightforward to deifne the same in mf_adapt instead of providing
% these data in the worksheet, for example in one line:
%
% WEL=[ones(9,1)*[1 1 8] (2:10)' [0 0 10 20 30 50 70 90 100]' ones(9,1)];
%
% We may leave empty matrices as shown here in mf_adapt of leave them out
% altogether.

% WEL=[]; GHB=[]; DRN=[]; RIV=[]; CHD=[]; PNTSRC=[];

% In the current example, the DRNs are defined in the worksheet.

% Again, mflab will only use what it needs.
