% Dutch top system consisting of a shallow Holocene aquifer underlain by a
% semi-pervious layer on top of a larger regional aquifer. The shallow
% aquiifer is intersected by may dithces in which the surface water level
% is maintained to manage the groundwater system.
% We simulate two parallel ditches in this shallow top system underlain by
% the confining bed and the regional system. THe ditches may cut through
% the confining layer.
% Input is recharge and seepage from below through the confining bed from
% the regional aquifer. Extraction is evapotranspiration and leakage to the
% ditches. These dithches may at times also infiltrate.
% The specific yield will be made dependent on the the water table
% elevation in the shallow aquifer.
% When this water table intersects ground surface, surface runoff occurs.
% There is no limit to the infiltration capacity of the soil.
% Evapotranspiration may be made depdendent on the extincktion depth.
% We may apply various tricks to influence the behavior of the boundary
% conditions such as more than one drainage level wiithin a single cell.
% The smallest conceivable model consists of only a single cell.
% It is sufficient to simuulate only half the strip of land between the
% ditch and its axis of symmetry at the center.
clear variables; close all;

basename='DutchTop';

% to ccycle through the list of timesteps if larger than the number of
% allowed or specified stress periods use mf_start.m. This initializes
% counters necessary


%% Read meteo time series

fid=fopen('PE-92-00.txt','r');
fgets(fid); A=fscanf(fid,'%d-%d-%d %f %f',[5,Inf])';
fclose(fid);

tne=[datenum(A(:,3),A(:,2),A(:,1)) A(:,[4 5])/1000]; % [t P N] % to mm/d

fprintf('Length of time series = %d\n',length(tne(:,1)));
%plot(tne(:,1),tne(:,2),'b',tne(:,1),tne(:,3),'r');

clear A;

%% mf2k has 1000 stress period limit, use mf2005 instead, see nam sheet

RECH=ones(1,1,length(tne(:,1))); RECH(1,1,:)=tne(:,2)-tne(:,3);
%EVTR=ones(1,1,length(tne(:,1))); EVTR(1,1,:)=tne(:,3);

%% Get data from GGOR basis data spreadsheet (should be database)

%[parnams,parvals]=getExcelData('GGOR_basisdata_Noorderpark','','Horizontal');
%[parnams,parvals]=getExcelData('GGOR_basisdata_test','','Horizontal');

[parvals parnams]=dbfread('GGOR_basisdata');

GP = [parvals{:,strmatchi('PEIL',parnams,'exact')}]';  GP=round(100*GP)/100; % average phreatic head
ZP = [parvals{:,strmatchi('ZOMERPEIL' ,parnams  )}]';  ZP=round(100*ZP)/100; % summer ditch level
WP = [parvals{:,strmatchi('WINTERPEIL',parnams  )}]';  WP=round(100*WP)/100; % winter dithc level
Z1 = GP-6;                                             % bottom of top layer (default thickness 6 m)
L  = [parvals{:,strmatchi('L', parnams,'exact'  )}]';  L(L<20)=20; L =round(L);          % width of parcels
XC = [parvals{:,strmatchi('XC',parnams,'exact'  )}]';  XC=round(XC);         % parcel center x
YC = [parvals{:,strmatchi('YC',parnams,'exact'  )}]';  YC=round(YC);         % parcel center y
K  = [parvals{:,strmatchi('K' ,parnams,'exact'  )}]';  K =round(10*K)/10;    % parcel hor k
MU = [parvals{:,strmatchi('MU',parnams,'exact'  )}]';  MU=round(100*MU)/100; % parcel Sy
C  = [parvals{:,strmatchi('C' ,parnams,'exact'  )}]';  C(C>1000)=1000; C =round(C);          % parcel resistance (k_vert/d)
q  = [parvals{:,strmatchi('KWEL_GEM',parnams    )}]';  q(isnan(q))=0;        % upward seepage
phi= [parvals{:,strmatchi('WVP',parnams,'exact' )}]'; phi(abs(phi)>10)=0; phi=round(phi*1000)/1000; % head in regional aquifer
AHN= [parvals{:,strmatchi('AHN_MED',parnams     )}]'; AHN(abs(AHN)>10)=0; AHN=round(AHN*1000)/1000; % ground elevations
AREA=[parvals{:,strmatchi('AREA',parnams     )}]';

%% Catagorizing works but the number of cross section is not sufficiently

% binsL   = binRange(L,10);     catL  =sum(L  *ones(size(binsL ( 1:end-1)))<ones(size(L )) *binsL ( 2:end),2); % parcel width

