hydro-topmodel 0.1.0

//! TOPMODEL(Beven & Kirkby 1979;1995 Fortran 版)——忠实分布式实现。
//!
//! Clean-room 从 R `topmodel` 包 v0.7.5 的 C 数值核心(`core_topmodel.c`/`get_f.c`/
//! `param_init.c`/`misc.c`)移植:算法(Beven 公式)不可版权,C 源仅作行为 oracle。
//! 11 参数 + TWI 直方图 + 河道延迟;单位 mm / hour(与 R 参考一致,便于 bit-exact 校验)。
//! 交叉验证:`tests/topmodel_crosscheck.rs`(待加)vs R `topmodel` on huagrahuma 数据集。
//!
//! 关键方程(eq 号同 Beven HESS 2021 / R 源注释):
//! - 基流退水 `qs = qss·exp(-S_mean/m)`(eq 6.33),`qss = exp(lnTe + ln(dt) - λ)`
//! - 局部亏缺 `S[j] = S_mean + m·(λ - atb[j])`(eq 18.8)
//! - 非饱和排水 `qv = Suz/(S·td)·dt`(eq 6.26)
//! - 下渗 Green-Ampt/Morel-Seytoux(Newton-Raphson,`get_f`)
//! - 路由:河道延迟 `Qt[k] += qt·Ad[j]`,k = i + j + ndelay

use hydro_core::{Forcing, HydroModel};
use serde::{Deserialize, Serialize};

const ZERO: f64 = 0.0000001;

// ── 流域几何(TWI 直方图 + 河道延迟)──

/// TWI 直方图:`atb` = ln(a/tanβ) 值(**降序**),`Aatb_r` = 累积面积比例(到该类下限)。
#[derive(Clone, Debug, Serialize, Deserialize, Default)]
pub struct TopidxHistogram {
    pub atb: Vec<f64>,
    pub Aatb_r: Vec<f64>,
}

/// 河道延迟:`d` = 河道距离(增序),`Ad_r` = 对应累积面积比例。
#[derive(Clone, Debug, Serialize, Deserialize, Default)]
pub struct ChannelDelay {
    pub d: Vec<f64>,
    pub Ad_r: Vec<f64>,
}

/// TOPMODEL 11 参数 + 流域几何。
#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct TopmodelParams {
    pub qs0: f64,    // 初始基流 [m]
    pub lnTe: f64,   // areal avg ln(T0) [m²/h]
    pub m: f64,      // 饱和亏缺衰减 [m]
    pub Sr0: f64,    // 初始根区亏缺 [m]
    pub Srmax: f64,  // 最大根区亏缺 [m]
    pub td: f64,     // 非饱和带时间延迟 [h/m]
    pub vch: f64,    // 河道波速 [m/h]
    pub vr: f64,     // 坡面波速 [m/h]
    pub K0: f64,     // 饱和导水率 [m/h]
    pub CD: f64,     // 毛管驱动 [m]
    pub dt: f64,     // 时步 [h]
    #[serde(default)]
    pub topidx: Option<TopidxHistogram>,
    #[serde(default)]
    pub channel: Option<ChannelDelay>,
    #[serde(default = "default_area")]
    pub area_km2: f64,
}

fn default_area() -> f64 { 1000.0 }

impl Default for TopmodelParams {
    fn default() -> Self {
        // 默认:单类 TWI 直方图(集总退化,供平台多模型对比;真实运行需提供 topidx)。
        Self {
            qs0: 0.001, lnTe: 5.0, m: 0.01, Sr0: 0.0, Srmax: 0.05,
            td: 30.0, vch: 100.0, vr: 100.0, K0: 3.0, CD: 1.0, dt: 1.0,
            topidx: Some(TopidxHistogram { atb: vec![0.0], Aatb_r: vec![1.0] }),
            channel: None,
            area_km2: 1000.0,
        }
    }
}

