% A matlab version of the simple model of:
% Galanti, E. and E. Tziperman, 2000. On ENSO's phase locking to
% the seasonal cycle in the fast SST, fast wave, and mixed mode
% regimes. Journal of the Atmospheric Sciences. 57, 2936-2950.  

clear all; close all;

% initialize h and T
%-------------------
delta_t=0.1;
i_start=40/delta_t;
i_end=1000*12/delta_t;
h(1:i_start-1)=1.e-4;
h(i_start:i_end)=0.0;
T(1:i_end)=0.0;
mu(1:i_end)=0.0;

% model parameters
%-----------------

Robert=0.25;
ru=1000;
b0=8.1E10;
W0=45.;
H1=75.;
gamma=0.75;
rE=0.9;
rW=0.75;
em=.033333333;
et=.25;
Co=7.E6;
dY=2.;
dT=.5;
var_b0=0.1;
month_max_b0=5;
mu_0=1.02;
L=1.5E7;
betta=5.96E-5;
Lo=sqrt(Co/betta);
Yn=dY*Lo;
tau1=L*Yn*Yn*betta/(Co*Co);
tau2=L/Co;
A=exp(-.05*dY*dY)*(1+.1*dY*dY)/(Yn*Yn);
em21=exp(    -tau2*em);
em22=exp(-0.5*tau2*em);
em11=exp(    -tau1*em);
em12=exp(-0.5*tau1*em);
tau_K=tau2;
tau_R=tau1;

% The parameters for the T_sub calculation
% ----------------------------------------

T1=28.;
T2=-40.;
b1=1./(80.);
b2=1./(33.);
H0=47.;
kappa=0.5*( b1*T1/(cosh(b1*H0)*cosh(b1*H0)) -b2*T2/(cosh(b2*H0)*cosh(b2*H0)) );
b_plus=T1*(1.-tanh(b1*H0));
b_minus=T2*(1.-tanh(b2*H0));
a_plus=1;
a_minus=1;
h_plus=b_plus*(a_plus-1)/a_plus;
h_minus=-b_minus*(a_minus-1)/a_minus;

% Solve the equation
% ------------------
for ii=i_start:i_end
  hE_coeff_a=0.472279;
  hE_coeff_b=-7.25359;
  hE_coeff_c=10.0873;

  % calculate h:
  mu(ii)=mu_0*(1.0+var_b0*cos(2*pi*(ii*delta_t-month_max_b0)/12));
  h(ii)= +em21*rW*em11*rE*h(ii-round((tau_K+tau_R)/delta_t))...
      -em22*rW*dT*em11*tau1*b0/(betta*ru)*A*mu(ii-round((tau_K+tau_R/2)/delta_t))*T(ii-round((tau_K+tau_R/2)/delta_t)) ...
      +em12*b0/(ru*Co)*dT*tau2*mu(ii-round((tau_K/2)/delta_t))*T(ii-round((tau_K/2)/delta_t));
  
  % hyperbolic function for Tsub (Munnich,Cane & Zebiak 1991):
  hh=kappa*h(ii-1);
  if(hh < h_minus)
    Tsub=-b_minus+b_minus/a_minus*(tanh(a_minus/b_minus *(hh-h_minus))+1);
  elseif (hh >=h_minus & hh < h_plus)
    Tsub=hh;
  else
    Tsub=b_plus+b_plus/a_plus*(tanh(a_plus/b_plus*(hh-h_plus))-1);
  end
  % advance T in time:
  T(ii)=T(ii-2)+2*delta_t*(-et*T(ii-1)-gamma*W0/H1*(T(ii-1)-Tsub));
  T(ii-1)=T(ii-1)+Robert*(T(ii)-2*T(ii-1)+T(ii-2));
end

time=[delta_t/12:delta_t/12:delta_t*i_end/12];
is=i_end-round(50*12/delta_t);
plot(time(is:i_end)-time(is),T(is:i_end))
title('Temperature time series')
xlabel('years')

figure
comet(h(is+1:1/delta_t:i_end),h(is-12/delta_t+1:1/delta_t:i_end-12/delta_t))
title('phase space trajectory')

figure
h=plot(h(i_start-12/delta_t:12/delta_t:i_end-12/delta_t),h(i_start:12/delta_t:i_end),'.');
set(h,'markersize',2)
title('phase space attractor')