function [Dprecip,Dimp,Dperv,Drunon,Drunoff,Devap,Drecharge,Dudrain,Dstore,tctpmax,tstponded,inot] = recarga_us(nprecip,Sd,avet,rundays,outpre,Summary,Record,dmin,Srz,St,Rsoilnumber,Ssoilnumber,Nsoilnumber,Krz,Kst,Kcz,Urate,CN,EIA,Arg,TribA);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%)
%  BASIC PROGRAM INFORMATION                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                           
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  FILE NAME:
%  recarga_us.m
%
%  GENERAL DESCRIPTION:
%  This program computes the values of the water balance terms (inflow,
%  infiltration, runoff, evapotranspiration, andrecharge) for a raingarden 
%  or area of lawn receiving runoff from animpervious surface, such as a 
%  roof.  
%  This version of the program takes English inputs, calculates in 
%  Metric and then outputs the results to Englishm,k.;
%
%
%  *MULTILAYER GREEN&AMPT!!!
%
%  ASSUMPTIONS:
%  1. Conditions corresponding to a summer rainy season (no snow).  
%  2. Vertical, one-dimensional infiltration in a uniform, partially 
%     saturated soil.
%  3. hwf (wetting front suction) is constant, and equal to root zone soil value
%     [this could be modified: follow wf and update hwf in time, but
%      it has been shown GA is not very sensitive to hwf parameter]
%
%  CALCULATION PROCEDURES:
%  1. Infiltration is calculated according to Green-Ampt.
%  2. Evaporation/evapotranspiration calculated using pan evaporation
%  records
%
%  REFERENCES:
%  1. van Genuchten, 1980  (soil characteristic functions)
%
%  REQUIRED FUNCTION M-FILES:
%  1. storageout.m
%  2. stormcheck.m, eventinfo.m
%  3. iga.m
%  4. relcond.m
%
%  HOW TO USE THE PROGRAM:
%  For each model run, the user will want to adjust the parameters and
%  input file in order to obtain the information he desires.  Adjust the
%  code at the appropriate locations prior to each model run.  Make sure
%  all units are consistent (see comments within code).  The input file
%  must be a text file containing two tab-delimited columns, time [h] and
%  precipitation depth [cm].  The program assumes the first row of each 
%  column is the column name, and eliminates it from calculations.  The
%  second must be 0 hours, O precipitation depth.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   UNIT CONVERSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
format long

C1 = 10.763910416709722573556258223628;  %m2 to ft2
C2 = 10.763910416709722573556258223628/43560;  %m2 to acre
C3 = 1/2.54;  %cm to in.
C4 = .0008107131937899125401845273267677; %m3 to acre-ft
%
B1 = 1/10.763910416709722573556258223628; %ft2 to m2
B2 = 4046.856422399999900285457752064;  %acre to m2
B3 = 2.54;  %in to cm
B4 = 1233.4818375475199544104; %acre-ft to m3
%Convert inputs from English to Metric
avet = avet*B3; dmin= dmin*B3;Srz=Srz*B3;St=St*B3;Krz=Krz*B3;Kst=Kst*B3;Kcz=Kcz*B3;Urate=Urate*B3;Arg=Arg*B1;TribA=TribA*B2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  PARAMETERS                                                             
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TRIBUTARY BASIN PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=9.8;		% acceleration of gravity [m/s2]
Sd_imp = Sd;           % Depression storage depth of impervious surface [cm]
Si_perv = Sd;          % Interception storage depth of pervious surface [cm]
Si_rg = Sd;            % Interception storage depth of pervious surface [cm]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EVAPORATION PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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]
Cpan = .75;            % {DBA} Pan Evaporation correction coefficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SOIL SENSITIVITY PARAMETERS:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
soilsdata=dlmread('soilsDNR.txt','\t',[0 1 10 4]);
%soilnames=textread('soilslist.txt','%s','delimiter','\n','whitespace','');
%STexture = soilnames(Ssoilnumber,1);
%RTexture = soilnames(Rsoilnumber,1);
NSATsz = soilsdata(Ssoilnumber,1);  %Porosity of control zone
NRESsz = soilsdata(Ssoilnumber,1)-soilsdata(Ssoilnumber,2); %Residual water content, NSAT - NEFF
NSATrz = soilsdata(Rsoilnumber,1);  %Porosity of control zone
NRESrz = soilsdata(Rsoilnumber,1)-soilsdata(Rsoilnumber,2); %Residual water content, NSAT - NEFF
AA = [0.01 0.01 0.01];				% alpha [cm-1]
Ss=[1e-6 1e-6 1e-6];   % specific storage [cm-1]
he=-1./AA;   % air entry suction [cm]
b = 3.1;                         %{DBA} Pore size distribution index, from C&N
rzwilt = 0.1;                       %{DBA} Root zone wilting point, sandy loam from C&N
nwiltfactor=1.21; %{DBA} for estimation of near wilting point
nsfactor=0.95;  % for tracking near saturated conditions in the root zone
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODELING PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEMPORAL:

time=rundays*24*60;			% min in data file
dt=15;						% TIME STEP [min]
stime=time/60;				% total time of simulation [h]
dt=dt/60;					% TIME STEP [h]

Scz = 10;             %Control Zone Thickness
N  = [1.6 1.378 1.242];
h0=[-2.4246e-004 -4.7992e-007 -1.2314e-010];
%NTexture = soilnames(Nsoilnumber,1);
NSATcz = soilsdata(Nsoilnumber,1);  %Porosity of control zone
NREScz = soilsdata(Nsoilnumber,1)-soilsdata(Nsoilnumber,2); %Residual water content, NSAT - NEFF
if St==0    %Allow for no storage by dividing control zone into two equal layers
    St=5;
    Scz=5;
    Kst = Kcz;
    NSATsz = NSATcz;
    NRESsz = NREScz;
    N  = [1.6 1.242 1.242];
    h0=[-2.4246e-004 -1.2314e-010 -1.2314e-010];
end
hwf=3*(2+3*(N(1)-1))/(1+3*(N(1)-1))*he(1)/2;%11;% wetting front suction of upper layer, Rawls&Brakensiek
NSAT = [NSATrz NSATsz NSATcz];				% water content at saturation
NRES = [NRESrz NRESsz NREScz];				% residual water content
KSAT = [Krz Kst Kcz];	% sat. hydraulic conductivity [cm/h]
thick = [Srz St Scz];		% thickness of RZ, TZ, CZ layers [cm]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)); 
%set init. moisture for initial pressure head
    %theta_init(i)=NRES(i);  %For use to begin simulation at field capacity
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sd_perv = 2.54*(1000/CN-10);     % Depression storage depth of pervious surface CN=80 [cm]
Rb = Arg/TribA;  %Ratio of Bioretention to Total Area
Ri = EIA/100/Rb;  %Ratio of impervious area to facility area
Rp = (1-EIA/100)/Rb; %Ratio of pervious area to facility area

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  UNDERDRAIN SIZING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Umax = (Urate/100)*(Arg/3600); %Max. Underdrain flow [m3/sec]
if Umax == 0
    udiam=0;
