function ds=l95ode(t,s,params)
%L95ODE The full two-tier lorenz 95 system
% ydot=l95ode(t,s,params)

% Leutbecher M, 2010, Predictability and ensemble forecasting with
% Lorenz-95 systems, lecture notes, ECMWF meteorological training
% course on Predictability, diagnostics and seasonal
% forecasting, http://www.ecmwf.int/.

if isempty(params.gufun) % the full model

  F = params.F;
  h = params.h;
  c = params.c;
  b = params.b;
  Fy = params.Fy;

  J = params.J;
  K = params.K;

  x = s(1:K);
  y = s(K+1:end);
  dx = zeros(K,1);
  dy = zeros(J*K,1);

  % 1) the slow part
  dx(1) = -x(end-1)*x(end)+x(end)*x(2)-x(1)+F-h*c/b*sum(y(1:J));
  dx(2) = -x(end)*x(1)+x(1)*x(3)-x(2)+F-h*c/b*sum(y(J+1:2*J));

  for k=3:K-1
    dx(k) = -x(k-2)*x(k-1)+x(k-1)*x(k+1)-x(k)+F-h*c/b*sum(y(J*(k-1)+1:J*k));
  end
  dx(k+1) = -x(end-2)*x(end-1)+x(end-1)*x(1)-x(end)+F-h*c/b*sum(y(J*k+1:J*(k+1)));

  % 2) the fast part
  dy(1) = -c*b*y(2)*(y(3)-y(end))-c*y(1)+c/b*Fy+h*c/b*x(1);
  for j=2:J*K-2
    dy(j) = -c*b*y(j+1)*(y(j+2)-y(j-1))-c*y(j)+c/b*Fy+h*c/b*x(1+floor((j-1)/J));
  end

  dy(end-1) = -c*b*y(end)*(y(1)-y(end-2))-c*y(end-1)+c/b*Fy+h*c/b*x(end);
  dy(end) = -c*b*y(1)*(y(2)-y(end-1))-c*y(end)+c/b*Fy+h*c/b*x(end);

  ds=[dx;dy];

else % the parameterized simpler model

  F = params.F;
  K =length(s);

  x = s;
  dx=zeros(length(x),1);

  gu = feval(params.gufun,x,params.gupars);  % deterministic parameterization
  eta = params.eta;   % stochastic forcing

  dx(1) = -x(end-1)*x(end)+x(end)*x(2)-x(1)+F-gu(1)+eta(1);
  dx(2) = -x(end)*x(1)+x(1)*x(3)-x(2)+F-gu(2)+eta(2);

  for k=3:K-1
    dx(k) = -x(k-2)*x(k-1)+x(k-1)*x(k+1)-x(k)+F-gu(k)+eta(k);
  end
  dx(k+1) = -x(end-2)*x(end-1)+x(end-1)*x(1)-x(end)+F-gu(k+1)+eta(k+1);
  ds=dx;

end

Published with MATLAB® 7.12