% This code run the analysis presented in 

% "Respiration regimes in rivers: Partitioning source‐specific respiration from metabolism time series" 
% by Bertuzzo et al., Limonology and Oceanography, 2022.
% DOI: 10.1002/lno.12207

% Run the code and select the site when prompted. The code run the model and
% reproduces figure 2 of the paper

% The code has been tested on Matlab R2021b

clear all   %clear all variables in the workspace
close all   %close all open figures

%% *************************************************************************
% Select experimental site from prompt
%***************************************************************************

for iter=1:100
    site= input('Select experimental site: 1) Black Earth; 2) Potomac; 3) Medina; 4) East Canyon : ');
    if not(site==1 || site==2 || site==3 || site==4 )
        disp('Please digit an integer number among 1, 2, 3 or 4');
    else
        break
    end
end

%% *************************************************************************
% Preliminary operations
%***************************************************************************

% Define colormap to be used
col=get(groot,'defaultAxesColorOrder');
col=col( [ 3 5 1 2],:);

%% *************************************************************************
% Load data of the site
%***************************************************************************

load(['site_',num2str(site),'_preprocessed_data']);  %load preprocessed data for the site
load(['site_',num2str(site),'_posterior_parameter']) %load parameter posterior distribution

%% *************************************************************************
% Plot panels of Figure 2 of the paper 
%***************************************************************************
[~,ER_output,ER_contribution]=model_partition(quantile(Posterior_sample,0.5),E); %run model with median parameter

alpha_face=0.8; %Set transparency 

figure('color','w','Position',[50 50 800 600]) % Create figure

subplot(2,1,1);
h=area(D.date,ER_output,'facecolor',[0.8 0.8 0.8]);
set(h,'displayname','Model')
hold on
plot(D.date,-D.ER,'ko','markerfacecolor','w','markersize',3)
h=plot(D.date,-D.ER,'ko','markerfacecolor','k','markerfacecolor','w','markersize',3,'linewidth',0.75);
set(h,'displayname','Data')
box off
ylabel('ER [gO_2 m^{-2} d^{-1}]','fontsize',12)
title(D.site_name,'fontweight','normal','fontsize',14)
l=legend('show');
set(l,'box','off','fontsize',12)

subplot(2,1,2);
hold on
hh=fill(([D.date; flipud(D.date)]),[zeros(size(ER_output));flipud(ER_contribution(:,1))],col(1,:),'edgecolor','n','facealpha',alpha_face);
set(hh,'displayname','Baseline')
hh=fill(([D.date; flipud(D.date)]),[ER_contribution(:,1);flipud(ER_contribution(:,1)+ER_contribution(:,4))],col(4,:),'edgecolor','n','facealpha',alpha_face);
set(hh,'displayname','Litter Fall')
hh=fill(([D.date ;flipud(D.date)]),[ER_contribution(:,1)+ER_contribution(:,4);flipud(ER_contribution(:,1)+ER_contribution(:,4)+ER_contribution(:,2))],col(2,:),'edgecolor','n','facealpha',alpha_face);
set(hh,'displayname','GPP')
hh=fill(([D.date; flipud(D.date)]),[ER_contribution(:,1)+ER_contribution(:,4)+ER_contribution(:,2);flipud(ER_contribution(:,1)+ER_contribution(:,4)+ER_contribution(:,2)+ER_contribution(:,3))],col(3,:),'edgecolor','n','facealpha',alpha_face);
set(hh,'displayname','Q')
hh=plot((D.date),ER_output,'-k');
set(hh,'HandleVisibility','off')
hh=legend('show');
set(hh,'box','off','fontsize',12)
box off
ylabel('ER [gO_2 m^{-2} d^{-1}]','fontsize',12)