end
hmax = (dmin + Srz)/100;  %Max. Head on underdrain [m]
Co = 0.62; g=9.81; pi=3.141592654; %constants
udiam = 2*sqrt(Umax/(Co*pi*sqrt(2*g*hmax)));  %underdrain diameter [m]
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  PERVIOUS SURFACE              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
hp_0 = 0;               % Ponded depth at impervious surface initialized [cm]
A_perv = Arg;           % Area of pervious surface (raingarden or lawn area) [m2]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  STORM SIMULATION              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Storm_gap = 6;  % Time with no rain defining gap between storms for impervious areas [h]
t=0;
% dt = 5/60;      % Time step used in model (use divisor of data interval) [h]
% stime=4;        % [h], end of simulation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  SOIL                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% theta_init=NSAT;  % Initial soil water moisture for Green-Ampt
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  INTIIALIZE VARIABLES AND COMPUTE VALUES                                
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  INFILTRATION                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

theta_RZ = theta_init(1);        % Soil moisture of root zone initialized
%theta_RZ = .2;         % Use to match event model
theta_TZ=theta_init(2); theta_CZ=theta_init(3);
curve_area = 0;             % Area under storm event intensity curve initialized [cm]
ponding = 0;                % Ponding indicator initialized to false [0 or 1]
F = 0;                      % Cumulative infiltration depth initialized [cm]
Fprev=0;
ttp = 0;                    % Total time ponded initialized [h]
tp = 0;                     % time of ponding for storm [h]
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  WATER BALANCE                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f = 0;          % Infiltration rate initialized [cm/h]
flim = 0;       % Used when saturated conditions impose limits on f.  [cm/h]
Fcum=0;         % Cumulative infiltration depth [cm]
hp=hp_0;
roff = 0;       % Runoff rate initialized [cm/h]
perc1=0; perc2=0;perc3(:,1)=0;       % Percolation rate initialized [cm/h]
if Record == 1
    WBSum(1,:)=[0 0 0 0 0 0 0 0 0 0 0];                       %Water Budget Summary initialized
end
Inflow(:,1)=0; udrain(:,1)= 0; infilt(:,1)= 0; Runoff(:,1)= 0; 
EvapoT(:,1)= 0; Recharge(:,1)=0;
x= 1; s = 1;
perc1lim=0; perc2lim=0; perc3lim=0; % Used when saturated condtns impose limits on perc. [cm/h]
V_inp = 0;      % Inflow volume initialized [m^3]
V_inf = 0;      % Infiltrated volume initialized [m^3]
V_roff = 0;     % Runoff volume initialized [m^3]
V_rech = 0;     % Recharge volume initialized [m^3]
V_ud = 0;       % Underdrain volume initialized [m^3]
V_et = 0;       % Evapotranspiration volume initialized [m^3]
V_pond = 0;     % Ponding volume initialized [m^3]  %TAKE OUT
tctpmax=0; tstponded=0; tctponded=0;    % initial counters for ponded time
tctoverflow=0;  %initial counters for overflow time
% track near saturated total time (>nsfactor):
tstrznsat=0; tsttznsat=0; tctrznsatmax=0; tcttznsatmax=0; tctrznsat=0; 
tcttznsat=0;
tstrzwilt=0; tctrzwilt=0; tctrzwiltmax=0; %{DBA}
inst=0; nstimes(1,:)=[0 0]; % counter for vector 'nstimes' that will store nsat times for simulation
inpt=0; ptimes(1,:)=[0 0]; % counter for vector 'ptimes' that will store ponded times for simulation
inot=0; oftimes(1,:)=[0 0]; prevof=0; % counter for vector 'oftimes' that will store overflow times for simulation
inwt=0; wtimes(1,:)=[0 0]; % {DBA} counter for vector 'wtimes' that will store wilted times for simulation
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  DATA RETRIEVAL AND PREPROCESSING                                       
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 [cm].  The first row of each
% column should contain the title of the column, as the model will not look
% at that row.
%

qs=nprecip(1,2);    %first value for surface flux [cm h-1]
es=nprecip(1,3);    %first value for surface flux
pdfile = nprecip; %dlmread('lawntest.txt','\t',1,0); % Read data file [h cm]
%
%
len = size(pdfile,1); % Calculate number of data points %PUT BACK IN
if stime>pdfile(len,1);stime=pdfile(len,1);end; % Set max simulation time as final time in data set [h] % PUT BACK IN
%
%
% Linearly interpolate the data set using the model time step.  This will
% produce a new two-column array with values of time and precipitation 
% depth at each multiple of the time step.
%
tk = 0:dt:stime; tk = tk'; % Create a column vector with multiples of dt
intensity = interp1(pdfile(:,1),pdfile(:,2),tk,'*linear'); % Intensity vector
intfile = [tk intensity]; % Interpolated time/intensity array [h cm/h]
Storm_gap = Storm_gap/dt; % Convert storm gap from [h] to # time steps
%
% Take the interpolated time/intensity array and determine which time steps
% correspond to storm events.  This is done using the stormcheck function.
% The function accepts a two-column time/intensity array and the storm gap.
% It adds a third column indicating whether the time step is part of a storm
% event [0 or 1] and a fourth column indicating the length of the event,
% shown at the first time step in an event.
%
datafile = stormcheck(intfile,Storm_gap); % Obtain storm information
Dprecip = sum(nprecip(:,2));
%
% Remove the depression storage in impervious surface from runon:
%
irain = [tk dt*datafile(:,2)*Ri];    % Compute total ppt on impervious area
evap = Evap_imp*dt;                     % Compute the evaporated depth per time step [cm]
%runon = storageout(irain,Sd_imp,evap);  % Remove water stored in impervious
runon1 = storageout(irain,Sd_imp*Ri,evap*Ri,datafile);  % Remove water stored in impervious
runon1(:,2) = runon1(:,2) / dt; 
Vimp = sum(runon1(:,2))* dt;
if Ri == 0;
    Dimp = 0;
else 
    Dimp = Vimp/Ri;
end
%
% Remove the interception storage in pervious surface from rain:
% 
rgrain = [tk dt*datafile(:,2)];          % Compute the total ppt on pervious area
rain = storageout(rgrain,Si_rg,evap,datafile);  % Remove water stored in vegetation
rain (:,2) = rain(:,2) / dt;
% Calc runoff from pervious areas using NRCS Method
Gap_P = 24;  % 24 hour gap for pervious areas
Gap_P = Gap_P/dt; % Convert storm gap from [h] to # time steps
pdatafile = stormcheck(intfile,Gap_P); % Obtain storm information
prain=[tk dt*pdatafile(:,2)];
if Rp == 0;
    runon2 = storageout(prain,0,evap,pdatafile);
else
    runon2 = scsII(prain,CN,evap,pdatafile);
end
runon2(:,2) = Rp*runon2(:,2) / dt;
Vperv = sum(runon2(:,2)) *  dt;
if Rp == 0;
    Dperv=0;
else
    Dperv = Vperv/Rp;
end
Dbasin=(Dprecip-Dimp)*(EIA)+(Dprecip-Dperv)*(1-EIA);
%
% Correct input file
%
runon(:,2) = runon1(:,2) + runon2(:,2);
datafile(:,2)=rain(:,2)+runon(:,2); % Sum effects of impervious & pervious
%
% Create vectors for convenient use of key values
%
inp = datafile(:,2);            % [cm/h]


indstorm = datafile(:,3);       % 0/1
ttstorm = datafile(:,4) * dt;
%
% Create output file for storing results of simulation
% delete and open file for output, print simulation parameters:
% outfile='lawntest_out.txt';

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MAIN SIMULATION LOOP                                                   
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
t=0; n=1;
% This indicates that the following loop begins at the BEGINNING of the first
% time step.
tcpu=cputime;    % starts counting CPU time
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while t < stime % Check that simulation time not passed
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHECK IF STORM STARTS TO GET INTO LOOP:
    if indstorm(n)==1
        %(n==1 & indstorm(n)==1) | (indstorm(n)==1 & indstorm(n-1)==0)
        test2 = indstorm(n); %TAKE OUT
        test3 = indstorm(n-1); %TAKE OUT
        % reset counters before beginning storm loop:
        te = 0; ponding = 0;
        F = 0; curve_area = 0;
        thetainiRZ = theta_RZ; thetainiTZ=theta_TZ; 
