function [logPosterior,ER_output,ER_contribution]=model_partition(PAR,E)

% This function run the partition model and compute the Log posterior 

% INPUT:
% "PAR": Vector of parameters
% "E": structure with model data and other information

% OUTPUT
% "logPosterior": Log of posterior distribution
% "ER_output": ER time series simulated by the model
% "ER_contribution": Time series (along the rows) of ER contributions as simulated by the model. 
                    % column 1: baseline
                    % column 2: GPP
                    % column 3: Q
                    % column 4: LF


% Add n_years_spin_up spin up to data
n_years_spin_up=2;
T_withSpinUp=[E.D.T(end-365*n_years_spin_up+1:end);E.D.T];
GPP_withSpinUp=[E.D.GPP(end-365*n_years_spin_up+1:end);E.D.GPP];
LF_withSpinUp=[E.D.LF(end-365*n_years_spin_up+1:end);E.D.LF];
Q_withSpinUp=[E.D.Q(end-365*n_years_spin_up+1:end);E.D.Q];

% Assign Parameters, see Table 2 of the paper 
ER_B_20=PAR(1);
theta_B=PAR(2);

beta_GPP=PAR(3);
theta_GPP=PAR(4);
tau_GPP=PAR(5);

beta_Q=PAR(6);
theta_Q=PAR(7);
tau_Q=PAR(8);
alpha_Q=PAR(9);

beta_LF=PAR(10);
theta_LF=PAR(11);
tau_LF=PAR(12);

rho=PAR(13);
sigma=PAR(14);

% Preallocation of variables
n_data=length(T_withSpinUp);  % Number of data points
ER_output=zeros(n_data,1);       % ER time series simulated by the model

% Compute Baseline contribution 
ER_B=ER_B_20*theta_B.^(T_withSpinUp-20);
ER_output=ER_output+ER_B; % Add Baseline contribution to total ER

% Compute GPP contribution 
length_lagfunction=ceil(tau_GPP*6); % Length of exponential lag function (6 times the avergae value to ensure proper normalization) 
lagfunction_GPP=diff(1-exp(-(0:length_lagfunction)'/tau_GPP))/(1-exp(-length_lagfunction/tau_GPP)); % Compute lag function
GPP_conv=conv(GPP_withSpinUp,lagfunction_GPP); % Compute convolution integral
ER_GPP=beta_GPP*theta_GPP.^(T_withSpinUp-20).*GPP_conv(1:n_data); % Compute GPP contribution
ER_output=ER_output+ER_GPP; % Add GPP contribution to total ER

% Compute Q (discharge) contribution
length_lagfunction=ceil(tau_Q*6); % Length of exponential lag function (6 times the avergae value to ensure proper normalization) 
trans_Q=diff(1-exp(-(0:length_lagfunction)'/tau_Q))/(1-exp(-length_lagfunction/tau_Q)); % Compute lag function
Q_conv=conv(Q_withSpinUp.^alpha_Q,trans_Q); % Compute convolution integral
ER_Q=beta_Q*theta_Q.^(T_withSpinUp-20).*Q_conv(1:n_data); % Compute Q contribution
ER_output=ER_output+ER_Q; % Add GPP contribution to total ER

% Compute LF (Litter Fall) contribution
length_lagfunction=ceil(tau_LF*6); % Length of exponential lag function (6 times the avergae value to ensure proper normalization)
trans_LF=diff(1-exp(-(0:length_lagfunction)'/tau_LF))/(1-exp(-length_lagfunction/tau_LF)); % Compute lag function
LF_conv=conv(LF_withSpinUp,trans_LF); % Compute convolution integral
ER_LF=beta_LF*theta_LF.^(T_withSpinUp-20).*LF_conv(1:n_data); % Compute LF contribution
ER_output=ER_output+ER_LF; % Add GPP contribution to total ER

% Cut spin up years form model output
ER_B=ER_B(365*n_years_spin_up+1:end);
ER_GPP=ER_GPP(365*n_years_spin_up+1:end);
ER_Q=ER_Q(365*n_years_spin_up+1:end);
ER_LF=ER_LF(365*n_years_spin_up+1:end);
ER_output=ER_output(365*n_years_spin_up+1:end);

% Create output variables
ER_contribution=[ER_B ER_GPP ER_Q ER_LF];

% Compute Log Posterior with autocorrelated error
Residuals=-E.D.ER-ER_output;
n_measurement=length(Residuals);
Residuals_2 = Residuals(2:n_measurement)-rho*Residuals(1:n_measurement-1);
logPosterior=...% Likelihood
            -(n_measurement/2)*log(2*pi)-(n_measurement/2)*log( sigma^2 / (1-rho^2)) - ...
            (1/2)*(1-rho^2)* ( Residuals(1)/sigma )^2 - (1/2) * sum((Residuals_2./sigma).^2 )...
            -4*0.5*log(2*pi*E.sigma_theta^2)...% (following) Prior part 
            -(theta_B-E.theta_mean)^2/(2*E.sigma_theta^2)...
            -(theta_GPP-E.theta_mean)^2/(2*E.sigma_theta^2)...
            -(theta_Q-E.theta_mean)^2/(2*E.sigma_theta^2)...
            -(theta_LF-E.theta_mean)^2/(2*E.sigma_theta^2);



