function [Tmm,hmm]=reosc_model(T0,h0,a11m,a12m,a21m,a22m,amp_t,amp_h,wwe,tdim)

% numerical parameters
mdim  = 12*tdim;                                 % last month
dt    = 40*60;                                    % length of one time step
sec   = 30*24*60*60;                              % seconds per month
tsm   = sec/dt;                                   % time steps per month

% model parameters in sec
a11=a11m/sec; 
a12=a12m/sec; 
a21=a21m/sec; 
a22=a22m/sec; 

% westerly wind burst scale
whscl = -0.008/0.021/2; % the relative scale for h eq. = a22/a12/2  eq. [3] yu et. al 2016


% noise parameter per time step
a1_noise  = exp(-dt/(3*24*3600));
amp_scl   = sqrt((2.5/sec)^2-(a1_noise*2.5/sec)^2);

% set initial values
T = T0;      % initial temperature anomaly
h = h0;      % initial thermocline depth anomaly
xiT = 0.0; 
xih = 0.0;

% TIME LOOP to integrate model
for t = 1:mdim                                      % loop over months
    m = mod(t-1,12)+1;                              % index of month
    for mo = 1:tsm                                  % loop over all points in time of the current month
        xiT = a1_noise*xiT+amp_scl*amp_t(m)*randn()+wwe(t); 
        xih = a1_noise*xih+amp_scl*amp_h(m)*randn()+whscl*wwe(t); 
        T   = T+(a11(m)*T+a12(m)*h+xiT).*dt; 
        h   = h+(a21(m)*T+a22(m)*h+xih).*dt; 
        Tx(mo) = T; hx(mo) = h;                    % vector of monthly results
    end
    Tmm(t) = mean(Tx);                             % compute mean of monthly temperatures
    hmm(t) = mean(hx);                             % compute mean of monthly temperatures
end

return