thetainiCZ=theta_CZ; % thetaini is initial soil moisture for event
        T = ttstorm(n); % Check length of storm, store as T
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        while (indstorm(n)==1 | ponding==1) & t <= stime   % INFILTRATION/SUPPLY LOOP%%%%%
            %(te < T | ponding==1) & t <= stime   % INFILTRATION/SUPPLY LOOP%%%%%
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %DETERMINE f (INFILTRATION RATE)
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            test = inp(n); %TAKE OUT
            Fprev=F;    % store previous value of storm cumulative infiltration for accounting
            if ponding == 0 % Rainfall infiltration case
                % Keep track of area under water input (infiltration) curve
                if n>1; curve_area = curve_area + (inp(n) + inp(n-1)) / 2 * dt; end; % trapezoidal rule
                f = inp(n); F = F + f * dt;     % update storm infiltration rate & depth
                
                
    
                
                
                if inp(n) > KSAT(1);              % only ponding if i>ksat
                    % Compute F assuming Green-Ampt > 0
 FGA = KSAT(1)*abs(hwf)*(NSAT(1) - thetainiRZ)/(inp(n) - KSAT(1));
                    if curve_area > FGA
                        ponding = 1;
                        tp = te;            % relative time (t is absolute)
                        Fp = curve_area;    % will overestimate tp and Fp
                    end
                end
            end
            if ponding == 1 % Ponded infiltration case
                A = (abs(hwf) + hp) * (NSAT(1) - thetainiRZ);        % Infilt. impetus
                [time,F1]=ode45(@iga,[t,t+dt],F,[],KSAT(1),A);    % Solve for F diff eqn
                F = F1(size(F1,1));                               % Retrieve final value
                f = iga(t,F,KSAT(1),A);                           % Solve for f
            end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Check for limits on infiltration
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %The infil. rate cannot be > input rate plus ponding.  Therefore impose limit below.
            f = min(inp(n)+hp/dt,f);
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %When any combination of the soil layers is saturated, infiltration is limited by
            %the smallest KSAT of the saturated layers.  The limit is assigned in the section 
            %that addresses soil layer saturation and taken into account here.
            if flim >0
                f = min(flim,f);
            end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %DETERMINE PERCOLATION RATES
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            perc1 = KSAT(1) * relcond(theta_RZ,AA(1),N(1),NRES(1),NSAT(1)); 
            perc2 = KSAT(2) * relcond(theta_TZ,AA(2),N(2),NRES(2),NSAT(2));
            perc3(n,1) = KSAT(3) * relcond(theta_CZ,AA(3),N(3),NRES(3),NSAT(3));
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Check for limits on percolation.
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %The percolation rate out of a soil layer cannot be greater than the KSAT of that layer.
            %It does not matter whether the layer is saturated or unsaturated.
            if perc1 > KSAT(1)
                perc1 = KSAT(1);
            end
            if perc2 > KSAT(2)
                perc2 = KSAT(2);
            end
            if perc3(n,1) > KSAT(3)
                perc3(n,1) = KSAT(3);
            end
                
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %When any combination of the soil layers is saturated, percolation rates are
            %limited by the smallest KSAT of the saturated layers.  Limits are assigned
            %in the section that addresses soil layer saturation and taken into account here.
            if perc1lim> 0
                perc1 = min(perc1lim,perc1);
            end
            if perc2lim > 0
                perc2 = min(perc2lim,perc2);
            end
            if perc3lim > 0
                perc3(n,1) = min(perc3lim,perc3(n,1));
            end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Determine Underdrain discharge rate.  Assume underdrain
            % discharge occurs only when Storage zone is saturated.  Assume
            % underdrain limited only by outgoing pipe size (orifice
            % equation) not by perforation in underdrain itself.
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if udiam > 0
                if theta_TZ==NSAT(2)
                    hu = (Srz*(theta_RZ-NRES(1))/(NSAT(1)-NRES(1)))/100;
                    if hp > 0
                        hu = (hp + Srz)/100;
                    end
                else
                    hu = 0;
                end
                uarea = pi*(udiam/2)^2;
                uflowmax = (Co*uarea*(2*g*hu)^.5)/Arg*3600*100;  %underdrain flow [cm/hr]
                uflow = min(uflowmax,perc1-perc2);   
                uflow = max(uflow,0);
            else
                uflow = 0; uflowmax = 0;
            end

            
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %When soil moisture decreases to residual levels, percolation equals zero.  But percolation
            %function makes percolation essentially equal zero.  So don't need to include special
            %limits.
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Update SOIL MOISTURE values and keep old values:
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            theta_RZold=theta_RZ; theta_TZold=theta_TZ; 
theta_CZold=theta_CZ;
            theta_RZ = theta_RZ + (f - perc1) * dt / thick(1); % no et when ponded
            theta_TZ = theta_TZ + (perc1- uflow - perc2) * dt / thick(2);
            theta_CZ = theta_CZ + (perc2 - perc3(n,1)) * dt / thick(3);
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Checks for limits on soil moisture:
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % 1. Lower: 
            %The soil moisture content of the layers cannot be less than the residual soil moisture.
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %if theta_RZ <= NRES(1)
                %deficit = NRES(1) - theta_RZ;
                %theta_RZ = NRES(1);
                %perc1 = perc1 - deficit*thick(1)/dt; 
                %theta_TZ = theta_TZold + (perc1-perc2)*dt/thick(2);
                %end
            if theta_TZ <= NRES(2)
                deficit = NRES(2) - theta_TZ;
                theta_TZ = NRES(2);
                perc2 = perc2 - deficit*thick(2)/dt; 
                theta_CZ = theta_CZold + (perc2-perc3(n,1))*dt/thick(3);
            end
            if theta_CZ <= NRES(3)
                deficit = NRES(3) - theta_CZ;
                theta_CZ = NRES(3);
                perc3(n,1) = perc3(n,1) - deficit*thick(3)/dt;
            end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % 2. Upper:
            %When any combination of the soil layers is saturated, infiltration and percolation
            %are limited by the smallest KSAT of the saturated layers. If the layers become
            %saturated in the current time step, then this section assigns appropriate limits 
            %for the next time step.  
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            flim=0; 
            perc1lim=0; perc2lim=0; perc3lim=0;
            if theta_RZ >= NSAT(1)  %only RZ saturated
                leftover = theta_RZ-NSAT(1);
                xswater = leftover*thick(1); %[cm]
                theta_RZ = NSAT(1);
                f = f - xswater/dt; %[cm/h]
                flim = KSAT(1); %[cm/h]
                perc1lim = KSAT(1); %[cm/h]
                if theta_TZ >= NSAT(2)  %RZ and TZ saturated
                    leftover = theta_TZ-NSAT(2);
                    xswater = leftover*thick(2); %[cm]
                    theta_TZ = NSAT(2);
                    %perc1 = perc1 - xswater/dt;
                    %perc1 used to calc. theta_RZ and theta_TZ but already know those values so
                    %not necessary to keep track of perc1.
                    f = f-xswater/dt; %[cm/h]
                    flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                    perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                    perc2lim = KSAT(2); %[cm/h]
                    if theta_CZ >= NSAT(3) %RZ, TZ, and CZ saturated
                        leftover = theta_CZ - NSAT(3);
                        xswater = leftover*thick(3); %[cm]
                        theta_CZ = NSAT(3);
                        %perc2 = perc2 - xswater/dt; [cm/h]
                        %perc1 = perc1 - xswater/dt; [cm/h]
                        f = f-xswater/dt; %[cm/h]
                        flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                        flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                        perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                        perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h]
                        perc2lim = min(KSAT(2),KSAT(3)); %[cm/h]
                        perc3lim = KSAT(3); %[cm/h]
                    end
                elseif theta_CZ >= NSAT(3)  %RZ and CZ saturated, TZ being filled
                    leftover = theta_CZ-NSAT(3);
                    xswater = leftover*thick(3); %[cm]
                    theta_CZ = NSAT(3);
                    theta_TZ = theta_TZ + xswater/thick(2);                  
                    %perc2 = perc2 - xswater/dt;  %[cm/h] no need to track perc2
                    %flim = KSAT(1); [cm/h]  value assigned above
                    %perc1lim = KSAT(1); [cm/h]  value assigned above
                    perc2lim = KSAT(3); %[cm/h] %
                    perc3lim = KSAT(3); %[cm/h]
                    if theta_TZ >= NSAT(2);  %RZ, TZ, and CZ saturated
                        leftover = theta_TZ-NSAT(2);
                        xswater = leftover*thick(2); %[cm]
                        theta_TZ = NSAT(2);
                        %perc1 = perc1 - xswater/dt; [cm/h] no need to track perc1
                        f = f - xswater/dt; %[cm/h]
                        flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                        flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                        perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                        perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h]
                        perc2lim = min(KSAT(2),KSAT(3)); %[cm/h] 
                    end
                end
            elseif theta_TZ >= NSAT(2)  %only TZ saturated, RZ being filled
                leftover = theta_TZ - NSAT(2);
                xswater = leftover*thick(2);  %[cm]
                theta_TZ = NSAT(2);
                theta_RZ = theta_RZ + xswater/thick(1);
                %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                perc1lim = KSAT(2)+uflowmax; %[cm/h] 
                perc2lim = KSAT(2); %[cm/h]
                if theta_RZ >= NSAT(1)%RZ and TZ saturated
                    leftover = theta_RZ - NSAT(1);
                    xswater = leftover*thick(1); %[cm]
                    theta_RZ = NSAT(1);  
                    f = f - xswater/dt; %[cm/h]
                    flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                    perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h] 
                    %perc2lim = KSAT(2); value assigned above
                elseif theta_CZ >= NSAT(3) %TZ and CZ saturated, RZ being filled
                    leftover = theta_CZ - NSAT(3);
                    xswater = leftover*thick(3); %[cm]
                    theta_CZ = NSAT(3);
                    theta_RZ = theta_RZ + xswater/thick(1);
                    %perc2 = perc2 - xswater/dt; [cm/h] no need to track perc2
                    %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                    perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                    perc2lim = min(KSAT(2),KSAT(3)); %[cm/h]
                    perc3lim = KSAT(3); %[cm/h]
                    if theta_RZ >=NSAT(1) %RZ,TZ, and CZ saturated
                        leftover = theta_RZ - NSAT(1);
                        xswater = leftover*thick(1); %[cm]
                        theta_RZ=NSAT(1);  
                        f = f - xswater/dt; %[cm/h] 
                        flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                        flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                        perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]  
                        perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h] 
                        %perc2lim = min(KSAT(2),KSAT(3)); % [cm/h]  value assigned above
                    end    
                end    
            elseif theta_CZ >= NSAT(3)  %only CZ saturated, TZ being filled
                leftover = theta_CZ - NSAT(3);
                xswater = leftover*thick(3);  %[cm]
                theta_CZ = NSAT(3);
                theta_TZ = theta_TZ + xswater/thick(2);               
                %perc2 = perc2 - xswater/dt; %[cm/h] no need to track perc2
                perc2lim = KSAT(3); %[cm/h]  
                perc3lim = KSAT(3); %[cm/h]
                if theta_TZ >=NSAT(2) %TZ and CZ saturated, RZ being filled
                    leftover = theta_TZ - NSAT(2); 
                    xswater = leftover*thick(2); %[cm]
                    theta_TZ = NSAT(2);  
                    theta_RZ = theta_RZ + xswater/thick(1);                  
                    %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                    perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                    perc2lim = min(KSAT(2),KSAT(3)); %[cm/h] 
                    if theta_RZ >= NSAT(1) %RZ, TZ, and CZ saturated
                        leftover = theta_RZ - NSAT(1);
                        xswater = leftover*thick(1); %[cm]
                        theta_RZ = NSAT(1);  
                        f = f - xswater/dt; %[cm/h]
                        flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                        flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                        perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h] 
                        perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h] 
                        %perc2lim = min(KSAT(2),KSAT(3)); [cm/h]  value assigned above
                    end
                end
            end
            
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %SURFACE WATER MASS BALANCE CALCULATIONS:
            hpold=hp;   % keep previous value of surface ponding depth
            hp = hp + dt * (inp(n) - f); % Potential (Note: no evaporation)
            
            if hp <= 0   % check for no ponding, drainage in the time step
                hp = 0;
                ponding=0;  % end of ponding as hp goes to zero
            else
                ponding=1;  % because hp>0
            end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if hp>dmin   % means water can spill from depression storage
                roff = (hp-dmin)/dt;
                hp = hp - dt * roff;
            else
                roff = 0;
            end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %if hp <= 0%Always 0!!!
            %    ttp = ttp + (te-tp);
            %    %roff = 0;
            %end;%             fprintf('V_inp= %f\n',V_inp); %TAKE OUT%             fprintf('hp= %f\n',hp); %TAKE OUT%                           fprintf('time= %f\n',t); %TAKE OUT
            %            fprintf('n= %f\n',n); % TAKE OUT%                         fprintf('Theta_RZ= %f\n',theta_RZ); %[TAKE THIS OUT]
            %             fprintf('Theta_TZ= %f\n',theta_TZ); %[TAKE THIS OUT]
            %fprintf('Theta_CZ= %f\n',theta_CZ); %[TAKE THIS OUT]
            %             fprintf('f= %f\n',f); %TAKE OUT
            %             fprintf('inp(n)= %f\n',inp(n)); %TAKE OUT
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % UPDATE CUMULATIVE WATER BUDGET TERMS:
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            V_inp = V_inp + inp(n) / 100 * A_perv * dt;
            V_roff = V_roff + roff / 100 * A_perv * dt;     % Runoff volume [m^3]
            V_inf = V_inf + f / 100 * A_perv * dt;          % Infiltrated volume [m^3]
            V_ud = V_ud + uflow / 100 * A_perv * dt;              % Underdrain volume [m^3]
            %V_rech = V_rech + perc3(n,1) / 100 * A_perv * dt;%             V_pond = V_pond + hp/ 100 * A_perv * dt;  %TAKE OUT
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Keep track of saturated soil and ponding indicators:
            %
            % Update values for average water content in root zone and transmission zone:
            if theta_RZold > nsfactor*NSAT(1)    % near saturated in past time step, either the beginning or intermediate
                if theta_RZ > nsfactor*NSAT(1)   % near sat. now too, so keep adding:
                    tstrznsat=tstrznsat+dt; % track total simulation near sat. time
                    tctrznsat=tctrznsat+dt; % continuous near sat. time
                else      % means check if maximum value in simulation so far:
                    if tctrznsat>tctrznsatmax; tctrznsatmax=tctrznsat; 
