% mf_setup  --- generic setup of MODFLOW model, reads data from spreadsheet
close all

P1='D:\Tijs\FinalModel\Model'
P2='D:\Tijs\FinalModel\processScripts'
path(path,P1); path(path,P2);

basename='testMODEL' 
demname='groundElev'

load([basename,'.mat']);

XLSF=['D:\Tijs\FinalModel\model\',basename,'.xls']

%% test
if 0
dx=(1:30); Nx=length(dx);
dy=(1:30);  Ny=length(dy);

top=zeros(length(dy),length(dx));
bodem1=top-10;
bodem2=top-20;
bodem3=top-40;

dx=1:10;
dy=1:10;

startheads=top-0.2; 
% end of test ================
end 
%% SIMULATION INFO
[Snum,Stxt]=xlsread(XLSF,'SIMULATION');
Spackage =Stxt(1,:); 
Sheaders =Stxt(2,:);
Stxt     =Stxt(3:end,:);
%% STRESS PERIOD INFO
[Pnum,Ptxt]=xlsread(XLSF,'PERIODS');
Ppackage =Ptxt(1,:); 
Pheaders =Ptxt(2,:);
Ptxt     =Ptxt(3:end,:);

%% LAYER INFO
[Lnum,Ltxt]=xlsread(XLSF,'LAYERS');
Lpackage =Ltxt(2,:); 
Lheaders =Ltxt(3,:);
Ltxt     =Ltxt(4:end,:);

%% NAM Namefile INFO
[Nnum,Ntxt,nam]=xlsread(XLSF,'NAMFILE');
writeNAM(basename,nam)

%% DISCRETIZATION struct buildup
dis.NROW=length(dy);  % default matlab variable
dis.NCOL=length(dx);  % default matlab variable
dis.NLAY=size(Lnum,1);
dis.NPER=size(Pnum,1);

dis.DELR   = dx;  % column sizes
dis.DELC   = dy;  % row sizes

dis.LAYCBD = Lnum(:,strmatchi('LAYCBD',Lheaders));

%% Layer elevations

%eval(['load ',demname]);  % off in test

eval(['dis.TOP=' demname,';']);

dis.BOTM =cell(dis.NLAY,1);
for iL=1:dis.NLAY
    varname=Ltxt{iL,strmatchi('BOTM',Lheaders)};
    %eval(['load ',varname]);
    eval(['dis.BOTM{iL} =',varname,';']);
end

%% STRESS PERIODS
dis.PERLEN = Pnum(:,strmatchi('PERLEN',Pheaders));
dis.NSTP   = Pnum(:,strmatchi('NSTP',  Pheaders));
dis.TSMULT = Pnum(:,strmatchi('TSMULT',Pheaders));
dis.SSTr   = Ptxt(:,strmatchi('SS/Tr' ,Pheaders));

writeDIS(basename,dis); % write the modflow DIS file

transientflag=0;
for i=1:length(dis.SSTr)
    if strcmp(upper(dis.SSTr{i}),'TR')
        transientflag=1;
    end
end
%% bas-file
bas.IBOUND=cell(dis.NLAY,1);
for iL=1:dis.NLAY
    varname=Ltxt{iL,strmatchi('IBOUND',Lheaders)};
    %eval(['load ',varname]);                        % voor test uit
    eval(['bas.IBOUND{iL} =',varname,';']);
end

bas.STRT=cell(dis.NLAY,1);
for iL=1:dis.NLAY
    varname=Ltxt{iL,strmatchi('STRT',Lheaders)};
    %eval(['load ',varname]);                        % voor test uit
    eval(['bas.STRT{iL} =',varname,';']);
end

writeBAS6(basename,bas);

%% LPF-file  (layer property flow package)
% unit to write cell-by-cell flows

lpf.Options=cell(1); % see documentation page 8.29

% per simulation
lpf.ILPFCB=Snum(1,strmatchi('ILPFCB',Sheaders));
lpf.HDRY  =Snum(1,strmatchi('HDRY'  ,Sheaders));
lpf.NPLPF =Snum(1,strmatchi('NPLPF' ,Sheaders));

% for all layers
lpf.LAYTYP=Lnum(:,strmatchi('LAYTYP',Lheaders));
lpf.LAYAVG=Lnum(:,strmatchi('LAYAVG',Lheaders));
lpf.CHANI =Lnum(:,strmatchi('CHANI' ,Lheaders));
lpf.LAYVKA=Lnum(:,strmatchi('LAYVKA',Lheaders));
lpf.LAYWET=Lnum(:,strmatchi('LAYWET',Lheaders));
if any(lpf.LAYWET)
    lpf.WETFCT=0.5;
    lpf.IWETIT=1;
    lpf.IHDWET=0;   % 0: use: h=BOT+WETFCT (hm-BOT), else use: h=BOT+WETFACT(THRESHOLD)
end

HK    =Lnum(:,strmatchi('HK'  ,Lheaders));     % kx
CHANI  =Lnum(:,strmatchi('HANI',Lheaders));     % kx
VKA   =Lnum(:,strmatchi('VKA' ,Lheaders));   % kz
SS    =Lnum(:,strmatchi('SS'  ,Lheaders));   % kz
Sy    =Lnum(:,strmatchi('Sy'  ,Lheaders));   % kz
WETDRY=Lnum(:,strmatchi('WETDRY'  ,Lheaders));   % kz

lpf.LAYCBD=dis.LAYCBD;     % need LAYBCD here too, separate aquitard  below current layer ?
lpf.HK  =cell(dis.NLAY,1);
lpf.HANI=cell(dis.NLAY,1);
lpf.VKA =cell(dis.NLAY,1);
lpf.SS  =cell(dis.NLAY,1);
lpf.Sy  =cell(dis.NLAY,1);
lpf.VKCB=cell(dis.NLAY,1);
lpf.WETDRY=cell(dis.NLAY,1);
for iL=1:dis.NLAY
    lpf.HK{iL} =   HK(iL)  *ones(dis.NROW,dis.NCOL);
    lpf.VKA{iL}=  VKA(iL)  *ones(dis.NROW,dis.NCOL);
    if CHANI(iL)<=0
        lpf.HANI{iL}=CHANI(iL) *ones(dis.NROW,dis.NCOL);
    end
    lpf.SS{iL}=   SS(iL)   *ones(dis.NROW,dis.NCOL);
    if lpf.LAYTYP(iL)~=0
        lpf.Sy{iL}=Sy(iL)  *ones(dis.NROW,dis.NCOL);
    end
    if dis.LAYCBD~=0
        lpf.LAYCBD{iL}=VKCB(iL)  *ones(dis.NROW,dis.NCOL);  % kz of aquitard
    end
    if lpf.LAYTYP(iL)~=0 && lpf.LAYWET(iL)~=0
        lpf.WETDRY{iL}=WETDRY(iL)* ones(dis.NROW,dis.NCOL);
    end
end

writeLPF(basename,lpf,transientflag);

%% rch-file

rch.NROW=dis.NROW;
rch.NCOL=dis.NCOL;

%1
%2
rch.NRCHOP =Snum(1,strmatchi('NRCHOP',Sheaders));
rch.IRCHCB =Snum(1,strmatchi('IRCHCB',Sheaders));

%3
%4
%5
rch.INRECH =Pnum(:,strmatchi('INRECH',Pheaders));
rch.INIRCH =Pnum(:,strmatchi('INIRCH',Pheaders));
%6
rch.RECH   =Pnum(:,strmatchi('RECH'  ,Pheaders));
%7
%8
rch.IRCH   =Pnum(:,strmatchi('IRCH'  ,Pheaders));

writeRCH(basename,rch) 

%% EVT-file

evt.NROW=dis.NROW;
evt.NCOL=dis.NCOL;

%1
%2
evt.NEVTOP =Snum(1,strmatchi('NEVTOP',Sheaders));
evt.IEVTCB =Snum(1,strmatchi('IEVTCB',Sheaders));

%3
%4
%5
evt.INSURF =Pnum(:,strmatchi('INSURF',Pheaders));
evt.INEVTR =Pnum(:,strmatchi('INEVTR',Pheaders));
evt.INEXDP =Pnum(:,strmatchi('INEXDP',Pheaders));
evt.INIEVT =Pnum(:,strmatchi('INIEVT',Pheaders));

%6
evt.SURF =cell(dis.NPER,1);
for iL=1:dis.NPER
varname=Ptxt{iL,strmatchi('SURF',Pheaders)};
    %eval(['load ',varname]);
    eval(['evt.SURF{iL} =',varname,';']);
end

%7
evt.EVTR =Pnum(:,strmatchi('EVTR',Pheaders));
%8
%9
evt.EXDP =Pnum(:,strmatchi('EXDP',Pheaders));
%10
evt.IEVT =Pnum(:,strmatchi('IEVT',Pheaders));
 
writeEVT(basename,evt) 

%% DRN-file
[drn.values,drn.headers]=xlsread(XLSF,'drains');
drn.IDRNCB =Snum(1,strmatchi('IDRNCB',Sheaders));
drn.ITMP   =Pnum(:,strmatchi('ITMP'  ,Pheaders));

% the ones(size(drains,1)) adds period number 1
drn.values=[drn.values;[ones(size(drains,1),1), drains]];

writeDRN(basename,drn)

%% LAK-file

lak.NLAKES =Snum(:,strmatchi('NLAKES',Sheaders));
lak.ILKCB  =Snum(:,strmatchi('ILKCB' ,Sheaders));
lak.THETA  =Snum(:,strmatchi('THETA' ,Sheaders));
lak.SSNCR  =Snum(:,strmatchi('SSNCR' ,Sheaders));
lak.NSSITR =Snum(:,strmatchi('NSSITR',Sheaders));
lak.STAGES =Snum(:,strmatchi('STAGES',Sheaders));
lak.SSMN   =Snum(:,strmatchi('SSMN'  ,Sheaders));
lak.SSMX   =Snum(:,strmatchi('SSMX'  ,Sheaders));
lak.NSLMS  =Snum(:,strmatchi('NSLMS' ,Sheaders));

lak.ITMPLK =Pnum(:,strmatchi('ITMPLK',Pheaders));
lak.ITMP1  =Pnum(:,strmatchi('ITMP1' ,Pheaders));
lak.LWRT   =Pnum(:,strmatchi('LWRT'  ,Pheaders));

lak.LKARR  =Pnum(:,strmatchi('LKARR' ,Pheaders));
lak.BDLKNC =Pnum(:,strmatchi('BDLKNC',Pheaders));
lak.NSLMS  =Pnum(:,strmatchi('NSLMS' ,Pheaders));
lak.PRCPLK =Pnum(:,strmatchi('PRCPLK',Pheaders));
lak.EVAPLK =Pnum(:,strmatchi('EVAPLK',Pheaders));
lak.RNF    =Pnum(:,strmatchi('RNF'   ,Pheaders));
lak.WTHDRW =Pnum(:,strmatchi('WTHDRW',Pheaders));
lak.SSMN   =Pnum(:,strmatchi('SSMN'  ,Pheaders));
lak.SSMX   =Pnum(:,strmatchi('SSMX'  ,Pheaders));

writeLAK(basename,lak,transientflag)
%% OC-file
oc.IHEDFM=Snum(:,strmatchi('IHEDFM',Sheaders));
oc.CHEDFM=Snum(:,strmatchi('CHEDFM',Sheaders));
oc.IHEDUN=Snum(:,strmatchi('IHEDUN',Sheaders));
oc.IDDNFM=Snum(:,strmatchi('IDDNFM',Sheaders));
oc.CDDNFM=Snum(:,strmatchi('CDDNFM',Sheaders));
oc.IDDUN =Snum(:,strmatchi('IDDUN'  ,Sheaders));

oc.INCODE =Pnum(:,strmatchi('INCODE',Pheaders));
oc.IHDDFL =Pnum(:,strmatchi('IHDDFL',Pheaders));
oc.IBUDFL =Pnum(:,strmatchi('IBUDFL',Pheaders));
oc.ICBCFL =Pnum(:,strmatchi('ICBCFL',Pheaders));

oc.Hdpr =Lnum(:,strmatchi('Hdpr',Lheaders));
oc.Ddpr =Lnum(:,strmatchi('Ddpr',Lheaders));
oc.Hdsv =Lnum(:,strmatchi('Hdsv',Lheaders));
oc.Ddsv =Lnum(:,strmatchi('Ddsv',Lheaders));

writeOC(basename,oc)

%% pcg-file
[pcg.values,pcg.headers]=xlsread(XLSF,'PCG');
writePCG(basename,pcg);

%% run programma
%eval(['dos(''D:\Tijs\FinalModel\Model\mf2k.exe ',basename,'.nam'')']);
eval(['dos(''',basename,'.nam'')']);

%% Read output

heads=readDat([basename,'.hds']);
bud=readBud([basename,'.bgt']);
contour(xm,ym,heads.values{1});
