% Riley (1946) phytoplankton model

function P = riley(p, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1)

P = NaN(24,1);
P(1) = 3.4;
Ph = p*I./(k.*z1).*(1 - exp(-k.*z1)).*OneMinusN.*OneMinusV;
R = R0*exp(r*T);
G = g*Z;

for t = 2:length(Z)
    P(t) = exp(15*(Ph(t) - R(t) - G(t)) + log(P(t-1)));
end

end