end; % track max cont. near sat. time
                    inst=inst+1; nstimes(inst,:)=[t tctrznsat];      % keep all times stored in vector [t_end tct]
                    tctrznsat=0;  % reset event counter
                end
            end
            %  {DBA}
            if theta_RZold < nwiltfactor*rzwilt
                if theta_RZ < nwiltfactor*rzwilt
                    tstrzwilt=tstrzwilt+dt; % track total simulation at wilting point
                    tctrzwilt=tctrzwilt+dt; % continuous wilting time
                else % check if max value in simulation
                    if tctrzwilt>tctrzwiltmax; tctrzwiltmax=tctrzwilt; 
end; % track max wilting time
                    inwt=inwt+1; wtimes(inwt,:)=[t tctrzwilt];    % keep track of all wilting times
                    tctrzwilt=0; %  reset counter
                end
            end
            % end {DBA} edits
            if theta_TZold > nsfactor*NSAT(2)    % near saturated in past time step, either the beginning or intermediate
                if theta_TZ > nsfactor*NSAT(2)   % near sat. now too, so keep adding:
                    tsttznsat=tsttznsat+dt; % track total simulation near sat. time
                    tcttznsat=tcttznsat+dt; % continuous near sat. time
                else      % means check if maximum value in simulation so far:
                    if tcttznsat>tcttznsatmax; tcttznsatmax=tcttznsat; 
