% Function to loops over the EBM and AMOC parts of the AMOC-emulator for all forcing scenarios

function [moutput,Toce,Soce,tatm,tbnd]=runmodel_coupledsystem(Par)

load('constants.mat');

% Use GCM-based regional temperatures to compute 'total atmospheric effect'
Total_atmos_effect

% Load GCM-based total atmospheric effect time series
load('atmeff.mat')

k=Par(1);
thau=Par(2);
F1=Par(3);
F2=Par(4);
V4=Par(5);
h1=Par(6);
h2=Par(7);

% Initialize arrays
moutput = zeros(1,(eq_length+trans_length)/ebm_res*amoc_res);
tatm=zeros((eq_length+trans_length)/ebm_res,jmt);
tbnd=zeros((eq_length+trans_length)/ebm_res,jmt);
Toce = zeros((eq_length+trans_length)/ebm_res,4);
Soce = zeros((eq_length+trans_length)/ebm_res,4);
EBM_time=1;
AMOC_time=1;

%% Run AMOC-emulator to equilibrium
if (eq_length > 0)
% Initialize arrays
    tatmnew=zeros(1,jmt);
    tbndnew=zeros(1,jmt);
    Tocenew=zeros(1,jmt-3);
    Socenew=zeros(1,jmt-3);     
    tatmold=zeros(1,jmt);
    tbndold=zeros(1,jmt);
    Toceold=zeros(1,jmt-3);
    Soceold=zeros(1,jmt-3);    

    for time=1:eq_length % Run AMOC-emulator for eq_length timesteps
        [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Toceold,thaui,atmeff(EBM_time,:)); % Run EBM using ocean temperatures from previous timestep
        if (EBM_time/(ebm_res/amoc_res)==floor(EBM_time/(ebm_res/amoc_res))) % If EBM_time/ebm_res/amoc_res is integer run AMOC
            [moutput(AMOC_time),Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbndnew,k,thau,F1,F2,V4,h1,h2); % Run AMOC using updated atmospheric boundary temperatures
            Toceold=Tocenew;
            Soceold=Socenew;
            AMOC_time=AMOC_time+1; % Increase AMOC timestep
        end
        tatmold=tatmnew;
        tbndold=tbndnew;
        EBM_time=EBM_time+1;
        if (EBM_time/ebm_res==floor(EBM_time/ebm_res)) % Save yearly EBM output values
            tatm(EBM_time/ebm_res,:)=tatmold;
            tbnd(EBM_time/ebm_res,:)=tbndold;
            Toce(EBM_time/ebm_res,:)=Toceold;
            Soce(EBM_time/ebm_res,:)=Soceold;            
        end
    end
    % Save temperature and salinity fields as equilibrium values
    tatmeq=tatmold;
    save('tatmeq.mat','tatmeq')
    tbndeq=tbndold;
    save('tbndeq.mat','tbndeq')
    Toceeq=Toceold;
    save('Toceeq.mat','Toceeq')
    Soceeq=Soceold;
    save('Soceeq.mat','Soceeq')
end

%% Run coupled EBM AMOC transient for historical period
if (trans_length_hist > 0)
    AMOC_time=max(AMOC_time,1);
    EBM_time=max(EBM_time,1);    
    % Load temperature and saliniyt fields as initial values
    load('tatmeq.mat','tatmeq')
    tatmold=tatmeq;
    load('tbndeq.mat','tbndeq')
    tbndold=tbndeq;
    load('Toceeq.mat','Toceeq')
    Toceold=Toceeq;
    load('Soceeq.mat','Soceeq')
    Soceold=Soceeq;

    for time=(1+eq_length):(eq_length+trans_length_hist*ebm_res) % Run AMOC-emulator for trans_length_hist*ebm_res timesteps
        [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Toceold,thaui,atmeff(EBM_time,:)); % Run EBM using ocean temperatures from previous timestep
        if (EBM_time/(ebm_res/amoc_res)==floor(EBM_time/(ebm_res/amoc_res))) % If EBM_time/ebm_res/amoc_res is integer run AMOC
            [moutput(AMOC_time),Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbndnew,k,thau,F1,F2,V4,h1,h2); % Run AMOC using updated atmospheric boundary temperatures
            Toceold=Tocenew;
            Soceold=Socenew;
            AMOC_time=AMOC_time+1; % Increase AMOC timestep
        end
        tatmold=tatmnew;
        tbndold=tbndnew;
        EBM_time=EBM_time+1;
        if (EBM_time/ebm_res==floor(EBM_time/ebm_res)) % Save yearly EBM output values
            EBM_time/ebm_res;
            tatm(EBM_time/ebm_res,:)=tatmold;
            tbnd(EBM_time/ebm_res,:)=tbndold;
            Toce(EBM_time/ebm_res,:)=Toceold;
            Soce(EBM_time/ebm_res,:)=Soceold;            
        end
    end
    % Save temperature and salinity fields as historical (2006) values
    tatmhist=tatmold;
    save('tatmhist.mat','tatmhist')
    tbndhist=tbndold;
    save('tbndhist.mat','tbndhist')
    Tocehist=Toceold;
    save('Tocehist.mat','Tocehist')
    Socehist=Soceold;
    save('Socehist.mat','Socehist')