%% Set the physical model parameters

P=[];
for i=1:length(GP)
    % coordinates and depth
    P(i).dBaseNr=i;           % identify this Xsection in the dbase file
    P(i).w   = L(i);          % half parcel size between water divide and center of ditch
    P(i).area= AREA(i);       % repersentative area of parcel
 
    P(i).x  = XC(i);          % parcel center x
    P(i).y  = YC(i);          % parcel center y
    P(i).D1 = 6;              % default thickness of top layer
    P(i).D2 = 25;             % default D of regional aquifer [m];
    P(i).z0 = AHN(i);         % parcel ground elevation
    P(i).z1 = P(i).z0-P(i).D1;  % default bottom elevation of top layer [mNAP]
    P(i).z2 = P(i).z1-P(i).D2;  % default bottom elevation of regional aquifer [mNAP];
    % heads
    P(i).h_mean   = GP(i);    % average phreatic head [NAP]
    P(i).h_summer = ZP(i);    % summer phreatic head [NAP]
    P(i).h_winter = WP(i);    % winter phreatic heat [NAP]
    P(i).phi      = phi(i);         % head in regional aqufier [NAP]
    % seepage
    P(i).q = q(i);            % vertical (upward) seepage [m/d]
    % yHraulic parameters
    P(i).hk1 = K(i);          % hor hudraulic conductivity of top layer [m/d]
    P(i).c   = C(i);          % vert resistance of top layer [d];
    P(i).vk1 = P(i).D1/P(i).c/2; % vert hydraulic conductiity of first layer [m/d]
    P(i).sy  = MU(i);         % specific yield of first layer [ ];
    P(i).ss  = 1e-5;          % default elastic storativity all layers [1/m]
    P(i).hk2 = 30;            % default k of regional aquifer [m/d];
    P(i).vk2 = 30/3;          % default vertical cond regional aqufier
    % Ditch
    P(i).wd    = 2;           % default width of ditch [m]
    P(i).cd    = 0.001;       % default entry resistance of ditch [d]
    P(i).dd    = 0.6;         % default ditch depth [m]
    P(i).omega = P(i).wd+2*P(i).dd;  % default circumference of ditch [m]
    % Drains
    P(i).cdr   = 0.1;         % default drain resistance [d]
    P(i).zdr   = P(i).z0;     % default drain depth [m]
end

MINAREA=1e4;            % [m2] minimum parcel area taken into account

P=P([P.area]>=MINAREA); % Only use parcels > 1 hectare

P=P(1:end);             % Use only these cross section

clear parnams parvals GP ZP WP L XC YC K MU C q phi AHN Z1  % clean up unnecessary variables


%% General data
dx=5;

xGr=[-0.001 0:dx:round(max([P.w]))]; % [m] width of model (max length of any parcel)
yGr=0:1:length(P);                % each x-section gets its own row (CHANI=0)
zGr=[0;-1;-2;-3];                 % default zGr values, will be adjusted to actual parcel
[xGr,yGr,xm,ym,Dx,Dy,Nx,Ny]=modelsize(xGr,yGr);

Z=NaN(Ny,Nx,3);
for iy=1:length(P)
    Z(iy,:,1)=P(iy).z0;
    Z(iy,:,2)=P(iy).z1;
    Z(iy,:,3)=P(iy).z2;
end

DZ=abs(diff(Z,1,3)); Nz=size(DZ,3);

[DX,DY]=meshgrid(Dx,Dy);


%% None number come useful for section from arrays layer on
iDRN  =1;    % zone of drains
iGHB  =2;    % zone of ditch
iWVP2 =3;    % zone of regional aquifer
iWEL  =4;    % zone where to dose seepage

%% Build IBOUND use zone number to identify zones
IBOUND=ones(Ny,Nx,Nz);
IBOUND(:,:,1)=iDRN;    % zone with drains (ground surface)
IBOUND(:,1,1)=iGHB;    % location of ditch
IBOUND(:,:,2)=iWVP2;   % second aquifer
IBOUND(:,1,2)=iWEL;    % location where seepage is added

% Because the width of the cross section varies and we only model one half
% of each cross section, some cells near the water divide will be inactive
% for smaller sections. This is effected by setting IBOUND to zero for
% those cells

% cells with xm>w of the section are made inactive
for iy=1:Ny, IBOUND(iy,xm>P(iy).w+eps,:)=0; end  % cells beyond water divide inactive

%% Build model arrays

HK    =ones(size(IBOUND));
VK    =ones(size(IBOUND));
SY    =ones(size(IBOUND));
SS    =ones(size(IBOUND));
STRTHD=ones(size(IBOUND));