end; % track max cont. near sat. time
                    tcttznsat=0;  % reset event counter
                end
            end
            % Update values for ponding times:
            if hpold >= dmin  % Track overflow events
                if hp >= dmin
                    tctoverflow=tctoverflow+dt;
                else
                    if t > prevof+24; inot=inot+1; oftimes(inot,:)=[t tctoverflow]; tctoverflow=0; prevof=t;end; % keep all times stored in vector
                end
            end  % Track overflow events
            if hpold > 0    % ponded in past time step, either the beginning or intermediate
                if hp > 0   % ponded now too, so keep adding:
                    tstponded=tstponded+dt; % track total ponded time
                    tctponded=tctponded+dt; % continuous ponded time
                else      % means ponding event ended, so check if maximum value in simulation so far:
                    if tctponded>tctpmax; tctpmax=tctponded; end; % track max cont. ponded t
                    inpt=inpt+1; ptimes(inpt,:)=[t tctponded];      % keep all times stored in vector [t_end tct]
                    tctponded=0;  % reset event counter
                end
            end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %RECORD WATER BALANCE TERMS
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Inflow(s,1)=inp(n)/4;
        udrain(s,1)= uflow/4;
        infilt(s,1)= f/4;
        Runoff(s,1)= roff/4;
        EvapoT(s,1)= 0;
        Recharge(s,1)=perc3(n,1)/4;
        s=s+1;
        if Record ==1 
            if mod(t,1)==0
                WBSum(x,:)=[t sum(Inflow(:,1)) hp sum(infilt(:,1)) sum(Runoff(:,1)) sum(udrain(:,1)) sum(Recharge(:,1)) sum(EvapoT(:,1)) theta_RZ theta_TZ theta_CZ];
                Inflow(:,1)=0;
                udrain(:,1)= 0;
                infilt(:,1)= 0;
                Runoff(:,1)= 0;
                EvapoT(:,1)= 0;
                Recharge(:,1)=0;
                x=x+1;
                s=1;
            end
        end

        
            te = te + dt;
            Fcum=Fcum+(F-Fprev);
            %qst(n,:)=[t Fcum f hp theta_RZ theta_TZ theta_CZ perc1 perc2 perc3];   % value to record
            t = t + dt;
            %set(handles.Progress,'String',num2str(Pcompl));
            n = n + 1;
            if t>stime, break, end
               
            %fprintf('n= %-3.0f| ',n);fprintf('In Storm\n',t);  % {DBA} screen print, omit if want to speed up program execution
        end % OF STORM/SUPPLY LOOP
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    else    % FREE DRAINAGE CASE (BETWEEN STORMS):
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        f = 0;      % no infiltration during hiatus periods
        roff = 0;   % no runoff either
        test = inp(n); %TAKE OUT
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %DETERMINE PERCOLATION RATES
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        perc1 = KSAT(1) * relcond(theta_RZ,AA(1),N(1),NRES(1),NSAT(1)); 
        perc2 = KSAT(2) * relcond(theta_TZ,AA(2),N(2),NRES(2),NSAT(2));
        perc3(n,1) = KSAT(3) * relcond(theta_CZ,AA(3),N(3),NRES(3),NSAT(3));
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Check for limits on percolation
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
        %The percolation rate out of a soil layer cannot be greater than the KSAT of that layer.
        %It does not matter whether the layer is saturated or unsaturated.
        if perc1 > KSAT(1)
            perc1 = KSAT(1);
        end
        if perc2 > KSAT(2)
            perc2 = KSAT(2);
        end
        if perc3(n,1) > KSAT(3)
            perc3(n,1) = KSAT(3);
        end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Determine Underdrain discharge rate.  Assume underdrain
        % discharge occurs only when Storage zone is saturated.  Assume
        % underdrain limited only by outgoing pipe size (orifice
        % equation) not by perforation in underdrain itself.
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if udiam > 0
            if theta_TZ==NSAT(2)
                hu = (Srz*(theta_RZ-NRES(1))/(NSAT(1)-NRES(1)))/100;
                if hp > 0
                    hu = (hp + Srz)/100;
                end
            else
                hu = 0;
            end
            uarea = pi*(udiam/2)^2;
            uflowmax = (Co*uarea*(2*g*hu)^.5/Arg)*3600*100;  %underdrain flow [cm/hr]
        else
            uflowmax = 0;
        end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %When any combination of the soil layers is saturated, percolation rates are
        %limited by the smallest KSAT of the saturated layers.  The limits are assigned
        %in the section that addresses soil layer saturation and taken into account here.
        if perc1lim> 0
            perc1 = min(perc1lim,perc1);
            perc1 = max(perc1,0);
        end
        if perc2lim > 0
            perc2 = min(perc2lim,perc2);
            perc2 = max(perc2,0);
        end
        if perc3lim > 0
            perc3(n,1) = min(perc3lim,perc3(n,1));
            perc3(n,1) = max(perc3(n,1),0);
        end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Limit uflowmax to min. of actual flow calculated per orifice equation
        %and max. available from percolation
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        uflow = min(uflowmax,perc1-perc2);
        uflow = max(uflow,0);
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %When soil moisture decreases to residual levels, percolation equals zero.  But percolation
        %function makes percolation essentially equal zero.  So don't need to include special
        %limits.
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Set EVAPORATION RATE
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if theta_RZ <= rzwilt,
            aw = 0;
        elseif theta_RZ >= NRES(1),
            aw = 1;
        else
            aw = (theta_RZ - rzwilt)/(NRES(1)-rzwilt);
        end
        Pevap = interp1(pdfile(:,1),pdfile(:,3),t,'*linear'); %{DBA} Pot. Evaporation
        et = (1-(1+1.3*aw).^-b)*Pevap;  
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Update SOIL MOISTURE values
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        theta_RZold=theta_RZ; theta_TZold=theta_TZ; theta_CZold=theta_CZ;
        theta_RZ = theta_RZ + (-perc1 - et)*dt/thick(1);
        theta_TZ = theta_TZ + (perc1 - uflow - perc2) * dt / thick(2);
        theta_CZ = theta_CZ + (perc2 - perc3(n,1)) * dt / thick(3);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Checks for limits on soil moisture:
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 2 Lower limits:
        %The soil moisture content of the layers cannot be less than the residual soil moisture.
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %if theta_RZ <= NRES(1)
            %deficit = NRES(1) - theta_RZ;
            %theta_RZ = NRES(1);
            %perc1 = perc1 - deficit*thick(1)/dt; 
            %perc1 = max(perc1,0);
            %theta_TZ = theta_TZold + (perc1-perc2)*dt/thick(2);
            %end
        if theta_TZ <= NRES(2)
            deficit = NRES(2) - theta_TZ;
            theta_TZ = NRES(2);
            perc2 = perc2 - deficit*thick(2)/dt;
            perc2 = max(perc2,0);
            theta_CZ = theta_CZold + (perc2-perc3(n,1))*dt/thick(3);
        end
        if theta_CZ <= NRES(3)
            deficit = NRES(3) - theta_CZ;
            theta_CZ = NRES(3);
            perc3(n,1) = perc3(n,1) - deficit*thick(3)/dt;
            perc3(n,1) = max(perc3(n,1),0);
        end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 2. Upper limits:
        %When any combination of the soil layers is saturated, infiltration and percolation
        %are limited by the smallest KSAT of the saturated layers. If the layers become
        %saturated in the current time step, then this section assigns appropriate limits 
        %for the next time step.  Since free drainage, and infiltration equals zero, do not
        %consider the case of the root zone being saturated.
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        perc1lim=0; perc2lim=0; perc3lim=0;
        if theta_TZ >= NSAT(2)  % TZ saturated, RZ being filled
            leftover = theta_TZ-NSAT(2);
            xswater = leftover*thick(2); %[cm]
            theta_TZ = NSAT(2);
            theta_RZ = theta_RZ + xswater/thick(1);        
            %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
            perc1lim = KSAT(2)+uflowmax; %[cm/h] 
            perc2lim = KSAT(2); %[cm/h]%             if theta_RZ >=NSAT(1) %RZ and TZ saturated%                 theta_RZ = NSAT(1);%                 perc1lim = min(KSAT(1),KSAT(2)); %[cm/h] 
