%% Inverse Mediterranean Box Model to calculate African Monsoon Runoff for S5 
% Amies, J. D., Rohling, E. J., Grant, K. M., Rodríguez-Sanz, L. & Marino, 
% G. (2019) ‘Quantification of African monsoon runoff during last 
% interglacial sapropel S5’. Paleoceanography and Paleoclimatology, 34. 
%
% Based on previous Mediterranean box models by Rohling et al. (2004, 2014)
%
% This model uses inputs of G.ruber(w) d18O, sea level and both alkenone and clumped isotope SST
%
%% African monsoon runoff calculation setup
% Set number of iterations for probabilistic assesment at lines 14-15
% Input data files at lines 16-36
% Switch HS11 Atlantic fw input on/off at lines 184-193
%% Set number of iterations for probabilistic assessment
N=10000;
%% Input data files
% Load sea level 
% Age in ka (col 1), RSL (col 2), 1SD error (col 3)
input=load('S5_Inputs_age_SL_alk_clump.txt');
% Input montecarlo datasets of eMed d180ruber 
load('eMed_compiled_d18O_ruber_MC_output_20190416.mat');
d18Oruber=y1;
load('eMed_compiled_d18O_ruber_MC_output_age_20190416.mat');
d18Oruber_age=t; 
% Input montecarlo datasets of eMed alkenone SST
load('eMed_alk_MC_output_20190416.mat');
alk=y1;
load('eMed_alk_MC_output_age_20190416.mat');
alk_age=t;
% Input montecarlo datasets of eMed clumped SST
load('eMed_clumped_MC_output.mat');
clump=y2;
load('eMed_clumped_MC_output_age.mat');  
clump_age=t2;
% load d18Osw for run with constrained Atlantic inflow over HS11
d18Osw=load('976_d18Osw_Skinner_Schackleton_2006.txt');
%% Interpolate input data to same age scale and interval
d=0.025;
min1 = 116.7;
max1 = 134.6;
t = min1: d : max1;
t = t';
in = interp1(input(:,1),input,t,'linear');
d18Oruber=interp1(d18Oruber_age,d18Oruber,t,'linear');
alk=interp1(alk_age,alk,t,'linear');
clump=interp1(clump_age,clump,t,'linear');
d18Osw=interp1(d18Osw(:,1),d18Osw,t,'linear');
% Normalise raw d18O inputs to present day using average for last 3 kyr
O_ru=0.42;
d18Oruber_norm = d18Oruber-O_ru;
% Normalise using correction for dsmlc to correspond to 0m at present day 
% (derived from original box model dsml to sea level correlations)
O_smlc=0.186;
d18Oruber_norm = d18Oruber_norm+O_smlc;

% calculate input dataset deviations from median
d18Oruber_med=prctile(d18Oruber_norm,50,2);
d18Oruber_dist=d18Oruber_norm-d18Oruber_med;
alk_med=prctile(alk,50,2);
alk_dist=alk-alk_med;
clump_med=prctile(clump,50,2);
clump_dist=clump-clump_med;

%% Set 'modern values used'
v=7.5;          % Wind speed (ms-1) 1945-1990 average (Garrett et al., 1993)
P=1.25*(10^12); % Precipitation volume (m^3) (Garrett et al., 1993)

% Set percentage of (random) uncertainty of Black Sea net freshwater flux for summer
sinc=100;

%% Set 'common numbers'
C=1.15*(10^-3);     % Exchange coefficient (for calculating E using latent heat flux calculations). (Garrett et al., 1993; Wells, 1995).
rho=1.2;            % Density of air (kg m^-3) at mean sea level pressure (1012 mbar)
Area=2.5*(10^12);   % Area of basin (2.5*(10^12) m^2)
e_Area=1.65*(10^12);% Area of eastern basin
w_Area=0.85*(10^12);% Area of western basin
p=1.012*(10^5);     % Total pressure (Pa) at mean sea level