/// 完整运行输出(逐时步序列,mm;供 oracle 逐变量比对)。
#[derive(Clone, Debug, Default)]
pub struct TopmodelOutput {
    pub Qt: Vec<f64>,      // 出流(路由后)[mm/时步]
    pub qs: Vec<f64>,      // 基流
    pub qo: Vec<f64>,      // 地表产流(含 fex)
    pub S_mean: Vec<f64>,  // 平均饱和亏缺
    pub f: Vec<f64>,       // 下渗
    pub fex: Vec<f64>,     // 超渗
    pub Ea: Vec<f64>,      // 实际蒸散发(面平均)
}

// ── 下渗:Green-Ampt/Morel-Seytoux(Newton-Raphson;跨时步状态)──

/// 下渗计算器(对应 get_f.c 的 static 变量 cumf/f_/pt/cnst/ponding)。
#[derive(Clone)]
struct Infiltration {
    cumf: f64,
    f_: f64,
    pt: f64,
    cnst: f64,
    ponding: bool,
}

impl Infiltration {
    fn new() -> Self {
        Self { cumf: 0.0, f_: 0.0, pt: 0.0, cnst: 0.0, ponding: false }
    }
    fn reset(&mut self) {
        self.cumf = 0.0; self.f_ = 0.0; self.pt = 0.0; self.cnst = 0.0; self.ponding = false;
    }