HK(:,:,1)    =[P.hk1]'*ones(size(xm));  % kh of top aquifer
HK(:,:,2)    =[P.hk2]'*ones(size(xm));  % kh of bottom aquifer
VK(:,:,1)    =[P.vk1]'*ones(size(xm));  % kz of top aquifer
VK(:,:,2)    =[P.vk2]'*ones(size(xm));  % kz of bottom aquifer
SS(:,:,1)    =[P.ss]' *ones(size(xm));  % SS of top aquifer
SS(:,:,2)    =[P.ss]' *ones(size(xm));  % SS of bottom aquifer
SY(:,:,1)    =[P.sy]' *ones(size(xm));  % SY of top aquifer
SY(:,:,2)    =[P.sy]' *ones(size(xm));  % SY of bottom aquifer
STRTHD(:,:,1)=[P.h_mean]' *ones(size(xm));  % start head of top aquifer
StRTHD(:,:,2)=[P.phi]'*ones(size(xm));  % start head of bottom aquifer

%% Drains at ground surface to compute surface runoff

CDRN=NaN(size(IBOUND));              % allocate memory for drain conductance
CDRN(:,:,1)=(Dy*Dx)./([P.cdr]'*ones(size(Dx)));  % comptue drain conductance
CDRN(IBOUND==0)=NaN;                 % remove inactive cells
IDRN=find(~isnan(CDRN));             % get global incices of drain locations

hDRN=NaN(size(IBOUND));              % allocate memory
hDRN(:,:,1)=[P.z0]'*ones(size(Dx));  % [Ny*Nx] drain head is equal to z0, ground surface

%% Head boundaries in the left-most cell to simulate the ditch

CGHB=0.5*[P.omega]'./[P.cd]'; % conductance of ditch (half circumfernce/resistance)
hGHB=[P.h_mean]';             % GHB heads in the left ditch only, whill change during stress periods
IGHB=find(IBOUND==iGHB);      % global indices of ditch

%% Wells to simulate seepage from bottom aquifer or vice versa

% The specific upward positive flux is given. We compute the total flux by 
% by multiplying this areal flux withe the area of the model
% This fulx is given and constant for all stress periods

qWEL=sum(([P.q]'*ones(size(Dx))).*(Dy*Dx).*(IBOUND(:,:,1)>0),2);
IWEL=find(IBOUND==iWEL);     % Global indices of well

%% using the global incices compute the Layer Row Column indices necesary
%  to specify the boundary conditions in MODFLOE throuth the DRN, GHB and
%  WEL packages

LRCdrn=cellIndices(IDRN,size(IBOUND),'LRC');
LRCghb=cellIndices(IGHB,size(IBOUND),'LRC');
LRCwel=cellIndices(IWEL,size(IBOUND),'LRC');

uDRN=ones(size(IDRN)); % a column of ones for easyier multiplying below
uGHB=ones(size(IGHB)); % same
uWEL=ones(size(IWEL)); % same

%% specify the arrays for the mentioned boundary conditions so that
% mfLab can generate the input files for MODFLOW

[pernams,pervals]=getPeriods(basename);
NPER=size(pervals,1); % number of stress periods defined in worksheet PER

[Y,M]=datevec(tne(:,1)); summer=(M>4 & M<10);


iPer=1;

if summer(iPer), h=[P.h_summer]; else h=[P.h_winter]; end  % NP heads

DRN=[iPer*uDRN LRCdrn hDRN(IDRN) CDRN(IDRN) ];
GHB=[iPer*uGHB LRCghb h' CGHB ];
WEL=[iPer*uWEL LRCwel qWEL];
for iPer=2:NPER
    if summer(iPer)==summer(iPer-1), % use previous
        GHB=[GHB; -iPer, ones(1,5)];
    else
        if summer(iPer), h=[P.h_summer]; else h=[P.h_winter]; end
        GHB=[GHB;  iPer*uGHB LRCghb h' CGHB]; % same in as first period
    end
    WEL=[WEL; -iPer ones(1,4)]; % same as in first period
    DRN=[DRN; -iPer ones(1,5)]; % same as in first period
end

%% What to do after mf_setup has finished and the model has run ?

%AFTERMFSETUP='mf_analyze';  % uncomment this line if you want to run
% mf_analyze immediately after mf_setup.  This has the risk that if modflow
% fails, mf_analyze uses the old .hds and .bgt files as it has no clue on
% whether modflow finished normally

save underneath P tne % P and tne may be needed by mf_analyze

