clear;
close all;
%
% input parameter file for use by rgarden.m, MUST be run first!!!!!!!!
%
% IMPORTANT: assumes rain data is in [mm] given at HOURLY intervals!!!
%
% betapw=-0.04; % Pan-Wierenga transform parameter, between [-.01 -0.1] cm-1
%%%%%%%%%%%%%%%%%%%%%%%%%%
% RAIN GARDEN PARAMETERS %
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Storm_gap = 2;  % Time with no rain defining gap between storms [h]
g=9.8;		% acceleration of gravity [m/s2]
Sd=.25;		% depression storage in pervious area [cm]
Sd_imp = 0.2;           % Depression storage depth of impervious surface [cm]
Si_rg = 0.2;          % Interception storage depth of pervious surface [cm]
Sd_perv = 6.4;            % Depression storage depth of pervious surface CN=80 [cm]
Arg=100;		% area of rain garden [m2]
% R=10;		% ratio of impervious to pervious area
% Aimp=Arg*R; % impervious area connected to rain garden [m2]
% dmin=15;		% min depth for weir [cm]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SOIL PARAMETERS          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Soil parameter values for van Genuchten, M. T. 1980.
% SOIL LAYERS (Maidment 1992, Chapter 5)
%
POR = [0.41 0.41 0.39];				% porosity
NSAT = [0.41 0.41 0.39];				% water content at saturation
NRES = [0.02 0.041 0.075];				% residual water content
wiltpoint=0.033;
rzwilt = 0.1;                       %{DBA} Root zone wilting point, sandy loam from C&N
rzfc= 0.21;                         %{DBA} Root zone field capacity, sandy loam from C&N
b = 3.1;                         %{DBA} Pore size distribution index, from C&N
AA = [0.01 0.01 0.01];				% alpha [cm-1]
N  = [1.6 1.378 1.242];
% KSAT = [10.16 15 1.0];				% sat. hydraulic conductivity [cm/h]
% thick = [30 90 10];		% thickness of RZ, TZ, CZ layers [cm]
he=-1./AA;   % air entry suction [cm]
hwf=3*(2+3*(N(1)-1))/(1+3*(N(1)-1))*he(1)/2;%11;% wetting front suction of upper layer, Rawls&Brakensiek
Ss=[1e-6 1e-6 1e-6];   % specific storage [cm-1]
h0=[-2.4246e-004 -4.7992e-007 -1.2314e-010];
nsfactor=0.95;  % for tracking near saturated conditions in the root zone
nwiltfactor=1.21; %{DBA} for estimation of near wilting point
%
%%%%%%%%%%%%%%%%%%%%%%%%%
% SET UP DISCRETIZATION %
%%%%%%%%%%%%%%%%%%%%%%%%%
%
G=1;		% gravity factor for Richards equation
w=1.0;		% weighting for Crank-Nicholson (1: fully implicit)
Kavetype=1; % hydraulic conductivity averaging (=1 arithmetic, =2 geometric, =3 harmonic)
%
%
% TEMPORAL:
% --------
time=8759*60;					% min in 6 months April 1 - Sep 30
dt=15;						% TIME STEP [min]
dtdata=1;                 % [h]
adaptimestep=1;     % flag to tell to use adaptive time step
dtmin=1e-6; dtmax=dtdata;   % [h]
nt=round(time/dt);				% total number of temporal steps
% stime=time/60;				% total time of simulation [h]
% dt=dt/60;					% TIME STEP [h]
%
% SPATIAL:
% -------
% totdepth=thick(1,1)+thick(1,2)+thick(1,3);	% spatial extent [cm]
dx=2;						% SPATIAL STEP [cm]
% m=totdepth./dx+1;	            % total number of spatial nodes
% x=0:dx:totdepth; x=x';	    	% cm, x downwards, i nodes downwards
% z=flipdim(x,1);					% z upwards
% rf=thick(1,1)./dx+1;ti=rf+1;	% boundary root zone - transm zone
% tf=(thick(1,1)+thick(1,2))/dx; ci=tf+1;	% TZ - confining layer
%
% CONVERGENCE CRITERIA:
% --------------------
rlim=10;        % limit of picard iterations acceptable at a time level [10 50] in HYDRUS
toleranceth=0.0001;	% for theta absolute tolerance, iteration convergence (mass balance)
tah=0.5;           % h absolute tolerance [cm] for time iteration
tahp=0.1*tah;          % h absolute tolerance [cm] for mod Picard iteration, <tah prudent!
trh=0.05;     % h relative tolerance for time iteration
trhp=0.1*trh;    % h relative tolerance for mod Picard iteration, <trh prudent e.g. 0.01*trh!
sf=0.8;
fracmin=0.8;            % heuristic time stepping parameter, [0.3 0.9] (note 0.8 = 1/1.25)
fracmax=1.5;           % id., [1.1 1.5] in HYDRUS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIAL CONDITIONS          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hin=-100;					% [cm]
% for i=1:3
%     theta_init(i)=THETAH(hin,AA(i),N(i),NSAT(i),NRES(i),Ss(i),h0(i));
% end
% hinit=ones(m,1);
% hinit=hin*hinit;			% initial value of pressure heads [cm]
% Sinit=zeros(m,1);			% initial value of ET
%
%%%%%%%%%%%%%%%%%%%%%%%
% EVAPOTRANSPIRATION: %
%%%%%%%%%%%%%%%%%%%%%%%
etrun=1;		% flag to indicate to rg.m whether ET will be modeled
% Sold=Sinit;             % initial value for transpiration
avet=2.3/10;	% average daily ET [cm/day]
hatm=-500000;   % [cm], equivalent to ~70% RH for Seattle - annual
% {DBA} Qevap=avet/24/100*Arg;  % total evaporation (soil+plants) [m3/h]
Evap_imp = 10/10/24;    % Average rate of evaporation from impervious surface [cm/h]
Evap_perv = avet/24;  % Average rate of evaporation from pervious surface [cm/h]
ET_avg = avet/24;     % Average evapotranspiration [cm/h]
Cpan = .75;            % {DBA} Pan Evaporation correction coefficient
%altitude=300;				% a.m.s.l. [m]
%beta=0.5;					% canopy attenuation coefficient (0.5 for grass)
%LAI=2;						% leaf area index (3-4 for grass, 2 for prairie)
%tau=exp(-beta*LAI);			% canopy transmissivity
%alpha=1.26;					% Priestley & Taylor coefficient
%metfile=dlmread('met.txt','\t',1,0);
%temp=interp1(metfile(:,1),metfile(:,2),t,'linear');
%rad=interp1(metfile(:,1),metfile(:,3),t,'linear');
%ep=EPPT(temp,rad,altitude);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set top boundary condition as head (TYPE 1) [cm]          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%tbc=1;					% head - Dirichlet
%hf=dlmread('tbch.txt','\t',1,0);
%ht=interp1(hf(:,1),hf(:,2),t,'nearest');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set top boundary condition as flux (TYPE 2) [cm/h]        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hsbc=0; %tells it's not a fixed head top BBCC
tbc=2;					% flux - Neuman
swrun=1;		% =0 means no raingarden surface water balance modeled
rgrun=1;        % =0 means no enclosure in raingarden i.e. no ponding and weir
if rgrun==0; dmin=Sd; end;          % no rg condition, depression storage minimum ponding
% hs=hin;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  RAIN DATA RETRIEVAL                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Read data file from working directory.  The file must be a tab-delimited
% text file with two columns.  The first column must be time [h] and the 
% second column must be precipitation depth [Mm].  The first row of each
% column should contain the title of the column, as the model will not look
% at that row.
%
% nprecip=dlmread('rainmad95h.txt','\t',[1 0 ceil(stime)+1 1]);
% nprecip(:,2)=nprecip(:,2)./10;%convert from [mm h-1] to [cm h-1]
% qs=nprecip(1,2);    %first value for surface flux [cm h-1]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set bottom boundary condition as pressure head (TYPE 1) [cm]  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bbc=1; hb=hin; hinit(m)=hb;				% head - Dirichlet
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set bottom boundary condition as unit gradient (gravity flux) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%bbc=3;
%
%%%%%%%%%%%%%%%%
% OUTPUT FILE: %
%%%%%%%%%%%%%%%%
% outfile='rg95madwiR10.txt';