%                 %perc2lim = KSAT(2); [cm/h]  value assigned above
%             elseif theta_CZ >= NSAT(3) %TZ and CZ saturated, RZ being filled
            if theta_CZ >= NSAT(3) %TZ and CZ saturated, RZ being filled
                leftover = theta_CZ - NSAT(3);
                xswater = leftover*thick(3); %[cm]
                theta_CZ = NSAT(3);
                theta_RZ = theta_RZ + xswater/thick(1);
                %perc2 = perc2 - xswater/dt; [cm/h] no need to track perc2
                %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc2
                perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                perc2lim = min(KSAT(2),KSAT(3)); %[cm/h]
                perc3lim = KSAT(3); %[cm/h]%                 if theta_RZ >=NSAT(1)%                     theta_RZ = NSAT(1);  %RZ, TZ, CZ saturated
%                     perc1lim = min(KSAT(1),KSAT(2)); %[cm/h] 
%                     perc1lim = min(KSAT(3),perc1lim); %[cm/h] 
%                     %perc2lim = min(KSAT(2),KSAT(3)); [cm/h]  value assigned above
%                 end
            end    
        elseif theta_CZ >= NSAT(3)  %only CZ saturated, TZ being filled
            leftover = theta_CZ - NSAT(3);
            xswater = leftover*thick(3); %[cm]
            theta_CZ = NSAT(3);
            theta_TZ = theta_TZ + xswater/thick(2); %WHAT IF THICK(2)=0??????????????????????
            %perc2 = perc2 - leftover*thick(3)/dt; %[cm/h] no need to track perc2
            perc2lim = KSAT(3); %[cm/h] 
            perc3lim = KSAT(3); %[cm/h]
            if theta_TZ >= NSAT(2) %TZ and CZ saturated, RZ being filled
                leftover = theta_TZ-NSAT(2);
                xswater = leftover*thick(2); %[cm]
                theta_TZ = NSAT(2);  
                theta_RZ = theta_RZ + xswater/thick(1);
                %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                perc2lim = min(KSAT(2),KSAT(3)); %[cm/h] 
                %perc3lim = KSAT(3); [cm/h]  value assigned above
%                 if theta_RZ >= NSAT(1) %RZ, TZ, CZ saturated
%                     theta_RZ = NSAT(1);
%                     perc1lim = min(KSAT(1),KSAT(2)); %[cm/h] 
%                     perc1lim = min(KSAT(3),perc1lim); %[cm/h] 
%                     % perc2lim = min(KSAT(2),KSAT(3)); [cm/h]  value assigned above
%                     %perc3lim = KSAT(3); [cm/h]  value assigned above
%                 end
            end
        end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % UPDATE CUMULATIVE WATER BUDGET TERMS:
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        V_inp = V_inp + inp(n) / 100 * A_perv * dt;
        V_roff = V_roff + roff / 100 * A_perv * dt; % Runoff volume [m^3]
        V_inf = V_inf + f / 100 * A_perv * dt; % Only exfiltration in this case [m^3]
        %V_rech = V_rech + perc3(n,1) / 100 * A_perv * dt;
        V_et = V_et + et / 100 * A_perv * dt; %TAKE OUT
        V_ud = V_ud + uflow / 100 * A_perv * dt;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Keep track of saturated soil and ponding indicators:
        %
        % Update values for average water content in root zone and transmission zone:
        if theta_RZold > nsfactor*NSAT(1)    % near saturated in past time step, either the beginning or intermediate
            if theta_RZ > nsfactor*NSAT(1)   % near sat. now too, so keep adding:
                tstrznsat=tstrznsat+dt; % track total simulation near sat. time
                tctrznsat=tctrznsat+dt; % continuous near sat. time
            else      % means check if maximum value in simulation so far:
                if tctrznsat>tctrznsatmax; tctrznsatmax=tctrznsat; end; 
% track max cont. near sat. time
                inst=inst+1; nstimes(inst,:)=[t tctrznsat];      % keep all times stored in vector [t_end tct]
                tctrznsat=0;  % reset event counter
            end
        end
        %  {DBA}
        if theta_RZold < nwiltfactor*rzwilt
            if theta_RZ < nwiltfactor*rzwilt
                tstrzwilt=tstrzwilt+dt; % track total simulation at wilting point
                tctrzwilt=tctrzwilt+dt; % continuous wilting time
            else % check if max value in simulation
                if tctrzwilt>tctrzwiltmax; tctrzwiltmax=tctrzwilt; end; 
% track max wilting time
                inwt=inwt+1; wtimes(inwt,:)=[t tctrzwilt];    % keep track of all wilting times
                tctrzwilt=0; %  reset counter
            end
        end
        % end {DBA} edits
        if theta_TZold > nsfactor*NSAT(2)    % near saturated in past time step, either the beginning or intermediate
            if theta_TZ > nsfactor*NSAT(2)   % near sat. now too, so keep adding:
                tsttznsat=tsttznsat+dt; % track total simulation near sat. time
                tcttznsat=tcttznsat+dt; % continuous near sat. time
            else      % means check if maximum value in simulation so far:
                if tcttznsat>tcttznsatmax; tcttznsatmax=tcttznsat; end; 
