function [s, L, EVT] = soilMoistureDynamics(rainfall, dt, s0, sw, sstar, Ks, b, ETsw, ETmax, n, Zr)
% SOIL_MOISTURE_DYNAMICS Bilancio di massa P - ET - L = Delta(s)*n*Zr per ogni timestep.
%
% Convenzione temporale: rainfall(i) cade nell'intervallo (i-1, i] e produce
% l'aggiornamento s(i-1) -> s(i), con L(i) ed EVT(i) riferiti allo stesso intervallo.
% s(1) e' lo stato INIZIALE (prima di qualunque pioggia), quindi L(1)=EVT(1)=0
% per definizione (nessun intervallo precede il primo dato).
%
% Input
% rainfall [mm]  = pioggia per timestep (rainfall(1) e' lo stato iniziale, vedi sopra)
% dt     [h]     = passo temporale
% s0     [-]     = soil moisture iniziale
% sw     [-]     = soil moisture al wilting point
% sstar  [-]     = soil moisture per ET massima
% Ks     [mm/h]  = conducibilita' idraulica satura
% b      [-]     = esponente Clapp-Hornberger
% ETsw   [mm/h]  = ET potenziale a sw (scalare o vettore length(rainfall))
% ETmax  [mm/h]  = ET potenziale a sstar (scalare o vettore length(rainfall))
% n      [-]     = porosita'
% Zr     [mm]    = profondita' radici
%
% Output
% s      [-]     = soil moisture per timestep (s(1) = s0)
% L      [mm]    = leaching per timestep (L(1) = 0)
% EVT    [mm]    = evapotraspirazione per timestep (EVT(1) = 0)

% disp('Soil moisture dynamics');

P = rainfall(:)';
ntot = length(P);

s   = NaN(1, ntot);
L   = zeros(1, ntot);
EVT = zeros(1, ntot);
s(1) = s0;

for i = 2:ntot
    if isscalar(ETmax)
        etMax_i = ETmax;
    else
        etMax_i = ETmax(i);
    end
    if isscalar(ETsw)
        etSw_i = ETsw;
    else
        etSw_i = ETsw(i);
    end

    [s(i), L(i), EVT(i)] = step_rk4(s(i-1), dt, P(i), n, Zr, sstar, sw, etMax_i, etSw_i, Ks, b);
end

%% Diagnostics
time_sstar  = sum(s >= sstar) / length(s) * 100;
time_sw     = sum(s <= sw) / length(s) * 100;
time_middle = 100 - time_sw - time_sstar;

Ptot    = sum(P);
Ltot    = sum(L);
ETtot   = sum(EVT);
storage = (s(end) - s(1)) * n * Zr;

balance_error = Ptot - ETtot - Ltot - storage;

fprintf('\n');
fprintf('SOIL MOISTURE DIAGNOSTICS\n');
fprintf('-------------------------\n');
fprintf('Time s <= sw:        %.2f %%\n', time_sw);
fprintf('Time sw < s < s*:    %.2f %%\n', time_middle);
fprintf('Time s >= s*:        %.2f %%\n', time_sstar);
fprintf('\n');
fprintf('Rainfall:            %.2f mm\n', Ptot);
fprintf('Evapotranspiration:  %.2f mm (%.1f %%)\n', ETtot, ETtot / Ptot * 100);
fprintf('Leaching:            %.2f mm (%.1f %%)\n', Ltot, Ltot / Ptot * 100);
fprintf('Storage change:      %.2f mm (%.1f %%)\n', storage, storage / Ptot * 100);
fprintf('\n');
fprintf('Mass balance error:  %.6e mm\n', balance_error);

end