    /// 对应 get_f(t, R, C, K0, m, dt)。R = 雨强 [m/h],t = 累计时间 [h]。
    fn get_f(&mut self, t: f64, r: f64, c: f64, k0: f64, m: f64, dt: f64) -> f64 {
        const TOLERANCE: f64 = 0.00001;
        const MAXITER: usize = 2000;
        const NTERMS: usize = 10;

        if t / dt == 1.0 {
            self.reset();
        }
        if r <= 0.0 {
            self.cumf = 0.0; self.ponding = false; self.f_ = 0.0; self.pt = 0.0;
            return 0.0;
        }

        let mut f1 = 0.0;
        if !self.ponding {
            if self.cumf > 0.0 {
                f1 = self.cumf;
                let mut r2 = -k0 / m * (c + f1) / (1.0 - (f1 / m).exp());
                if r > r2 {
                    self.f_ = self.cumf;
                    self.pt = t - dt;
                    self.ponding = true;
                    // goto cont1
                    self.cnst = Self::compute_cnst(self.f_, c, m);
                    self.f_ += r * (t - self.pt) / 2.0;
                    // fall through to Newton-Raphson below
                } else {
                    // no ponding path continues below
                    let f2 = self.cumf + r * dt;
                    r2 = -k0 / m * (c + f2) / (1.0 - (f2 / m).exp());
                    if f2 == 0.0 || r < r2 {
                        let f = r;
                        self.cumf += f * dt;
                        self.ponding = false;
                        return f;
                    }
                    self.f_ = self.cumf + r2 * dt;
                    let mut f2m = f2;
                    let mut f1m = f1;
                    let mut i = 0;
                    while i < MAXITER {
                        r2 = -k0 / m * (c + self.f_) / (1.0 - (self.f_ / m).exp());
                        let diff;
                        if r2 > r {
                            f1m = self.f_;
                            self.f_ = (self.f_ + f2m) / 2.0;
                            diff = self.f_ - f1m;
                        } else {
                            f2m = self.f_;
                            self.f_ = (self.f_ + f1m) / 2.0;
                            diff = self.f_ - f2m;
                        }
                        if diff.abs() < TOLERANCE { break; }
                        i += 1;
                    }
                    if i == MAXITER { return -9999.0; }
                    self.pt = t - dt + (self.f_ - self.cumf) / r;
                    if self.pt > t {
                        let f = r;
                        self.cumf += f * dt;
                        self.ponding = false;
                        return f;
                    }
                    self.cnst = Self::compute_cnst(self.f_, c, m);
                    self.f_ += r * (t - self.pt) / 2.0;
                    self.ponding = true;
                }
            } else {
                // cumf == 0 path
                let f2 = self.cumf + r * dt;
                let r2 = -k0 / m * (c + f2) / (1.0 - (f2 / m).exp());
                if f2 == 0.0 || r < r2 {
                    let f = r;
                    self.cumf += f * dt;
                    self.ponding = false;
                    return f;
                }
                self.f_ = self.cumf + r2 * dt;
                let mut f2m = f2;
                let mut f1m = f1;
                let mut i = 0;
                while i < MAXITER {
                    let r2 = -k0 / m * (c + self.f_) / (1.0 - (self.f_ / m).exp());
                    let diff;
                    if r2 > r {
                        f1m = self.f_;
                        self.f_ = (self.f_ + f2m) / 2.0;
                        diff = self.f_ - f1m;
                    } else {
                        f2m = self.f_;
                        self.f_ = (self.f_ + f1m) / 2.0;
                        diff = self.f_ - f2m;
                    }
                    if diff.abs() < TOLERANCE { break; }
                    i += 1;
                }
                if i == MAXITER { return -9999.0; }
                self.pt = t - dt + (self.f_ - self.cumf) / r;
                if self.pt > t {
                    let f = r;
                    self.cumf += f * dt;
                    self.ponding = false;
                    return f;
                }
                self.cnst = Self::compute_cnst(self.f_, c, m);
                self.f_ += r * (t - self.pt) / 2.0;
                self.ponding = true;
            }
        }

        // Newton-Raphson(ponding 发生后)
        let mut i = 0;
        while i < MAXITER {
            let fc = self.f_ + c;
            let mut sum = 0.0;
            let mut factorial = 1.0;
            for j in 1..=NTERMS {
                factorial *= j as f64;
                sum += (fc / m).powi(j as i32) / (j as f64 * factorial);
            }
            let g1 = -((fc.ln() - (fc.ln() + sum) / (c / m).exp() - self.cnst) / (k0 / m)) - (t - self.pt);
            let g2 = ((self.f_ / m).exp() - 1.0) / (fc * k0 / m);
            let diff = -g1 / g2;
            self.f_ += diff;
            if diff.abs() < TOLERANCE { break; }
            i += 1;
        }
        if i == MAXITER { return -9999.0; }

        if self.f_ - self.cumf < r * dt {
            let f = (self.f_ - self.cumf) / dt;
            self.cumf = self.f_;
            self.f_ += f * dt;
            f
        } else {
            let f = r;
            self.cumf += f * dt;
            self.ponding = false;
            self.pt = 0.0;
            f
        }
    }
}

impl Infiltration {
    fn compute_cnst(f_: f64, c: f64, m: f64) -> f64 {
        const NTERMS: usize = 10;
        let fc = f_ + c;
        let mut cnst = 0.0;
        let mut factorial = 1.0;
        for j in 1..=NTERMS {
            factorial *= j as f64;
            cnst += (fc / m).powi(j as i32) / (j as f64 * factorial);
        }
        fc.ln() - (fc.ln() + cnst) / (c / m).exp()
    }
}

// ── 派生量:λ(地形指数面平均)+ Ad(河道延迟面积分布)──

/// λ = areal integral of ln(a/tanβ) = Σ Aatb_r[i]·(atb[i]+atb[i-1])/2(对应 get_lambda)。
fn get_lambda(atb: &[f64], aatb_r: &[f64]) -> f64 {
    let n = atb.len().min(aatb_r.len());
    let mut ret = 0.0;
    for i in 1..n {
        ret += aatb_r[i] * (atb[i] + atb[i - 1]) / 2.0;
    }
    ret
}