end

%% Run coupled EBM AMOC transient for rcp45_nomelt
if (trans_length_rcp45nomelt > 0)
    AMOC_time=max(AMOC_time,1);
    EBM_time=max(EBM_time,1);    
    % Load temperature and salinity fields for historical (2006) period
    load('tatmhist.mat','tatmhist')
    tatmold=tatmhist;
    load('tbndhist.mat','tbndhist')
    tbndold=tbndhist;
    load('Tocehist.mat','Tocehist')
    Toceold=Tocehist;
    load('Socehist.mat','Socehist')
    Soceold=Socehist;

    for time=(1+eq_length+trans_length_hist*ebm_res):(eq_length+(trans_length_hist+trans_length_rcp45nomelt)*ebm_res) % Run AMOC-emulator for trans_length_rcp45nomelt*ebm_res timesteps
        [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Toceold,thaui,atmeff(EBM_time,:)); % Run EBM using ocean temperatures from previous timestep
        if (EBM_time/(ebm_res/amoc_res)==floor(EBM_time/(ebm_res/amoc_res))) % If EBM_time/ebm_res/amoc_res is integer run AMOC
            [moutput(AMOC_time),Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbndnew,k,thau,F1,F2,V4,h1,h2); % Run AMOC using updated atmospheric boundary temperatures
            Toceold=Tocenew;
            Soceold=Socenew;
            AMOC_time=AMOC_time+1; % Increase AMOC timestep
        end
        tatmold=tatmnew;
        tbndold=tbndnew;
        EBM_time=EBM_time+1;
        if (EBM_time/ebm_res==floor(EBM_time/ebm_res)) % Save yearly EBM output values
            EBM_time/ebm_res;
            tatm(EBM_time/ebm_res,:)=tatmold;
            tbnd(EBM_time/ebm_res,:)=tbndold;
            Toce(EBM_time/ebm_res,:)=Toceold;
            Soce(EBM_time/ebm_res,:)=Soceold;            
        end
    end
end

%% Run coupled EBM AMOC transient for rcp45gGISmelt
if (trans_length_rcp45gGISmelt > 0)
    AMOC_time=max(AMOC_time,1);
    EBM_time=max(EBM_time,1);    
    % Load temperature and salinity fields for historical (2006) period
    load('tatmhist.mat','tatmhist')
    tatmold=tatmhist;
    load('tbndhist.mat','tbndhist')
    tbndold=tbndhist;
    load('Tocehist.mat','Tocehist')
    Toceold=Tocehist;
    load('Socehist.mat','Socehist')
    Soceold=Socehist;
    
    for time=(1+eq_length+(trans_length_hist+trans_length_rcp45nomelt)*ebm_res):(eq_length+(trans_length_hist+trans_length_rcp45nomelt+trans_length_rcp45gGISmelt)*ebm_res) % Run AMOC-emulator for trans_length_rcp45gGISmelt*ebm_res timesteps
        [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Toceold,thaui,atmeff(EBM_time,:)); % Run EBM using ocean temperatures from previous timestep
        if (EBM_time/(ebm_res/amoc_res)==floor(EBM_time/(ebm_res/amoc_res))) % If EBM_time/ebm_res/amoc_res is integer run AMOC
            [moutput(AMOC_time),Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbndnew,k,thau,F1,F2,V4,h1,h2); % Run AMOC using updated atmospheric boundary temperatures
            Toceold=Tocenew;
            Soceold=Socenew;
            AMOC_time=AMOC_time+1; % Increase AMOC timestep
        end
        tatmold=tatmnew;
        tbndold=tbndnew;
        EBM_time=EBM_time+1;
        if (EBM_time/ebm_res==floor(EBM_time/ebm_res)) % Save yearly EBM output values
            EBM_time/ebm_res;
            tatm(EBM_time/ebm_res,:)=tatmold;
            tbnd(EBM_time/ebm_res,:)=tbndold;
            Toce(EBM_time/ebm_res,:)=Toceold;
            Soce(EBM_time/ebm_res,:)=Soceold;            
        end
    end