% track max cont. near sat. time
                tcttznsat=0;  % reset event counter
            end
        end
        %         V_et = V_et + ET_avg / 100 * A_perv * dt;
        %         Stg = V_inf-V_rech-V_et; %TAKE OUT  %Change in water stored in the soil [m^3]
        DSSW1=(theta_RZ - theta_init(1))*thick(1)/100*A_perv; %TAKE OUT
        DSSW2=(theta_TZ - theta_init(2))*thick(2)/100*A_perv; %TAKE OUT
        DSSW3=(theta_CZ - theta_init(3))*thick(3)/100*A_perv; %TAKE OUT
        %         DSSW=DSSW1+DSSW2+DSSW3; %TAKE OUT
        %         if Stg < DSSW %TAKE OUT
        %             fprintf('n= %e\n', n); %TAKE OUT
        % %         fprintf('Stg = %e\n',Stg); %TAKE OUT
        %             fprintf('DSSW= %e\n',DSSW); %TAKE OUT
        %         elseif Stg > DSSW %TAKE OUT
        %             fprintf('n= %e\n', n); %TAKE OUT
        % %         fprintf('Stg = %e\n',Stg); %TAKE OUT
        %             fprintf('DSSW= %e\n',DSSW); %TAKE OUT
        %         end %TAKE OUT
        %         fprintf(' %4.0f',n);
%                           fprintf('time= %f\n',t); %TAKE OUT
        %         fprintf('  %3.3f',theta_RZ);
        %         fprintf(' %e',V_inp-V_roff);
        %         fprintf('  %4.2f\n',hp);
%         qst(n,:)=[t Fcum f hp theta_RZ theta_TZ theta_CZ perc1 perc2 perc3];   % value to record
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %RECORD WATER BALANCE TERMS
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Inflow(s,1)=inp(n)/4;
        udrain(s,1)= uflow/4;
        infilt(s,1)= f/4;
        Runoff(s,1)= roff/4;
        EvapoT(s,1)= et/4;
        Recharge(s,1)=perc3(n,1)/4;
        s=s+1;
        if Record == 1
            if mod(t,1)==0
                WBSum(x,:)=[t sum(Inflow(:,1)) hp sum(infilt(:,1)) 
sum(Runoff(:,1)) sum(udrain(:,1)) sum(Recharge(:,1)) sum(EvapoT(:,1)) 
theta_RZ theta_TZ theta_CZ];
                Inflow(:,1)=0;
                udrain(:,1)= 0;
                infilt(:,1)= 0;
                Runoff(:,1)= 0;
                EvapoT(:,1)= 0;
                Recharge(:,1)=0;
                x=x+1;
                s=1;
            end
        end

        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        t = t + dt;
        %fprintf('n= %-3.0f| ',n);fprintf('free draining\n',t);  % {DBA} screen print, omit if want to speed up program execution
       Pcompl = 100*n/(4*stime);
        %set(handles.Progress,'String',num2str(Pcompl));
        n = n + 1;
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end     % END OF SIMULATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHECK FOR LAST VALUES OF NEAR SATURATION AND PONDING TIMES:
if tctrznsat>tctrznsatmax; tctrznsatmax=tctrznsat; end; % track max cont. near sat. time
if sum(nstimes(:,2))<tstrznsat; inst=inst+1; nstimes(inst,:)=[t tctrznsat]; end;    % keep all times stored in vector
if tcttznsat>tcttznsatmax; tcttznsatmax=tcttznsat; end; % track max cont. near sat. time
if tctponded>tctpmax; tctpmax=tctponded; end; % track max cont. ponded t
if sum(nstimes(:,2))<tstrznsat; inst=inst+1; nstimes(inst,:)=[t tctrznsat]; end;    % keep all times stored in vector
%CHECK MASS BALANCE:
DSSW1=(theta_RZ - theta_init(1))*thick(1)/100*A_perv;
DSSW2=(theta_TZ - theta_init(2))*thick(2)/100*A_perv;
DSSW3=(theta_CZ - theta_init(3))*thick(3)/100*A_perv;
DSSW=DSSW1+DSSW2+DSSW3;
unitSM=DSSW/A_perv;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaporated water from Priestley & Taylor:
%if etrun==1
%   Qevap=Arg*tau*ep(n)/100;   % potential soil evaporation
%else
%   Qevap=0;
%end
% DO SOMETHING TO WARN IF BEGINNING WITH SUPERSATURATED CONDITIONS
% COULD PUT SOMETHING IN THE STORM TIME LOOP TO DETECT WHETHER THERE
%    IS A TIME DURING EVENT WHEN PONDING = 0, AND IF SO, TO START OVER AND 
%    MAKE THE STORM GAP SMALLER.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  SUMMARIZE OUTPUT DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dVsfce=(hp-hp_0)/100*A_perv; % sfce storage change
Vinsoil=V_inp-V_roff-dVsfce;
infdiff=(V_inf-Vinsoil)/Vinsoil*100;
Vr=V_inf-V_et-DSSW-V_ud;   % check on Vrecharge by mass balance
V_rech = sum(perc3(:,1))/100*A_perv*dt; %{DBA}
gwrdiff=(V_rech-Vr)/Vr*100;
gwravgdepth=0.5*(V_rech+Vr)/((1+Ri+Rp)*A_perv)*100; % [cm]
unitgwr=0.5*(V_rech+Vr)/A_perv; %{DBA} [m3/m2] recharge per unit RGarea%unitgwr=V_rech/A_perv;
dVsoil=V_inf-V_rech-V_et;
unitet=V_et/A_perv; %{DBA} [m3/m2] evaporation per unit RG area
unitron=V_inp/A_perv; %{DBA} [m3/m2] runon per unit RG area
unitroff=V_roff/A_perv; %{DBA} [m3/m2] runoff per unit RG area
unitud=V_ud/A_perv; %{DBA} [m3/m2] underdrain vol per unit RG area
Drunon=(100*Arg*unitron)/(Arg+TribA);
Drunoff=(100*Arg*unitroff)/(Arg+TribA);
Drecharge=(100*Arg*unitgwr)/(Arg+TribA);
Devap=(100*Arg*unitet)/(Arg+TribA);
Dudrain=(100*Arg*unitud)/(Arg+TribA);
Dstore=(100*Arg*unitSM)/(Arg+TribA);
dVs=V_inp-V_roff-V_inf;      % check on sfce storage change by mass balance
ttcpu=cputime-tcpu;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  CREATE OUTPUT FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Summary == 1;
    outfile=strcat(outpre,'.txt');
    fid=fopen(outfile,'w'); fclose(fid);	% delete previous version & open
	fid=fopen(outfile,'at');
	fprintf(fid,'RECARGA 2.X\n');
	fprintf(fid,'Date:\t %s\n',datestr(clock));
	fprintf(fid,'Output File:\t %s\n',outfile);
	fprintf(fid,'Number of time steps=\t %-4.0f\n',n);
	fprintf(fid,'CPU elapsed time (s)=\t %-10.0f\n\n',ttcpu);
	fprintf(fid,'INPUT TERMS\n');
	fprintf(fid,'%%Impervious=\t %g\n',EIA);
	fprintf(fid,'Pervious CN=\t %g\n',CN);
	fprintf(fid,'Facility Area (sf)=\t %g\n',Arg*C1);
	fprintf(fid,'Trib. Area (AC)=\t %g\n',TribA*C2);
	fprintf(fid,'RATIOimp2rg=\t %g\n',Ri);
	fprintf(fid,'RATIOperv2rg=\t %g\n\n',Rp);
	fprintf(fid,'LAYER       \t DEPTH (in)  \t Ksat (in/hr)\t TEXTURE     \n');
	fprintf(fid,'Depression  \t     %g      \t     ---      \t   ---       \n',dmin*C3);
	fprintf(fid,'Root Zone   \t     %g      \t',thick(1)*C3);
	fprintf(fid,'    %4.2f    \t',KSAT(1)*C3);
	fprintf(fid,'    %g      \n',Rsoilnumber);
	fprintf(fid,'Storage     \t     %g      \t',thick(2)*C3);
	fprintf(fid,'    %4.2f     \t',KSAT(2)*C3);
	fprintf(fid,'    %g      \n',Ssoilnumber);
	fprintf(fid,'Native      \t     ---     \t     %4.2f   \t',KSAT(3)*C3);
	fprintf(fid,'    %g      \n',Nsoilnumber);
	fprintf(fid,'\nInit. Cond.: suction head (in)=\t %5.0f\n',hin*C3);
	fprintf(fid,'Max. Underdrain flow (in/hr)=\t %5.3f\n',Umax*C3);
	fprintf(fid,'Underdrain diam (in)=\t %5.3f\n',udiam*100*C3);
	
	fprintf(fid,'\nWATER BALANCE TERMS\n\n');
	fprintf(fid,'           \t Volume(ac-ft) \t Depth (in)\t %%of Inflow \n');
	fprintf(fid,'Inflow     \t %10.2f \t',V_inp*C4);
	fprintf(fid,' %9.2f \t',Drunon*C3);
	fprintf(fid,' %9.2f \n',100*Drunon/Drunon);
	fprintf(fid,'Runoff     \t %10.2f \t',V_roff*C4);
	fprintf(fid,' %9.2f \t',Drunoff*C3);
	fprintf(fid,' %9.2f \n',100*Drunoff/Drunon);
	fprintf(fid,'Recharge   \t %10.2f \t',V_rech*C4);
	fprintf(fid,' %9.2f \t',Drecharge*C3);
	fprintf(fid,' %9.2f \n',100*Drecharge/Drunon);
	fprintf(fid,'Evaporation\t %10.2f \t',V_et*C4);
	fprintf(fid,' %9.2f \t',Devap*C3);
	fprintf(fid,' %9.2f \n',100*Devap/Drunon);
	fprintf(fid,'Underdrain \t %10.2f \t',V_ud*C4);
	fprintf(fid,' %9.2f \t',Dudrain*C3);
	fprintf(fid,' %9.2f \n',100*Dudrain/Drunon);
	fprintf(fid,'Storage    \t %10.2f \t',DSSW*C4);
	fprintf(fid,' %9.2f \t',Dstore*C3);
	fprintf(fid,' %9.2f \n',100*Dstore/Drunon);
	fprintf(fid,'Infiltrated\t %10.2f \n',V_inf*C4);
	fprintf(fid,'Precip.    \t            \t %9.2f\n',Dprecip*C3);
	fprintf(fid,'Imp. Inflow\t %10.2f \t',Vimp*C4);
	fprintf(fid,' %9.2f \t',Dimp*C3);
	fprintf(fid,' %9.2f \n',100*Vimp/V_inp);
	fprintf(fid,'Perv.Inflow\t %10.2f \t',Vperv*C4);
	fprintf(fid,' %9.2f \t',Dperv*C3);
	fprintf(fid,' %9.2f \n',100*Vperv/V_inp);
	fprintf(fid,'Basin Stayon(in)\t\t %9.2f\n',Dbasin*C3);
    fprintf(fid,'Total Stayon(in)\t\t %9.2f\n\n',(Dprecip-Dudrain-Drunoff)*C3);
	fprintf(fid,'Mass Balance Checks:\n');
	fprintf(fid,'Vinsoil=input-roff-dvh\n');
	fprintf(fid,'Vinsoil (ac-ft)=\t %e\n',Vinsoil*C4);
	fprintf(fid,'Vinf difference (%%)=\t %4.1f\n',infdiff);
	fprintf(fid,'Vrecharg1=inf-et-dsws\n');
	fprintf(fid,'Vrecharg1 (ac-ft)=\t %e\n',Vr*C4);
	fprintf(fid,'Vrech difference (%%)=\t %4.1f\n',gwrdiff);
	fprintf(fid,'dVsoil=inf-rech-et\n');
	fprintf(fid,'dVsoil=\t %e\n',dVsoil*C4);
	fprintf(fid,'dVsfce (ac-ft)=\t %e\n',dVsfce*C4);
	fprintf(fid,'dVs2=input-roff-inf\n');
	fprintf(fid,'dVs2 (ac-ft)= %e\n',dVs*C4);
	%
	fprintf(fid,'\nPLANT SURVIVABILITY TERMS\n');
	fprintf(fid,'total time ponded (h)=\t %5.1f\n',tstponded);
	fprintf(fid,'max time ponded (h)=\t %5.1f\n',tctpmax);
	fprintf(fid,'total time RZ saturated (h)=\t %5.1f\n',tstrznsat);
	fprintf(fid,'total time RZ at wilting point (h)=\t %5.1f\n', tstrzwilt);  % {DBA}
	fprintf(fid,'total time TZ saturated (h)=\t %5.1f\n',tsttznsat);
	fprintf(fid,'max time RZ saturated (h)=\t %5.1f\n',tctrznsatmax);
	fprintf(fid,'max time RZ at wilting point (h)=\t %5.1f\n',tctrzwiltmax);  %{DBA}
	fprintf(fid,'max time TZ saturated (h)=\t %5.1f\n',tcttznsatmax);
	
	fprintf (fid,'near saturated times for simulation=\t %6.0f\n', inst); 