/// 河道延迟 → (tch, ndelay, nreach, Ad)。对应 get_Ad。
fn compute_ad(d: &[f64], ad_r: &[f64], vch_dt: f64, vr_dt: f64) -> (Vec<f64>, usize, usize, Vec<f64>) {
    let nch = d.len().min(ad_r.len());
    if nch == 0 {
        return (Vec::new(), 0, 0, Vec::new());
    }
    let mut tch = vec![0.0; nch];
    tch[0] = d[0] / vch_dt;
    for i in 1..nch {
        tch[i] = tch[0] + (d[i] - d[0]) / vr_dt;
    }
    let mut nreach = tch[nch - 1] as usize;
    if (nreach as f64) < tch[nch - 1] { nreach += 1; }
    let ndelay = tch[0] as usize;
    nreach = nreach.saturating_sub(ndelay);
    if nreach == 0 {
        return (tch, ndelay, 0, Vec::new());
    }
    let mut ad = vec![0.0; nreach];
    for i in 0..nreach {
        let t = (ndelay + i + 1) as f64;
        if t > tch[nch - 1] {
            ad[i] = 1.0;
        } else {
            for j in 1..nch {
                if t <= tch[j] {
                    ad[i] = ad_r[j - 1] + (ad_r[j] - ad_r[j - 1]) * (t - tch[j - 1]) / (tch[j] - tch[j - 1]);
                    break;
                }
            }
        }
    }
    // 差分:累积 → 单步
    let mut a1 = ad[0];
    for i in 1..nreach {
        let a2 = ad[i];
        ad[i] = a2 - a1;
        a1 = a2;
    }
    (tch, ndelay, nreach, ad)
}

// ── 派生状态(对应 param_init)──

struct Derived {
    lambda: f64,
    qss: f64,
    ndelay: usize,
    nreach: usize,
    ad: Vec<f64>,
    nidxclass: usize,
    atb: Vec<f64>,
    aatb_r: Vec<f64>,
}

impl Derived {
    fn from_params(p: &TopmodelParams) -> Self {
        let (atb, aatb_r) = match &p.topidx {
            Some(h) if !h.atb.is_empty() => (h.atb.clone(), h.Aatb_r.clone()),
            _ => (vec![0.0], vec![1.0]),
        };
        let nidxclass = atb.len();
        let lambda = get_lambda(&atb, &aatb_r);
        let ln_te_dt = p.lnTe + p.dt.ln();
        let qss = (ln_te_dt - lambda).exp();
        let (ndelay, nreach, ad) = match &p.channel {
            Some(ch) if !ch.d.is_empty() => {
                let vch_dt = p.vch * p.dt;
                let vr_dt = p.vr * p.dt;
                let (_tch, nd, nr, ad) = compute_ad(&ch.d, &ch.Ad_r, vch_dt, vr_dt);
                (nd, nr, ad)
            }
            _ => (0, 0, Vec::new()),
        };
        Self { lambda, qss, ndelay, nreach, ad, nidxclass, atb, aatb_r }
    }
}

// ── 忠实批量核心:run_topmodel_full(供 oracle 比对)──

