%%

clear all; clc;

dt = 0.01;             % timestep (days)
tmax = 365*10;     % integration length (days)
nt = tmax / dt;
alpha = 1/365;
q = 1/1825;
tau = 1/60;
alpha_T = 0.2;      
alpha_S = 0.8;

%epsilon = -0.1;

%%% Prescribed:
T_i     = 15;
T_b     = -5;
T_atm   = -10;

S_i     = -1;
S_b     = 2;

FS      = 3*q;


rho_b = -alpha_T*T_b + alpha_S*S_b;

S = nan(nt,1);
T = nan(nt,1);
conv = nan(nt,1);
nonconv = nan(nt,1);

% initial values
S(1)=0;
T(1)=-5;

Teq = nan(nt,1);
Seq = nan(nt,1);

for i=1:nt
    
    rho = -alpha_T*T(i) + alpha_S*S(i);
    
    if rho_b - rho > 0
        delta_rho = 0;
        nonconv(i) = 1;
    else
        delta_rho = 1;
        conv(i) = 1;
    end
    
    kappa = q + tau*delta_rho;
    Phi_S = FS  + q*S_i + tau*delta_rho*S_b;
    
    Phi_T = alpha*T_atm+q*T_i+tau*delta_rho*T_b;

%    Teq(i) = Phi_T / (alpha+kappa);%(q+alpha+tau*delta_rho);
%    Seq(i) = Phi_S / (kappa);%(q+tau*delta_rho);
    
    
    dSdt = Phi_S - kappa*S(i);
    
    dTdt = Phi_T - (alpha+kappa)*T(i);

    if i<nt
        S(i+1) = S(i) + dSdt*dt;
        T(i+1) = T(i) + dTdt*dt;
    end
end
    Phi_S = FS + q*S_i + tau*S_b;
    Phi_T = alpha*T_atm+q*T_i+tau*T_b;
    kappa = q + tau;
    Teq(1) = Phi_T / (alpha+kappa);%(q+alpha+tau*delta_rho);
    Seq(1) = Phi_S / (kappa);%(q+tau*delta_rho);
    
    Phi_S = FS  + q*S_i;
    Phi_T = alpha*T_atm+q*T_i;
    Teq(2) = Phi_T / (alpha+q);%(q+alpha+tau*delta_rho);
    Seq(2) = Phi_S / (q);%(q+tau*delta_rho);

rho = -alpha_T*T+alpha_S*S;
Tvec = -8:8;
Svec = -3:3;
%[Tgrid,Sgrid] = meshgrid(Tvec,Svec);
%rho_grid = (-alpha_T*Tgrid + alpha_S*Sgrid) - rho_b;

time = [dt:dt:tmax]/365;
figure(1);
subplot(2,2,1:2)
    plot(time,-alpha_T*T(1:nt),'.',time,alpha_S*S(1:nt),'.',time,rho-rho_b,'.')
    grid on;
    legend('-\alpha_T T','\alpha_S S','\rho - \rho_b','FontSize',14)
    xlabel('Time (yr)','FontSize',14)
subplot(2,2,3)
    hold off; plot(Teq,Seq,'*',Teq.*nonconv,Seq.*nonconv,'*')
    legend('Convecting','Nonconvecting','Location','NorthWest')
    hold on; plot(T(1:nt).*conv,S(1:nt).*conv,T(1:nt).*nonconv,S(1:nt).*nonconv)
    xlabel('Temperature ','FontSize',14)
    ylabel('Salinity ','FontSize',14)
    grid on;
subplot(2,2,4)
%    C = [-1:0.05:1];
    hold off; plot(T(0.5*nt:nt).*conv(0.5*nt:nt),S(0.5*nt:nt).*conv(0.5*nt:nt),T(0.5*nt:nt).*nonconv(0.5*nt:nt),S(0.5*nt:nt).*nonconv(0.5*nt:nt),'LineWidth',2)
%    hold on; axis tight; [C,h]=contour(Tvec,Svec,rho_grid,C,'k');
    xlabel('Temperature','FontSize',14)
    ylabel('Salinity','FontSize',14)
%    clabel(C,h);



Phi_S = alpha_S*(FS+q*(S_i-S_b))
Phi_T = -alpha_T*(alpha*(T_atm-T_b)+q*(T_i-T_b));
kappa_tau = (q+tau) / (q+alpha+tau);
kappa_zero = q / (q+alpha);


fprintf('must be greater than %.5f\n',-kappa_tau*Phi_T)
fprintf('must be less than %.5f\n',-kappa_zero*Phi_T)