end

%% Run coupled EBM AMOC transient for rcp85_nomelt
if (trans_length_rcp85nomelt > 0)
    AMOC_time=max(AMOC_time,1);
    EBM_time=max(EBM_time,1);    
    % Load temperature and salinity fields for historical (2006) period
    load('tatmhist.mat','tatmhist')
    tatmold=tatmhist;
    load('tbndhist.mat','tbndhist')
    tbndold=tbndhist;
    load('Tocehist.mat','Tocehist')
    Toceold=Tocehist;
    load('Socehist.mat','Socehist')
    Soceold=Socehist;

    for time=(1+eq_length+(trans_length_hist+trans_length_rcp45nomelt+trans_length_rcp45gGISmelt+trans_length_rcp45rGISmelt)*ebm_res):(eq_length+(trans_length_hist+trans_length_rcp45nomelt+trans_length_rcp45gGISmelt+trans_length_rcp45rGISmelt+trans_length_rcp85nomelt)*ebm_res) % Run AMOC-emulator for trans_length_rcp85nomelt*ebm_res timesteps
        [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Toceold,thaui,atmeff(EBM_time,:)); % Run EBM using ocean temperatures from previous timestep
        if (EBM_time/(ebm_res/amoc_res)==floor(EBM_time/(ebm_res/amoc_res))) % If EBM_time/ebm_res/amoc_res is integer run AMOC
            [moutput(AMOC_time),Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbndnew,k,thau,F1,F2,V4,h1,h2); % Run AMOC using updated atmospheric boundary temperatures
            Toceold=Tocenew;
            Soceold=Socenew;
            AMOC_time=AMOC_time+1; % Increase AMOC timestep
        end
        tatmold=tatmnew;
        tbndold=tbndnew;
        EBM_time=EBM_time+1;
        if (EBM_time/ebm_res==floor(EBM_time/ebm_res)) % Save yearly EBM output values
            EBM_time/ebm_res;
            tatm(EBM_time/ebm_res,:)=tatmold;
            tbnd(EBM_time/ebm_res,:)=tbndold;
            Toce(EBM_time/ebm_res,:)=Toceold;
            Soce(EBM_time/ebm_res,:)=Soceold;            
        end
    end
end

%% Run coupled EBM AMOC transient for rcp85gGISmelt
if (trans_length_rcp85gGISmelt > 0)
    AMOC_time=max(AMOC_time,1);
    EBM_time=max(EBM_time,1);
    % Load temperature and salinity fields for historical (2006) period
    load('tatmhist.mat','tatmhist')
    tatmold=tatmhist;
    load('tbndhist.mat','tbndhist')
    tbndold=tbndhist;
    load('Tocehist.mat','Tocehist')
    Toceold=Tocehist;
    load('Socehist.mat','Socehist')
    Soceold=Socehist;

    for time=(1+eq_length+(trans_length_hist+trans_length_rcp45nomelt+trans_length_rcp45gGISmelt+trans_length_rcp45rGISmelt+trans_length_rcp85nomelt)*ebm_res):(eq_length+(trans_length_hist+trans_length_rcp45nomelt+trans_length_rcp45gGISmelt+trans_length_rcp45rGISmelt+trans_length_rcp85nomelt+trans_length_rcp85gGISmelt)*ebm_res) % Run AMOC-emulator for trans_length_rcp85gGISmelt*ebm_res timesteps
        [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Toceold,thaui,atmeff(EBM_time,:)); % Run EBM using ocean temperatures from previous timestep
        if (EBM_time/(ebm_res/amoc_res)==floor(EBM_time/(ebm_res/amoc_res))) % If EBM_time/ebm_res/amoc_res is integer run AMOC
            [moutput(AMOC_time),Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbndnew,k,thau,F1,F2,V4,h1,h2); % Run AMOC using updated atmospheric boundary temperatures
            Toceold=Tocenew;
            Soceold=Socenew;
            AMOC_time=AMOC_time+1; % Increase AMOC timestep
        end
        tatmold=tatmnew;
        tbndold=tbndnew;
        EBM_time=EBM_time+1;
        if (EBM_time/ebm_res==floor(EBM_time/ebm_res)) % Save yearly EBM output values
            EBM_time/ebm_res;
            tatm(EBM_time/ebm_res,:)=tatmold;
            tbnd(EBM_time/ebm_res,:)=tbndold;
            Toce(EBM_time/ebm_res,:)=Toceold;
            Soce(EBM_time/ebm_res,:)=Soceold;            
        end
    end
end

end