/// 完整运行 TOPMODEL(逐时步循环,对应 R 的 run_topmodel 调用序列)。
/// 输入:参数 + 雨量(mm/时步)+ ETp(mm/时步);输出:全序列(mm)。
pub fn run_topmodel_full(p: &TopmodelParams, rain: &[f64], etp: &[f64]) -> TopmodelOutput {
    let n = rain.len().max(etp.len());
    let d = Derived::from_params(p);
    let qs0_dt = p.qs0 * p.dt;

    let mut out = TopmodelOutput {
        Qt: vec![0.0; n], qs: vec![0.0; n], qo: vec![0.0; n],
        S_mean: vec![0.0; n], f: vec![0.0; n], fex: vec![0.0; n], Ea: vec![0.0; n],
    };
    let mut srz = vec![p.Sr0; d.nidxclass]; // 根区亏缺
    let mut suz = vec![0.0; d.nidxclass];   // 非饱和蓄水
    let mut infl = Infiltration::new();

    // S_mean[0] = -m·log(qs0_dt/qss)(对应 param_init)
    out.S_mean[0] = if d.qss > 0.0 && qs0_dt > 0.0 {
        -p.m * (qs0_dt / d.qss).ln()
    } else { 0.0 };

    // Qt 初值(对应 param_init):Qt[0..ndelay]=qs0;reach 段退水 qs0·(1-ΣAd)
    for i in 0..n.min(d.ndelay) {
        out.Qt[i] = qs0_dt;
    }
    {
        let mut a = 0.0;
        for i in 0..d.nreach {
            a += d.ad[i];
            let k = d.ndelay + i;
            if k < n { out.Qt[k] = qs0_dt * (1.0 - a); }
        }
    }

    for i in 0..n {
        let r = rain.get(i).copied().unwrap_or(0.0);
        let e = etp.get(i).copied().unwrap_or(0.0);

        // 下渗 + 超渗
        let t = (i as f64 + 1.0) * p.dt;
        let mut fi = p.dt * infl.get_f(t, r / p.dt, p.CD, p.K0, p.m, p.dt);
        if fi < 0.0 { fi = r; }
        out.f[i] = fi;
        out.fex[i] = r - fi;

        // 基流
        out.qs[i] = d.qss * (-out.S_mean[i] / p.m).exp();

        let mut qo_total = 0.0;
        let mut qv_total = 0.0;
        let mut ea_total = 0.0;
        let mut ex_prev = 0.0;

        for j in 0..d.nidxclass {
            let aatb_local = (d.aatb_r[j]
                + if j < d.nidxclass - 1 { d.aatb_r[j + 1] } else { 0.0 }) / 2.0;

            // 局部亏缺 eq 18.8
            let mut s = out.S_mean[i] + p.m * (d.lambda - d.atb[j]);
            if s < 0.0 { s = 0.0; }

            // 根区:雨进根区(亏减少),溢出补非饱和
            srz[j] -= fi;
            if srz[j] < 0.0 {
                suz[j] -= srz[j];
                srz[j] = 0.0;
            }

            // 饱和溢流 ex(Suz 超过 S)
            let mut ex = 0.0;
            if suz[j] > s {
                ex = suz[j] - s;
                suz[j] = s;
            }

            // 非饱和排水 eq 6.26
            let mut qv = 0.0;
            if s > 0.0 {
                qv = suz[j] / (s * p.td) * p.dt;
                if qv > suz[j] { qv = suz[j]; }
                suz[j] -= qv;
                if suz[j] < ZERO { suz[j] = 0.0; }
                qv *= aatb_local;
            }
            qv_total += qv;

            // ET eq 6.27(从根区抽取,亏缺增加)
            let mut ea = 0.0;
            if e > 0.0 {
                ea = e * (1.0 - srz[j] / p.Srmax);
                if ea > p.Srmax - srz[j] { ea = p.Srmax - srz[j]; }
            }
            srz[j] += ea;
            ea_total += aatb_local * ea;

            // 饱和面产流 qo(相邻类 ex;分支1用原始 Aatb_r[j],分支2用 local)
            if j > 0 {
                let qo = if ex > 0.0 {
                    d.aatb_r[j] * (ex_prev + ex) / 2.0
                } else if ex_prev > 0.0 {
                    aatb_local * ex_prev / (ex_prev - ex) * ex_prev / 2.0
                } else { 0.0 };
                qo_total += qo;
            }
            ex_prev = ex;
        }

        out.qo[i] = qo_total + out.fex[i];
        let qt = out.qo[i] + out.qs[i];

        // S_mean 更新(qs 出、qv 入)
        out.S_mean[i] += out.qs[i] - qv_total;
        if i + 1 < n { out.S_mean[i + 1] = out.S_mean[i]; }
        out.Ea[i] = ea_total;

        // 路由:Qt[k] += qt·Ad[j]
        if d.nreach > 0 {
            for j in 0..d.nreach {
                let k = i + j + d.ndelay;
                if k > n - 1 { break; }
                out.Qt[k] += qt * d.ad[j];
            }
        } else {
            // 无河道延迟:即时出流
            out.Qt[i] += qt;
        }
    }
    out
}