fprintf(fid,'%6.1f\t %6.1f\n',nstimes');
	fprintf (fid,'wilting point times for simulation=\t %6.0f\n',inwt); 
fprintf(fid,'%6.1f\t %6.1f\n',wtimes'); % {DBA}
	fprintf (fid,'ponded times for simulation=\t %6.0f\n', inpt); 
fprintf(fid,'%6.1f\t %6.1f\n',ptimes');
	fprintf (fid,'overflow times for simulation=\t %6.0f\n',inot); 
fprintf(fid,'%6.1f\t %6.1f\n',oftimes');
	fclose(fid);
end
if Record == 1;
    WBSum(:,2)=WBSum(:,2)*C3; WBSum(:,3)=WBSum(:,3)*C3; 
WBSum(:,4)=WBSum(:,4)*C3; WBSum(:,5)=WBSum(:,5)*C3;
    
WBSum(:,6)=WBSum(:,6)*C3;WBSum(:,7)=WBSum(:,7)*C3;WBSum(:,8)=WBSum(:,8)*C3;
    WBfile = strcat(outpre,'_WB.txt');    
	fid=fopen(WBfile,'w'); fclose(fid);	% delete previous version & open
	fid=fopen(WBfile,'at');
	fprintf(fid,'Time(hr)\t Runon(in)\t Ponding(in)\t Infil(in)\t Runoff(in)\t Drain(in)\t Recharge(in)\t ET(in)\t ThetaRZ\t ThetaSZ\t ThetaCZ\n');
	fprintf(fid,'%6.3f\t %6.3f\t 6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %.3f\n',WBSum');
	fclose(fid);
end

%  Convert output to English
Dprecip=Dprecip*C3;Dimp=Dimp*C3;Dperv=Dperv*C3;Drunon=Drunon*C3;Drunoff=Drunoff*C3;
Devap=Devap*C3;Drecharge=Drecharge*C3;Dudrain=Dudrain*C3;Dstore =Dstore*C3;
end
end