% Function that simulated AMOC part of AMOC-emulator following Zickfeld et al.
function [mnew,Tocenew,Socenew] = runmodel_amoc(Toceold,Soceold,AMOC_time,tbnd,k,thau,F1t,F2t,V4,h1,h2)

%% Load constants
persistent amoc_con tempGISmelt

if AMOC_time==1 % Load constants and files in first timestep
    load('constants.mat');
    load('bandweight.mat');
    load('Forcing_data/Tglob_ref_gcm.mat');
    amoc_con=zeros(11,5);
    amoc_con(1:6,1)=[thaui;yr2sec;rhowat;Tglob_ref;dt_amoc;S0];
    amoc_con(7,1:3)=lambdai;
    amoc_con(8:10,1:4)=[V;Teq;Seq];
    amoc_con(9,:)=bandweight;
end

Ts=tbnd(3:5); % Use 'tbnd' (atmospheric surface temperatures) instead of T-start values as in Zickfeld et al.

lambda=amoc_con(7,1:3)/amoc_con(1,1)*thau; % Correct lambda values for changes in thau
amoc_con(8,4)=V4; % Use correct volume of oceanic box 4.

%% GISmelt
if AMOC_time==1 % Load GISmelt forcing file in first timestep
    load('Forcing_data/GISmelt_gcm.mat');
    tempGISmelt=GISmelt;
end
GISmelt = tempGISmelt(AMOC_time); % Get GIS melt forcing value for current AMOC timestep

%% Calculate total freshwater fluxes
dTglob=sum(tbnd(2:6).*amoc_con(9,:))/sum(amoc_con(9,:))-amoc_con(4,1); % Calculate dTglob (global temperature change) value for EBM_time-step
F1=(F1t+h1*dTglob)*10^6*amoc_con(2,1); % Calculate F1 from base value and temperature dependent part using linear relation F1 and global temperature change
F2=(F2t+h2*dTglob)*10^6*amoc_con(2,1); % Calculate F2 from base value and temperature dependent part using linear relation F2 and global temperature change
F2GIS=GISmelt*10^6*amoc_con(2,1);

%% Run AMOC model
if (AMOC_time>1 || eq_length==0)
    [drho1]=density(Toceold(1),Soceold(1));
    [drho3]=density(Toceold(3),Soceold(3));
    m = max(100,k/amoc_con(3,1)*(drho3-drho1)); % Do not allow (nearly) negative AMOC strength values
    Tocenew(1)=Toceold(1)+(m/amoc_con(8,1)*(Toceold(4)-Toceold(1))+lambda(1)*(Ts(1)-Toceold(1)))*amoc_con(5,1);
    Tocenew(2)=Toceold(2)+(m/amoc_con(8,2)*(Toceold(1)-Toceold(2))+lambda(2)*(Ts(2)-Toceold(2)))*amoc_con(5,1);
    Tocenew(3)=Toceold(3)+(m/amoc_con(8,3)*(Toceold(2)-Toceold(3))+lambda(3)*(Ts(3)-Toceold(3)))*amoc_con(5,1);
    Tocenew(4)=Toceold(4)+(m/amoc_con(8,4)*(Toceold(3)-Toceold(4)))*amoc_con(5,1);
    Socenew(1)= Soceold(1)+(m/amoc_con(8,1)*(Soceold(4)-Soceold(1))+amoc_con(6,1)*F1/amoc_con(8,1))*amoc_con(5,1);
    Socenew(2)= Soceold(2)+(m/amoc_con(8,2)*(Soceold(1)-Soceold(2))-amoc_con(6,1)*(F1-F2)/amoc_con(8,2))*amoc_con(5,1);
    Socenew(3)= Soceold(3)+(m/amoc_con(8,3)*(Soceold(2)-Soceold(3))-amoc_con(6,1)*(F2+F2GIS)/amoc_con(8,3))*amoc_con(5,1);
    Socenew(4)= Soceold(4)+(m/amoc_con(8,4)*(Soceold(3)-Soceold(4)))*amoc_con(5,1);
else
    % Use initial ocean temperature and salinity values at AMOC_time=1
    [drho1]=density(amoc_con(9,1),amoc_con(10,1));
    [drho3]=density(amoc_con(9,3),amoc_con(10,3));
    m = max(100,k/amoc_con(5,1)*(drho3-drho1)); % Do not allow (nearly) negative AMOC strength values
    Tocenew(1)=amoc_con(9,1)+(m/amoc_con(8,1)*(amoc_con(9,4)-amoc_con(9,1))+lambda(1)*(Ts(1)-amoc_con(9,1)))*amoc_con(5,1);
    Tocenew(2)=amoc_con(9,2)+(m/amoc_con(8,2)*(amoc_con(9,1)-amoc_con(9,2))+lambda(2)*(Ts(2)-amoc_con(9,2)))*amoc_con(5,1);
    Tocenew(3)=amoc_con(9,3)+(m/amoc_con(8,3)*(amoc_con(9,2)-amoc_con(9,3))+lambda(3)*(Ts(3)-amoc_con(9,3)))*amoc_con(5,1);
    Tocenew(4)=amoc_con(9,4)+(m/amoc_con(8,4)*(amoc_con(9,3)-amoc_con(9,4)))*amoc_con(5,1);
    Socenew(1)=amoc_con(10,1)+(m/amoc_con(8,1)*(amoc_con(10,4)-amoc_con(10,1))+amoc_con(6,1)*F1/amoc_con(8,1))*amoc_con(5,1);
    Socenew(2)=amoc_con(10,2)+(m/amoc_con(8,2)*(amoc_con(10,1)-amoc_con(10,2))-amoc_con(6,1)*(F1-F2)/amoc_con(8,2))*amoc_con(5,1);
    Socenew(3)=amoc_con(10,3)+(m/amoc_con(8,3)*(amoc_con(10,2)-amoc_con(10,3))+amoc_con(6,1)*F2/amoc_con(8,3))*amoc_con(5,1);
    Socenew(4)=amoc_con(10,4)+(m/amoc_con(8,4)*(amoc_con(10,3)-amoc_con(10,4)))*amoc_con(5,1);
end
mnew=m;
end