// ── HydroModel trait 包装(平台集成;逐时步)──

pub struct TopmodelModel {
    params: TopmodelParams,
    derived: Derived,
    srz: Vec<f64>,
    suz: Vec<f64>,
    infl: Infiltration,
    s_mean: f64,
    qt_buf: Vec<f64>, // 路由缓冲(Qt[i])
    step_i: usize,
    discharge_m3s: f64,
}

fn mm_to_m3s(mm: f64, area_km2: f64, dt_h: f64) -> f64 {
    mm * area_km2 / (dt_h * 3.6)
}

impl HydroModel for TopmodelModel {
    type Params = TopmodelParams;

    fn new(params: Self::Params) -> Self {
        let derived = Derived::from_params(&params);
        let qs0_dt = params.qs0 * params.dt;
        let s_mean0 = if derived.qss > 0.0 && qs0_dt > 0.0 {
            -params.m * (qs0_dt / derived.qss).ln()
        } else { 0.0 };
        Self {
            srz: vec![params.Sr0; derived.nidxclass],
            suz: vec![0.0; derived.nidxclass],
            infl: Infiltration::new(),
            s_mean: s_mean0,
            derived,
            params,
            qt_buf: Vec::new(),
            step_i: 0,
            discharge_m3s: 0.0,
        }
    }

    fn step(&mut self, f: &Forcing, dt_h: f64) {
        let p = &self.params;
        let i = self.step_i;
        let r = f.p_mm.max(0.0);
        let e = f.pet_mm.max(0.0);
        let nidx = self.derived.nidxclass;

        let t = (i as f64 + 1.0) * p.dt;
        let mut fi = p.dt * self.infl.get_f(t, r / p.dt, p.CD, p.K0, p.m, p.dt);
        if fi < 0.0 { fi = r; }
        let fex = r - fi;

        let qs = self.derived.qss * (-self.s_mean / p.m).exp();

        let mut qo_total = 0.0;
        let mut qv_total = 0.0;
        let mut ea_total = 0.0;
        let mut ex_prev = 0.0;

        for j in 0..nidx {
            let aatb_local = (self.derived.aatb_r[j]
                + if j < nidx - 1 { self.derived.aatb_r[j + 1] } else { 0.0 }) / 2.0;
            let mut s = self.s_mean + p.m * (self.derived.lambda - self.derived.atb[j]);
            if s < 0.0 { s = 0.0; }

            self.srz[j] -= fi;
            if self.srz[j] < 0.0 {
                self.suz[j] -= self.srz[j];
                self.srz[j] = 0.0;
            }
            let mut ex = 0.0;
            if self.suz[j] > s {
                ex = self.suz[j] - s;
                self.suz[j] = s;
            }
            let mut qv = 0.0;
            if s > 0.0 {
                qv = self.suz[j] / (s * p.td) * p.dt;
                if qv > self.suz[j] { qv = self.suz[j]; }
                self.suz[j] -= qv;
                if self.suz[j] < ZERO { self.suz[j] = 0.0; }
                qv *= aatb_local;
            }
            qv_total += qv;

            let mut ea = 0.0;
            if e > 0.0 {
                ea = e * (1.0 - self.srz[j] / p.Srmax);
                if ea > p.Srmax - self.srz[j] { ea = p.Srmax - self.srz[j]; }
            }
            self.srz[j] += ea;
            ea_total += aatb_local * ea;

            // 饱和面产流 qo(相邻类 ex)
            if j > 0 {
                let qo = if ex > 0.0 {
                    self.derived.aatb_r[j] * (ex_prev + ex) / 2.0
                } else if ex_prev > 0.0 {
                    aatb_local * ex_prev / (ex_prev - ex) * ex_prev / 2.0
                } else { 0.0 };
                qo_total += qo;
            }
            ex_prev = ex;
        }

        qo_total += fex;
        let qt = qo_total + qs;
        self.s_mean += qs - qv_total;

        // 路由到 Qt 缓冲
        if self.derived.nreach > 0 {
            for j in 0..self.derived.nreach {
                let k = i + j + self.derived.ndelay;
                if k >= self.qt_buf.len() { self.qt_buf.resize(k + 1, 0.0); }
                self.qt_buf[k] += qt * self.derived.ad[j];
            }
        } else {
            if i >= self.qt_buf.len() { self.qt_buf.resize(i + 1, 0.0); }
            self.qt_buf[i] += qt;
        }

        let qt_i = self.qt_buf.get(i).copied().unwrap_or(0.0);
        self.discharge_m3s = mm_to_m3s(qt_i, p.area_km2, dt_h);
        self.step_i += 1;
    }

