% Example see USGS Modflow 2000 manual, Open-File Report 00-92
% TO 090806 091129

% McDonald and Harbaugh, 1988, appendix D (Modular Three-Dimensional Finite
%-Difference Groundwater Flow Model, Book 6, Chapter A1 USGS)
% The sampe problem is intended to illustrate input and output form the program.
% There are three layers, as shown in the accompanying illustration, which
% are separated fro each other by confining layers. Each layer is a square
% 75000 feet on a side and is divided by a grid into 15 rows and 15 columns
% which form  squares 5000 feet on a side. Flow within te confining layers
% is not simulated, but the efects of the confining layers on tlow betwen
% the active alyers are incoporated in the vertical leakance (Vcont) terms.
% Flow nto the system is infiltraiton form precipitation; flow out of the
% system is to buriekd drain tubes, discharing wells, and a lake which is
% represented by a constant-head bounadry.
%
% The drains and wells are specified in the worksheets in the workbook <<basename>>.xls
% where <<basename>> in this example is simply "ex1"
%
% The first layer is a water table lager with k=0.001 ft/s and bottom at
% z=-150 ft.
%The second and third layers (aquifers) are confined with transmissivities
%of 0.01 and 0.02 ft^2/s respectively. The leakances of the two confining
%layers are 2e-5 and 1e-5 /s respectively.
%Recharge is 3e-8 ft/s.
% The head at the left side (column number 1) is constant and equal to zero
% in the two uppermost layers.
% There are 15 wells as specified in the worksheet WEL.
% there are 8 drains as specified in the worksheet DRN.
% Although the elevations of the lower layers are not necessary for the
% computation and have not been provided originally by McDonald and
% Harbaugh, they have been specified in the manual for MDLFOW2000 (mf2k) by
% Hill and others. These values are used here.

%% 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;

%% The model name. Every model has a "basename" and all associated files
% will be of the form basename.ext where ext depends on the type of info
% specified in the nam file (basename.nam)

basename='ex1';

%% Specify grid line coordinates
xGr=0:5000:75000;  % feet
yGr=0:5000:75000;  % feet

%% 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. This
% tremendously faciltates merging grids. modelcare further computes the
% coordinates xm and ym of the cell centers, the cell sizes (DELX=DELC
% DELY=DELR) and finally returns the size of the grid (in terms of number of
% columns and rows). There is also a function modelsize3 to do the same in
% 3D. See help modelsize3. Lastly, modelsize puts the x-coordinates in
% ascending order and the y (and in 3D also the z) coordinates in
% descending order, so that the first line (and top plane) in the matrix/array
% of grid values always corresponds with the highest y and or z coordinate,
% which implies that a listing of the matrix has the same orentation as the
% coordinates themselves.

%% 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.

z=[200 -150 -200 -250 -350 -450];  % 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

% this was advanced ! Alternative
% Z=NaN(NROW,NCOL,numel(z));
% for iz=1:length(z), Z(:,:,iz)=z(iz); end

%% IBOUND array
NLAY=3;

% 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
IBOUND(:,1,1:2)=-1;             % then change left column layers 1 & 2 into fixed heads

%% STRTHD

STRTHD=zeros(size(IBOUND));  % all zeros, matrix same size as IBOUND

%% TRAN and VCONT

% The model has one water table (phreatic) aquifer and two confined
% aqufiers. For the water table aquifer we have to provide the horizontal
% conductivity (KH alternative HK) and for the two confined aquifers we have
% to provide their transmissivities TRAN. Notice that the variable names in
% mflab are the same as those in the MODFLOW manual to prevent confusion.
% Firther, because we use the BCF package, we have to provide the VCONT of
% the two aquitards:

% We need as many layers as there are water table aquifers (1 in this case)
% Generally, we provide mflab with a matix of the size of the entire grid.
% mflab picks out the layers that it needs (which hare defined as water
% table layers). As a consequence it is sufficient when the number of layers
% provided equals the highest layer number that is convertible to a water
% table aquifer. In this case we need only one layer:
HY =1e-3*ones(NROW,NCOL);     % hor confuctivity top water table layer ft/s