%% Set 'modern values'
Sinp=36.2;          % Present-day inflow salinity in Strait of Gibraltar
dinp=1;             % Present-day inflow oxygen isotope ratio of water in Strait of Gibraltar 
Xp=1.3*(10^12);     % Excess of evaporation over freshwater input (m^3). Determined using modern X average (~0.535) multiplied by Area (see Rohling 1999)
Rp=0.45*(10^12);    % Runoff (m^3). (Garrett et al. 1993)
Ap=23*(10^12);      % Volume of present-day Atlantic flow through Strait of Gibraltar   
Bp=0.2*(10^12);     % Black Sea Input (m^3) before damming (Tolmazin 1985)
dB=-9;              % mean d18O for B (from Swart 1991, in Rohling 1999)
dR=-6;              % mean d18O for R, non-monsoon times (Rohling 1999) - excess monsoon runoff has monsoon d18O value
Np=8.4*(10^10);     % Runoff from Nile (pre-Aswan damming) (m^3) (Rohling 2015)
Np_hf=1.5*(10^11);  % Runoff from Nile high flood year (m^3) (Rohling 2015)

%% Box model
Results=zeros(size(in,1),206); % Preallocate array for results
% Interval loop
    for f=1:size(in,1)    % loop to run for each age step
        out=zeros(N,42);   % Preallocate array for loop output
        % Perform calculations N times, which gives the matrix of results over which probability intervals are determined
        for a=1:N   
            dXR=-10+randn;
            SL=in(f,2)+((in(f,3))*randn);
            Ts=alk_med(f,:)+4+emprand(alk_dist(f,:));
            Tw=alk_med(f,:)-2+emprand(alk_dist(f,:));
            Tmon=clump_med(f,:)+emprand(clump_dist(f,:));
            zwmlin=150;
            dflc=d18Oruber_med(f,:)+emprand(d18Oruber_dist(f,:));
            zsml=30+(5/3)*randn; % depth of summer mixed layer (zsml) to 30±5m where 5 is a 3 sigma range (Nykjaer, 2009)
            % Calculate DT and zUSML (as reasoned in manuscript) 
            DT=Tmon-Ts;
            if DT <= 0
                zusml=zsml;
            else
                zusml=(-4.2545+0.2485*randn)*DT+30; % randomiser covers range of zUSML:DT relationships in experiment Tmon-Ts(+4)
            end
            sl=-SL; % set sign convention for sea-level (model uses lowering as a + value)
            r=0.7+(0.05/3)*randn;   % relative humidity (Rohling, 1999) and so that modern freshwater balance and salinities agree with observations
            V=7.5+sl/120+(1/3)*randn;   % wind speed to modern annual mean (Garrett et al., 1993), with shift to stronger (winter-type) values at glacial times, and with a 3 sigma uncertainty of 1 m/s
            % set isotopic values for dR and dB to modern values with 3 sigma ranges over ±1ppt
            dRw=dR-1+(1/3)*randn;   % d18O runoff winter
            dRs=dR+1+(1/3)*randn;   % d18O runoff summer
            dRmon=dR+1+(1/3)*randn; % d18O runoff monsoon
            dPw=dRw;                % d18O precipitation winter
            dPs=dRs;                % d18O precipitation summer
            dPmon=dRw;              % d18O precipitation monsoon
            dBw=dB+(1/3)*randn;     % d18O Black Sea fw input winter
            dBs=dB+(1/3)*randn;     % d18O Black Sea fw input summer
            dBmon=dB+(1/3)*randn;   % d18O Black Sea fw input monsoon
            % set R and B fluxes to modern values with uncertainties over 3 sigma ranges of ±10%
            if sl<=80
                B=(1+(0.1/3)*randn)*Bp;
            else
                B=0;    % B is set to zero when sea level is 80 m or more below present
            end
            DTw=1.5+(1/3)*randn;                %Winter air-sea T difference (Rohling 1999)
            DTs=0.5+(1+(1/3)*randn)*(sl/120);   %Summer air-sea T difference (Rohling 1999), Artale et al. (2002) say should be very small in strong wind conditions so DTs becomes equal to DTw at glacial times (locking in winter mode).
            DTmon=DTs;
            Tssth=Tw;
            Tint=Tw-(1+(1/3)*randn);    % Intermediate waters approx 1C cooler (Rohling et al. 2004)
            Tatlantic=Tint+(1/3)*randn; % Set Atlantic inflow T so that there is no heat gain/loss in Mediterranean beyond a 3 sigma = ±1 deg C difference
            Tas=Ts-DTs;
            Taw=Tw-DTw;
            Tamon=Tmon-DTmon;
            % Calculate SUMMER 4 MONTH evaporation, following Abbott and Tabony (1985)
            Ls=(2500.84-2.34*Ts)*10^3; % Latent heat of vaporization for summer
            ess=10^2*exp(55.17-6803*(Ts+273.15)^(-1)-5.07*log(Ts+273.15)); % Summer saturation vapour pressure at SST
            eas=10^2*exp(55.17-6803*(Tas+273.15)^(-1)-5.07*log(Tas+273.15)); % Summer saturation vapour pressure at T 10m above sea-level
            qss=(ess/(p-ess))*(18.0153/28.965); % Summer saturation mixing ratio at SST (Wells 1986).
            qas=(eas/(p-eas))*(18.0153/28.965); % Summer saturation mixing ratio at air T 10m above sea-level (Wells 1986).
            Es=(rho*Ls*C*V*(qss-r*qas)*Area*1.26*10^(-2))/3; % Calculate Es (latent heat flux eqn) in m3/y over summer 4 months
            Ps=(0.4+(0.1/3)*randn)*Es; % maintain modern E:P proportion (Garrett et al., 1993)
            Rs=(1+(0.1/3)*randn)*Rp*(Ps/P); % change amount of runoff in proportion to change in amount of precipitation
            % Calculate MONSOON 2 MONTH evaporation
            Lmon=(2500.84-2.34*Tmon)*10^3;
            esmon=10^2*exp(55.17-6803*(Tmon+273.15)^(-1)-5.07*log(Tmon+273.15));
            eamon=10^2*exp(55.17-6803*(Tamon+273.15)^(-1)-5.07*log(Tamon+273.15));
            qsmon=(esmon/(p-esmon))*(18.0153/28.965);
            qamon=(eamon/(p-eamon))*(18.0153/28.965);
            Emon=(rho*Lmon*C*V*(qsmon-r*qamon)*Area*1.26*10^(-2))/6; % Calculate Emon in m3/y over monsoon 2 months
            Pmon=(0.4+(0.1/3)*randn)*Emon;
            Rmon=(1+(0.1/3)*randn)*Rp*(Pmon/P);
            % Calculate WINTER 6 MONTH evaporation
            Lw=(2500.84-2.34*Tw)*10^3;
            esw=10^2*exp(55.17-6803*(Tw+273.15)^(-1)-5.07*log(Tw+273.15));
            eaw=10^2*exp(55.17-6803*(Taw+273.15)^(-1)-5.07*log(Taw+273.15));
            qsw=(esw/(p-esw))*(18.0153/28.965);
            qaw=(eaw/(p-eaw))*(18.0153/28.965);
            Ew=(rho*Lw*C*V*(qsw-r*qaw)*Area*1.26*10^(-2))/2;    % Calculate Ew in m3/y over winter 6 months
            Pw=(0.4+(0.1/3)*randn)*Ew;
            Rw=(1+(0.1/3)*randn)*Rp*(Pw/P);
            % calculate freshwater budget
            Bs=(Bp/2+(sinc/100)*(B-Bp))*(4/6);
            Bmon=(Bp/2+(sinc/100)*(B-Bp))*(2/6);
            Bw= Bp/2+((100-sinc)/100)*(B-Bp);
            % Fractionation factors for oxygen isotopes at calcification:
            alphacalcTs=exp((2.78/(Ts+273.15)^2)*10^3-3.39*10^(-3));
            alphacalcTmon=exp((2.78/(Tmon+273.15)^2)*10^3-3.39*10^(-3));
            alphacalcTw=exp((2.78/(Tw+273.15)^2)*10^3-3.39*10^(-3));
            alphacalcTint=exp((2.78/(Tint+273.15)^2)*10^3-3.39*10^(-3));
            alphacalcTatlantic=exp((2.78/(Tatlantic+273.15)^2)*10^3-3.39*10^(-3));
            dfl=((dflc-10^3*log(alphacalcTmon))*1.03092)+30.92;
            zssth=70+5*randn; % assumed depth of zssth (shown not to affect runoff results)
            % Calculate excess of evaporation over freshwater inputs
            Xs=Es-Ps-Rs-Bs;
            Xw=Ew-Pw-Rw-Bw;
            Sinsap=Sinp+(sl/120); % past inflow salinity value, as function of present inflow values plus general glacial enrichment
            %% Switch isotopically depleted HS11 Atlantic inflow ON/OFF by commenting out relevant section below
            % OFF: No isotopically depleted HS11 Atlantic inflow
            dinsap=dinp+sl*0.009; % past inflow oxygen isotope value, as function of present inflow values plus general glacial enrichment