    fn discharge(&self) -> f64 { self.discharge_m3s }
    fn state(&self) -> serde_json::Value {
        serde_json::json!({ "S_mean": self.s_mean, "step": self.step_i })
    }
    fn reset(&mut self) {
        let p = self.params.clone();
        *self = Self::new(p);
    }
    fn name(&self) -> &'static str { "TOPMODEL" }
    fn params(&self) -> &Self::Params { &self.params }
    fn params_mut(&mut self) -> &mut Self::Params { &mut self.params }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn lambda_single_class() {
        // 单类 atb=[0], Aatb_r=[1] → λ=0(get_lambda 从 i=1 起,无项)
        assert!(get_lambda(&[0.0], &[1.0]).abs() < 1e-12);
    }

    #[test]
    fn lambda_two_class() {
        // atb=[2,0], Aatb_r=[0.3,1.0] → λ = Aatb_r[1]*(atb[1]+atb[0])/2 = 1.0*(0+2)/2 = 1.0
        let l = get_lambda(&[2.0, 0.0], &[0.3, 1.0]);
        assert!((l - 1.0).abs() < 1e-12, "λ={}", l);
    }

    #[test]
    fn rain_produces_flow() {
        let mut p = TopmodelParams::default();
        p.m = 0.01; p.lnTe = (5.0_f64).ln(); p.qs0 = 0.01; p.Srmax = 0.05; p.Sr0 = 0.0;
        p.K0 = 0.0; p.CD = 0.0; p.dt = 1.0;
        let mut m = TopmodelModel::new(p);
        for _ in 0..10 { m.step(&Forcing { p_mm: 30.0, pet_mm: 0.0, t_c: 20.0 }, 1.0); }
        assert!(m.discharge() > 0.0, "rain should produce flow: {}", m.discharge());
    }

    #[test]
    fn no_mass_creation_constant_rain() {
        let params = TopmodelParams::default();
        let area = params.area_km2;
        let mut m = TopmodelModel::new(params);
        let n = 200;
        let mut sum_q_mm = 0.0;
        for _ in 0..n {
            m.step(&Forcing { p_mm: 5.0, pet_mm: 0.0, t_c: 20.0 }, 1.0);
            sum_q_mm += m.discharge() * 1.0 * 3.6 / area;
        }
        let sum_rain = 5.0 * n as f64;
        assert!(sum_q_mm <= sum_rain * 1.05, "出流>降雨(质量凭空):{:.1}>{:.1}", sum_q_mm, sum_rain);
    }

    #[test]
    fn dyn_dispatch() {
        let mut m: Box<dyn hydro_core::DynHydroModel> = Box::new(TopmodelModel::new(TopmodelParams::default()));
        m.step(&Forcing { p_mm: 50.0, pet_mm: 1.0, t_c: 20.0 }, 1.0);
        assert!(m.discharge() >= 0.0);
        assert_eq!(m.name(), "TOPMODEL");
    }
}