% This script loads existing output from the 1D estuary model 
% and calculates residence times and fluxes according to different methods.

% Clear workspace first.
clear

spday = 24*3600;   % seconds per day
sphr  = 3600;      % seconds per hour

% Ask which case to load.
disp(' ')
disp('Which case do you want to load?')
disp(' ')
filelist = dir('1Dsim*.mat');
for i=1:length(filelist)
    disp([num2str(i) ') ' filelist(i).name]) 
end
disp(' ')
fno = input('Enter number: ');
disp(' ')

% Load chosen file
load(filelist(fno).name)

% Get the name of the tracer
if ~isempty(strfind(filelist(fno).name,'DIC')), xname = 'DIC'; end
if ~isempty(strfind(filelist(fno).name,'SF6')), xname = 'SF6'; end
if ~isempty(strfind(filelist(fno).name,'dye')), xname = 'dye'; end

% Calculate basic parameters
dx = sim.x(2) - sim.x(1);
dt = sim.t(2) - sim.t(1);

% Check if this is the small (50 km) or large (100 km) domain
lg_domain = ~isempty(strfind(filelist(fno).name,'100km')); 

% Define spatial indices for the estuary
if lg_domain
    % find spatial indices for 0 to 80 km
    ind_estuary = find(sim.x<=80000);
else
    % find spatial indices for 15 to 30 km
    %ind_estuary = find((sim.x>=15000) & (sim.x<=30000));
    % find spatial indices for 0 to 30 km
    ind_estuary = find(sim.x<=30000);
end


% Calculate residence times:

% 1. water volume / discharge
% 0 to 30 km (or 0 to 80 km for the large domain)

if lg_domain
    R1 = 80000 / sim.net_discharge;
else
    R1 = 30000 / sim.net_discharge;
end
disp([' Residence time 1 (volume / discharge): ' num2str(R1/spday) ' d'])
disp(' ')

% 2. decay rate
    
% time series of tracer inventory (over ind_estuary) from daily transects
% first, calculate daily inventory in region of interest 
% variable: sim.tracer.dist, sum along the 2nd dimension
x_inv = sum(sim.tracer.dist(:,ind_estuary)*dx,2);
time = sim.tracer.dist_time/spday; % time in days

% plot daily inventory on linear and log scales
figure
plot(time,x_inv,'.-','Linewidth',1,'MarkerSize',6)
xlabel('time (d)')
ylabel([xname ' inventory (mmol)'])
set(gca,'FontSize',16)

figure
plot(time,log(x_inv),'.-','Linewidth',1,'MarkerSize',6)
xlabel('time (d)')
ylabel(['log(' xname ' inventory)'])
set(gca,'FontSize',16)

calc_fit = input([' Would you like to calculate a linear fit for log(' xname ')? [yes = 1, no = 0] ']);
if calc_fit

    disp(' Which time interval would you like to use for fitting the decay rate?')
    time_start = input(' Enter first day and last day. First day: ');
    time_end = input(' Last day: ');
    % determine corresponding indices of the time vector 
    ind_start = find(round(time) == time_start);
    ind_end = find(round(time) == time_end);
    ind = ind_start:ind_end;

    % linear fit in log space
    [p,s] = polyfit(time(ind),log(x_inv(ind)),1);
    hold on
    plot(time(ind),polyval(p,time(ind)),'k-')

    disp(' ')
    disp([' Residence time 2 (decay rate, i.e. inverse slope of the black line): ' num2str(-1/p(1)) ' d'])
    disp(' ')

    write_data = input(' Would you like to write the data (time and log(x)) to a text file of offline calculations? [yes = 1, no = 0] ');
    if write_data

       fn = 'data.txt';
       fid = fopen(fn,'w');
       fprintf(fid,['time, log(' xname ')\n']);
       fprintf(fid,'%.2f %.2f\n',[time(ind) log(x_inv(ind))]');
       fclose(fid); 
       disp([' Data was written to file: ' fn]) 

    end
end



% Calculate fluxes:
% 3 ways to determine fluxes

% 1. Inventory / residence time
% What?s needed for this calculation:
% Inventory of DIC
% Residence times calculated above


% 2. Tidal covariance (concentration and current velocity)
% What?s needed for this calculation:
% Time series of [DIC] at mouth and middle of estuary
% Time series of current velocities

% 3. Concentration x net flow
% What?s needed for this calculation:
% Average of time series of [DIC] at mouth and middle of estuary
% Net discharge from model


% 4. Ed Boyle?s method
% TBD
% 