% We need only the transmissivities of layers 2 and 3. However, we provide
% mflab with a matrix  having the size of the entire grid. However, only
% the layers that MODFLOW actually needs (in this case layers 2 and 3) have
% to have the actual numeric values. Starting with a matrix consisting of
% NaNs is effective to prevent putting the numbers in the wrong layer.
TRAN=NaN(NROW,NCOL,NLAY);     % use NaN or os when not used
TRAN(:,:,2)=1e-2;             % transmissivity of second aquifer [ft^2/s]
TRAN(:,:,3)=2e-2;             % transmissivity of aqufier 3

% The number of VCONT layers cannot be larger than NLAY-1. As before, the
% VCONT array needs not have more layers than the number of the lowest
% layer having an aquitard beneath it.
VCONT=NaN(NROW,NCOL,NLAY-1);  % use NaN or so if no VCONT beloow layer
VCONT(:,:,1)=2e-8;            % the Leakance (thickness/conductivity of aquitard
VCONT(:,:,2)=1e-8;            % VCONT of second aquitard

%% RCH (RECH)

% Recharge can be define both in the workbook in worksheet PER or
% separately in mf_adapt like is done here.
% The rule is as follows, if the INRECH value in the stress period is zero, 
% then mflab expects the values to be defined in mf_adapt (here) in the variable RECH.
% Else it takes the RECH value from the PER worksheet. But the worksheet
% can only provide one value per stress period for the entire model.
% Therefore, providing cell by cell values have to be done in mf_adapt.
%
% To obtain maximum flexibility in a general fashion, RECH defined in mf_adapt
% has to be defined as a cell array, with one cell per stress period. Each
% cell can either hold an array of size (NROW,NCOL) with the cell-by-cell
% values for that stress period, or a single value as a constant for the
% entire model.
RECH={3.0e-8*ones(NROW,NCOL,1)};  % recharge in first stress period as cell array

% The rule by which mflab interprets the input is as follows:
% For any stress period in the PER worksheet where INRECH>0, it uses the
% the spreadsheet value for the entire  model whether or not RCH is defined
% in mf_lab. For any stress period where INRECH=0 the mf_lab defined values
% are used.

% The same method is applied in similar circumstances like EVTR and other
% parameters required by the EVT package.
%
% If no recharge is desired, switch the RCH package off in the NAM
% worksheet

%% Storage coefficients (only if transient)
% not used now, they are only needed for transient flow situations
% You may provide them, but mflab will only use them when necessary, i.e.
% when at least one stress period is transient.
%SY  =0.25 *ones(NROW,NCOL,NLAY);
%SS  =0.001*ones(NROW,NCOL,NLAY);

%% 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.

%% Stess periods

% stress periods are defined entirely in the PER worksheet. In this example
% the necessary stress period definition looks like (see worksheet)

%          DIS/BTN  DIS/BTN   DIS/BTN   DIS
%   IPER   PERLEN   NSTP      TSMULT    Transient
%      1     86400         1        1           0

% The first line shows the packages that need the specific stress period
% data. The second line shows the actual stress period parameters. The
% third and following lines define the actual stress period parameter
% values. IPER is the stress period number, PERLEN its length, NSTP the
% number of steps within this stress period. TSMULT the mutiplier applied
% to the next time step length in this stress period and Transient defines
% whether this stress period is steady state or transient. Use 0 for
% steady state and any non-zero value of transient. This replaces SS and TR
% in the MODLFOW manual because numeric values allow much easier handling.
% Only if any stress period is transient will mflab read SS and or SY.

%% Layers

% Layer properties are defined in the LAY worksheet. In this case this
% worksheet looks like this:
%
%         DIS     BCF   BCF/LPF  BCF/LPF  BCF/LPF  BCF/LPF  LPF    LPF
%LAYER    LAYCBD  TPRY  LAYAVG   LAYCON   LAYWET   WETDRY   CHANI  LAYVKA
%    1         1     1        0        1        1        1      1       0
%    2         1     1        0        0        1        1      1       0
%    3         0     1        0        0        1        1      1       0
%
% the top line show the packages that need the layer perameter. The second
% line shows the actual parameter and the 3 and following lines show the
% parameter values.
% In mflab the LAYER numbers are obligatory. On the other hand, you need
% only provide the layers that are uninque. In this case 2, would suffice
% as layers 2 and following have the same type properties (not
% conductivities). Notice tha LAYCBD of the third (in this case last layer)
% is 0 as there is no aquitards below the last layer. However as this is
% always the case, mflab takes care of that automatically, so that
% definition of the third layer is not necessary in this case.
%
% In general, using the layer numbers mflab assumes that missing layers
% have the same properties as the last layer for which properties were
% defined. So if layers 1 2 3 8 20 are defined in a 50 layer model, then
% layers 3-7 have the same type properties as layer 2, layers 9-19 as layer
% 8 and layers 21-50 as layer 20. Note however, that hydraulic properties,
% are all defined in mf_lab and will generally differe from layer to layer
% despite all layers haveing the same type properties defined in the LAY
% worksheet.
%

%% Transfer data to mf_setup

% finally the defined data are saved for layer retrieval. This allows
% cleaning the workspace and retrieving the data in a clean environment as
% is usual when running mf_analyze.

% this saving is done in mf_setup immediately after having run mf_adapt
% (this file).

% mf_setup saves only the data necessary for mfmab in the <<basename>>.mat file
% then clears the workspace, retrieves the mat file into the clean
% workspace and starts writing input files for modflow, mt3dms, seawat or
% whatever models and or packages are wanted according to the switches in
% teh NAM worksheet.



%% Running mflab by typing mf_setup
% This is what shows in the command window of MATLAB after typing mf_setup
%
%
% mf_setup
% ======= Start construction of your model ============
% Your current HOME directory is
% Z:\tolsthoorn On My Mac\GRWMODELS\mflab\examples\mf2k\ex1
% Running model construction script mf_adapt..
% ... mf_adapt finished, model has been constructed.
% Saving data to <<ex1.mat>>
% save([basename,'.mat'],'xGr','yGr','Z','IBOUND','STRTHD','DELR','DELC','KH',
%          'TRAN','VCONT','RECH','WEL','DRN','RIV','GHB','CHD','PNTSRC');
% .... Saved !
% Cleanup workspace ....... workspace cleaned and <<ex1.mat>> reloaded.
% Defining paths to your excecutables
% 
% Starting generation of input files ....
% Preparing name file struct.
% Basename current model is 'ex1'
% Getting nam file data from ex1.xls
% # MODFLOW etc. writeNAM 29-Nov-2009 21:54:46
% Getting simulation parameters from worksheet MFLOW    in ex1.xls
% Getting stress period data from    worksheet PER      in ex1.xls
% Getting layer info from            worksheet LAY      in ex1.xls
% Generating basic struct
% # MODFLOW writeBAS6 29-Nov-2009 21:54:49
% Generating discretization struct
% # MODFLOW writeDIS 29-Nov-2009 21:54:49
% Generating BCF struct and file
% # MODFLOW writeBCF 29-Nov-2009 21:54:49
% Generating RCH struct and file
% # MATLAB writeRCH, 29-Nov-2009 21:54:49
% Generating WEL struct and file
% # MATLAB writeBCN for WEL,  29-Nov-2009 21:54:51
% Generating DRN struct and file
% # MATLAB writeBCN for DRN,  29-Nov-2009 21:54:52
% Generating OC struct and file
% # MATLAB writeOC 29-Nov-2009 21:54:52
% Generating PCG struct and file
% # MATLAB writePCG 29-Nov-2009 21:54:52
% .... finished generation of input files, all files closed! Time= 29-Nov-2009 21:54:52!
% 
% Elapsed time is 7.535516 seconds.
%  
%                                  MODFLOW-2000 
%    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL 
%                             Version 1.18.01 06/20/2008                       
%  
% Using NAME file: mf2k.nam                                                                                                                                                                                                 
% Run start date and time (yyyy/mm/dd hh:mm:ss): 2009/11/29 21:54:52 
%  
% Run end date and time (yyyy/mm/dd hh:mm:ss): 2009/11/29 21:54:52 
% Elapsed run time:  0.172 Seconds 
%  
% Normal termination of MODFLOW-2000 
%
%% The next step would be to run mf_analyze
% load mf_analyze in the Matlab editor and run it block by block or in its
% entirety.