%             % ON: Isotopically depleted HS11 Atlantic inflow (Skinner & Shackleton,2006) to Med over H11 (Marino et al., 2015)
%             if in(f,1)<=135+(1/3)*randn && in(f,1)>=130+(2/3)*randn
%                 dinsap=d18Osw(f,2)+(0.12/3)*randn; % past inflow oxygen isotope value over H11, according to Skinner and Schakletone 2006
%             else
%                 dinsap=dinp+sl*0.009; % past inflow oxygen isotope value, as function of present inflow values plus general glacial enrichment
%             end
            %%
            phi=1-0.5*(sl/120); % exchange reduction factor PHI, as function of sea level(Rohling 1991a,1994,1999; Rohling & Bryden,1994)
            Qinp=Ap; % present day inflow from Atlantic
            % calculate the fractionation factors and residence times
            Tcw=Tw-(5+(1/3)*randn);
            Tcs=Ts-(5+(1/3)*randn);
            Tcmon=Tmon-(5+(1/3)*randn);
            % water oxygen isotope fractionation factors at evaporation:
            alphas=exp((1.137/(Ts+273.15)^2)*10^3-0.4156/(Ts+273.15)-2.0667*10^(-3));
            alphaw=exp((1.137/(Tw+273.15)^2)*10^3-0.4156/(Tw+273.15)-2.0667*10^(-3));
            alphamon=exp((1.137/(Tmon+273.15)^2)*10^3-0.4156/(Tmon+273.15)-2.0667*10^(-3));
            dEeqs=-log(alphas)*10^3;
            dEeqw=-log(alphaw)*10^3;
            dEeqmon=-log(alphamon)*10^3;
            % water oxygen isotope fractionation factor at condensation:
            alphacs=exp((1.137/(Tcs+273.15)^2)*10^3-0.4156/(Tcs+273.15)-2.0667*10^(-3));
            alphacw=exp((1.137/(Tcw+273.15)^2)*10^3-0.4156/(Tcw+273.15)-2.0667*10^(-3));
            alphacmon=exp((1.137/(Tcmon+273.15)^2)*10^3-0.4156/(Tcmon+273.15)-2.0667*10^(-3));
            datms=dRs-log(alphacs)*10^3;
            datmw=dRw-log(alphacw)*10^3;
            datmmon=dRmon-log(alphacmon)*10^3;
            Zs=qas/qss;
            Zw=qaw/qsw;
            Zmon=qamon/qsmon;
            x=zsml/(zssth+zsml); % fraction of summer mixed-layer depth relative to winter mixed-layer depth
            % volumes of different water mass layers (boxes):
            Volssth=Area*zssth;
            Volsml=Area*zsml;
            Volfl=Area*zusml;
            Volbfl=Area*(zsml-zusml);
            % T (residence time for sml, based on pre-inversion model runs)
            T=11;
            dinit=dinsap; % initial d18O - make equal to past inflow value here, then will change with every iteration in loop
            Sinit=Sinsap; % initial S - make equal to past inflow value here, then will change with every iteration in loop
            %% Box-loop calculations
            % loop set to iterate k=25 times to reflect number of years in each time step (steady state is reached in this time)
            for k=0:25
                dEs=dinit+dEeqs;
                Ssml=(((1/T)*Sinsap+((T-1)/T)*Sinit)*Volsml)/(Volsml-Xs);
                dsml=(((1/T)*dinsap+((T-1)/T)*dinit)*Volsml+Rs*dRs+Bs*dBs+Ps*dPs-Es*dEs)/(Volsml-Xs);
                for m=1:50
                    dEs=((((1+dsml*10^(-3))/alphas-Zs*r-datms*10^(-3)*Zs*r)/(1.0142*(1-Zs*r)))-1)*10^3;
                    dsml=(((1/T)*dinsap+((T-1)/T)*dinit)*Volsml+Rs*dRs+Bs*dBs+Ps*dPs-Es*dEs)/(Volsml-Xs);
                end
                Sssth=(((1/T)*Sinsap+((T-1)/T)*Sinit)*Volssth)/Volssth;
                dssth=(((1/T)*dinsap+((T-1)/T)*dinit)*Volssth)/Volssth;
                %% Monsoon box calculations
                dbfl=dsml;
                dEfl=dsml+dEeqmon;
                Sbfl=Ssml;
                XRs= ((Volfl*dsml+Rmon*dRmon+Bmon*dBmon+Pmon*dPmon-Emon*dEfl)-(dfl*Volfl-dfl*Emon+dfl*Pmon+dfl*Rmon+dfl*Bmon))/(dfl-dXR);
                Xfl=Emon-Pmon-Rmon-Bmon-XRs;
                Sfl=(Volfl*Ssml)/(Volfl-Xfl);
                Qfl=Volfl-Xfl;
                Qbfl=Volbfl;
                Qsml=Volsml-Xs;
                Qssth=Volssth;
                dEw=(zsml*dsml+zssth*dssth)/(zsml+zssth)+dEeqw;
                Swml=(Qfl*Sfl+Qbfl*Sbfl+Qssth*Sssth)/(Qfl+Qbfl+Qssth-Xw);
                dwml=(Qfl*dfl+Qbfl*dbfl+Qssth*dssth+Rw*dRw+Bw*dBw+Pw*dPw-Ew*dEw)/(Qfl+Qbfl+Qssth-Xw);
                for h=1:50
                    dEw=((((1+dwml*10^(-3))/alphaw-Zw*r-datmw*10^(-3)*Zw*r)/(1.0142*(1-Zw*r)))-1)*10^3;
                    dwml=(Qfl*dfl+Qbfl*dbfl+Qssth*dssth+Rw*dRw+Bw*dBw+Pw*dPw-Ew*dEw)/(Qsml+Qssth-Xw);
                end
                Sinit=Swml;
                dinit=dwml;
                dsmlc=10^3*log(alphacalcTs)+(dsml-30.92)/1.03092;
                dssthc=10^3*log(alphacalcTw)+(dssth-30.92)/1.03092;
                dwmlc=10^3*log(alphacalcTw)+(dwml-30.92)/1.03092;
                dintc=10^3*log(alphacalcTint)+(dwml-30.92)/1.03092;
                dac=10^3*log(alphacalcTatlantic)+(dinsap-30.92)/1.03092;
                zwml=zsml+zssth;
            end
            % Correct XRs for area of eMed
            XRs=(XRs/Area)*e_Area;
            %% Runoff calculations - as reasoned in manuscript
            if DT <= 0 
                R=Rs+Rw+Rmon;
            else
                R=(XRs/(sinc/100))+(Rs+Rw+Rmon);
            end
            % Calculate monsoon intensification compared to total Med runoff
            M=R/Rp;
            % Calculate volume of African runoff (S5 total runoff - total Rp + Nile contribution)
            RAf=R-((Rs+Rw+Rmon)-Np);
            % monsoon intensification based on Nile high flood year
            MAf=RAf/Np_hf;
        % Other calculations       
        Xsap=Xs+Xfl+Xw; % total excess of evaporation over year
        gamma=Xsap/Xp;  % exchange reduction factor GAMMA as function of buoyancy loss after Rohling (1991b,1994,1999; Rohling & Bryden,1994)
        omega=gamma^(1/3);
        Qinsap=omega*phi*Qinp;  % exchange relative to present
        Qoutsap=Qinsap-Xsap;
        %calculate pycnocline shoaling, after Rohling & Bryden (1994)
        alpha=(((Ap/Xp)*phi*omega-gamma)/(Ap/Xp-1))*(Sinp/Sinsap);
        outa=[-sl,Sfl,Ssml,Sssth,Swml,dfl,dsml,dssth,dwml,dflc,dsmlc,dssthc,dwmlc,dintc,dac,Sinsap,dinsap,zsml,zssth,Xsap,Ew,Es,Emon,zwml,Xw,Xs,Xfl,zusml,M,R,DT,Pw,Ps,Pmon,Bw,Bs,Bmon,Rw,Rs,Rmon,RAf,MAf];
        out(a,:)=outa;
        end
        subout=out(1:N,:);    % take results for each f time step iteration
        % take median and percentiles of each property output by model
        Sfldata=[median(subout(:,2)),prctile(subout(:,2),16),prctile(subout(:,2),2.5),prctile(subout(:,2),84),prctile(subout(:,2),97.5)];
        Ssmldata=[median(subout(:,3)),prctile(subout(:,3),16),prctile(subout(:,3),2.5),prctile(subout(:,3),84),prctile(subout(:,3),97.5)];
        Sssthdata=[median(subout(:,4)),prctile(subout(:,4),16),prctile(subout(:,4),2.5),prctile(subout(:,4),84),prctile(subout(:,4),97.5)];
        Swmldata=[median(subout(:,5)),prctile(subout(:,5),16),prctile(subout(:,5),2.5),prctile(subout(:,5),84),prctile(subout(:,5),97.5)];
        dfldata=[median(subout(:,6)),prctile(subout(:,6),16),prctile(subout(:,6),2.5),prctile(subout(:,6),84),prctile(subout(:,6),97.5)];
        dsmldata=[median(subout(:,7)),prctile(subout(:,7),16),prctile(subout(:,7),2.5),prctile(subout(:,7),84),prctile(subout(:,7),97.5)];
        dssthdata=[median(subout(:,8)),prctile(subout(:,8),16),prctile(subout(:,8),2.5),prctile(subout(:,8),84),prctile(subout(:,8),97.5)];
        dwmldata=[median(subout(:,9)),prctile(subout(:,9),16),prctile(subout(:,9),2.5),prctile(subout(:,9),84),prctile(subout(:,9),97.5)];
        dflcdata=[median(subout(:,10)),prctile(subout(:,10),16),prctile(subout(:,10),2.5),prctile(subout(:,10),84),prctile(subout(:,10),97.5)];
        dsmlcdata=[median(subout(:,11)),prctile(subout(:,11),16),prctile(subout(:,11),2.5),prctile(subout(:,11),84),prctile(subout(:,11),97.5)];
        dssthcdata=[median(subout(:,12)),prctile(subout(:,12),16),prctile(subout(:,12),2.5),prctile(subout(:,12),84),prctile(subout(:,12),97.5)];
        dwmlcdata=[median(subout(:,13)),prctile(subout(:,13),16),prctile(subout(:,13),2.5),prctile(subout(:,13),84),prctile(subout(:,13),97.5)];
        dintcdata=[median(subout(:,14)),prctile(subout(:,14),16),prctile(subout(:,14),2.5),prctile(subout(:,14),84),prctile(subout(:,14),97.5)];
        dacdata=[median(subout(:,15)),prctile(subout(:,15),16),prctile(subout(:,15),2.5),prctile(subout(:,15),84),prctile(subout(:,15),97.5)];
        Sinsapdata=[median(subout(:,16)),prctile(subout(:,16),16),prctile(subout(:,16),2.5),prctile(subout(:,16),84),prctile(subout(:,16),97.5)];
        dinsapdata=[median(subout(:,17)),prctile(subout(:,17),16),prctile(subout(:,17),2.5),prctile(subout(:,17),84),prctile(subout(:,17),97.5)];
        zsmldata=[median(subout(:,18)),prctile(subout(:,18),16),prctile(subout(:,18),2.5),prctile(subout(:,18),84),prctile(subout(:,18),97.5)];
        zssthdata=[median(subout(:,19)),prctile(subout(:,19),16),prctile(subout(:,19),2.5),prctile(subout(:,19),84),prctile(subout(:,19),97.5)];
        Xsapdata=[median(subout(:,20))/Area,(prctile(subout(:,20),16))/Area,(prctile(subout(:,20),2.5))/Area,(prctile(subout(:,20),84))/Area,(prctile(subout(:,20),97.5))/Area];
        Ewdata=[median(subout(:,21))/Area,(prctile(subout(:,21),16))/Area,(prctile(subout(:,21),2.5))/Area,(prctile(subout(:,21),84))/Area,(prctile(subout(:,21),97.5))/Area];
        Esdata=[median(subout(:,22))/Area,(prctile(subout(:,22),16))/Area,(prctile(subout(:,22),2.5))/Area,(prctile(subout(:,22),84))/Area,(prctile(subout(:,22),97.5))/Area];
        Emondata=[median(subout(:,23))/Area,(prctile(subout(:,23),16))/Area,(prctile(subout(:,23),2.5))/Area,(prctile(subout(:,23),84))/Area,(prctile(subout(:,23),97.5))/Area];
        zwmldata=[median(subout(:,24)),prctile(subout(:,24),16),prctile(subout(:,24),2.5),prctile(subout(:,24),84),prctile(subout(:,24),97.5)];
        Xwdata=[median(subout(:,25))/Area,(prctile(subout(:,25),16))/Area,(prctile(subout(:,25),2.5))/Area,(prctile(subout(:,25),84))/Area,(prctile(subout(:,25),97.5))/Area];
        Xsdata=[median(subout(:,26))/Area,(prctile(subout(:,26),16))/Area,(prctile(subout(:,26),2.5))/Area,(prctile(subout(:,26),84))/Area,(prctile(subout(:,26),97.5))/Area];
        Xfldata=[median(subout(:,27))/Area,(prctile(subout(:,27),16))/Area,(prctile(subout(:,27),2.5))/Area,(prctile(subout(:,27),84))/Area,(prctile(subout(:,27),97.5))/Area];   
        zusmldata=[median(subout(:,28)),prctile(subout(:,28),16),prctile(subout(:,28),2.5),prctile(subout(:,28),84),prctile(subout(:,28),97.5)];
        Mdata=[median(subout(:,29)),prctile(subout(:,29),16),prctile(subout(:,29),2.5),prctile(subout(:,29),84),prctile(subout(:,29),97.5)];
        Rdata=[median(subout(:,30)),prctile(subout(:,30),16),prctile(subout(:,30),2.5),prctile(subout(:,30),84),prctile(subout(:,30),97.5)];
        DTdata=[median(subout(:,31)),prctile(subout(:,31),16),prctile(subout(:,31),2.5),prctile(subout(:,31),84),prctile(subout(:,31),97.5)];
        Pwdata=[median(subout(:,32))/Area,(prctile(subout(:,32),16))/Area,(prctile(subout(:,32),2.5))/Area,(prctile(subout(:,32),84))/Area,(prctile(subout(:,32),97.5))/Area];
        Psdata=[median(subout(:,33))/Area,(prctile(subout(:,33),16))/Area,(prctile(subout(:,33),2.5))/Area,(prctile(subout(:,33),84))/Area,(prctile(subout(:,33),97.5))/Area];
        Pmondata=[median(subout(:,34))/Area,(prctile(subout(:,34),16))/Area,(prctile(subout(:,34),2.5))/Area,(prctile(subout(:,34),84))/Area,(prctile(subout(:,34),97.5))/Area];
        Bwdata=[median(subout(:,35))/Area,(prctile(subout(:,35),16))/Area,(prctile(subout(:,35),2.5))/Area,(prctile(subout(:,35),84))/Area,(prctile(subout(:,35),97.5))/Area];
        Bsdata=[median(subout(:,36))/Area,(prctile(subout(:,36),16))/Area,(prctile(subout(:,36),2.5))/Area,(prctile(subout(:,36),84))/Area,(prctile(subout(:,36),97.5))/Area];
        Bmondata=[median(subout(:,37))/Area,(prctile(subout(:,37),16))/Area,(prctile(subout(:,37),2.5))/Area,(prctile(subout(:,37),84))/Area,(prctile(subout(:,37),97.5))/Area];
        Rwdata=[median(subout(:,38))/Area,(prctile(subout(:,38),16))/Area,(prctile(subout(:,38),2.5))/Area,(prctile(subout(:,38),84))/Area,(prctile(subout(:,38),97.5))/Area];
        Rsdata=[median(subout(:,39))/Area,(prctile(subout(:,39),16))/Area,(prctile(subout(:,39),2.5))/Area,(prctile(subout(:,39),84))/Area,(prctile(subout(:,39),97.5))/Area];
        Rmondata=[median(subout(:,40))/Area,(prctile(subout(:,40),16))/Area,(prctile(subout(:,40),2.5))/Area,(prctile(subout(:,40),84))/Area,(prctile(subout(:,40),97.5))/Area];
        RAfdata=[median(subout(:,41)),prctile(subout(:,41),16),prctile(subout(:,41),2.5),prctile(subout(:,41),84),prctile(subout(:,41),97.5)];
        MAfdata=[median(subout(:,42)),prctile(subout(:,42),16),prctile(subout(:,42),2.5),prctile(subout(:,42),84),prctile(subout(:,42),97.5)];
        resa=[-sl,Sfldata,Ssmldata,Sssthdata,Swmldata,dfldata,dsmldata,dssthdata,dwmldata,dflcdata,dsmlcdata,dssthcdata,dwmlcdata,dintcdata,dacdata,Xsapdata,Ewdata,Esdata,Emondata,zwmldata,Xwdata,Xsdata,Xfldata,zusmldata,zsmldata,zssthdata,Mdata,Rdata,DTdata,Sinsapdata,dinsapdata,Pwdata,Psdata,Pmondata,Bwdata,Bsdata,Bmondata,Rwdata,Rsdata,Rmondata,RAfdata,MAfdata];
        Results(f,:)=resa;% compile results matrix from each age step iteration
    end