kongdd · May 29, 2026 07:46
diff --git a/index.html b/index.html
 <!DOCTYPE html>
 <html lang="en">
 <head>
 <meta charset="UTF-8">
 <meta name="viewport" content="width=device-width, initial-scale=1.0">
 <title>BEPS.jl — Wiki</title>
 <script src="https://cdn.jsdelivr.net/npm/marked@11.0.0/marked.min.js"></script>
 <script src="https://cdn.jsdelivr.net/npm/mermaid@11/dist/mermaid.min.js"></script>
 <style>

 *{margin:0;padding:0;box-sizing:border-box}
 :root{
  --bg:#ffffff;--sidebar-bg:#f8f9fb;--border:#e5e7eb;
  --text:#1e293b;--text-muted:#64748b;--primary:#2563eb;
  --primary-soft:#eff6ff;--hover:#f1f5f9;--code-bg:#f1f5f9;
  --radius:8px;--shadow:0 1px 3px rgba(0,0,0,.08);
 }
 body{font-family:-apple-system,BlinkMacSystemFont,'Segoe UI',Roboto,sans-serif;
  line-height:1.65;color:var(--text);background:var(--bg)}

 .layout{display:flex;min-height:100vh}
 .sidebar{width:280px;background:var(--sidebar-bg);border-right:1px solid var(--border);
  position:fixed;top:0;left:0;bottom:0;overflow-y:auto;padding:24px 16px;
  display:flex;flex-direction:column;z-index:10}
 .content{margin-left:280px;flex:1;padding:48px 64px;max-width:960px}

 .sidebar-header{margin-bottom:20px;padding-bottom:16px;border-bottom:1px solid var(--border)}
 .sidebar-title{font-size:16px;font-weight:700;color:var(--text);display:flex;align-items:center;gap:8px}
 .sidebar-title svg{flex-shrink:0}
 .sidebar-meta{font-size:11px;color:var(--text-muted);margin-top:6px}
 .nav-section{margin-bottom:2px}
 .nav-item{display:block;padding:7px 12px;border-radius:var(--radius);cursor:pointer;
  font-size:13px;color:var(--text);text-decoration:none;transition:all .15s;
  white-space:nowrap;overflow:hidden;text-overflow:ellipsis}
 .nav-item:hover{background:var(--hover)}
 .nav-item.active{background:var(--primary-soft);color:var(--primary);font-weight:600}
 .nav-item.overview{font-weight:600;margin-bottom:4px}
 .nav-children{padding-left:14px;border-left:1px solid var(--border);margin-left:12px}
 .nav-group-label{font-size:11px;font-weight:600;color:var(--text-muted);
  text-transform:uppercase;letter-spacing:.5px;padding:12px 12px 4px;user-select:none}
 .sidebar-footer{margin-top:auto;padding-top:16px;border-top:1px solid var(--border);
  font-size:11px;color:var(--text-muted);text-align:center}

 .content h1{font-size:28px;font-weight:700;margin-bottom:8px;line-height:1.3}
 .content h2{font-size:22px;font-weight:600;margin:32px 0 12px;padding-bottom:6px;border-bottom:1px solid var(--border)}
 .content h3{font-size:17px;font-weight:600;margin:24px 0 8px}
 .content h4{font-size:15px;font-weight:600;margin:20px 0 6px}
 .content p{margin:12px 0}
 .content ul,.content ol{margin:12px 0 12px 24px}
 .content li{margin:4px 0}
 .content a{color:var(--primary);text-decoration:none}
 .content a:hover{text-decoration:underline}
 .content blockquote{border-left:3px solid var(--primary);padding:8px 16px;margin:16px 0;
  background:var(--primary-soft);border-radius:0 var(--radius) var(--radius) 0;
  color:var(--text-muted);font-size:14px}
 .content code{font-family:'SF Mono',Consolas,'Courier New',monospace;font-size:13px;
  background:var(--code-bg);padding:2px 6px;border-radius:4px}
 .content pre{background:#1e293b;color:#e2e8f0;border-radius:var(--radius);padding:16px;
  overflow-x:auto;margin:16px 0}
 .content pre code{background:none;padding:0;font-size:13px;line-height:1.6;color:inherit}
 .content table{border-collapse:collapse;width:100%;margin:16px 0}
 .content th,.content td{border:1px solid var(--border);padding:8px 12px;text-align:left;font-size:14px}
 .content th{background:var(--sidebar-bg);font-weight:600}
 .content img{max-width:100%;border-radius:var(--radius)}
 .content hr{border:none;border-top:1px solid var(--border);margin:32px 0}
 .content .mermaid{margin:20px 0;text-align:center}

 .menu-toggle{display:none;position:fixed;top:12px;left:12px;z-index:20;
  background:var(--bg);border:1px solid var(--border);border-radius:var(--radius);
  padding:8px 12px;cursor:pointer;font-size:18px;box-shadow:var(--shadow)}
 @media(max-width:768px){
  .sidebar{transform:translateX(-100%);transition:transform .2s}
  .sidebar.open{transform:translateX(0);box-shadow:2px 0 12px rgba(0,0,0,.1)}
  .content{margin-left:0;padding:24px 20px;padding-top:56px}
  .menu-toggle{display:block}
 }
 .empty-state{text-align:center;padding:80px 20px;color:var(--text-muted)}
 .empty-state h2{font-size:20px;margin-bottom:8px;border:none}

 </style>
 </head>
 <body>
 <button class="menu-toggle" id="menu-toggle" aria-label="Toggle menu">&#9776;</button>
 <div class="layout">
 <nav class="sidebar" id="sidebar">
 <div class="sidebar-header">
 <div class="sidebar-title">
 <svg width="18" height="18" viewBox="0 0 24 24" fill="none" stroke="currentColor" stroke-width="2"><path d="M2 3h6a4 4 0 014 4v14a3 3 0 00-3-3H2z"/><path d="M22 3h-6a4 4 0 00-4 4v14a3 3 0 013-3h7z"/></svg>
 BEPS.jl
 </div>
 <div class="sidebar-meta" id="meta-info"></div>
 </div>
 <div id="nav-tree"></div>
 <div class="sidebar-footer">Generated by GitNexus</div>
 </nav>
 <main class="content" id="content">
 <div class="empty-state"><h2>Loading…</h2></div>
 </main>
 </div>
 <script>
 var PAGES = {"c-library-bindings":"# C Library Bindings\n\n# C Library Bindings (`src/clang/`)\n\nThis module provides Julia wrappers around `libbeps`, the compiled C library implementing the original BEPS physics. It enables the C execution path (`inter_prg_c`) while sharing data structures and dispatch interfaces with the pure-Julia path.\n\n## Module Structure\n\nFour files compose the module:\n\n| File | Role |\n|------|------|\n| `BEPS_c.jl` | Module entry point; top-level wrappers (`inter_prg_c`, `s_coszs`, `lai2`, rainfall stages, parameter readers) |\n| `SOIL_c.jl` | Soil-specific C wrappers; overloads Julia soil functions for `Soil_c` dispatch |\n| `module.jl` | Per-process physics wrappers: radiation, photosynthesis, aerodynamics, evaporation, sensible heat, surface temperature |\n| `snowpack_stage.jl` | Three-stage snowpack wrappers with ergonomic overloads for `Layer3`/`Layer2` types |\n| `struct_SOIL.jl` | `Soil_c` struct definition — C-layout-compatible mutable struct |\n\n## The `Soil_c` Struct\n\n`Soil_c` is the central state object for the C path, mirroring the C `struct Soil` memory layout exactly. It holds all 5-layer soil state as `NTuple{10,Cdouble}` fields, which are stack-allocated and `isbitstype`, making `Ref(p)` safe for `ccall`.\n\n```julia\n@with_kw mutable struct Soil_c <: AbstractSoil{Float64,5}\n  n_layer     ::Cint   = Cint(5)\n  θ           ::NT10   = nzero(10)   # volumetric water content per layer\n  Tsoil_c     ::NT10   = nzero(10)   # current soil temperature\n  K_sat       ::NT10   = nzero(10)   # saturated hydraulic conductivity\n  # ... ~40 fields total\nend\n```\n\nAll array fields use `NT10 = NTuple{10,Cdouble}` regardless of `n_layer`. The extra slots beyond the active layer count are unused padding that matches the C struct layout. **Do not change field types or order** — the C library reads memory offsets directly.\n\n`Soil_c` is constructed during model setup via `_init_soil` → `Init_Soil_Parameters` → `Init_Soil_T_θ!`.\n\n## Calling Convention\n\nEvery wrapper follows the same pattern:\n\n1. Allocate output scalars with `init_dbl()` (returns `Ref{Cdouble}`)\n2. Call `ccall((:symbol, libbeps), Cvoid, (input_types..., Ptr_types...), inputs..., refs...)`\n3. Dereference outputs with `ref[]` and return them\n\nStruct arguments passed by value (e.g., `Leaf`) go directly into `ccall` signature. Struct arguments passed by pointer use `Ref(struct)`. Mutable outputs always use `Ptr{Cdouble}`.\n\n```julia\n# Typical pattern from photosynthesis_c\nGs_w   = init_dbl()\naphoto = init_dbl()\nci     = init_dbl()\nccall((:photosynthesis, libbeps), Cvoid, (..., Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),\n      ..., Gs_w, aphoto, ci)\nGs_w[], aphoto[], ci[]\n```\n\n## Dual Dispatch via Method Overloading (`SOIL_c.jl`)\n\n`SOIL_c.jl` doesn't define new functions — it extends Julia functions imported from `BEPS` with a `Soil_c` method. This means `SoilPhysics` code calls `UpdateHeatFlux(p, ...)` and dispatches to the C implementation when `p::Soil_c`, or to the Julia implementation when `p::SoilModel`.\n\nFunctions overloaded this way include:\n\n- `UpdateRootFraction!(p::Soil_c)` → `SoilRootFraction`\n- `Init_Soil_Parameters(p::Soil_c, ...)` → `Init_Soil_Parameters`\n- `Init_Soil_T_θ!(p::Soil_c, ...)` → `Init_Soil_Status`\n- `soil_water_factor_v2(p::Soil_c)` → `soil_water_factor_v2`\n- `Root_Water_Uptake(p::Soil_c, ...)` → `Soil_Water_Uptake`\n- `UpdateThermal_Cv`, `UpdateThermal_κ`, `Update_ice_ratio`\n- `UpdateHeatFlux`, `Update_Tsoil_c`, `Update_G`, `UpdateSoilMoisture`\n\nNote that C function names sometimes differ from Julia names (e.g., `SoilRootFraction` vs `UpdateRootFraction!`, `Update_Cs` vs `UpdateThermal_Cv`). The Julia names are the stable interface.\n\n## Snowpack Overloads (`snowpack_stage.jl`)\n\nEach snowpack stage has two overloads:\n\n- **Low-level**: accepts raw `TypeRef` (`Ref` or `RefArray`) arguments — matches C signature directly\n- **High-level**: accepts `Layer3{Float64}` / `Layer2{Float64}` mutable structs, wraps fields in `Ref`, calls the low-level version, then writes results back into the struct fields\n\n```julia\n# Low-level\nsnowpack_stage1(Tair, prcp, mass_snow_o_last, ..., mass_snow_o::TypeRef, ...)\n\n# High-level — preferred for integration code\nz_snow = snowpack_stage1(Tair, prcp, lai_o, lai_u, clumping,\n                         mass_snow_pre, mass_snow, perc_snow, area_snow,\n                         z_snow, ρ_snow, albedo_v_snow, albedo_n_snow)\n```\n\nThe high-level overload returns scalar outputs (e.g., `z_snow`) and mutates the `Layer3`/`Layer2` structs in place.\n\n## Physics Wrappers (`module.jl`)\n\nThese wrap individual C physics routines callable for testing or as C-path alternatives in `inter_prg_c`. The main wrappers:\n\n| Julia name | C symbol | Returns |\n|---|---|---|\n| `photosynthesis_c` | `photosynthesis` | `Gs_w, aphoto, ci` |\n| `netRadiation_c` | `netRadiation` | `netRad_o, netRad_u, netRad_g` |\n| `evaporation_soil_c` | `evaporation_soil` | `evapo_soil, evapo_water_g, evapo_snow_g` |\n| `evaporation_canopy_c` | `evaporation_canopy` | `evapo_water_o, evapo_water_u, evapo_snow_o, evapo_snow_u` |\n| `aerodynamic_conductance_c` | `aerodynamic_conductance` | `ra_o, ra_u, ra_g, G_o_a, G_o_b, G_u_a, G_u_b` |\n| `transpiration_c` | `transpiration` | `trans_o, trans_u` |\n| `sensible_heat_c` | `sensible_heat` | `SH_o, SH_u, SH_g` |\n| `surface_temperature_c` | `surface_temperature` | `heat_flux, T_ground, T_any0, T_soil0, T_snow, T_snow1, T_snow2` |\n| `Leaf_Temperature_c` | `Leaf_Temperature` | `Cdouble` (scalar) |\n| `Leaf_Temperatures_c` | `Leaf_Temperatures` | mutates `Tc::Leaf` |\n\n`aerodynamic_conductance_c` uses `Ref(0.0)` for outputs rather than `init_dbl()` — functionally equivalent but inconsistent with the rest of the module.\n\n## Top-Level Entry Point (`inter_prg_c`)\n\n`inter_prg_c` is the C-path equivalent of `inter_prg_jl`. It computes `CosZs` via `s_coszs` (the one Julia call made before entering C), then dispatches the full hourly timestep to `(:inter_prg_c, libbeps)`.\n\n```julia\nfunction inter_prg_c(jday, rstep, lon, lat, lai, Ω,\n                     meteo, parameter, var_o, var_n, soilp::Soil_c,\n                     mid_res, mid_ET, cache; debug=false)\n  CosZs = s_coszs(jday, rstep, lat, lon)\n  ccall((:inter_prg_c, libbeps), Cvoid, (...), ..., Ref(soilp), Ref(mid_res), Ref(mid_ET))\nend\n```\n\n`parameter` and `var_o`/`var_n` are `Vector{Float64}` passed as `Ptr{Cdouble}` — these carry vegetation parameters and per-timestep state that the C library manages internally.\n\n## Data Flow: Setup to Execution\n\n```mermaid\ngraph TD\n  setup_model --> _init_soil\n  _init_soil --> Soil_c\n  _init_soil --> Init_Soil_Parameters\n  _init_soil --> Init_Soil_T_θ!\n  inter_prg_c --> s_coszs\n  inter_prg_c --> libbeps_inter_prg_c[\"(:inter_prg_c, libbeps)\"]\n  SoilPhysics --> Root_Water_Uptake\n  SoilPhysics --> Update_Tsoil_c\n  SoilPhysics --> Update_G\n```\n\n## What Is and Isn't Exported\n\nThe module exports a minimal surface:\n\n- `inter_prg_c` — the C timestep entry point\n- `Soil_c` — the soil state struct\n- `rainfall_stage1`, `rainfall_stage2` — rainfall interception stages\n- `snowpack_stage1`, `snowpack_stage2`, `snowpack_stage3`\n\nIndividual physics wrappers (`photosynthesis_c`, `netRadiation_c`, etc.) are exported from `module.jl` but not re-exported from `BEPS_c.jl`. They are available as `clang.photosynthesis_c` and are used directly in tests.\n\nSeveral functions are commented out (`soilresp`, `rainfall_stage3`, `plantresp`) — their C implementations exist in `libbeps` but are not currently used via the Julia interface.\n\n## Adding a New C Wrapper\n\n1. Identify the C function signature in `beps.h`\n2. Add the wrapper to the appropriate file (`module.jl` for physics, `SOIL_c.jl` for soil ops, `BEPS_c.jl` for orchestration)\n3. Allocate output refs with `init_dbl()`; pass input structs with `Ref(struct)`\n4. If the function operates on `Soil_c`, extend the existing Julia function rather than defining a new name\n5. Add to the `export` list in `BEPS_c.jl` only if it's part of the public C-path interface","chinaflux-case-studies":"# ChinaFlux Case Studies\n\n# ChinaFlux Case Studies (`Project_ChinaFlux/`)\n\nThis module drives multi-site calibration of BEPS against Chinese flux tower observations. It loads hourly meteorological forcing and daily flux/LAI/soil observations for ~12 forest sites, runs the forward model, evaluates goodness of fit, optimizes a small parameter subset via SCEUA, then saves results for R-based analysis.\n\n## Workflow Overview\n\n```mermaid\nflowchart TD\n    A[st_flux12.csv<br/>site metadata] --> B[RunModel per site]\n    F[FORCING<br/>hourly, all sites] --> C[LoadData]\n    D[FluxLAISoil_Daily<br/>per site] --> C\n    C --> B\n    B --> E[simulate<br/>hourly forward model]\n    E --> G[BEPS_GOF<br/>pre-optim metrics]\n    G --> H[optim<br/>SCEUA parameter search]\n    H --> I[goodness<br/>post-optim metrics]\n    I --> J[JLD2 output<br/>gof / gof_opt / data_sim]\n```\n\n## Data Sources\n\nAll input data lives in `ChinaFlux2026/` (path resolved via `path_mnt`):\n\n| File | Contents |\n|---|---|\n| `data/st_flux12.csv` | Site table: `lon`, `lat`, `VegType`, `SoilType`, `z_Uz`, `z_overstory`, `z_SM`, `z_TS` |\n| `data-raw/BEPS/Forcing_Met_Hourly_BEPS_Forest_sp12_hourly_v*.csv` | Stacked hourly forcing for all sites, one `site` column key |\n| `data-raw/Daily/BEPS/{SITE}_FluxLAISoil_Daily_v*.csv` | Per-site daily observations: GPP, ET, Hs, LAI, SM layers, TS layers |\n\nThe site `MF_乔灌混交林_燕山` is excluded (`SITES_bad`) due to data issues.\n\n## `LoadData(SITE)` — Per-site data preparation\n\nReturns `(dates_local, forcing::MetSeries, lai, FluxALL)`.\n\nSteps performed:\n1. Reads per-site daily file; renames standard columns (`LAI_glass_G005 → lai`, `GPP → GPP_obs`, etc.)\n2. Calls `normalize_flux_obs!(FluxALL)` to unit-standardize observations\n3. Slices the global `FORCING` DataFrame to this site; renames canopy-level columns to BEPS names (`Ta_canopy → Tair`, `RH_canopy → RH`, `WS_canopy → Uz`)\n4. Truncates forcing to match daily obs length (warns if mismatch exceeds zero rows)\n5. Calls `sanitize_forcing!(d_forcing)` — gap-fills via `_sanitize_column!` → `_interp_invalid!`\n6. Constructs `MetSeries(; ntime, Rs, Rln_in, Tair, RH, Uz, Prcp)`\n\n**Time convention:** timestamps in the daily file are local CST (UTC+8). `RunModel` converts to UTC before passing to `simulate`:\n```julia\ndates_UTC = dates_local .- Hour(8)\n```\n\n## `RunModel(SITE; ...)` — Forward run + calibration\n\n```julia\nRunModel(SITE; maxn=1000, outdir=\"Project_ChinaFlux/OUTPUT\",\n         goal=:NSE, goal_multiplier=-1)\n```\n\nSkips if the output JLD2 already exists (`isfile(fout) && return`).\n\n**Step 1 — Site configuration**\n\nReads site metadata from `st_full` and initializes the model:\n```julia\nmodel = ParamBEPS(VegType, SoilType)\nmodel.veg.z_wind = z_Uz\nmodel.veg.z_canopy_o = z_overstory\nstate = InitState0(model, forcing)\n```\n\nSoil depth arrays are parsed from comma-separated centimeter strings in the site table and converted to meters:\n```julia\ndepths_SM = map(x -> parse(Int, x), split(st.z_SM, \",\")) ./ 100\n```\n\n**Step 2 — Forward run and pre-optimization GOF**\n\n```julia\ndf_fluxes, df_ET, states, caches = simulate(forcing, lai, dates_UTC; ps=model, state, lon, lat)\ngof, data_sim, data_obs = BEPS_GOF(df_fluxes, states, dates_local, FluxALL; depths_SM, depths_TS)\n```\n\n`BEPS_GOF` (in `src/Optim.jl`) aligns simulated hourly output with daily observations via `_align_daily_data` and `_daily_obs_dates`, then computes NSE/KGE/R² for Flux (GPP, ET, Hs), SM layers, and TS layers.\n\n**Step 3 — Parameter optimization**\n\nFive parameters are calibrated:\n\n| Path | Description | Default | Bounds |\n|---|---|---|---|\n| `[:r_drainage]` | Drainage fraction | 0.5 | (0.2, 0.7) |\n| `[:veg, :Ω]` | Clumping index | 0.85 | (0.3, 1.0) |\n| `[:veg, :g1_w]` | Stomatal slope (MED/BBL) | 8.0 | (1.0, 20.0) |\n| `[:veg, :g0_w]` | Minimum stomatal conductance | 0.0175 | (0.001, 0.1) |\n| `[:veg, :VCmax25]` | Max carboxylation rate at 25°C | 29.0 | (5.0, 200.0) |\n\nThe `paths` list uses `ModelParams.jl` accessor syntax — each entry is a key path into `ParamBEPS`, allowing `optim` to set/restore values without knowing the struct layout:\n```julia\nopts = [\n    (; path=[:r_drainage],     name=:r_drainage, value=model.r_drainage),\n    (; path=[:veg, :Ω],        name=:Ω,          value=model.veg.Ω),\n    ...\n] |> DataFrame\n```\n\n`optim(model, forcing, lai, dates_UTC; paths, maxn, kw_loss...)` runs SCEUA with `maxn` function evaluations. `goal=:NSE` with `goal_multiplier=-1` minimizes negative NSE (i.e., maximizes NSE).\n\n**Step 4 — Output**\n\nResults saved to `{outdir}/BEPS_{SITE}.jld2`:\n```julia\njldsave(fout; gof_opt, gof, theta_opt, data_sim, data_obs)\n```\n\n| Key | Contents |\n|---|---|\n| `gof` | Pre-optimization metrics (Dict with keys `Flux`, `SM`, `TS`) |\n| `gof_opt` | Post-optimization metrics (same structure) |\n| `theta_opt` | Optimized parameter vector |\n| `data_sim` | Simulated daily time series |\n| `data_obs` | Observed daily time series (aligned) |\n\n## Batch execution\n\nSites are processed sequentially with error isolation:\n```julia\nfor SITE in SITES\n    try\n        RunModel(SITE; maxn=1000, outdir=\"Project_ChinaFlux/OUTPUT/NSE\",\n                 goal=:NSE, goal_multiplier=-1)\n    catch ex\n        @error \"Error processing site $SITE: $ex\"\n    end\nend\n```\n\nThe `isfile(fout) && return` guard at the top of `RunModel` makes this idempotent — re-running skips completed sites.\n\n## Visualization (`main_vis.jl`)\n\n`plot_result(sim, obs)` takes the `data_sim` and `data_obs` DataFrames and produces three PNG files:\n\n- `Figure1_Fluxes.png` — GPP, ET, Hs time series (3-panel)\n- `Figure1_SM.png` — all `SM*` columns vs observations\n- `Figure1_TS.png` — all `TS*` columns vs observations\n\n`select_vars` is a local helper that extracts columns by prefix into a matrix for multi-layer SM/TS plots.\n\n## R Post-processing\n\nAfter Julia runs complete, two R scripts analyze the saved JLD2 results using `JuliaCall`:\n\n**`Figure1_GOF.Rmd`** — scatter plots of pre- vs post-optimization KGE/NSE/R² for each variable and site. The `check_var` helper converts Julia `Symbol` arrays to R strings before binding.\n\n**`Table1_GOF.R`** — exports three Excel workbooks:\n- `Table1_GOF_V4.xlsx` — Flux (GPP + ET) GOF by site × index\n- `Table1_SM_GOF.xlsx` — Soil moisture GOF (commented out by default)\n- `Table1_TS_GOF.xlsx` — Soil temperature GOF\n\nThe standard reshape pattern throughout both files:\n```r\ngof_prev  <- map(lst, \\(l) l$gof$Flux)     %>% melt_list(\"site\")\ngof_optim <- map(lst, \\(l) l$gof_opt$Flux) %>% melt_list(\"site\")\ngof <- list(\"Original\"=gof_prev, \"Optimized\"=gof_optim) %>% melt_list(\"type\")\n# → wide format: site × var × index with Original / Optimized columns\n```\n\n## Adding a New Site\n\n1. Add a row to `st_flux12.csv` with the site's metadata fields.\n2. Ensure a matching row exists in the stacked forcing CSV (keyed by `site`).\n3. Create `{SITE}_FluxLAISoil_Daily_v*.csv` in `data-raw/Daily/BEPS/`.\n4. Re-run `case01_ChinaFlux.jl` — the new site will be picked up automatically.\n\nThe `SITES_bad` exclusion list at the top of `case01_ChinaFlux.jl` is the right place to skip problematic sites without modifying source data.","data-types-state":"# Data Types & State\n\n# Data Types & State (`src/DataType/`)\n\nThis module defines every struct, constant, and initialization routine that the rest of BEPS.jl operates on. It enforces a clean JAX-style separation: **parameters** (`ParamBEPS`) are immutable after setup, **state** (`StateBEPS`) is mutable and advances each timestep, and **output** (`Flux`, `ETFlux`) accumulates per-hour results.\n\n---\n\n## Type Hierarchy\n\n```mermaid\ngraph TD\n    ParamBEPS --> HydraulicProfile\n    ParamBEPS --> ThermalProfile\n    ParamBEPS --> ParamVeg\n    StateBEPS --> SnowLand\n    StateBEPS --> TriSolver\n    LeafCache --> PhotoConsts\n    LeafCache --> AeroConsts\n    LeafCache --> Leaf\n```\n\n---\n\n## Parameters: `ParamBEPS`\n\n`ParamBEPS{FT,N,H,T}` is the top-level parameter container. The type parameters encode the floating-point precision (`FT`), the number of soil layers (`N`), and the concrete hydraulic and thermal profile types (`H`, `T`). Because `N` is a type parameter — not a runtime field — dispatches on layer count are zero-cost.\n\n```julia\nParamBEPS{FT,N,H,T} <: AbstractSoilModel{FT,N}\n  dz         ::Vector{FT}      # layer thicknesses [m], default: (0.05, 0.10, 0.20, 0.40, 1.25)\n  r_drainage ::FT              # surface runoff rate, bounded (0.2, 0.7)\n  ψ_min      ::FT              # water stress onset [cm]; 1000 for broadleaf forest, else 3300\n  alpha      ::FT              # stress shape parameter (He 2017 JGR-B Eq. 4)\n  hydraulic  ::H               # Campbell retention + Ksat layers\n  thermal    ::T               # thermal conductivity + heat capacity layers\n  veg        ::ParamVeg{FT}    # LAI, albedo, Vcmax25, stomatal slope, clumping…\n```\n\n### Constructors\n\nThe preferred entry point is the vegetation/soil-type constructor:\n\n```julia\nps = ParamBEPS(\"evergreen_needleleaf\", \"silty_clay_loam\")         # N=5, FT=Float64\nps = ParamBEPS(\"deciduous_broadleaf\",  \"sandy_loam\"; N=5, FT=Float32)\n```\n\nThis calls `InitParam_Veg` and `InitParam_Soil` from `Params/Param_Init.jl` to populate lookup-table values, then auto-selects `ψ_min`/`alpha` based on `veg.is_bforest`.\n\nThe lower-level constructor accepts explicit `HydraulicProfile`/`ThermalProfile` objects for testing or custom parameterizations:\n\n```julia\nParamBEPS{FT,N}(hydraulic, thermal; dz, kwargs...)\n```\n\n### Property forwarding\n\nTwo `Base.getproperty` overloads make `HydraulicProfile` and `ThermalProfile` transparent at the `ParamBEPS` level:\n\n- `ps.K_sat` → `ps.hydraulic.profile.Ksat`\n- `ps.b`, `ps.θ_sat`, `ps.ψ_sat` → forwarded through `hydraulic.profile`\n- `ps.r_root_decay` → forwarded to `ps.veg.r_root_decay`\n- `ps.N` → `nlayer(ps)` (type-level, not a field)\n\n---\n\n## State: `StateBEPS`\n\n`StateBEPS <: AbstractSoil{Float64,5}` holds everything that evolves over time. Fields are grouped by role:\n\n| Group | Fields | Notes |\n|---|---|---|\n| Snow temperatures | `Tsnow_c`, `Tsnow_p` (both `SnowLand`) | current/previous snowpack temperatures |\n| Water stores | `z_water`, `z_snow`, `r_rain_g` | depths [m] and ground rainfall rate [m/s] |\n| Soil moisture | `θ`, `θ_prev`, `ψ`, `ψ_prev`, `ψ_next` | volumetric and matric potential [negative cm] |\n| Soil temperature | `Tsoil_c`, `Tsoil_p` | current and previous [°C] |\n| Water stress | `f_soilwater`, `f_water`, `f_stress`, `w_norm`, `f_temp` | scalar and per-layer factors |\n| Hydraulics | `K`, `K₊ₕ`, `Kmid`, `Kavg`, `r_waterflow` | conductivities [cm h⁻¹] and flow |\n| Heat | `Cv`, `κ`, `G`, `ETi` | volume heat capacity, thermal conductivity, fluxes |\n| Bonan-Q0 solver | `Δz_cm`, `Δz₊ₕ_cm`, `∂θ∂ψ`, `tri` | tridiagonal moisture solver workspace |\n| Root distribution | `f_root` | fraction per layer, recomputed at setup |\n\n### Field classification\n\nTwo compile-time tuples classify fields for time-series export:\n\n```julia\nVARS_SCALAR  # AbstractFloat fields eligible for scalar export\nVARS_VECTOR  # Vector{Float64} fields eligible for matrix export\n             # excludes geometry (dz, Δz_cm, Δz₊ₕ_cm) and solver scratch (ψ_prev, ψ_next, ∂θ∂ψ)\n```\n\n`split_vars(vars)` partitions a user-supplied list against these tuples and returns `(Val{sf}, Val{vf})` used by `save_state!`.\n\n### `SnowLand`\n\nA lightweight `FieldVector{5,FT}` holding the five surface/snowpack temperatures that evolve across hourly steps. The `clamp!(land, Tair)` method bounds all temperatures within ±2 °C of air temperature to prevent runaway numerical drift.\n\n### `state_hydraulic`\n\n`state_hydraulic(st::StateBEPS)` returns a `NamedTuple` that maps `StateBEPS` fields to the names expected by the Bonan-Q0 solver in `SoilDiffEqs.jl`. This decouples field naming in `StateBEPS` from the solver's generic interface — the solver sees `Δz` while `StateBEPS` stores `dz`.\n\n---\n\n## Legacy C-Compatible State: `Soil`\n\n`Soil <: AbstractSoil{Float64,5}` predates `StateBEPS` and exists for C library compatibility. It fuses parameters and state into a single flat struct mirroring the C `SOIL` layout. All vectors are pre-allocated at length 10 (`MAX_LAYERS`).\n\nTwo functions bridge the two representations:\n\n- `Params2Soil!(soil, params)` — writes `ParamBEPS` fields into a `Soil` (used before calling `libbeps`)\n- `Soil2Params!(params, soil)` — reads a `Soil` back into `ParamBEPS` and updates layer caches\n\n`StateBEPS(soil)` constructs a state from a `Soil` by copying the state-tagged fields. `State2Soil!(soil, st)` writes them back using `@pack!`.\n\n---\n\n## Initialization: `setup`\n\n`setup` is the single recommended entry point for constructing a ready-to-run model:\n\n```julia\nstate, ps = setup(\"evergreen_needleleaf\", \"silty_clay_loam\")\nstate, ps = setup(\"evergreen_needleleaf\", \"silty_clay_loam\";\n  Ta=15.0, Tsoil=12.0, θ0=0.35, z_snow=0.0,\n  r_drainage=0.5, r_root_decay=0.95)\nstate, ps = setup(ps)   # reinitialize state from existing ParamBEPS\n```\n\nInternally, `setup` calls `_init_state(ps, Tsoil, Ta, θ0, z_snow)`, which:\n\n1. Allocates a blank `StateBEPS`\n2. Copies `dz` from `ps` and computes `Δz₊ₕ_cm` (inter-layer midpoint distances for the finite-difference solver)\n3. Calls `UpdateRootFraction!(st, ps)` to distribute root fraction across layers\n4. Calls `Init_Soil_T_θ!(st, ...)` to set initial soil temperature and moisture profiles\n\n`setup_model` is the legacy entry point, kept for backward compatibility with C-path tests. It returns `(soil, state, ps)` and dispatches on `version=\"julia\"` vs `version=\"c\"`.\n\n---\n\n## Forcing: `Met` and `MetSeries`\n\n`Met` holds a single hourly timestep of atmospheric forcing:\n\n| Field | Units | Notes |\n|---|---|---|\n| `Rs` | W m⁻² | inward shortwave |\n| `Rln_in` | W m⁻² | inward longwave; `NaN` triggers model estimation |\n| `Tair` | °C | 2 m air temperature |\n| `RH` | % | relative humidity |\n| `Prcp` | mm h⁻¹ | precipitation |\n| `Uz` | m s⁻¹ | wind speed at measurement height |\n\n`MetSeries` is the time-series counterpart, generated by `@DefFluxSeries MetSeries = Met`. `fill_met!(met, forcing, i)` copies timestep `i` from `MetSeries` into a single `Met` struct, which is then passed to `inter_prg_jl`. This avoids per-timestep allocation.\n\n---\n\n## Output: `Flux` and `ETFlux`\n\n`Flux` accumulates the primary hourly outputs written to the result DataFrame:\n\n- Carbon: `gpp_o_sunlit`, `gpp_u_sunlit`, `gpp_o_shaded`, `gpp_u_shaded` [μmol m⁻² s⁻¹]; `GPP`, `NPP`, `NEP` [gC m⁻²]\n- Energy: `SH`, `LH`, `Net_Rad` [W m⁻²]\n- Water: `Trans`, `Evap` [mm]; `z_water`, `z_snow` [m]; `ρ_snow` [kg m⁻³]\n\n`ETFlux` carries the detailed ET partitioning at instantaneous rates [kg m⁻² s⁻¹]. `update_ET!(etflux, flux, Ta)` integrates these rates over the timestep (`step = 3600 s`) and fills the summary fields in `flux`.\n\nBoth types extend `AbstractFlux`. `@DefFluxSeries FluxSeries = Flux` auto-generates the corresponding time-series containers via the macro system (see below).\n\n---\n\n## Per-Leaf Cache: `LeafCache`\n\n`LeafCache` holds intermediate variables for a single leaf computation. All `Leaf`-typed fields represent overstory/understory pairs. Key groupings:\n\n- **Conductances:** `Gs_old/new`, `Gc`, `Gh`, `Gw`, `Gw_wet` — stomatal and boundary-layer conductances at various reference points\n- **Temperatures:** `Tc_old/new`, `Cs_old/new`, `Ci_old/new` — leaf surface/intercellular CO₂ and temperature iteration\n- **Radiation:** `Rn`, `Rns`, `Rnl` — net, shortwave, longwave radiation components\n- **Precomputed constants:** `pc::PhotoConsts` and `ac::AeroConsts` — temperature-dependent photosynthesis constants and aerodynamic geometry, recomputed when inputs change rather than each inner iteration\n\n`ALL_VARS_CACHE` lists the exportable `Leaf` fields. `CacheSeries` and `save_cache!` follow the same `Val{VAR}`-dispatch pattern as `StateSeries`.\n\n### `PhotoConsts` and `AeroConsts`\n\n`PhotoConsts{T}` caches the five Farquhar kinetic terms (Γ*, K, Rd/Jmax/Vcmax temperature factors) for a given leaf temperature. Call `PhotoConsts!(pc, T_leaf_K)` to update in place.\n\n`AeroConsts{T}` caches aerodynamic geometry terms (friction velocity, Monin-Obukhov coefficient, boundary-layer resistances, attenuation factors) that are expensive to compute but change only with canopy structure and wind, not with the canopy energy iteration. `AeroConsts!(ac, ...)` recomputes these; `ra_updateH(...)` then applies the stability correction in the inner loop using only the current sensible heat flux.\n\n---\n\n## Time-Series Output: `StateSeries` and `CacheSeries`\n\nBoth containers use `Val{VARS}`-dispatched `@generated` functions so the field-copy loop is unrolled at compile time — no runtime iteration over field names.\n\n```julia\n# Allocate\nsf, vf = split_vars([:z_snow, :θ, :Tsoil_c])\nstates = StateSeries(sf, vf, N, n_time)\n\n# Fill each timestep (called inside simulate's inner loop)\nsave_state!(states, state, t, sf, vf)\n\n# Index a single timestep back\nstates[42]   # → Dict of scalar and vector values at t=42\n```\n\n`StateSeries` is a type alias for `NamedTuple{(:scalars,:vectors), ...}` — scalars become `Vector{Float64}`, vectors become `Matrix{Float64}(n_time × n_layer)`.\n\n---\n\n## Macro System\n\n`@DefFluxSeries` and `@DefStateSeries` generate time-series structs from single-timestep structs. Given `Flux`, `@DefFluxSeries FluxSeries = Flux` produces a `FluxSeries{FT} <: AbstractFluxSeries{FT}` where every non-vector field becomes `Vector{FT}`. The `setindex!` method on `AbstractSeries` is `@generated` to inline per-field assignment with zero runtime overhead.\n\n`@DefState` and `@DefFlux` are lightweight builders for simple structs with all-`FT` fields, used in specialized submodules.\n\n---\n\n## Connections to the Rest of the Codebase\n\n- `inter_prg_jl` (`src/inter_prg.jl`) receives `(met, ps, state, flux, etflux, cache)` and mutates `state`, `flux`, and `etflux` in place each hour.\n- `UpdateSoilMoisture!` (`src/SoilPhysics/UpdateSoilMoisture.jl`) reads the Bonan-Q0 solver fields from `StateBEPS` via `state_hydraulic(st)`.\n- `simulate` (`src/beps_modern.jl`) calls `split_vars` / `split_cache_vars` to build `Val`-typed export specs, then calls `save_state!` and `save_cache!` in the inner loop.\n- The C path (`src/clang/`) uses `Soil` and `Soil_c` directly; `Params2Soil!` / `State2Soil!` bridge the two worlds before and after each C call.\n- `interp_depths(SM, z_obs)` re-interpolates simulated soil moisture from the five model layer centers to arbitrary observation depths, useful when comparing against in-situ sensors.","hydrology-snowpack":"# Hydrology & Snowpack\n\n# Hydrology & Snowpack Module\n\n**Files:** `src/rainfall_stage.jl`, `src/snowpack.jl`\n\nThis module tracks the movement of water through the canopy-atmosphere interface at each hourly timestep. It handles two parallel processes — liquid rainfall interception and solid snowfall interception — using a three-layer cascade (overstory → understory → ground). Each process is split into numbered stages that bracket the evaporation/sublimation calculations in `inter_prg_jl`.\n\n---\n\n## Canopy Interception Architecture\n\nBoth rainfall and snow use the same physical model: Beer-Lambert canopy interception, a finite storage capacity proportional to LAI, and a layer cascade where throughfall from one layer becomes input to the next.\n\n```mermaid\nflowchart TD\n    P[\"Precipitation (prcp, m s-1)\"] --> OS[\"Overstory<br/>m_water/snow.o\"]\n    OS -- \"throughfall\" --> US[\"Understory<br/>m_water/snow.u\"]\n    US -- \"throughfall\" --> GR[\"Ground<br/>z_water / m_snow.g\"]\n    OS -- \"evap/sublim\" --> ATM[\"Atmosphere\"]\n    US -- \"evap/sublim\" --> ATM\n```\n\n**Interception fraction** for each layer follows Beer-Lambert:\n```\nτ = 1 − exp(−lai · Ω)\n```\nwhere `Ω` is the clumping index. `τ` is the fraction of incoming precipitation captured by that layer's canopy.\n\n**Storage capacity:**\n- Water: `mMax_water = 0.1 · lai` [kg m⁻²]\n- Snow: `massMax_snow = 0.1 · lai` [kg m⁻²]\n\n---\n\n## Rainfall Interception\n\n### `water_change(m_water_pre, prcp, lai, Ω; kstep)`\n\nCore computation for a single canopy layer. Converts precipitation flux to mass increment, adds it to the pre-existing stored water, clamps to capacity, and returns the updated storage, fractional fill, and mass increment (which drives throughfall).\n\n```julia\n# prcp: [m s-1]  →  mass added: prcp * kstep * ρ_w * τ  [kg m-2]\nΔm_water_o = max(m_water_o - m_water_pre, 0.0)\n```\n\nOnly the increase above the previous timestep's value is considered new interception (i.e., negative changes don't propagate as throughfall).\n\n### `rainfall_stage1_jl(Tair, prcp, frac_water, m_water, m_water_pre, lai_o, lai_u, Ω)`\n\n**Called before evaporation** in the timestep loop.\n\nPartitions precipitation through the canopy cascade:\n1. If `Tair ≤ 0`, all precipitation is treated as snow — `prcp` is zeroed here (snowpack handles it).\n2. Overstory intercepts using `water_change`; the mass increment `Δm_water_o` is withheld from the layer below.\n3. Understory intercepts from the remaining throughfall.\n4. Returns `prcp_g` — the ground-level precipitation rate [m s⁻¹] after both layers have absorbed their share.\n\nMutates `m_water` and `frac_water` (both `Layer2{Float64}`) in place.\n\n### `rainfall_stage2_jl(evapo_water_o, evapo_water_u, mass_water)`\n\n**Called after evaporation** is computed.\n\nDrains canopy water stores by the evaporated amount, floored at zero:\n```julia\nmass_water.o = max(mass_water.o - evapo_water_o * kstep, 0.0)\n```\n`evapo_water_o` is a flux [kg m⁻² s⁻¹]; multiplying by `kstep` converts to [kg m⁻²] mass removed.\n\n---\n\n## Snowpack\n\nThe snowpack module mirrors the rainfall structure but adds a third complication: snow on the ground undergoes densification, melting, and freezing, requiring three separate stages.\n\n### `snow_change(m_snow_pre, snowrate, kstep, ρ_new_snow, lai, Ω)`\n\nCanopy snow interception for one layer. Analogous to `water_change` but operates on snow mass and also returns `area_snow` (fractional canopy area covered by snow), which affects radiation calculations:\n```julia\nareaMax_snow = 0.01 * lai   # max area fraction\narea_snow = perc_snow * areaMax_snow\n```\n\n### Stage 1 — Accumulation: `snowpack_stage1_jl(...)`\n\n**Called before sublimation.**\n\nProcesses snowfall when `Tair < 0`. New snow density is computed empirically from air temperature:\n```julia\nρ_new = 67.9 + 51.3 * exp(Tair / 2.6)   # [kg m-3]\n```\nPrecipitation (rain rate) is converted to snow depth rate:\n```julia\nsnowrate = prcp * ρ_w / ρ_new   # [m s-1]\n```\n\nThe cascade then mirrors rainfall interception: overstory → understory → ground. The ground snow layer uses mass-weighted density mixing to track the bulk snowpack density `ρ_snow`:\n```julia\nρ_snow[] = (ρ_snow[] * z_snow + ρ_new * δ_zs) / (z_snow + δ_zs)\n```\nWhen no new snow falls, density evolves by exponential compaction toward 250 kg m⁻³.\n\n**Snow cover fraction** for the ground layer uses a 50 cm threshold:\n```julia\nfrac_snow.g = min(m_snow.g / (0.05 * ρ_snow[]), 1.0)\n```\n\n**Snow albedo** decays exponentially toward a minimum when snow is aging, and resets to maximum when fresh snow falls:\n```julia\n# Aging:\nalbedo_v_snow[] = (albedo_v_snow[] - 0.70) * exp(-0.005 * kstep / 3600) + 0.70\n# Fresh:\nalbedo_v_snow[] = 0.94\n```\n\nReturns `z_snow` (ground snow depth, capped at 10 m).\n\n### Stage 2 — Sublimation: `snowpack_stage2_jl(evapo_snow_o, evapo_snow_u, m_snow)`\n\n**Called after sublimation fluxes are computed.**\n\nReduces canopy snow mass by sublimation, floored at zero. Does not touch ground snow (sublimation from the ground surface is handled separately).\n\n### Stage 3 — Melt/Freeze: `snowpack_stage3_jl(Tair, Tsnow, Tsnow_last, ρ_snow, z_snow, z_water, m_snow)`\n\n**Called at the end of the snowpack update sequence.**\n\nConverts between ground snow (`m_snow.g`, `z_snow`) and surface liquid water (`z_water`) depending on snowpack temperature.\n\nThe energy available for phase change is computed by `cal_melt`:\n```julia\nE = z_snow * ρ_snow * (Tsnow - 0) * cp_ice   # [J m-2]\nΔm = E / λ_fusion                              # [kg m-2]\n```\nPositive `Δm` → melt; negative → freeze.\n\nThree depth regimes apply different melt schemes:\n\n| Depth (`z_snow`) | Melt scheme |\n|---|---|\n| ≤ 0.02 m | Temperature-index: `ms_melt = Tair · 0.0075 · kstep/3600 · 0.3` |\n| 0.02–0.05 m | Full energy balance: `ms_melt = min(Δm, ms_sup)` |\n| > 0.05 m | 2% of energy balance: `ms_melt = min(0.02·Δm, ms_sup)` (shallow active layer approximation) |\n\nThe 2% factor for deep snowpacks reflects that only the surface layer (~2 cm of a thick pack) exchanges heat on a sub-daily timescale. The code comments note this is a known approximation and the correct form should scale with `max(zs_sup * 0.02, 0.02)`.\n\nAfter melt/freeze, `z_snow` and `z_water` are updated consistently:\n```julia\nz_snow  = zs_sup  + (mw_frozen - ms_melt) / ρ_snow\nz_water = z_water + (ms_melt  - mw_frozen) / ρ_w\n```\n\n---\n\n## Integration with `inter_prg_jl`\n\nThe stages are called in this order within each hourly timestep:\n\n```\nrainfall_stage1_jl   ← interception before evap (determines prcp_g)\nsnowpack_stage1_jl   ← accumulation before sublimation\n  ↓ [evaporation/sublimation computed]\nrainfall_stage2_jl   ← drain canopy water by evaporation\nsnowpack_stage2_jl   ← drain canopy snow by sublimation\n  ↓ [soil heat / Tsnow update]\nsnowpack_stage3_jl   ← melt/freeze ground snow → z_water\n  ↓\nUpdateSoilMoisture!  ← z_water enters soil infiltration\n```\n\n`prcp_g` from `rainfall_stage1_jl` feeds directly into the soil moisture update as the ground-level water input. `z_water` produced by `snowpack_stage3_jl` is also routed there as meltwater.\n\n---\n\n## Key Types\n\n| Type | Description |\n|---|---|\n| `Layer2{Float64}` | Two-field struct: `.o` (overstory), `.u` (understory) |\n| `Layer3{Float64}` | Three-field struct: `.o`, `.u`, `.g` (ground) |\n| `Ref{Float64}` | Used for scalar state updated in-place (`ρ_snow`, `albedo_v_snow`, etc.) |\n\nUnits throughout: mass in [kg m⁻²], depth in [m], fluxes in [m s⁻¹] or [kg m⁻² s⁻¹], time in [s] with `kstep = 360` s default.\n\n---\n\n## Known Issues / TODOs\n\n- **Deep snowpack melt (`snowpack_stage3_jl`):** The 2% active-layer approximation is hardcoded. The correct expression `max(zs_sup * 0.02, 0.02)` is noted in a comment but not yet applied.\n- **New snow density clamp:** The `ρ_new` clamp to `[50, 200]` kg m⁻³ is present in the formula comment but commented out in the code — at very warm subfreezing temperatures the empirical formula can produce anomalously high densities.","model-core-api":"# Model Core & API\n\n# Model Core & API\n\nThe model core consists of three layers: the module entry point (`BEPS.jl`), the simulation API (`simulate` / `beps_main`), and the hourly time-step engine (`inter_prg_jl`). Together they form the complete pipeline from forcing data to flux output.\n\n---\n\n## Module Entry Point (`BEPS.jl`)\n\n`BEPS.jl` handles three responsibilities at load time:\n\n1. **Dependency wiring** — re-exports `Serialization`, `DelimitedFiles`, and `ModelParams`; imports `DataFrame` without re-exporting it.\n2. **Library resolution** — resolves `libbeps` at module initialization using a platform-aware filename (`libbeps-{os}-{arch}.{ext}`). It checks `deps/` first; if absent, falls back to `LazyArtifacts`.\n3. **Include order** — `SPAC` and `DataType` are included before `BEPS_modules.jl` because the physics modules depend on those types.\n\n`path_proj(f...)` is a convenience helper that anchors paths to the package root (one level above `src/`). Tests use it extensively to load forcing data.\n\n---\n\n## Simulation APIs\n\n### `simulate` — Modern API\n\n```julia\ndf_flux, df_ET, states, caches = simulate(forcing, lai, dates;\n    ps::ParamBEPS, state::StateBEPS,\n    lon=120.0, lat=20.0,\n    SM_obs=nothing, TS_obs=nothing,\n    VARS_STATE=DEFAULT_STATE_EXPORT,\n    VARS_CACHE=DEFAULT_CACHE_EXPORT)\n```\n\n`simulate` is the primary entry point. Key design decisions:\n\n- **State isolation**: `state = deepcopy(state)` at entry — the caller's `StateBEPS` is never mutated.\n- **SM/Tsoil prescription**: When `SM_obs` or `TS_obs` are provided (e.g., for calibration runs where observed soil state should drive the simulation), the corresponding `state` fields are overwritten each timestep *before* calling `inter_prg_jl`. The flags `fix_sm` and `fix_Tsoil` are forwarded to the sub-hourly loop to suppress the corresponding update equations.\n- **Output series**: Four collectors are pre-allocated — `FluxSeries`, `ETSeries`, `StateSeries`, and `CacheSeries`. `StateSeries` and `CacheSeries` save only the variables listed in `VARS_STATE` / `VARS_CACHE`, keeping memory use proportional to what's actually needed.\n- **LAI indexing**: `lai` is daily; the hourly index maps to a day via `ceil(Int, i / 24)`.\n- **Hours are 0-based** in `simulate` (`hour.(dates)`) versus 1-based in `beps_main` (`hour.(dates) .+ 1`). This matters when passing `hour` into `inter_prg_jl`, which uses it for solar angle calculation.\n\n`beps_modern` is an alias for `simulate` kept for backward compatibility.\n\n### `beps_main` — Legacy API\n\n```julia\ndf_flux, df_ET, states = beps_main(forcing, lai, dates;\n    version=\"julia\",   # or \"c\"\n    VegType=\"default\", SoilType=\"silty_clay_loam\", ...)\n```\n\n`beps_main` exists specifically to test consistency between the Julia and C execution paths. It calls `setup_model` (which returns the legacy `soil` struct alongside `state` and `params`) and dispatches to either `inter_prg_jl` or `inter_prg_c` based on `version`. It does not return a `CacheSeries`. New code should use `simulate` instead.\n\n---\n\n## Hourly Time-Step Engine\n\n### `inter_prg_jl`\n\n`inter_prg_jl` runs a single model hour. Internally it runs a sub-hourly loop (`kloop = step / kstep`, default 10 iterations at 360 s each) over the full energy balance. The structure within that loop is:\n\n```mermaid\ngraph TD\n    A[snowpack_stage1] --> B[rainfall_stage1]\n    B --> C[soil_water_factor_v2]\n    C --> D[solve_canopy_energy_balance!]\n    D --> E[transpiration_jl]\n    E --> F[evaporation_canopy_jl]\n    F --> G[rainfall_stage2 / snowpack_stage2]\n    G --> H[evaporation_soil_jl]\n    H --> I[surface_temperature!]\n    I --> J[snowpack_stage3]\n    J --> K[sensible_heat_jl]\n    K --> L[UpdateHeatFlux]\n    L --> M[Root_Water_Uptake]\n    M --> N[UpdateSoilMoisture]\n```\n\nEach stage is a side-effect-free function operating on stack-allocated locals; only `state` fields are mutated. After the sub-hourly loop, accumulated results are packed into `mid_ET` and `mid_flux`.\n\n**GPP unit conversion** happens at the very end:\n```julia\nmid_flux.GPP = GPP * 12 * step * 1e-6  # [µmol m⁻² s⁻¹] → [gC m⁻²]\n```\n`step` is the full hourly interval in seconds (3600), not `kstep`.\n\n**Soil albedo** is computed each sub-step from `state.θ_prev[2]` (layer 2 moisture), linearly interpolating between `α_soil_dry` and `α_soil_sat`.\n\n### `solve_canopy_energy_balance!`\n\nThis function was extracted from `inter_prg_jl` to isolate the iterative inner loop. It converges canopy temperatures (`Tc_old` / `Tc_new`) to within `atol=0.02°C` across four leaf fractions (overstory sunlit/shaded, understory sunlit/shaded):\n\n1. Update aerodynamic resistances via `ra_updateH` (uses previous-iteration sensible heat `H_canopy_o` to estimate stability correction).\n2. Compute conductances `Gh`, `Gw_wet`, `Gw`, `Gc` from `Ga`, `Gb`, `Gs`.\n3. Recompute longwave radiation (`netRadiation_LW!`) — this is the only radiation term that depends on `Tc` and must re-run each iteration. Shortwave (`netRadiation_SW!`) is computed once before the loop.\n4. Run `photosynthesis` (daytime only) to update `Gs_new`.\n5. Solve leaf temperatures via `Leaf_Temperatures_jl`.\n6. Check convergence; if not converged after `maxn=10` iterations, fall back to `Tc = Tair`.\n\nAfter convergence, `GPP` is computed as `multiply!(GPP, Ac, LAI)` (element-wise product of net assimilation rate and leaf area index by fraction).\n\n---\n\n## Data Flow Summary\n\n```mermaid\ngraph LR\n    MS[MetSeries] -->|fill_met!| M[Met]\n    LAI[lai vector] --> IP[inter_prg_jl]\n    M --> IP\n    PS[ParamBEPS] --> IP\n    ST[StateBEPS] <-->|mutated| IP\n    IP --> FL[Flux]\n    IP --> ET[ETFlux]\n    FL --> FS[FluxSeries → DataFrame]\n    ET --> ES[ETSeries → DataFrame]\n    ST -->|save_state!| SS[StateSeries]\n    LC[LeafCache] -->|save_cache!| CS[CacheSeries]\n```\n\n`ParamBEPS` is read-only inside `inter_prg_jl`. `StateBEPS` is the only struct that carries state across timesteps; `LeafCache` is overwritten each sub-step and its end-of-hour snapshot is saved to `CacheSeries`.\n\n---\n\n## Prescribing Observed State\n\nThe `SM_obs` and `TS_obs` arguments to `simulate` support data-assimilation workflows where measured soil profiles should override model-computed values:\n\n```julia\n# SM_obs: nlayer × ntime matrix, rows = soil layers 1–5\nsimulate(forcing, lai, dates; ps, state,\n    SM_obs = observed_theta,   # skips UpdateSoilMoisture\n    TS_obs = observed_tsoil)   # skips Tsoil update in UpdateHeatFlux\n```\n\nThe prescription happens *before* `inter_prg_jl` is called for each timestep. `state.θ_prev` is set to the current `state.θ` first (preserving the previous-step value the physics needs), then `state.θ[1:5]` is overwritten from the observation matrix column. The `fix_sm` / `fix_Tsoil` flags are then passed into `inter_prg_jl` and forwarded to `UpdateSoilMoisture` / `UpdateHeatFlux` to suppress the update equations while still allowing dependent calculations (e.g., `ice_ratio` is still updated from the prescribed temperatures).\n\n---\n\n## Key Exports\n\n| Symbol | Source | Notes |\n|--------|--------|-------|\n| `simulate` | `beps_modern.jl` | Primary API |\n| `beps_modern` | `beps_modern.jl` | Alias for `simulate` |\n| `beps_main` | `beps_main.jl` | Legacy/testing only |\n| `path_proj` | `BEPS.jl` | Path helper relative to package root |\n| `inter_prg_jl` | `inter_prg.jl` | Not exported; called by both APIs |\n| `Soil_c` | `clang/BEPS_c.jl` | Re-exported for C-path users |","model-parameters":"# Model Parameters\n\n# Model Parameters Module\n\n`src/DataType/Params/`\n\nThis module defines all physical parameters for the BEPS model — vegetation characteristics, soil hydraulics, soil thermal properties, and photosynthesis kinetics. It is the single source of truth for model configuration, and is consumed by `setup()` to initialize both `ParamBEPS` (parameter container) and `StateBEPS` (initial state).\n\n---\n\n## Parameter Types\n\n### `ParamVeg{FT}`\n\nVegetation parameters with bounds metadata for optimization. Defined with `@bounds @with_kw`, every field carries a valid range `(lo, hi)` used by `filter_params` / `update!` from `ModelParams.jl`.\n\nKey fields:\n\n| Field | Description | Units |\n|---|---|---|\n| `LAI_max_o` / `LAI_max_u` | Peak LAI for overstory / understory | m² m⁻² |\n| `α_canopy_vis` / `α_canopy_nir` | Canopy albedo (visible / NIR) | — |\n| `α_soil_sat` / `α_soil_dry` | Soil albedo at wet and dry endpoints | — |\n| `Ω` | Clumping index | — |\n| `z_canopy_o` / `z_canopy_u` | Canopy height (overstory / understory) | m |\n| `g1_w` / `g0_w` | Ball-Berry stomatal slope and intercept | — / mol m⁻² s⁻¹ |\n| `VCmax25` | Maximum Rubisco carboxylation rate at 25°C | μmol m⁻² s⁻¹ |\n| `N_leaf` | Leaf nitrogen content | g m⁻² |\n| `slope_Vc` | Vcmax–N relationship slope | — |\n\nTwo boolean flags control canopy structure: `has_understory` (false for CRO and C4) and `is_bforest` (true for DBF/EBF, affects stomatal ψ thresholds).\n\n### `ParamPhoto_Farquhar{FT}`\n\nStandalone photosynthesis parameters for the Farquhar model, independent of `ParamVeg`. Used directly by the photosynthesis standalone path (`test-photosynthesis_standalone.jl`). Contains Vcmax25, Jmax25, temperature response energies (`evc`, `ejm`, `erd`), optimum temperatures, Michaelis-Menten constants (`kc25`, `ko25`, `tau25`), electron transport shape (`qalpha`, `theta2`), and Ball-Berry coefficients.\n\n`InitParam_Photo_Farquhar(VegType::Int; FT)` constructs this from hardcoded defaults; VegType is accepted for future per-PFT loading.\n\n### Soil Parameters\n\nSoil hydraulics and thermal properties are not single structs but composites delegated to `ModelParams.jl` types:\n\n- `HydraulicProfile{FT,N}` — wraps `CampbellLayers{FT,N}` (retention curve parameters) and `KvLayers{FT,N}` (vertical conductivity profile)\n- `ThermalProfile{FT,N}` — wraps `ThermalBaseLayers{FT,N}` (dry conductivity, bulk density, organic matter fraction)\n\n`N` is the number of soil layers (default 5, depths `[0.05, 0.10, 0.20, 0.40, 1.25]` m).\n\n---\n\n## Initialization Functions\n\n### `InitParam_Veg(VegType; FT) → ParamVeg{FT}`\n\nLoads from `data/ParamVeg.json` and `data/ParamGeneral.json`. The JSON key lookup goes through `veg_name(code) → veg_param_key(name)` — notably, CRO maps to `\"default\"` key in the JSON.\n\n```julia\nps_veg = InitParam_Veg(\"ENF\")         # by name\nps_veg = InitParam_Veg(1)             # by BEPS land-cover code\n```\n\nSupported VegTypes and their codes:\n\n| Name | Code | Notes |\n|---|---|---|\n| `ENF` | 1 | Evergreen needleleaf forest |\n| `DNF` | 2 | Deciduous needleleaf |\n| `DBF` | 6 | Deciduous broadleaf — `is_bforest=true` |\n| `EBF` | 9 | Evergreen broadleaf — `is_bforest=true` |\n| `Shrub-SH` | 13 | |\n| `CRO` | 25 | Cropland — no understory, loads `\"default\"` JSON key |\n| `C4` | 40 | C4 grass — no understory |\n| `default` | -1 | |\n\n### `InitParam_Soil(SoilType, N, FT) → (HydraulicProfile, ThermalProfile)`\n\nLooks up `SOIL_PARAMS` by name or integer index, interpolates five-layer tabulated values to `N` layers via `_fit_layers`, and constructs the `ModelParams.jl` composite types.\n\n**Unit conversion note:** `SOIL_PARAMS` stores `ψ_sat` as positive suction in meters (Campbell 1974 convention). `InitParam_Soil` converts to negative centimeters before constructing `HydraulicProfile`, which is what `soil_moisture_Q0!` expects.\n\n```julia\nhydraulic, thermal = InitParam_Soil(\"silty_clay_loam\", 5, Float64)\n```\n\nThe 11 supported soil textural classes, in order:\n\n`sand → loamy_sand → sandy_loam → loam → silty_loam → sandy_clay_loam → clay_loam → silty_clay_loam → sandy_clay → silty_clay → clay`\n\nOut-of-range integer codes fall back to `clay` (index 11).\n\n---\n\n## Lookup Tables (`GlobalData.jl`)\n\n`SOIL_PARAMS` is a compile-time constant NamedTuple array. Each entry holds per-layer arrays for `b` (Campbell exponent), `K_sat` (cm h⁻¹), and scalar `θ_sat`, `θ_res`, `κ_dry`.\n\n`find_SoilType` and `find_VegType` resolve string or integer identifiers to indices — used throughout `Param_Init.jl` and `ultilize.jl`.\n\n---\n\n## Initialization Flow\n\n```mermaid\ngraph TD\n    setup[\"setup(VegType, SoilType)\"] --> IV[\"InitParam_Veg()\"]\n    setup --> IS[\"InitParam_Soil()\"]\n    IV --> JSON[\"ParamVeg.json / ParamGeneral.json\"]\n    IS --> SOIL[\"SOIL_PARAMS constant\"]\n    IS --> ML[\"ModelParams.jl<br/>CampbellLayers / HydraulicProfile / ThermalProfile\"]\n    IV --> PV[\"ParamVeg{FT}\"]\n    IS --> HP[\"HydraulicProfile{FT,N}\"]\n    IS --> TP[\"ThermalProfile{FT,N}\"]\n    PV & HP & TP --> PB[\"ParamBEPS{FT,N,H,T}\"]\n```\n\n---\n\n## Optimization Interface\n\n`ParamVeg` and `ParamPhoto_Farquhar` expose their fields through `ModelParams.jl`:\n\n```julia\n# Get optimizable parameters as a flat vector\nx, names = parameters(ps.veg)\n\n# Update struct from flat vector\nupdate!(ps.veg, x)\n\n# Get bounds for optimizer\nlo, hi = split_bounds(bounds(ps.veg))\n\n# Filter to subset\nactive = filter_params(ps.veg, [:VCmax25, :g1_w, :Ω])\n```\n\nThe `@bounds` macro attaches `(lo, hi)` tuples to field metadata; `bounds()` extracts them in the same order as `parameters()`.\n\n---\n\n## Deprecated API\n\n`deprecated/` contains the legacy flat-vector representation used by the original C-binding path.\n\n- **`ReadParamVeg(lc) → Vector{Float64}`** — returns a 48-element parameter vector. Index mapping is fixed (e.g., `param[37]` = VCmax25). Still used internally by `inter_prg_c`.\n- **`theta2par!(par, theta)` / `par2theta!(theta, par)`** — convert between the flat vector and `ParamVeg`. These are the bridge between old and new APIs.\n- **`Init_Soil_Parameters(soil, VegType, SoilType, r_root_decay)`** — populates the mutable `Soil` struct used by the C path. Sets `ψ_min` and `alpha` based on broadleaf flag.\n- **`Init_Soil_var`** — wraps `Init_Soil_Parameters` + `Init_Soil_T_θ!` + `InitState!` for full legacy initialization.\n\nDo not use these in new Julia-path code. They exist to keep `inter_prg_c` working without changes.","optimization-calibration":"# Optimization & Calibration\n\n# Optimization & Calibration (`src/Optim.jl`)\n\nThis module provides parameter calibration for BEPS using the SCE-UA global optimizer. It wraps the forward model (`simulate`) with a parameter-update/restore cycle, aggregates simulation output to daily flux and state diagnostics, computes multi-variable goodness-of-fit scores, and exposes a scalar loss function that SCE-UA minimizes.\n\n## Call chain\n\n```mermaid\ngraph TD\n    optim --> _loss\n    _loss --> loss\n    loss --> goodness\n    goodness --> predict\n    goodness --> BEPS_GOF\n    predict --> InitState0\n    predict --> simulate\n    BEPS_GOF --> normalize_flux_obs!\n    BEPS_GOF --> _align_daily_data\n```\n\n`optim` is the entry point. Every layer below it is also callable independently, which makes forward evaluation and diagnostic plotting easy without running a full calibration.\n\n---\n\n## `InitState0`\n\n```julia\nInitState0(model::ParamBEPS{FT}, forcing::MetSeries{FT}) -> StateBEPS\n```\n\nBuilds a warm-start `StateBEPS` from the first air temperature in `forcing`. Soil temperature is set to `Tair[1]`; initial volumetric soil moisture is 80 % of saturation (`θ_sat[1] * 0.8`), chosen to stay above the wilting point without prescribing a measured value. Snow depth starts at zero.\n\nDelegates to `BEPS._init_state` (`src/DataType/setup.jl`).\n\n---\n\n## `predict`\n\n```julia\npredict(theta, model, forcing, lai, dates_UTC; paths, lon, lat)\n    -> (df_fluxes, df_ET, states, caches)\n```\n\nRuns one complete forward simulation under the parameter vector `theta`:\n\n1. `deepcopy(model)` — isolates the trial from the caller's model.\n2. Records `theta_prev` via `parameters(model; paths)` so the original values can be restored.\n3. `BEPS.update!(model, paths, theta)` — writes the trial parameters into the copy.\n4. `InitState0` — spins up initial state.\n5. `simulate(...)` — runs the full hourly time loop.\n6. `BEPS.update!(model, paths, theta_prev)` — restores the original; the deepcopy means this is a no-op for the caller's model but keeps internal bookkeeping clean.\n\n`paths` is a `ModelParams` path selector that maps parameter vector indices to struct fields (e.g. `ps.veg.Vcmax25`, `ps.hydraulic.b`). Only the parameters selected by `paths` are varied; everything else stays fixed.\n\n---\n\n## `goodness`\n\n```julia\ngoodness(theta, model, forcing, lai, dates_UTC;\n    paths, lon, lat, depths_SM, depths_TS, FluxDay,\n    ignored...) -> (gof, data_sim, data_obs)\n```\n\nCalls `predict`, converts timestamps from UTC to local time (`+8 h`), then delegates to `BEPS_GOF` for scoring.\n\nThe `ignored...` keyword absorbs extra kwargs passed from `loss`/`optim` (e.g. `goal`, `goal_multiplier`) so callers don't need separate keyword lists for the two layers.\n\nReturns:\n- `gof` — named tuple `(Flux, SM, TS)`, each a `DataFrame` with one row per variable.\n- `data_sim`, `data_obs` — aligned daily DataFrames with columns `date`, `GPP`, `ET`, `Hs`, `SM_*cm`, `TS_*cm`.\n\n---\n\n## `loss`\n\n```julia\nloss(theta, model, forcing, lai, dates_UTC;\n    goal=:NSE, goal_multiplier=-1, kw_loss...) -> Float64\n```\n\nCollapses `goodness` to a scalar:\n\n```\nmean(gof.Flux[1:2, goal]) * goal_multiplier\n```\n\nRows 1–2 of `gof.Flux` are GPP and ET; `Hs` (sensible heat) is excluded by default. `goal` selects the metric column (`:NSE`, `:KGE`, or any column produced by `ModelParams.GOF`). The default `goal_multiplier = -1` turns maximisation of NSE into minimisation for SCE-UA.\n\n---\n\n## `optim`\n\n```julia\noptim(model, forcing, lai, dates_UTC;\n    paths, maxn=200, kw_loss...) -> theta::Vector\n```\n\nWraps SCE-UA (`ModelParams.sceua`) around the loss function:\n\n1. Extracts parameter bounds (`params.bound`) and initial values (`params.value`) from `parameters(model; paths)`.\n2. Constructs a closure `_loss(theta)` that forwards all keyword arguments to `loss`.\n3. Calls `sceua(_loss, u0, lb, ub; maxn, verbose=true, parallel=true)`.\n4. Returns the best-found parameter vector `theta`.\n\n`maxn` is the total number of model evaluations. SCE-UA runs complex shuffling in parallel (`parallel=true`), so wall time scales roughly with `maxn / ncores` rather than linearly.\n\n---\n\n## `BEPS_GOF`\n\n```julia\nBEPS_GOF(df_fluxes, states, dates_hour, FluxALL;\n    depths_SM, depths_TS) -> (gof, data_sim, data_obs)\n```\n\nThe diagnostic core. Handles unit reconciliation, temporal aggregation, depth interpolation, and date alignment before computing scores.\n\n**Steps:**\n\n1. **Normalise observations** — `normalize_flux_obs!(FluxALL)` applies unit conversions in-place (defined in `check_forcing.jl`).\n2. **Extract observation columns** — GPP, ET, Hs, and dynamically-named SM/TS columns (`SM_5cm`, `TS_10cm`, etc.) constructed from `depths_SM`/`depths_TS` (units: m → column suffix in cm).\n3. **Daily aggregation of simulated fluxes** — `agg_daily(df_fluxes, dates_hour)` collapses the hourly `Flux` DataFrame.\n4. **Daily aggregation of soil profiles** — `agg_daily(SM_sim, dates_hour)` and `agg_daily(TS_sim, dates_hour)` aggregate the full `N_layer × N_hour` profile matrices from `states.vectors.θ` and `states.vectors.Tsoil_c`.\n5. **Depth interpolation** — `interp_depths(SM_day, depths_SM)` interpolates the multi-layer soil profile to the observation depths, aligning model layers with sensor positions.\n6. **Date alignment** — if simulation and observation date ranges differ (missing-data gaps, partial years), `_align_daily_data` trims both to their intersection.\n7. **GOF scoring** — `ModelParams.GOF` computes NSE, KGE, RMSE, Bias, R² for each variable. Results are collected into three DataFrames: `gof.Flux` (GPP, ET, Hs), `gof.SM` (one row per depth), `gof.TS` (one row per depth).\n\n### Variable naming convention\n\nDepth-indexed columns follow the pattern `SM_{d}cm` / `TS_{d}cm` where `d = Int(depth_m * 100)`. The same convention is used in both the input observation DataFrame and the output `data_sim`/`data_obs`, so column names align automatically.\n\n---\n\n## Integration with `case01_ChinaFlux.jl`\n\nThe external entry point `RunModel` (in `Project_ChinaFlux/case01_ChinaFlux.jl`) calls `BEPS_GOF` directly — bypassing `optim` — for post-hoc evaluation of a fixed parameter set. This is the primary diagnostic/plotting path. The call graph confirms:\n\n```\nRunModel → BEPS_GOF → _align_daily_data\nRunModel → BEPS_GOF → _daily_obs_dates\n```\n\n`_daily_obs_dates` reads the date column from `FluxALL` to determine whether an explicit observation calendar exists; if `FluxALL` has no date column (or it is `nothing`), simulation dates are used as the reference.\n\n---\n\n## Adding a new calibration target\n\nTo add a new flux variable (e.g. soil heat flux `G`):\n\n1. Ensure `df_fluxes` contains the column and `FluxALL` has the matching observation column.\n2. Add a row to the `gof_Flux` construction in `BEPS_GOF`:\n   ```julia\n   (; var=\"G\", GOF(data_obs.G, data_sim.G)...)\n   ```\n3. Update `normalize_flux_obs!` in `check_forcing.jl` if unit conversion is needed.\n4. Adjust the `loss` row-index slice (`gof.Flux[1:2, goal]`) to include or exclude the new variable from the scalar objective.","other-01":"# Other — 01_参数化方案\n\n# Parameterization Scheme (`01_参数化方案`)\n\nThis module defines the parameter type hierarchy for BEPS. Every tuneable physical constant lives here — with default values, calibration bounds, and physical units co-located in a single struct declaration.\n\n## Parameter Hierarchy\n\n```mermaid\ngraph TD\n  ParamBEPS --> ParamSoilHydraulicLayers\n  ParamBEPS --> ParamSoilThermalLayers\n  ParamBEPS --> ParamVeg\n  ParamSoilHydraulicLayers --> ParamSoilHydraulic[\"ParamSoilHydraulic × N\"]\n  ParamSoilThermalLayers --> ParamSoilThermal[\"ParamSoilThermal × N\"]\n```\n\n`ParamBEPS` is the single object passed through the model. The three sub-structs are logically independent and can be initialized, tuned, or serialized separately.\n\n---\n\n## `ParamSoilHydraulic{FT}`\n\nEncodes the **Campbell (1974) soil water retention curve** for one soil layer.\n\n| Field | Default | Bounds | Units | Description |\n|-------|---------|--------|-------|-------------|\n| `θ_vfc` | 0.30 | (0.10, 0.45) | — | Volumetric field capacity |\n| `θ_vwp` | 0.10 | (0.02, 0.30) | — | Volumetric wilting point |\n| `θ_sat` | 0.45 | (0.25, 0.70) | — | Porosity / volumetric saturation |\n| `K_sat` | 5.0 | (0.01, 50.0) | cm h⁻¹ | Saturated hydraulic conductivity |\n| `ψ_sat` | −0.5 | (−2.0, −0.01) | m | Matric potential at saturation |\n| `b` | 5.0 | (1.5, 15.0) | — | Campbell shape parameter |\n\n`ψ_sat` and `b` together define the Brooks-Corey / Campbell ψ(θ) and K(θ) curves used throughout `SoilPhysics/`. The retention relationship is ψ = ψ_sat · (θ/θ_sat)^(−b).\n\n---\n\n## `ParamSoilThermal{FT}`\n\nProvides layer-level thermal properties for the heat diffusion solver in `UpdateHeatFlux.jl`.\n\n| Field | Default | Bounds | Units | Description |\n|-------|---------|--------|-------|-------------|\n| `κ_dry` | 0.2 | (0.05, 0.5) | W m⁻¹ K⁻¹ | Dry soil thermal conductivity |\n| `ρ_soil` | 1300.0 | (800, 1800) | kg m⁻³ | Bulk density |\n| `V_SOM` | 0.02 | (0.0, 0.3) | — | Organic matter volume fraction |\n\nActual effective conductivity is computed at runtime from these properties combined with the current soil moisture state.\n\n---\n\n## `ParamVeg{FT}`\n\nAll vegetation parameters for the overstory/understory canopy model.\n\n**Structural flags** (not calibrated):\n- `has_understory` — enables the two-layer canopy path in `inter_prg_jl`\n- `is_bforest` — selects deciduous vs. evergreen phenology logic\n\n**Canopy optics** — `α_canopy_vis/nir` control shortwave absorption by leaves; `α_soil_sat/dry` bound the two-end soil albedo used in `netRadiation.jl`.\n\n**Stomatal conductance** — Ball-Berry model:  \ngs = g1_w · A · rH / Cs + g0_w  \n`g1_w` (slope) and `g0_w` (intercept for H₂O) are the primary calibration targets for transpiration.\n\n**Photosynthesis** — `VCmax25` (Rubisco capacity at 25 °C) drives the Farquhar model. Alternatively, it can be estimated from `N_leaf` and `slope_Vc` via the Vcmax–leaf N regression used in `VCmax.jl`.\n\n**Root distribution** — `r_root_decay` is the exponential decay constant controlling how root density decreases with depth, used to compute layer-wise root fractions in `StateBEPS`.\n\n---\n\n## `ParamBEPS{FT}` — Top-level Container\n\n```julia\n@bounds @with_kw_noshow mutable struct ParamBEPS{FT<:AbstractFloat}\n  N::Int = 5\n  dz::Vector{FT} = FT[0.05, 0.10, 0.20, 0.40, 1.25]\n  r_drainage::FT = 0.50 | (0.2, 0.7)\n  ψ_min::FT = 33.0               # stress onset [m]\n  alpha::FT = 0.4                # He 2017 JGR-B Eq. 4 stress shape\n  hydraulic::ParamSoilHydraulicLayers{FT}\n  thermal::ParamSoilThermalLayers{FT}\n  veg::ParamVeg{FT}\nend\n```\n\n**Soil discretization** — `N=5` layers with thicknesses `dz` summing to ~2 m. The default profile (5, 10, 20, 40, 125 cm) matches the BEPS V2023 discretization. Changing `N` requires matching changes to `dz` and triggers re-initialization of the `*Layers` sub-structs.\n\n**Soil water stress** — `ψ_min` is the matric potential threshold at which plant stress begins (~0.1–0.33 MPa). `alpha` controls the shape of the stress reduction function (He et al. 2017 JGR-B, Eq. 4). Both influence `soil_water_factor_v2.jl`.\n\n**`r_drainage`** — fraction of surface water that drains laterally per timestep. Acts as a simple parameterization of overland flow; a Manning-formula replacement is noted as a future option in the code.\n\n---\n\n## The `@bounds @with_kw` Pattern\n\nAll structs use macros from `ModelParams.jl`:\n\n```julia\nθ_sat::FT = FT(0.45) | (0.25, 0.70)   # default | (lower, upper)\n```\n\n- `@with_kw` enables keyword construction and copies: `ParamVeg(VCmax25=60.0)`\n- `@bounds` attaches the `(lower, upper)` tuple as metadata accessible via `bounds(ps.veg, :VCmax25)`, used by the SCE-UA optimizer in `Optim.jl`\n- All structs are `mutable` so parameters can be updated in-place during calibration loops\n\n---\n\n## Initialization in Practice\n\nThe normal entry point is `setup()` in `DataType/setup.jl`, which constructs `ParamBEPS` from lookup tables keyed by vegetation type and soil texture class:\n\n```julia\nstate, ps = setup(\"evergreen_needleleaf\", \"silty_clay_loam\")\n# ps.hydraulic, ps.thermal populated from GlobalData.jl lookup tables\n# ps.veg populated from Param_Init.jl / InitParam_Veg\n```\n\nDirect construction is used in calibration workflows where individual parameters are swept or optimized against flux tower observations.","other-agents-md":"# Other — AGENTS.md\n\n# AGENTS.md — AI Agent Guide for BEPS.jl\n\n`AGENTS.md` is the authoritative onboarding document for AI coding agents working on this repository. It captures everything a model needs to contribute correctly without reading every source file first: project identity, physics scope, file layout, naming conventions, coding patterns, and workflow rules.\n\n---\n\n## Purpose and Audience\n\nThis file is not user documentation — it is a structured briefing for AI agents (Claude, Copilot, etc.) that will read and modify source code. It complements `CLAUDE.md` (detailed technical reference) and `README.md` (user-facing quickstart). When an agent is given a task in this repository, `AGENTS.md` is the first file it should internalize.\n\n---\n\n## What the Document Covers\n\n### Project Identity (§1)\n\nEstablishes that BEPS.jl is a Julia port of the Boreal Ecosystem Productivity Simulator, running at hourly resolution and simulating:\n\n- Canopy photosynthesis (sunlit/shaded leaf separation, Farquhar model)\n- Radiation transfer (direct/diffuse, VIS/NIR)\n- Energy balance (SH, LE, net radiation)\n- Soil physics: multi-layer water (Richards equation) and heat (thermal diffusion)\n- Snowpack (three-stage scheme)\n- Hydrology: interception, infiltration, root water uptake\n\nPerformance context: the Julia implementation runs ~2.5× faster than the reference C library (`libbeps`), which remains the physical ground truth for validation.\n\n### Tech Stack (§2)\n\nKey dependencies an agent needs to recognize:\n\n| Package | Role in code |\n|---|---|\n| `Parameters` | `@with_kw` struct definitions throughout `DataType/` |\n| `UnPack` | `@unpack`/`@pack!` for field extraction in physics functions |\n| `StaticArrays` | Performance-critical small arrays |\n| `ModelParams` | Parameter struct handling (`ParamVeg`, `ParamSoilHydraulicLayers`) |\n| `ComponentArrays` | Structured parameter arrays |\n\n### File Layout (§3)\n\nThe tree maps every directory to its physical meaning. Key areas:\n\n- `src/DataType/` — all structs: `StateBEPS`, `ParamBEPS`, `Met`, `Flux`, `ETFlux`\n- `src/SoilPhysics/` — Richards equation, thermal diffusion, water stress\n- `src/SPAC/` — leaf-scale calculations, LAI partitioning, Vcmax\n- `src/inter_prg.jl` — the hourly integration loop that calls every physics module\n- `test/modules/` — per-module unit tests, each validated against C output\n\n### Naming Conventions (§5)\n\nThe conventions section is critical for generating physically correct variable names:\n\n| Pattern | Meaning | Examples |\n|---|---|---|\n| Greek letters | Physical quantities | `θ` (moisture), `ψ` (potential), `κ` (conductivity), `ρ` (density) |\n| `_c` suffix | Current timestep value | `Tsoil_c`, `Tsnow_c` |\n| `_p` suffix | Previous timestep value | `Tsoil_p` |\n| `_sat` / `_vwp` / `_vfc` | Saturation / wilting point / field capacity | `θ_sat`, `K_sat` |\n| `f_` prefix | Fraction or stress factor | `f_soilwater`, `f_root` |\n| `r_` prefix | Rate | `r_rain_g`, `r_drainage` |\n| `z_` prefix | Depth | `z_snow`, `z_water` |\n\n**Hard rule:** use `T` or `Tsoil` for temperature — never `temp`.\n\n### Integration Loop Logic (§6.3)\n\nThe sub-hourly loop in `inter_prg.jl` runs 10 iterations per hour (timestep `kstep = 360 s`). The documented call sequence is the correct ordering for all physics:\n\n```\nsnowpack_stage1 → rainfall_stage1 → soil_water_factor_v2\n→ canopy energy balance → transpiration → evaporation_canopy\n→ rainfall_stage2 → snowpack_stage2 → evaporation_soil\n→ surface_temperature! → snowpack_stage3 → sensible_heat\n→ UpdateHeatFlux → Root_Water_Uptake → UpdateSoilMoisture\n```\n\nAltering this ordering without understanding physical coupling will produce incorrect results.\n\n### Safety: Numerical Clamping (§5.3)\n\nAll agents must apply physical bounds when writing or modifying state updates:\n\n```julia\nclamp(θ, θ_vwp, θ_sat)        # soil moisture\nclamp(T, Ta - 25, Ta + 25)    # temperatures relative to air temp\nclamp(G, -200, 200)            # heat flux [W/m²]\nclamp(ice_ratio, 0.0, 1.0)    # frozen fraction\n```\n\nMissing clamps on state variables is a common source of physically impossible values propagating across timesteps.\n\n### Testing Philosophy (§7)\n\nThe primary validation method is Julia-vs-C comparison. Every physics module has a corresponding test that:\n\n1. Initialises identical state in both `StateBEPS` (Julia) and `Soil_c` (C struct layout)\n2. Runs both implementations for the same inputs\n3. Asserts field-by-field equality within tolerance (`tol=1e-7`)\n\nThis means the C library in `deps/BEPS.c` is the physical ground truth. When a Julia result diverges from C, the C result is correct unless there is an identified C bug.\n\n### Development Workflow (§10)\n\nWhen adding a new physics process, the prescribed sequence is:\n\n1. Define struct in `DataType/` with `@with_kw`\n2. Add initialization in `setup.jl` or `Param_Init.jl`\n3. Implement physics in a new module file\n4. Insert call at the correct position in `inter_prg.jl`\n5. Export from `BEPS_modules.jl`\n6. Write a C-comparison test in `test/modules/`\n7. Document in `docs/modules/`\n\nSkipping step 6 means the implementation has no physical validation.\n\n---\n\n## Relationship to Other Documentation\n\n```mermaid\ngraph TD\n    A[AGENTS.md] -->|links to| B[CLAUDE.md]\n    A -->|links to| C[README.md]\n    B --> D[src/ architecture detail]\n    B --> E[Type hierarchy]\n    C --> F[User quickstart]\n    A --> G[docs/modules/*.md]\n    G --> H[Per-module physics docs]\n```\n\n`AGENTS.md` is the entry point; `CLAUDE.md` is the deep reference an agent should follow up with for struct signatures and data flow detail.\n\n---\n\n## Key Invariants for Contributors\n\nThese constraints are stated in `AGENTS.md` and must not be violated:\n\n- **Language split:** Chinese for conversation and inline comments where helpful; English for all identifiers, docstrings, and commit messages.\n- **Physical clamping is mandatory** on all state writes — omitting it will cause silent numerical instability.\n- **C parity is the correctness criterion** for soil physics — if Julia and C diverge, investigate before changing either.\n- **`inter_prg.jl` call order is load-bearing** — the sequence encodes physical dependencies (e.g., soil moisture stress must precede canopy conductance).\n- **CI targets Windows** because `libbeps.dll` is Windows-only; Linux/macOS runs Julia-only mode and cannot run C-comparison tests.","other-artifacts-toml":"# Other — Artifacts.toml\n\n# Artifacts.toml — Prebuilt `libbeps` Binary Registry\n\n`Artifacts.toml` declares the prebuilt C library artifacts that Julia's `Pkg` artifact system uses to fetch and verify `libbeps` — the compiled C core of BEPS.jl — at install time.\n\n## What It Does\n\nRather than requiring users to compile the BEPS C model from source, `Artifacts.toml` encodes download URLs, content hashes, and platform selectors for pre-built `libbeps` shared library archives. Julia's `Pkg.Artifacts` machinery reads this file to select the correct binary for the current platform, download it once, cache it in the Julia depot, and make it available to the package at runtime.\n\nThe `lazy = true` flag on all entries means the artifact is **not** fetched during `Pkg.instantiate()` — it is only downloaded when the code actually requests it (i.e., when `libbeps` is first loaded).\n\n## Covered Platforms\n\n| OS | Architecture | Archive |\n|---|---|---|\n| Windows | x86\\_64 | `beps-windows-x86_64.tar.gz` |\n| Linux (glibc) | x86\\_64 | `beps-linux-x86_64.tar.gz` |\n| macOS | x86\\_64 (Intel) | `beps-macos-x86_64.tar.gz` |\n| macOS | aarch64 (Apple Silicon) | `beps-macos-arm64.tar.gz` |\n\nAll four archives come from the `v2026.05.03` release of [`CUG-hydro/BEPS.c`](https://github.com/CUG-hydro/BEPS.c).\n\n## Entry Structure\n\nEach `[[libbeps]]` block is a platform-specific artifact entry:\n\n```toml\n[[libbeps]]\narch = \"x86_64\"\ngit-tree-sha1 = \"...\"   # SHA1 of the unpacked directory tree (integrity check)\nlazy = true              # deferred download — only fetched on first use\nos = \"linux\"\nlibc = \"glibc\"           # optional: further narrows the match (Linux only)\n\n    [[libbeps.download]]\n    sha256 = \"...\"       # SHA256 of the .tar.gz archive itself\n    url = \"...\"          # HTTPS URL to download from\n```\n\n`git-tree-sha1` is the hash Julia uses to verify the *unpacked* tree after extraction. `sha256` covers the archive before extraction. Both must match for `Pkg.Artifacts` to accept the download.\n\n## Connection to the Rest of the Codebase\n\nThe artifact is consumed by `src/clang/BEPS_c.jl`, which loads the shared library at module init:\n\n```julia\n# src/clang/BEPS_c.jl (simplified)\nconst libbeps = LazyArtifacts.find_artifact_dir(\"libbeps\")\n```\n\nFrom there, `src/BEPS.jl` conditionally exposes the C execution path (`inter_prg_c`) alongside the pure-Julia path (`inter_prg_jl`). The top-level `beps_main(...; version=\"c\")` call routes through this layer.\n\n```mermaid\ngraph LR\n    A[Artifacts.toml] -->|depot cache| B[libbeps .so/.dll/.dylib]\n    B --> C[src/clang/BEPS_c.jl]\n    C --> D[inter_prg_c]\n    D --> E[beps_main version=c]\n```\n\n## Updating the Artifacts\n\nWhen a new `libbeps` release is cut from `CUG-hydro/BEPS.c`:\n\n1. Build the platform archives and publish them as GitHub Release assets.\n2. Compute `sha256` for each `.tar.gz` (`sha256sum` / `certutil`).\n3. Extract each archive and compute `git-tree-sha1` with `julia -e 'using Pkg.GitTools; println(bytes2hex(GitTools.tree_hash(\"<dir>\")))'`.\n4. Update the corresponding `[[libbeps]]` block with the new URL, `sha256`, and `git-tree-sha1`.\n\nThe `ArtifactUtils.jl` package (`bind_artifact!`) can automate steps 2–4.","other-case01-ipynb":"# Other — case01.ipynb\n\n# case01.ipynb — Single-Site BEPS Evaluation Notebook\n\nThis notebook runs the BEPS model against a single ChinaFlux forest site and evaluates simulated GPP and ET against eddy-covariance observations. It is the primary interactive entry point for per-site model validation and parameter tuning.\n\n## Purpose\n\nThe notebook exercises the full `simulate` pipeline for one site from the ChinaFlux 2026 forest dataset (12 sites, deciduous broadleaf forest class), computes daily goodness-of-fit statistics, and exposes intermediate state and canopy conductance time series for diagnostic plotting.\n\n## Workflow\n\n```mermaid\nflowchart LR\n    A[Hourly Met CSV] --> B[MetSeries]\n    C[Daily FluxLAI CSV] --> D[lai / GPP_obs / ET_obs]\n    B & D --> E[simulate]\n    F[ParamBEPS + setup] --> E\n    E --> G[df_fluxes / df_ET]\n    G --> H[agg_daily → GOF]\n    E --> I[states / caches]\n```\n\n## Input Data\n\nTwo CSV files are read from `data-raw/BEPS/` inside the ChinaFlux2026 package directory:\n\n| File | Frequency | Key columns |\n|---|---|---|\n| `Forcing_Met_Hourly_BEPS_Forest_sp12_hourly_v20260507.csv` | Hourly | `site`, `time`, `Ta_canopy`, `RH_canopy`, `WS_canopy`, `Rs`, `Rln_in`, `Prcp` |\n| `Forcing_FluxLAI_Daily_BEPS_Forest_sp12_v20260506.csv` | Daily | `site`, `GPP`, `ET`, `LAI_glass_G005` |\n\nBoth files cover all 12 forest sites; the notebook filters to a single site with `FORCING.site .== SITE` where `SITE = SITES[8]`.\n\n## Helper Functions\n\n### `agg_daily(df_fluxes, dates)`\n\nAggregates hourly model output to daily totals using `nansum`. Returns a named tuple:\n\n```julia\n(; GPP_sim, ET_sim, dates_day)\n```\n\n`ET_sim` is the sum of `Trans + Evap` columns from `df_fluxes`. The `apply(...; by=dates_day, fun)` call groups by `Date.(dates)`.\n\n### `parse_time(x::AbstractString)`\n\nParses ISO 8601 timestamps from the forcing CSV (`yyyy-mm-ddTHH:MM:SSZ` format) into `DateTime`. Applied as a broadcast over the `time` column before constructing `MetSeries`.\n\n## Model Configuration\n\nThe site used here (index 8) is a deciduous broadleaf forest (DBF) at Baotianman (宝天曼), 115.06°E, 26.74°N.\n\n```julia\nVegType = 1      # deciduous broadleaf forest\nSoilType = 8\nmodel = ParamBEPS(VegType, SoilType)\nmodel.veg.z_wind = 39.6      # measurement height (m)\nmodel.veg.VCmax25 = 56.25    # maximum carboxylation rate (µmol m-2 s-1)\nmodel.veg.g1_w   = 4.8       # Ball-Berry stomatal slope\nclumping = 0.58              # canopy clumping index\n```\n\nInitial conditions are set to field capacity (`θ_vfc[1]`) and the first-hour air temperature for soil temperature. `setup(model; Ta, Tsoil, θ0, z_snow)` returns an initialized `(state, model)` pair.\n\n## Running the Simulation\n\n```julia\ndf_fluxes, df_ET, states, caches = simulate(forcing, lai, dates_model;\n    ps=model, state,\n    VARS_STATE=VARS_STATE, VARS_CACHE=VARS_CACHE,\n    lon=115.06, lat=26.74, clumping)\n```\n\n`dates_model` is `dates .- Hour(8)` — local CST timestamps shifted to UTC before passing to the model.\n\n**Unit conversions applied to observations before comparison:**\n- `GPP_obs = -GPP_obs` — flips GEE sign convention (negative uptake) to GPP (positive)\n- `ET_obs = ET_obs * 1e6 / 86400` — converts MJ m⁻² day⁻¹ to W m⁻²\n\n## Outputs\n\n### Flux output (`df_fluxes`)\n\nHourly `DataFrame` with at minimum `GPP`, `Trans`, `Evap` columns. Aggregated to daily totals by `agg_daily` for GOF evaluation.\n\n### State series (`states`)\n\n`StateSeries` holding 10 variables at every hourly timestep across 113,952 steps (~13 years):\n\n| Variable | Shape | Description |\n|---|---|---|\n| `z_water` | (113952,) | Surface water depth |\n| `ρ_snow` | (113952,) | Snow density |\n| `z_snow` | (113952,) | Snow depth |\n| `r_rain_g` | (113952,) | Ground rain rate |\n| `f_soilwater` | (113952,) | Soil water stress factor |\n| `θ` | (5, 113952) | Volumetric soil moisture, 5 layers |\n| `Tsoil_c` | (5, 113952) | Soil temperature (°C), 5 layers |\n| `ETi` | (5, 113952) | Layer evapotranspiration |\n| `r_waterflow` | (5, 113952) | Inter-layer water flow |\n| `G` | (5, 113952) | Soil heat flux, 5 layers |\n\n### Cache series (`caches`)\n\nPer-leaf conductance and concentration variables, shape (113952, 4) — 4 leaf types (overstory sunlit/shaded, understory sunlit/shaded):\n\n`Ci_new`, `Cs_new`, `Tc_new`, `Gs_new`, `Gc`, `Gh`, `Gw`, `Gw_wet`\n\n### Goodness-of-fit\n\nDaily GPP and ET are evaluated with `GOF(obs, sim)` from `RTableTools`, which returns NSE, R², KGE, R, RMSE, MAE, bias, and bias_perc. The example run achieves GPP NSE = 0.84, ET NSE = 0.79 over ~13 years of data.\n\n## Extending to Other Sites\n\nChange `SITE = SITES[N]` and update `lon`, `lat`, `z_wind`, and the vegetation/soil type integers to match the target site. The rest of the pipeline is site-agnostic.","other-ch01-soilmoisture":"# Other — ch01_SoilMoisture\n\n# Soil Moisture Module\n\nThe soil moisture module implements a 1D, 5-layer Richards-equation solver coupled to the SPAC model. It tracks volumetric water content (θ), drives root water uptake, and feeds water stress signals back to photosynthesis and stomatal conductance. The solver runs as a substep loop inside each hourly `inter_prg_jl` call.\n\n---\n\n## Soil Column Layout\n\nThe column is discretized into five layers with fixed thicknesses:\n\n| Layer | `dz` (m) | Depth to bottom (m) |\n|-------|-----------|----------------------|\n| 1     | 0.05      | 0.05                 |\n| 2     | 0.10      | 0.15                 |\n| 3     | 0.20      | 0.35                 |\n| 4     | 0.40      | 0.75                 |\n| 5     | 1.25      | 2.00                 |\n\nAll per-layer arrays in `StateBEPS` (`θ`, `Tsoil_c`, `ice_ratio`) and `HydraulicProfile` (`K_sat`, `ψ_sat`, `b`, `θ_sat`, `θ_vwp`, `κ_dry`) are indexed 1–5 following this layout. Layer thicknesses `dz` live in `ParamBEPS.hydraulic`.\n\n---\n\n## State and Parameters\n\n**State** (in `StateBEPS`, updated each timestep):\n\n| Field | Description | Units |\n|-------|-------------|-------|\n| `θ[i]` | Volumetric water content, layer i | m³ m⁻³ |\n| `z_water` | Surface ponding depth | m |\n| `Tsoil_c[i]` | Soil temperature, layer i | °C |\n| `ice_ratio[i]` | Ice fraction, layer i | — |\n\n**Parameters** (in `ParamBEPS.hydraulic`, type `HydraulicProfile`):\n\n| Field | Description | Units |\n|-------|-------------|-------|\n| `K_sat[i]` | Saturated hydraulic conductivity | cm h⁻¹ |\n| `ψ_sat[i]` | Saturated matric potential (negative) | cm |\n| `b[i]` | Campbell pore-size distribution index | — |\n| `θ_sat[i]` | Saturated volumetric water content | m³ m⁻³ |\n| `θ_vwp[i]` | Permanent wilting-point water content | m³ m⁻³ |\n| `r_drainage` | Surface retention fraction | — |\n| `kstep` | Outer timestep | s |\n\nCampbell constitutive relations tie `ψ` and `K` to `θ`:\n\n```\nψ_i = ψ_sat[i] · (θ[i] / θ_sat[i])^(-b[i])\nK_i = K_sat[i] · (θ[i] / θ_sat[i])^(2b[i] + 3)\n```\n\n---\n\n## Per-Timestep Flow\n\n```mermaid\ngraph TD\n    A[inter_prg_jl] --> B[soil_water_factor_v2<br/>water stress f_ψ, f_T → w_i]\n    B --> C[UpdateSoilThermalConductivity + Update_Cs]\n    C --> D[surface_temperature!<br/>energy balance → T_surf, E_soil]\n    D --> E[UpdateHeatFlux<br/>heat diffusion across layers]\n    E --> F[UpdateSoilMoisture!<br/>Richards solver + root uptake]\n    F --> G[StateBEPS updated]\n```\n\n`UpdateSoilMoisture!` is called once per hour but internally runs an adaptive substep loop.\n\n---\n\n## Surface Infiltration and Runoff\n\nEntry point: `UpdateSoilMoisture!` in `src/SoilPhysics/UpdateSoilMoisture.jl`.\n\n**Maximum infiltration capacity** (Darcy, layer 1):\n\n```\nI_max = K_sat[1] · [1 - (θ_sat[1] - θ[1]) / dz1_cm · (ψ_sat[1] · b[1]) / θ_sat[1]]\n```\n\nBecause `ψ_sat < 0`, the second term is positive — a drier surface draws more water in.\n\n**Available supply rate** (from ponding + precipitation, converted to cm h⁻¹):\n\n```\nA = max(z_water / kstep + r_rain_g, 0) × 3600 × 100\n```\n\n**Actual infiltration:**\n\n```\nI = min(A, I_max)\n```\n\n**Surface water balance** — excess after infiltration splits between runoff and new ponding:\n\n```\nD_ex  = (A_SI - I_SI) · kstep          # [m]\nR_s   = D_ex · (1 - r_drainage)\nz_water(t+1) = D_ex · r_drainage\n```\n\n`r_drainage` ∈ [0.2, 0.7] (default 0.5) controls how much excess water lingers on the surface vs. immediately leaves as runoff. Flat terrain → higher values.\n\n---\n\n## Interlayer Hydraulic Conductivity\n\nThe interface conductivity between layers i and i+1 is computed in `UpdateSoilMoisture!` as:\n\n```\nK_half = (K[i]·ψ[i] + K[i+1]·ψ[i+1]) / (ψ[i] + ψ[i+1])\n       · (b[i] + b[i+1]) / (b[i] + b[i+1] + 6)\n```\n\nThe first factor is a ψ-weighted harmonic-style mean; the second is a Campbell nonlinearity correction. This formula treats the b values from each layer separately at their endpoints rather than pre-averaging them — the correction `b_eff / (b_eff + 6)` at a heterogeneous interface is a flux-weighted result, not an arithmetic mean of b.\n\n**Downward flux** (cm h⁻¹, positive = downward):\n\n```\nQ[i] = K_half · (1 - (ψ[i+1] - ψ[i]) / dz_half_cm)\n```\n\nGravity contributes the `+1` term; the capillary gradient term uses negative-convention pressure heads.\n\n**Flux cap** to prevent over-draining layer i+1 within one substep:\n\n```\nQ_max[i] = ((θ_sat[i+1] - θ[i+1]) · dz[i+1] / dt + S[i+1]) × 3600 × 100\nQ[i] ← min(Q[i], Q_max[i])\n```\n\n---\n\n## Moisture Update and Adaptive Substeps\n\nBoundary conditions: `Q[0] = I` (infiltration), `Q[N] = 0` (no-flow bottom).\n\nExplicit update per substep `δt`:\n\n```\nθ[1] += (I - Q[1] - S[1]) · δt / dz[1]\nθ[i] += (Q[i-1] - Q[i] - S[i]) · δt / dz[i],   i = 2..N\nθ[i] ← clamp(θ[i], θ_vwp[i], θ_sat[i])\n```\n\nψ and K are recalculated at the start of every substep. Substep length is selected based on the current maximum flux magnitude:\n\n| max |Q[i]| (cm h⁻¹) | δt (s) |\n|------------------------|--------|\n| > 3.6                  | 1      |\n| > 0.36                 | 30     |\n| otherwise              | 360    |\n\nThe final substep is clipped to the remaining time in the outer step.\n\n---\n\n## ET Sink Terms\n\n`S[i]` (m s⁻¹) is the total water extraction from layer i, computed before `UpdateSoilMoisture!` is called.\n\n```\nS[1] = (T / ρ_w) · w[1] + E_s / ρ_w\nS[i] = (T / ρ_w) · w[i],    i = 2..N\n```\n\n`T` = total canopy transpiration (overstory + understory), `E_s` = soil evaporation from `evaporation_soil.jl`. The per-layer weights `w[i]` come from `soil_water_factor_v2`.\n\n---\n\n## Root Water Uptake Weights\n\nImplemented in `src/SoilPhysics/soil_water_factor_v2.jl` and `UpdateRootFraction!`.\n\n**Root distribution** (exponential decay, `r_root_decay = β`):\n\n```\nf_root[1]   = 1 - β^Z[1]\nf_root[i]   = β^Z[i-1] - β^Z[i],   2 ≤ i < N\nf_root[N]   = β^Z[N-1]\nZ[i] = 100 · Σ dz[1..i]              # cumulative depth in cm\n```\n\n**Water stress** (`f_ψ`) and **temperature stress** (`f_T`):\n\n```\nχ_i = -ψ[i]   # suction head, positive\n\nf_ψ[i] = 1                                            if χ ≤ ψ_min\nf_ψ[i] = [1 + ((χ - ψ_min) / ψ_min)^α]^(-1)         if χ > ψ_min\n\nf_T[i] = 1 - exp(-0.02 · T_soil[i]²)   if T_soil > 0\nf_T[i] = 0                              if T_soil ≤ 0\n```\n\n`ψ_min` and `α` are PFT-specific:\n\n| PFT | ψ_min (cm) | α |\n|-----|-----------|---|\n| DBF, EBF | 1000 | 1.5 |\n| Others | 3300 | 0.4 |\n\nψ_min = 1000 cm ≈ 0.10 MPa and ψ_min = 3300 cm ≈ 0.33 MPa are empirical onset thresholds, not stomatal closure points (which occur around 0.66–2.75 MPa for trees).\n\n**Normalized weights:**\n\n```\na[i] = f_root[i] · f_ψ[i] · f_T[i]\nw[i] = a[i] / Σ a[j]\n```\n\nThese weights flow into `S[i]` above and also directly into the `f_water` field of `StateBEPS`, which gates stomatal conductance in `photosynthesis.jl`.\n\n---\n\n## Connection to the Rest of the Model\n\n| Caller / Downstream | What it receives |\n|---------------------|-----------------|\n| `inter_prg_jl` | Calls `soil_water_factor_v2` → `UpdateHeatFlux` → `UpdateSoilMoisture!` in sequence |\n| `photosynthesis.jl` | Reads `state.f_water` (derived from `w[i]`) to scale Vcmax and stomatal conductance |\n| `evaporation_soil.jl` | Provides `E_s` that enters `S[1]` |\n| `surface_temperature!` | Reads `θ[1]` and `Tsoil_c[1]` for surface energy balance; updates them after convergence |\n| `UpdateHeatFlux` | Uses `κ_dry`, `ice_ratio`, `θ` to compute thermal conductivity and heat capacity before running the heat diffusion step |\n| `simulate` (modern API) | Accepts `SM_obs` keyword to skip `UpdateSoilMoisture!` entirely and prescribe θ from observations |\n\n---\n\n## JAX-Style API\n\nAll physics functions have been adapted to accept `(st::StateBEPS, ps::ParamBEPS)` rather than the legacy `(soil::Soil)` signature, enabling gradient-based calibration and ensemble runs without mutating shared parameter structs. The old `Soil` struct and `Sync_Param_to_Soil!` are retained for C-path compatibility only. `surface_temperature!` is the one function not yet adapted (pending `TransientCache` refactor).","other-ch02-soiltemperature":"# Other — ch02_SoilTemperature\n\n# Soil Temperature Module (`ch02_SoilTemperature`)\n\n## Overview\n\nThis module simulates 1D vertical heat conduction through soil and snowpack. It solves for surface and snow-layer temperatures at each timestep, computes the upper boundary flux `G[1]` entering the soil column, propagates heat through soil layers, and updates the freeze-thaw state.\n\nThe computation chain at each sub-hourly step is:\n\n```\nUpdateThermal(κ, Cv) → surface_temperature! → G[1] → UpdateHeatFlux → Update_ice_ratio\n```\n\n---\n\n## Grid Layout\n\nThe soil column uses a node-centered grid with `N+2` indices:\n\n| Index | Layer | Temperature | Flux `G[i]` |\n|-------|-------|-------------|-------------|\n| 1 | Surface skin | `T_surf` | `G[1]` — upper boundary flux into soil |\n| 2 | Soil layer 1 | `T_soil[1]` | flux leaving layer 1 bottom |\n| … | … | … | … |\n| N+1 | Soil layer N | `T_soil[N]` | flux leaving layer N bottom |\n| N+2 | Deep boundary | `T_ann_mean` (virtual) | bottom boundary flux |\n\nThe default configuration uses `N=5` layers with thicknesses `[0.05, 0.10, 0.20, 0.40, 1.25] m`. Note the index offset: `soil.κ[1]` (first physical layer) maps to `κ[2]` in the working array passed to `surface_temperature_jl`.\n\n---\n\n## Thermal Parameters\n\n### Volumetric Heat Capacity `Cv[i]`\n\nThree-component mix (Chen 2007, Eq. 18):\n\n$$C_{v,i} = 2\\times10^6 \\frac{\\rho_i}{2650} + 10^6 \\theta_i [4.2(1-f_{ice,i}) + 2.09 f_{ice,i}] + 2.5\\times10^6 V_{SOM,i}$$\n\nwhere $\\rho_i$ is bulk density, $\\theta_i$ is volumetric water content, $f_{ice,i}$ is ice fraction, and $V_{SOM,i}$ is organic matter volume fraction.\n\n### Soil Thermal Conductivity `κ[i]`\n\n$$\\kappa_i = \\left(\\kappa_{dry,i}^{1-\\theta_{sat,i}} \\cdot k_{ice}^{1.2\\theta_i f_{ice,i}} \\cdot k_w^{\\theta_i(1-f_{ice,i})} - 0.15\\right) \\frac{\\theta_i}{\\theta_{sat,i}} + 0.15$$\n\nwith $k_{ice} = 2.1\\,\\text{W m}^{-1}\\text{K}^{-1}$, $k_w = 0.61\\,\\text{W m}^{-1}\\text{K}^{-1}$. Result is clamped to $\\geq 0.15$.\n\n### Snow Thermal Conductivity `κ_dry_snow` — `cal_κ_snow`\n\nJordan (1991) empirical fit:\n\n$$\\kappa_{snow} = 0.021 + 4.2\\times10^{-4}\\rho_{snow} + 2.2\\times10^{-9}\\rho_{snow}^3$$\n\n---\n\n## Surface Temperature Solvers\n\nTwo numerical methods are used depending on whether the layer is at the surface (needs to account for aerodynamic exchange) or an interior node.\n\n### `solve_imp` — Implicit Surface Temperature\n\nSolves the surface energy balance for a control volume of thickness $\\Delta z_M$:\n\n$$\\Delta M (T^* - T^t) = G_g \\frac{z_{rad}}{z_c} - \\rho_a c_p \\frac{T^* - T_{up}}{r_a} - \\eta_c \\kappa_{bot} \\frac{T^* - T_{bot}}{z_c}$$\n\nwhere $\\Delta M = C_v \\Delta z_M / \\Delta t$. Multiplying through by $r_a z_c$ and collecting $T^*$:\n\n$$T^* = \\frac{T^t \\cdot I + G_g r_a z_{rad} + \\rho_a c_p T_{up} z_c + \\eta_c r_a \\kappa_{bot} T_{bot}}{\\rho_a c_p z_c + \\eta_c r_a \\kappa_{bot} + I}, \\quad I = \\Delta M \\cdot r_a z_c$$\n\nKey parameters:\n- `z` ($z_c$): characteristic exchange distance passed to the function\n- `z_rad`: depth over which net radiation $G_g$ acts (may differ from $z_c$)\n- `η_c`: conductance scaling factor; for bare-soil case `η_c = 2.0` (half-layer distance to node center)\n- `T_up`: upper boundary (air or lower canopy temperature)\n- `T_bot`: lower boundary temperature (soil or lower snow node)\n\nResult is clamped to $[\\mu - 25,\\, \\mu + 25]$ °C relative to a reference temperature $\\mu$ to suppress nonphysical jumps.\n\n> **Important distinction**: $\\Delta z_M$ (storage thickness) and $z_c$ (exchange distance) are not always the same. In the deep-snow case, the top snow node stores heat over `dz_snow_s1 = 0.02 m` but exchanges with the layer below over `dz_snow_s12 = 0.04 m`.\n\n### `step_exp` — Explicit Interior Update\n\nFor interior snow nodes (no aerodynamic exchange):\n\n$$T^{t+1} = T^t + \\frac{F_{in} - F_{out}}{C_v \\Delta z} \\Delta t$$\n\n---\n\n## `surface_temperature_jl` / `surface_temperature!`\n\nThe main entry point. Dispatches to one of three snow-cover regimes based on `z_snow`.\n\n```mermaid\nflowchart TD\n    A[surface_temperature_jl] --> B{z_snow?}\n    B -- ≤ 0.02 m --> C[Case 1: No snow<br/>solve_imp for T_surf<br/>G1 = κ₁·ΔT/Δz]\n    B -- 0_02[\"0.02\"]–0.05 m --> D[Case 2: Patchy snow<br/>solve T_soil0, T_snow0<br/>weighted T_mix0, T_surf]\n    B -- > 0.05 m --> E[Case 3: Deep snow<br/>3-layer snow model<br/>T_snow0 implicit<br/>T_snow1,2 explicit]\n    C --> F[Return G1, T_surf, T_mix0, SnowLand]\n    D --> F\n    E --> F\n```\n\n### `SnowLand` State Variables\n\n| Variable | Symbol | Active in |\n|----------|--------|-----------|\n| `T_snow0` | $T_{snow,0}$ | Patchy snow, Deep snow — top snow node / mixed surface temp |\n| `T_snow1` | $T_{snow,1}$ | Deep snow only — middle snow node |\n| `T_snow2` | $T_{snow,2}$ | Deep snow only — bottom snow node |\n| `T_mix0` | $T_{mix,0}$ | All cases — snow-soil interface weighted temperature |\n| `T_surf` | $T_{surf}$ | All cases — upward longwave radiation temperature |\n\n---\n\n## Three Snow-Cover Cases\n\n### Case 1: No Snow / Negligible Snow (`z_snow ≤ 0.02 m`)\n\nSingle-layer model. `T_surf = T_mix0 = T_soil0 = T_snow0`.\n\nSurface temperature from `solve_imp` with `T_bot = T_soil1_last`, then:\n\n$$G_1 = \\kappa_1 \\frac{T_{surf} - T_{soil,1}^t}{\\Delta z_1 / 2}$$\n\nFlux clamped to $[-100, 100]\\,\\text{W m}^{-2}$.\n\n### Case 2: Thin / Patchy Snow (`0.02 < z_snow ≤ 0.05 m`)\n\nTwo-tile approach: bare soil fraction $(1-f_{snow})$ and snow-covered fraction $f_{snow}$.\n\n**Snow-soil interface temperature** `T_mix0` (weighted energy balance):\n\n$$T_{int} = \\frac{\\kappa_1 T_{soil,1}^t / \\Delta z_{soil} + \\kappa_{snow} T_{snow,0} / z_{snow} + \\Delta M_1 T_{mix,0}^t}{\\kappa_1 / \\Delta z_{soil} + \\kappa_{snow} / z_{snow} + \\Delta M_1}$$\n\nwhere $\\Delta M_1 = C_{v,1} \\times 0.02\\,\\text{m} / \\Delta t$ — an artificial 2 cm interface storage layer that gives the interface temperature temporal memory.\n\n> **Known issue (TS_Case02.md)**: The snow conductance term in the numerator should be $\\kappa_{snow} T_{snow} / z_{snow}$, but the current code omits the $z_{snow}$ divisor, producing a dimensional error. The denominator has the correct $\\kappa_{snow}/z_{snow}$ form.\n\nArea-weighted surface temperatures:\n\n$$T_{surf} = (1 - f_{snow}) T_{soil,0} + f_{snow} T_{snow,0}$$\n$$T_{mix,0} = (1 - f_{snow}) T_{soil,0} + f_{snow} T_{int}$$\n\nUpper boundary flux:\n\n$$G_1 = f_{snow} \\frac{\\kappa_{snow}(T_{snow,0} - T_{soil,1}^t)}{z_{snow} + \\Delta z_{soil}} + (1 - f_{snow}) \\frac{\\kappa_1 (T_{soil,0} - T_{soil,1}^t)}{\\Delta z_{soil}}$$\n\n### Case 3: Deep Snow (`z_snow > 0.05 m`)\n\nThree-node snowpack. `T_snow0` solved implicitly; `T_snow1`, `T_snow2`, and `T_mix0` stepped explicitly.\n\nLayer layout:\n\n| Node | Variable | Thickness |\n|------|----------|-----------|\n| Top | `T_snow0` | `dz_snow_s1 = 0.02 m` |\n| Mid | `T_snow1` | `dz_snow_s2 = 0.02 m` |\n| Bot | `T_snow2` | `z_snow − 0.04 m` |\n| Sub-snow soil | `T_mix0` | 0.02 m (virtual) |\n\nInternal snow fluxes:\n\n$$G_{snow,0} = \\kappa_{snow}(T_{snow,0} - T_{snow,1}^t) / 0.04$$\n$$G_{snow,1} = \\kappa_{snow}(T_{snow,1}^t - T_{snow,2}^t) / (z_{snow} - 0.02)$$\n$$G_{snow,2} = \\frac{T_{snow,2}^t - T_{mix,0}^t}{0.5(z_{snow}-0.04)/\\kappa_{snow} + 0.02/\\kappa_1}$$\n\nSub-snow soil surface update:\n\n$$T_{mix,0}^{t+1} = T_{mix,0}^t + (G_{snow,2} - G_1^t)\\Delta t / (C_{v,1} \\times 0.02)$$\n\n> **Known issue (TS_Case03.md)**: The conductance distances in the current implementation are approximately 2× too large. The fix requires changing `solve_imp`'s `z` argument from `dz_snow_s12 (0.04)` to `dz_snow_s1 (0.02)` for the top-layer solve, and correcting the L1→L2 and L2→soil flux distances accordingly.\n\n`T_surf = T_snow0` in this regime.\n\n---\n\n## Internal Soil Heat Conduction — `UpdateHeatFlux`\n\nAfter `G[1]` is set, fluxes between soil layers use harmonic-mean thermal resistance:\n\n$$G_{i+1} = \\frac{2(T_{soil,i}^t - T_{soil,i+1}^t)}{\\Delta z_i/\\kappa_i + \\Delta z_{i+1}/\\kappa_{i+1}}, \\quad i = 1,\\ldots,N-1$$\n\nBottom boundary (connects to annual mean air temperature $T_{ann}$):\n\n$$G_{N+1} = \\kappa_N (T_{soil,N}^t - T_{ann}) / (D_f + \\Delta z_N/2)$$\n\nInternal fluxes are clamped to $[-200, 200]\\,\\text{W m}^{-2}$.\n\nTemperature update (explicit Euler, $S^H = 0$):\n\n$$T_{soil,i}^{t+1} = T_{soil,i}^t + (G_i - G_{i+1})\\Delta t / (C_{v,i} \\Delta z_i)$$\n\nResult clamped to $[-50, 50]$ °C. If `fix_Tsoil=true`, the temperature update is skipped but ice ratio still updates.\n\n---\n\n## Freeze-Thaw State — `Update_ice_ratio`\n\nFreezing (temperature crosses 0°C from above, $f_{ice} < 1$):\n\n$$Q_f = (0 - T^{t+1}) C_{v,i} \\Delta z_i$$\n$$f_{ice,i} \\leftarrow \\min\\!\\left(1,\\; f_{ice,i} + Q_f / (L_f \\rho_w \\theta_i \\Delta z_i)\\right), \\quad T^{t+1} \\leftarrow 0\\text{°C}$$\n\nMelting (temperature crosses 0°C from below, $f_{ice} > 0$):\n\n$$Q_m = T^{t+1} C_{v,i} \\Delta z_i$$\n$$f_{ice,i} \\leftarrow \\max\\!\\left(0,\\; f_{ice,i} - Q_m / (L_f \\rho_w \\theta_i \\Delta z_i)\\right), \\quad T^{t+1} \\leftarrow 0\\text{°C}$$\n\nAfter a moisture update changes $\\theta_i$, the ice fraction is rescaled to keep frozen water mass conserved:\n\n$$f_{ice,i} \\leftarrow \\min\\!\\left(1,\\; f_{ice,i} \\cdot \\theta_i^t / \\theta_i^{t+1}\\right)$$\n\nConstants: $L_f = 3.34\\times10^5\\,\\text{J kg}^{-1}$, $\\rho_w = 1000\\,\\text{kg m}^{-3}$.\n\n---\n\n## Coupling to Other Modules\n\n- **Soil moisture** (`UpdateSoilMoisture`): consumes `Tsoil_c` to compute freeze factor and root water uptake temperature stress; writes updated $\\theta_i$ which feeds back into `Cv`, `κ`, and `f_ice` scaling.\n- **Net radiation / aerodynamics**: `surface_temperature!` receives `Rn_g` (net ground radiation), `Gheat_g` (aerodynamic heat conductance), and ET fluxes to form $G_g$.\n- **Snow density** (`snow_density.jl`): provides `ρ_snow` used in `cal_κ_snow` and the `ΔM_snow` storage term.","other-ch03-surfacefluxes":"# Other — ch03_SurfaceFluxes\n\n# Surface Fluxes: Aerodynamic Conductance\n\nThe `ch03_SurfaceFluxes` module implements the aerodynamic transport pathways between the atmosphere, canopy layers, and soil surface. It provides friction velocity, stability corrections, and layer-wise resistance values that drive sensible heat, latent heat, and CO₂ flux calculations throughout the model.\n\nTwo implementations exist: `aerodynamic_conductance.jl` (V1, legacy) and `aerodynamic_conductance_V2.jl` (V2, current). V2 fixes several physical inconsistencies present in V1 and should be used for all new work.\n\n---\n\n## Physical Framework\n\n### Above-Canopy: Monin-Obukhov Similarity Theory\n\nWind and temperature profiles above the canopy follow MOST. The dimensionless wind gradient satisfies:\n\n```\nk(z-d)/u★ · ∂u/∂z = φₘ((z-d)/L)\n```\n\nIntegrating from the momentum roughness height `z₀ₘ` to observation height `z` yields the log-wind profile with stability correction:\n\n```\nu(z) = u★/k · [ln((z-d)/z₀ₘ) − Ψₘ((z-d)/L)]\n```\n\nFriction velocity is solved directly from this by rearranging for `u★`. Because `u★` depends on `L` (Obukhov length) and `L` depends on `u★`, the standard approach is:\n\n1. Neutral estimate: `u★⁰ = k·u / ln((z-d)/z₀ₘ)` (set Ψₘ = 0)\n2. Compute `L = −ρcₚT·u★³ / (kg·H)`\n3. Recompute `u★` using the full formula with `Ψₘ`\n\nThe aerodynamic resistance above the canopy is derived from the analogous heat transport equation:\n\n```\nra_o = 1/(k·u★) · [ln((z-d)/z₀ₕ) − Ψₕ((z-d)/L)]\n```\n\n**Critical distinction:** `u★` is computed using *momentum* roughness length `z₀ₘ ≈ 0.08h`. The aerodynamic resistance `ra_o` uses *thermal* roughness length `z₀ₕ ≈ 0.1·z₀ₘ`, i.e., `z₀ₕ ≈ 0.008h`. Conflating them is the primary error in V1 (see V1 vs V2 below).\n\n### Within-Canopy: K-Theory\n\nInside the canopy, turbulent diffusivity decays exponentially with depth. The eddy diffusivity profile is:\n\n```\nKₕ(z) = Kₕ(top) · exp(−γ · (1 − z/h))\n```\n\nThe aerodynamic resistance between any two heights `z₁` and `z₂` is:\n\n```\nra(z₁, z₂) = ∫[z₁→z₂] 1/Kₕ(z) dz\n            = h/(γ·Kₕ) · [exp(γ(1−z₁/h)) − exp(γ(1−z₂/h))]\n```\n\nThis integral is applied separately for each canopy layer.\n\n---\n\n## Resistance Network\n\nThe canopy is treated as a three-node resistance network connecting ground, understory air, and overstory air to the reference atmosphere:\n\n```mermaid\ngraph TB\n    ATM[\"Atmosphere (z_wind)\"]\n    AERO[\"ra_o\"]\n    OVR[\"Overstory air node\"]\n    RAU[\"ra_u\"]\n    UND[\"Understory air node\"]\n    RAG[\"ra_g (local)\"]\n    GND[\"Ground\"]\n\n    ATM --> AERO --> OVR --> RAU --> UND --> RAG --> GND\n```\n\nGround and each vegetation layer inject fluxes in *parallel* into their local air node. **The layer resistances `ra_g`, `ra_u`, `ra_o` must not be summed in series as a single pathway** for individual layer fluxes — doing so violates energy conservation. The exception is soil evaporation: `ra_g` passed to `evaporation_soil` is the *total* series resistance from ground to atmosphere (`ra_g_local + ra_u + ra_o`), giving `G_heat,g = 1/ra_g`.\n\n### Layer Definitions\n\n| Symbol | Description |\n|---|---|\n| `ra_o` | Above canopy: reference height → canopy top (MOST) |\n| `ra_u` | Canopy interior: understory top `hᵤ` → overstory top `h`, governed by `γ_{o,h}` |\n| `ra_g` (local) | Near-ground: surface → `hᵤ`, governed by `γ_{g,h}` |\n\nThe continuity condition at `hᵤ` ensures physical consistency:\n\n```\nKₕ,ᵤ = Kₕ,ₒ · exp(−γ_{o,h} · (1 − hᵤ/h))\n```\n\nThe near-ground layer then integrates from this base value using its own extinction coefficient.\n\n---\n\n## Extinction Coefficient Decoupling (V2)\n\nV2 uses four distinct extinction coefficients rather than a single shared `γ`:\n\n| Symbol | Governs |\n|---|---|\n| `γ_{o,m}` | Momentum (wind speed) decay through overstory |\n| `γ_{o,h}` | Heat/scalar turbulence decay through overstory |\n| `γ_{u,m}` | Momentum decay through understory |\n| `γ_{g,h}` | Heat/scalar decay in near-ground layer |\n\nWind speed attenuation (used to compute effective leaf wind speed and boundary layer resistance `rb`) uses momentum coefficients. Turbulent diffusivity profiles and aerodynamic resistances use heat/scalar coefficients.\n\n---\n\n## V1 vs V2: Key Differences\n\n### Monin-Obukhov Length\n\nV1 computes `L` with inverted numerator/denominator, giving units of m⁻¹ instead of m. The error cancels in the stability parameter `ξ = (z-d)/L` because V1 uses multiplication `ξ = (z-d)·L_V1` while V2 uses division `ξ = (z-d)/L_V2`. However the stability correction in `ra_o` then diverges:\n\n| | Formula |\n|---|---|\n| V1 | `ra_o = 1/(k·u★) · [ln((z-d)/z₀ₘ) + n·ξ]` |\n| V2 | `ra_o = 1/(k·u★) · [ln((z-d)/z₀ₕ) − Ψₕ(ξ)]` |\n\nV1 uses additive stability correction (incorrect) and `z₀ₘ` (wrong roughness scale). V2 uses the proper integral form `Ψₕ` and `z₀ₕ`.\n\n### Quantitative Impact\n\nUnder typical conditions (`h = 20 m`, `u = 3 m/s`, `H = 100 W/m²`):\n\n| Variable | V1 | V2 |\n|---|---|---|\n| `ra_o` [s/m] | 6.87 | 22.28 |\n| `Gₕ` [m/s] | 0.0276 | 0.0150 |\n| Canopy Tₛ [°C] | 22.3 | 24.6 |\n| SH [W/m²] | 64 | 70 |\n\nV2 produces systematically higher `ra_o` (by approximately `kB⁻¹/(k·u★) ≈ 5.75/u★` s/m), lower conductance, higher canopy temperatures, and higher sensible heat fluxes. The effect is largest under dry/low-transpiration conditions where canopy temperature is sensitive to `Gₕ`.\n\n---\n\n## Forcing Variable Reference Height\n\nAll meteorological inputs (`Uz`, `Tair`, `RH`) must be measured at the **same reference height** `z_wind`, set via `ps.veg.z_wind`. This is the height at which the log-wind profile is evaluated.\n\n**Consistency requirement:** `ra_o` describes transport from `z_wind` to canopy top. If `Tair` comes from a different height, the driving gradient `Tc − Tair` is physically inconsistent.\n\n### Height Validity\n\nThe log profile is only valid for `z_wind > d + z₀ₘ ≈ 0.9h`. For a 15 m canopy, `d ≈ 12 m` and `z₀ₘ ≈ 1.2 m`, so measurements below ~13.2 m lie inside the canopy and are invalid inputs.\n\n### Common Situations\n\n**Low vegetation (`h < 1.5 m`, `z_wind = 2 m`):** Standard 2 m station data works directly. Set `z_wind = 2`.\n\n**Tall forest + flux tower (`z_wind = 30 m`):** Use tower-top observations for all three variables. 2 m measurements are inside the canopy, where `RH` is elevated (underestimates VPD → lower transpiration) and daytime `Tair` is suppressed by shading (overestimates sensible heat).\n\n**Reanalysis data (ERA5):** Tair is at 2 m, wind at 10 m. Set `z_wind` to the wind speed height (10 m), or extrapolate wind to a consistent height using the log profile, and use co-located temperature.\n\n**Wind measured inside canopy:** If `z_wind < d`, the log profile fails. Two-step extrapolation is required: (1) use the exponential within-canopy profile to back-calculate `u(h)`, then (2) apply the log profile to extrapolate above canopy to a valid reference height.\n\n### Diagnostic Methods for V1/V2 Validation\n\nTwo inversion approaches allow comparison against eddy-covariance observations:\n\n**Method A — Direct resistance inversion (requires LST):**  \nIf canopy temperature and observed sensible heat flux are available:\n```\nr_H = ρₐ·cₚ·(Tc − Tair) / SH_obs\nra_o^EC = r_H − 0.5·rb_o\n```\nCompare `ra_o^EC` against V1/V2 predictions.\n\n**Method B — Penman-Monteith inversion (LE only):**  \nBack-calculate canopy conductance `Gc` from observed latent heat. Using V1 vs V2 `Ga` values, compare inferred `Gc` against modeled stomatal conductance to identify which version yields physiologically plausible values.\n\n---\n\n## Relevant Source Files\n\n| File | Contents |\n|---|---|\n| `src/aerodynamic_conductance.jl` | V1 implementation (legacy) |\n| `src/aerodynamic_conductance_V2.jl` | V2 implementation (current) |\n| `src/DataType/AeroConsts.jl:42` | `aero_exp_terms`, `ustar` and `ra_o` computation |\n| `src/DataType/Params/Params.jl:32` | `z_wind` parameter definition |\n| `src/DataType/Met.jl` | `Met` struct — note: `Tair` comment \"2m\" is a legacy annotation, actual requirement is measurement height |\n| `src/heat_H_and_LE.jl` | Consumes `ra_o`, `rb_o` to compute sensible and latent heat fluxes |\n| `src/evaporation_soil.jl` | Consumes total-column `ra_g` for soil evaporation |","other-ch04-photosynthesis":"# Other — ch04_Photosynthesis\n\n# ch04 — Photosynthesis\n\nBEPS.jl contains two independent Farquhar-Ball-Berry implementations. They share the same FvCB (1980) biochemical framework but differ in temperature responses, kinetic constants, electron transport, and solver strategy. Understanding both is necessary for any work that touches leaf-level carbon fluxes.\n\n---\n\n## Two Implementations\n\n| | `photosynthesis_jl` | `Photosynthesis` (standalone) |\n|---|---|---|\n| **File** | `src/photosynthesis.jl` | `src/standalone/Photosynthesis/` |\n| **Called from** | `inter_prg_jl` (main time-step loop) | Parameter optimization, single-leaf studies |\n| **Leaf temperature** | Iterated via energy balance (`T_l ≠ T_a`) | Assumed `T_l = T_a` |\n| **Leaf surface RH** | Back-calculated from latent heat flux LE | Atmospheric RH passed directly |\n| **Vcmax canopy gradient** | `VCmax()` in `src/SPAC/VCmax.jl` (N gradient) | Uniform Vcmax25 for sun/shade |\n| **Solver** | Analytical cubic/quadratic | Fixed-point iteration (≤15 steps) |\n| **Kinetic constants** | Harley & Baldocchi (1995) | Bernacchi et al. (2001) |\n| **T-response shape** | TBOLTZ bell (T_opt = 301 K ≈ 28°C) | Medlyn normalized Arrhenius (~15°C optimum) |\n\nUse `photosynthesis_jl` for full coupled simulation. Use the standalone module for optimization loops where you need a cheap, decoupled single-leaf model.\n\n---\n\n## Physics Overview\n\n### CO₂ Diffusion Path\n\nCO₂ moves from atmosphere to intercellular space through two resistances in series. All conductances in the solver are CO₂ molar conductances (mol m⁻² s⁻¹); water vapor conductances are divided by 1.6 before use.\n\n```\nc_s = c_a − A / g_b\nc_i = c_a − A (1/g_b + 1/g_s)\nA   = A_g − R_d\n```\n\nDark respiration is suppressed to 40% of its dark-adapted value under illumination (PPFD > 5 μmol m⁻² s⁻¹):\n\n```\nR_d = 0.4 × R_{d,25} × f_T(e_rd)        (light)\nR_d =       R_{d,25} × f_T(e_rd)        (dark)\nR_{d,25} = 0.004657 × Vcmax25\n```\n\n### Farquhar Model\n\nRubisco-limited (`W_c`) and light-limited (`W_j`) assimilation share a unified form:\n\n```\nA_g = a(c_i − Γ) / (e·c_i + b)\n```\n\n| Rate | a | e | b |\n|---|---|---|---|\n| W_c (Rubisco) | Vcmax | 1 | K_c(1 + O₂/K_o) |\n| W_j (light)   | J_x  | 4 | 8Γ |\n\nΓ is the CO₂ compensation point (no dark respiration): `Γ = 0.5·O₂/τ × 1000`.\n\nElectron transport rate (`photosynthesis_jl`):\n```\nJ_x = J_max × PPFD / (PPFD + 2.1·J_max)        # Chen 1999, hyperbolic saturation\nPPFD = 4.55 × 0.5 × R_sn\n```\n\nSucrose-export limitation (third limitation, Collatz):\n```\nA_sucrose = Vcmax/2 − R_d\n```\nFinal assimilation: `A = min(W_c − R_d, W_j − R_d, A_sucrose)`. In practice, `W_c` and `W_j` select the active limitation branch first; `A_sucrose` clips the result afterward.\n\n---\n\n## Temperature Responses\n\n### `TBOLTZ` — `photosynthesis_helper.jl`\n\nUsed for Vcmax and Jmax in `photosynthesis_jl`. Bell-shaped (Harley & Tenhunen 1991), peak at T_opt ≈ 28°C:\n\n```\nf_TBOLTZ(T_l) = rate × H_k·exp(e_a(T_l − T_opt)/(R·T_opt·T_l))\n                      / (H_k − e_a(1 − exp(H_k(T_l − T_opt)/(R·T_opt·T_l))))\n```\n\n`H_k = 200000` J mol⁻¹, `e_vc = e_jm = 55000` J mol⁻¹.\n\nNot normalized at 25°C — the function returns the absolute rate, not a dimensionless factor.\n\n### Medlyn Normalized Arrhenius — `Photosynthesis/temperature.jl:fTv`\n\nUsed in the standalone module. Strictly equals Vcmax25 at 25°C:\n\n```\nf_Medlyn(T_l) = Vcmax25\n              × exp((T_l − T_K25)·e_vc / (T_K25·R·T_l))                 # Arrhenius gain\n              × (1 + exp((S·T_K25 − H_d)/(T_K25·R)))\n                / (1 + exp((S·T_l − H_d)/(T_l·R)))                       # normalized deactivation\n```\n\n`H_d = 200000` J mol⁻¹, `S = 640` J mol⁻¹ K⁻¹, `e_vc = 30000` J mol⁻¹. Optimum ≈ 15°C — more appropriate for boreal conifers. At temperatures above 20°C this gives significantly lower Vcmax than TBOLTZ.\n\n### Arrhenius for K_c, K_o, τ, R_d\n\n```\nf_T(e_act) = exp((T_l − T_K25) × e_act / (T_K25 × R × T_l))\n```\n\nKinetic constants at 25°C differ between implementations (see comparison table in `photosynthesis_comparison.typ` §3). Key impact: the standalone module's `K = K_c(1 + O₂/K_o)` is ~1.8× larger, suppressing W_c at the same c_i.\n\n---\n\n## Ball-Berry Stomatal Conductance\n\n```\ng_{s,w} = g_{0,w} + g_{1,w} × RH_l × β_soil × A / c_s\ng_{s,c} = g_{s,w} / 1.6\n```\n\n`g_0` (minimum conductance) and `g_1` (slope) are vegetation parameters from `ParamVeg`. `β_soil` is the soil water stress factor (0–1) from `soil_water_factor_v2.jl`.\n\nIn `photosynthesis_jl`, leaf-surface RH is derived from the latent heat flux, coupling stomatal conductance back to the energy balance:\n\n```\nρ_v = (LE/λ)·r_v + ρ_{v,a}        # absolute humidity at leaf surface\ne   = ρ_v · T_l / 0.2165\nRH_l = 1 − (e_s − e) / e_s\n```\n\nThe standalone module uses atmospheric RH directly, skipping this coupling entirely.\n\n---\n\n## Coupled Solver\n\nThis is the core algorithmic piece of `photosynthesis_jl`.\n\n```mermaid\nflowchart TD\n    A[Inputs: T_l, R_sn, LE, gb_w, Vcmax25, β_soil, g0, g1, ci_old, ca] --> B[Temperature responses<br/>Kc Ko τ → Γ K Rd Vcmax Jmax]\n    B --> C[Compute Jx from PPFD]\n    C --> D[Estimate Wc Wj using ci_old]\n    D --> E{Wj > Rd AND Wc > Rd?}\n    E -->|yes| F[solve_cubic<br/>Ball-Berry coupled]\n    E -->|no| G[solve_quad<br/>gs = g0]\n    F --> H[An = min An, A_sucrose<br/>An = max An, 0]\n    G --> H\n    H --> I[Back-calculate cs, gs_w, ci]\n```\n\n### Quadratic Branch — `gs = g0`\n\nActivated when either W_j ≤ R_d or W_c ≤ R_d (dim light or very low temperature). With `g_s = g_0` fixed and `D = 1/g_b + 1/g_0`:\n\n```\nα₂·A² + β₂·A + γ₂ = 0\n\nα₂ = −e·D\nβ₂ = e·c_a + b + D·(a − e·R_d)\nγ₂ = R_d·(e·c_a + b) − a·(c_a − Γ)\n\nA = (−β₂ + √(β₂² − 4α₂γ₂)) / (2α₂)\n```\n\nSince `α₂ < 0` (parabola opens downward), this selects the physically meaningful smaller positive root.\n\n### Cubic Branch — Ball-Berry Coupled\n\nThe full coupled system (Farquhar + Ball-Berry + Fick) eliminates `c_i` and `c_s`, yielding a cubic in `A`. Define abbreviations:\n\n```\nM = a·Γ + b·R_d\nN = a − e·R_d\nθ = c·g_b − g_0          (c = g_1·RH_l·β_soil, CO₂ units)\nα = 1 + g_0/g_b − c\n```\n\nCubic coefficients after clearing denominators and collecting by power of A:\n\n```\np = (e·β + b·θ − a·α + e·R_d·α) / (e·α)\nq = (e·γ + b·γ/c_a − a·β + a·Γ·θ + e·R_d·β + R_d·b·θ) / (e·α)\nr = (−a·γ + a·Γ·γ/c_a + e·R_d·γ + R_d·b·γ/c_a) / (e·α)\n```\n\nwhere `α = 1 + g_0/g_b − c`, `β = −c_a·(g_b + 2g_0 − c·g_b)`, `γ = g_0·g_b·c_a²`.\n\nSolved with the Vieta trigonometric substitution (Press, *Numerical Recipes*):\n\n```\nQ   = (p² − 3q) / 9\nU   = (2p³ − 9pq + 27r) / 54\nψ   = arccos(clamp(U/√Q³, −1, 1))        # clamp guards against floating-point overrun\n\nA_k = −2√Q · cos((ψ + 2kπ)/3) − p/3,    k = 0, 1, −1\n```\n\nRoot selection (`findroot`):\n\n| Roots | Rule |\n|---|---|\n| All positive | Smallest positive root |\n| Two negative, one positive | The unique positive root |\n| One negative, two positive | Middle positive root |\n\n---\n\n## Canopy Scaling — `VCmax()` in `src/SPAC/VCmax.jl`\n\n`photosynthesis_jl` receives separate Vcmax values for sunlit and shaded leaves, computed by `VCmax()` using a vertical nitrogen gradient (Chen et al. 2012):\n\n```\nVcmax_sunlit  ∝ Vcmax25 × χ·N·K · (1 − e^{−(Kn+K)L}) / ((Kn+K)(1 − e^{−KL}))\nVcmax_shaded  ∝ Vcmax25 × χ·N · [e^{−Kn·L}/Kn − (1 − e^{−(Kn+K)L})/(Kn+K)]\n                                 / [L − (1 − e^{−KL})/K]\n```\n\n`K = 0.5·Ω/cos(θ_s)`, `K_n = 0.3`. At high LAI, shaded-leaf Vcmax can fall to 20–50% of sunlit. The standalone module skips this entirely — both layers use the same Vcmax25.\n\n---\n\n## Known Bugs Fixed in Standalone Module (2026-04-29)\n\nFive errors were corrected in `src/standalone/Photosynthesis/`:\n\n| # | Location | Error | Fix |\n|---|---|---|---|\n| 1 | `core.jl` — Γ* | `Γ* = 40.0/τ` → ~0.015 μmol mol⁻¹ | `Γ* = 105000/τ` → ~40.4 μmol mol⁻¹ |\n| 2 | `ParamPhoto.jl` — K_o,25 | `K_o,25 = 248000` (×1000 wrong) | `K_o,25 = 248` mmol mol⁻¹ |\n| 3 | `temperature.jl` — fTv | No normalization; Vcmax(25°C) ≈ 3724 | Medlyn normalization; Vcmax(25°C) = Vcmax25 exactly |\n| 4 | `core.jl` — TPU limit | `W_p = 0.167·Vcmax` | `W_p = 0.5·Vcmax` (= 3·TPU, TPU ≈ Vcmax/6) |\n| 5 | `photosynthesis.jl` — default g_b | `g_b = 0.01` mol m⁻² s⁻¹ (50× too small) | `g_b = 0.5` mol m⁻² s⁻¹ |\n\nBug 1 caused Γ* to be ~2700× too small, effectively eliminating the CO₂ compensation point. Bug 2 made the O₂/K_o term ~1000× too small, collapsing the apparent K_m. Both would have suppressed W_c unrealistically.\n\n---\n\n## Open Inconsistencies\n\nTwo residual divergences between the implementations warrant attention before unifying them:\n\n**Ball-Berry c_s simplification.** The standalone module passes `c_s = c_a` into the Ball-Berry equation rather than `c_s = c_a − A/g_b`. This underestimates CO₂ depletion at the leaf surface under high assimilation, causing `g_s` and `c_i` to be slightly overestimated. Fix: `c_s = c_a − An/gb` in the iterative loop.\n\n**Temperature response mismatch.** TBOLTZ peaks at ~28°C; the Medlyn scheme peaks at ~15°C. At 25–35°C (the main operating range for most flux tower sites), TBOLTZ gives 20–60% higher Vcmax. If the standalone module is used for temperate ecosystems, consider switching `e_vc` to 55000 J mol⁻¹ with T_opt = 301 K to match `photosynthesis_jl`.\n\n---\n\n## References\n\n1. Farquhar et al. (1980). *Planta* 149, 78–90.\n2. Ball, Woodrow & Berry (1987). *Progress in Photosynthesis Research* 4, 221–224.\n3. Harley & Baldocchi (1995). *Plant, Cell & Environment* 18(10), 1146–1156.\n4. Chen et al. (1999). *Ecological Modelling* 124, 99–119.\n5. Bernacchi et al. (2001). *Plant, Cell & Environment* 24(2), 253–259.\n6. Medlyn et al. (2002). *Plant, Cell & Environment* 25(9), 1167–1179.","other-claude-md":"# Other — CLAUDE.md\n\n# BEPS.jl — Developer Reference\n\nBEPS (Boreal Ecosystem Productivity Simulator) is a coupled soil-plant-atmosphere (SPAC) model that simulates carbon, water, and energy fluxes at hourly resolution. The package is written in Julia with an optional C execution path via a precompiled `libbeps` shared library.\n\n## Two Execution Paths\n\nEvery simulation runs through one of two code paths, selectable at runtime:\n\n| Path | Entry point | When to use |\n|------|-------------|-------------|\n| Pure Julia | `inter_prg_jl` | Development, debugging, custom modifications |\n| C binding | `inter_prg_c` | Validated legacy runs, bit-for-bit C reproducibility |\n\nBoth paths share the same forcing data and output format. The `version=\"c\"` keyword on `beps_main` or `inter_prg_c` directly switches to the C path.\n\n## Setup and Initialization\n\nAll runs start by constructing a parameter struct and a mutable state struct from ecosystem type strings:\n\n```julia\nstate, ps = setup(\"evergreen_needleleaf\", \"silty_clay_loam\")\n```\n\n`setup` (in `DataType/setup.jl`) reads vegetation and soil lookup tables (`GlobalData.jl`) via `InitParam_Veg`, `InitParam_Soil`, and `Init_Soil_T_θ!`, then assembles:\n\n- `ps::ParamBEPS{FT,N,H,T}` — immutable parameter container holding `hydraulic::HydraulicProfile`, `thermal::ThermalProfile`, and `veg::ParamVeg`\n- `state::StateBEPS` — mutable state holding per-layer soil moisture θ, soil temperature Tsoil, snowpack depth, root fractions, and water stress factors\n\nThe JAX-style convention is that `ps` is never mutated during a run; all evolving quantities live in `state`.\n\n## Simulation Loop\n\n```mermaid\ngraph TD\n    A[simulate] --> B[fill_met! — timestep i]\n    B --> C[inter_prg_jl — one hour]\n    C --> D[netRadiation]\n    C --> E[aerodynamic_conductance]\n    C --> F[surface_temperature / Leaf]\n    C --> G[photosynthesis]\n    C --> H[evaporation_canopy / evaporation_soil]\n    C --> I[rainfall_stage / snowpack]\n    C --> J[UpdateSoilMoisture! / UpdateHeatFlux!]\n    J --> K[Flux + ETFlux written]\n```\n\n`simulate` (in `beps_modern.jl`) drives the loop. For each hourly timestep it calls `fill_met!(met, forcing, i)` to copy the current row of `MetSeries` into the single-step `Met` struct, then dispatches to `inter_prg_jl`. Physics execute in fixed order: radiation partitioning → aerodynamic resistance → leaf/surface temperature iteration → Farquhar photosynthesis → ET → rainfall interception and snowpack → soil moisture and heat diffusion.\n\n### Prescribing State Variables\n\nWhen observations are available, pass them to bypass the internal solver:\n\n```julia\nsimulate(forcing, lai, dates; ps, state,\n  SM_obs=soil_moisture_matrix,   # skips UpdateSoilMoisture!\n  TS_obs=soil_temp_matrix)       # skips Tsoil update\n```\n\nThis is useful for diagnostic runs or model calibration where you want to evaluate flux responses under observed soil conditions.\n\n## Key Types\n\n**`ParamBEPS{FT,N,H,T}`** (`DataType/BEPS_Param.jl`) — top-level parameter container. Type parameters encode the float type (`FT`), number of soil layers (`N`), hydraulic profile type (`H`), and thermal profile type (`T`). Never mutated after `setup`.\n\n**`StateBEPS`** (`DataType/BEPS_State.jl`) — everything that evolves hour-to-hour: soil moisture profile, soil temperatures, snow depth, interception stores, root fractions, and water stress scalars. `SnowLand` is a nested struct within it for snowpack bookkeeping.\n\n**`Met` / `MetSeries`** (`DataType/Met.jl`) — `MetSeries` holds full-run arrays (Tair, SWin, LWin, RH, wind, precip). `fill_met!` projects one row into the scalar `Met` struct for the physics routines to consume.\n\n**`Flux` / `ETFlux`** (`DataType/OUTPUT.jl`) — `Flux` holds hourly totals (GPP, NPP, SH, LH, snow depth, water depth). `ETFlux` disaggregates evapotranspiration into overstory transpiration, understory transpiration, canopy interception evaporation, snow evaporation, and soil evaporation.\n\n**`LeafCache`** (`DataType/LeafCache.jl`) — scratch space for per-leaf temperature and conductance iterates. Allocated once and reused each timestep to avoid allocation pressure in the inner loop.\n\n**`StateSeries` / `CacheSeries`** (`DataType/StateSeries.jl`) — optional time-series recorders. Control which variables are saved via `VARS_STATE` and `VARS_CACHE` keyword arguments to `simulate`.\n\n## APIs\n\n### Modern API (`beps_modern.jl`)\n\n```julia\ndf_flux, df_ET, states, caches = simulate(forcing, lai, dates;\n  ps, state,\n  lon=120.0, lat=40.0,\n  SM_obs=nothing,\n  TS_obs=nothing,\n  VARS_STATE=[:θ, :Tsoil_c, :z_snow],\n  VARS_CACHE=[:Gs_o, :Gs_u])\n```\n\nReturns four values: a `DataFrame` of hourly fluxes, a `DataFrame` of ET partitioning, a `StateSeries`, and a `CacheSeries`. The preferred entry point for new work.\n\n### Legacy API (`beps_main.jl`)\n\n```julia\ndf_flux, df_ET, states = beps_main(forcing, lai, dates;\n  lon=120.0, lat=20.0,\n  VegType=\"default\", SoilType=\"silty_clay_loam\",\n  version=\"julia\")   # \"c\" invokes libbeps\n```\n\nWraps `setup` internally. Use `version=\"c\"` to invoke the C library path for validation against the original Fortran/C implementation.\n\n## Subsystem Modules\n\n| Subsystem | Files | Responsibility |\n|-----------|-------|----------------|\n| **SPAC helpers** | `SPAC/` | `fill_met!`, LAI partitioning, Vcmax from leaf N, snow density, vegetation/soil type lookup |\n| **Soil physics** | `SoilPhysics/` | `UpdateSoilMoisture!` (delegates to `SoilDiffEqs.jl`), `UpdateHeatFlux!`, water stress factor |\n| **Canopy** | `photosynthesis.jl`, `photosynthesis_helper.jl` | Farquhar A–Ci, Vcmax/Jmax temperature response |\n| **Aerodynamics** | `aerodynamic_conductance.jl`, `_V2.jl` | Aerodynamic resistance; V2 is the current default |\n| **Surface energy** | `heat_H_and_LE.jl`, `surface_temperature.jl` | Sensible/latent heat, leaf temperature iteration |\n| **Hydrology** | `evaporation_soil.jl`, `evaporation_canopy.jl`, `rainfall_stage.jl`, `snowpack.jl` | ET components, interception, snow accumulation/melt |\n| **C bindings** | `clang/` | `libbeps` loading, `Soil_c` struct, C snowpack wrappers |\n| **Data types** | `DataType/` | All structs, `@bounds` macro, `setup`, lookup tables |\n| **Standalone** | `standalone/Photosynthesis/` | Farquhar model without canopy coupling — useful for leaf-scale testing |\n\n## C Library Integration\n\n`libbeps` is loaded at module init from `deps/` or via `LazyArtifacts`. The `Soil_c` struct in `clang/SOIL_c.jl` mirrors the C memory layout exactly; do not reorder its fields. C wrappers are in `clang/BEPS_c.jl`. When adding new C-path functionality, match the field ordering to the C header and add a Julia-side wrapper that converts between `StateBEPS` and `Soil_c`.\n\n## Testing\n\n```bash\njulia --project test/runtests.jl            # full suite\njulia --project test/test-beps_modern.jl   # modern API end-to-end\njulia --project test/test-beps_main.jl     # C vs Julia consistency\njulia --project test/test-photosynthesis_standalone.jl\njulia --project test/test-soil_sm.jl       # soil moisture integration\n```\n\nThe C/Julia consistency test (`test-beps_main.jl`) is the regression guard for any changes touching `inter_prg_jl` or soil physics — if it diverges, the Julia path has drifted from the validated C reference.\n\n## Adding a New Physics Module\n\n1. Create `src/your_module.jl` and add it to `BEPS_modules.jl`.\n2. Accept `met::Met`, `ps::ParamBEPS`, `state::StateBEPS` as arguments — never accept raw arrays.\n3. Write outputs into `flux::Flux` or `etflux::ETFlux`; do not return new structs from inner physics functions.\n4. Wire the call into `inter_prg_jl` at the appropriate point in the radiation → aerodynamics → photosynthesis → hydrology → soil sequence.\n5. Run `gitnexus_impact` on any symbol you modify before committing, then `gitnexus_detect_changes()` before the final commit.","other-codecov-yml":"# Other — codecov.yml\n\n# `codecov.yml` — Coverage Reporting Configuration\n\nConfigures [Codecov](https://codecov.io) to exclude non-production source paths from coverage analysis. Without this file, Codecov would include debug utilities, deprecated code, and C bindings in coverage metrics, distorting the picture of how well the core Julia model is tested.\n\n## Ignored Paths\n\n| Pattern | Reason |\n|---|---|\n| `src/debug/*`, `*/debug/*` | Debugging utilities not part of the production model |\n| `src/*/deprecated`, `src/*/*/deprecated` | Legacy code preserved for reference but not actively maintained |\n| `src/DataType/Params/deprecated/*` | Deprecated parameter definitions |\n| `src/*/temp`, `src/*/backup/*` | Scratch and backup files |\n| `src/clang/*` | C library bindings — coverage is meaningless here since the logic lives in `libbeps`, not Julia |\n| `src/standalone/*` | Standalone sub-models (e.g., the Farquhar photosynthesis module) tested independently |\n\n## What Gets Measured\n\nEverything under `src/` that is not matched by the above patterns — the core simulation modules:\n\n- `inter_prg.jl` and the physics modules it orchestrates (radiation, aerodynamics, photosynthesis, ET, snow, soil)\n- `DataType/` structs and setup helpers\n- `SoilPhysics/` and `SPAC/` subsystems\n- The `simulate` and `beps_main` public APIs\n\n## Maintenance Notes\n\nWhen adding new top-level directories under `src/`, check whether they should be excluded. The convention is:\n\n- **Exclude** anything in a `deprecated/`, `temp/`, `backup/`, or `debug/` subdirectory.\n- **Include** all physics, data type, and solver code that the test suite exercises.\n\nThe `src/standalone/` exclusion is intentional: standalone modules have their own test files (`test-photosynthesis_standalone.jl`) and are typically benchmarked separately from the coupled model.","other-debug":"# Other — debug\n\n# `test/debug` — Longwave Radiation Debug Module\n\nThis directory contains development-time scratch code for validating the multi-layer longwave radiation (Rnl) scheme. It is **not part of the production model** — these files exist to prototype and cross-check the physics before the formulations are integrated into `src/`.\n\n## Purpose\n\nThe core problem being debugged is the net longwave radiation balance across a three-layer canopy stack:\n\n```\n  L_sky (atmosphere)\n     ↕\n  Overstory (o)   — LAI_o, ϵ_o, τ_o\n     ↕\n  Understory (u)  — LAI_u, ϵ_u, τ_u\n     ↕\n  Ground (g)      — ϵ_g\n```\n\nEach layer emits blackbody radiation, transmits a fraction downward/upward, reflects a fraction back, and absorbs the rest. The net radiation at each layer (Rnl_o, Rnl_u, Rnl_g) must balance: `Rnl_o + Rnl_u + Rnl_g ≈ L_sky − L_up_total`.\n\nThe two `debug_Rln` implementations explore two different formulations of this balance and are compared against each other.\n\n## Files\n\n### `debug_Rln.jl` — Direct flux decomposition\n\n`debug_Rln` computes net longwave radiation by explicitly tracking upward and downward flux streams and their inter-layer reflections:\n\n```julia\nRl_up_u   = L_u * (1.0 - τ_u) + L_g * τ_u        # upward from understory\nRl_down_o = L_sky * τ_o + L_o * (1.0 - τ_o)       # downward through overstory\nRl_down_u = Rl_down_o * τ_u + L_u * (1.0 - τ_u)   # downward through understory\n```\n\nReflection terms are then computed explicitly:\n\n```julia\nLf_u2o = Rl_down_o * (1.0 - τ_u) * ρ_u  # understory reflects upward to overstory\nLf_o2u = Rl_up_u   * (1.0 - τ_o) * ρ_o  # overstory reflects downward to understory\n```\n\nNet radiation at each layer absorbs incoming fluxes minus twice the layer's own emission (the factor of 2 accounts for both upward and downward faces of a canopy layer):\n\n```julia\nRnl_o = (ϵ_o * (L_sky + Rl_up_u + Lf_u2o) - 2 * L_o) * (1.0 - τ_o)\nRnl_u = (ϵ_u * (Rl_down_o + L_g + Lf_o2u + Lf_g2u) - 2 * L_u) * (1.0 - τ_u)\nRnl_g = (ϵ_g * (Rl_down_u + Lf_u2g) - L_g) * (1.0 - τ_u)\n```\n\nThe helper `cal_Rln_net(L_up, L_down, Tair, ϵ, τ)` captures the common pattern: `((L_up + L_down) * ϵ - 2 * Rl) * (1 - τ)`.\n\n### `debug_Rln2.jl` — Two-stream propagation formulation\n\nThis file reformulates the same physics using explicit two-stream propagators, following the structure of Bonan (2019) Eq. 14.122 more closely.\n\n`cal_Rln_down(Rdown_i2, Rup_i1, Tair, ϵ, τd)` computes the downward flux emerging from a layer given the incoming flux from above and the upward flux from below:\n\n```julia\nRdown_i2 * (τd + (1 - τd) * τl) +   # transmitted + leaf-transmitted\n  Rup_i1 * (1 - τd) * ρ +            # back-reflected upward flux\n  Rl_i1  * (1 - τd)                  # layer self-emission downward\n```\n\n`cal_Rln_up` mirrors this for the upward direction. The asymmetry between the two functions (`τl = 0` for down, `τl = 1 - ϵ` for up) reflects the directional treatment of leaf transmittance.\n\n`debug_Rln` in this file then chains these propagators top-down then bottom-up and reports both the propagator-based net radiation and a closed-form algebraic alternative (`Rnl_o`, `Rnl_u`, `Rnl_g` via the expanded formulas at the bottom of the function), so the two approaches can be compared at runtime with `@show`.\n\n### `debug-Rn.qmd` — Quarto scratch notebook\n\nA short derivation sketch for `cal_Rln_out`, exploring how a single layer's outgoing longwave is composed from transmission, reflection, and emission terms. This is the conceptual precursor to both `.jl` files.\n\n### `debug.md` — Annotation of production-candidate formulas\n\nContains a candidate production implementation with inline commentary explaining what each term represents physically, including known issues flagged with `# miss ϵ_u` and `# miss a ϵ_g`. These are deliberate markers identifying terms that were missing from an earlier version of the net radiation code.\n\n### `debug_snow.jl` — C vs. Julia consistency test\n\nUnrelated to the radiation work. Loads observed forcing, cools air temperature by 5 K to stress snowpack processes, runs both `beps_main(... version=\"julia\")` and `beps_main(... version=\"c\")`, and computes the percentage difference. The threshold comment `# 0.09%` records the current maximum observed divergence between the two backends.\n\n## Key Helpers\n\n| Function | File | What it computes |\n|---|---|---|\n| `cal_Rln(ϵ, T)` | `src/` (called externally) | Stefan-Boltzmann blackbody emission: `ϵ σ T⁴` |\n| `cal_Rln_net` | `debug_Rln.jl`, `debug_Rln2.jl` | Net Rnl for one layer given up/down fluxes |\n| `cal_Rln_out` | `debug_Rln.jl` | Outgoing flux from a layer (transmission + reflection + emission) |\n| `cal_Rln_down` / `cal_Rln_up` | `debug_Rln2.jl` | Two-stream propagators (downward / upward) |\n| `debug_Rln` | both `.jl` files | Top-level driver; returns `(Rnl_o, Rnl_u, Rnl_g)` |\n\n## Connection to Production Code\n\n`cal_Rln` is imported from `src/` (via `pow` in `src/SPAC/helper.jl`). The formulas developed here feed into `src/netRadiation.jl`, where the production longwave scheme lives. The commented-out blocks in `debug_Rln.jl` and the `# miss` annotations in `debug.md` document the iterative process of correcting the reflection and transmission terms before they were committed to the main source.\n\n`debug_snow.jl` uses `beps_main` and `path_proj` directly from the installed package, making it a lightweight regression harness rather than a unit test — it lives here rather than in `test/` proper because it requires the serialized forcing file `data/p1_meteo`.\n\n## Usage\n\nThese files are not included in `test/runtests.jl`. To run interactively:\n\n```julia\njulia --project\ninclude(\"test/debug/debug_Rln.jl\")\ndebug_Rln()  # returns (Rnl_o, Rnl_u, Rnl_g) for default parameters\n```\n\n```julia\n# Cross-check both formulations at identical inputs\ninclude(\"test/debug/debug_Rln2.jl\")\ndebug_Rln()  # prints both propagator result and algebraic result via @show\n```","other-deps":"# Other — deps\n\n# `deps/` — Native Library Artifacts\n\nThis directory manages the `libbeps` shared library that backs BEPS.jl's C execution path (`inter_prg_c`). It handles two concerns: building the library from source on Windows, and generating cross-platform `Artifacts.toml` entries so Julia's artifact system can fetch pre-built binaries at install time.\n\n## Files\n\n### `build.ps1` — Local Windows Build\n\n```powershell\ncd BEPS.c\nmake\ncp libbeps.dll ..\ncd ..\n```\n\nA minimal PowerShell build script for developers who need to compile `libbeps` locally on Windows. It expects the [CUG-hydro/BEPS.c](https://github.com/CUG-hydro/BEPS.c) source tree to be checked out as a sibling directory named `BEPS.c`, then copies the resulting DLL up one level into `deps/` where BEPS.jl's module loader can find it.\n\nUse this when you need a debug build, are developing against a local fork of BEPS.c, or the pre-built artifacts haven't been published yet for a new release.\n\n### `gen_artifacts.jl` — Artifact Registry Generator\n\nRegenerates `Artifacts.toml` for a new release of `libbeps`. Run it once when a new tagged release is published to GitHub:\n\n```bash\njulia deps/gen_artifacts.jl v2026.05.03\n# or, accepting the default tag:\njulia deps/gen_artifacts.jl\n```\n\nThe script iterates over four supported platforms:\n\n| Platform string | OS | Arch | Library filename |\n|---|---|---|---|\n| `windows-x86_64` | Windows | x86_64 | `libbeps-windows-x86_64.dll` |\n| `linux-x86_64` | Linux | x86_64 | `libbeps-linux-x86_64.so` |\n| `macos-x86_64` | macOS | x86_64 | `libbeps-macos-x86_64.dylib` |\n| `macos-arm64` | macOS | aarch64 | `libbeps-macos-arm64.dylib` |\n\nFor each platform it:\n1. Downloads `beps-{platform}.tar.gz` from the GitHub release\n2. Computes the SHA-256 of the tarball\n3. Unpacks into a Julia artifact tree and records its content hash\n4. Calls `bind_artifact!` with `lazy = true`, writing a platform-tagged entry into `Artifacts.toml`\n\nThe `lazy = true` flag means Julia will only download the artifact for the user's actual platform on first use — not all four at install time.\n\n## Relationship to the Rest of the Package\n\n`libbeps` is loaded at module init in `src/clang/BEPS_c.jl`. The C path is selected at runtime by passing `version=\"c\"` to `beps_main`, or by calling `inter_prg_c` directly. The `Artifacts.toml` generated here is what `LazyArtifacts` resolves to find the correct shared library for the current platform.\n\n```mermaid\nflowchart LR\n    A[\"Artifacts.toml\"] -->|resolved by LazyArtifacts| B[\"src/clang/BEPS_c.jl\"]\n    B --> C[\"inter_prg_c\"]\n    D[\"gen_artifacts.jl\"] -->|writes| A\n    E[\"build.ps1\"] -->|local override| F[\"deps/libbeps.dll\"]\n    F -->|fallback path| B\n```\n\n## Adding a New Platform\n\n1. Ensure the BEPS.c CI produces a `beps-{platform}.tar.gz` archive in the GitHub release containing a single shared library file.\n2. Add a matching named tuple to the `PLATFORMS` constant in `gen_artifacts.jl`.\n3. Re-run the script with the target tag to regenerate `Artifacts.toml`.\n4. Commit `Artifacts.toml` — this is the file users and CI actually consume.\n\n## Adding a New Release\n\n```bash\n# After tagging vYYYY.MM.DD in CUG-hydro/BEPS.c:\njulia deps/gen_artifacts.jl vYYYY.MM.DD\ngit add Artifacts.toml\ngit commit -m \"bump libbeps to vYYYY.MM.DD\"\n```","other-docs":"# Other — docs\n\n# BEPS.jl Documentation (`docs/`)\n\nThe `docs/` directory contains developer-oriented reference material for the BEPS.jl parameter management system, including a quickstart README and an interactive Jupyter notebook demonstrating struct construction and initialization.\n\n---\n\n## `docs/README.md` — Parameter Management Overview\n\nThis document is the entry point for understanding how BEPS.jl handles model parameters. It covers three concerns:\n\n**Parameter definition syntax.** Parameters are declared in structs using a `field::Type = default | (min, max)` pattern via the `@bounds` macro (defined in `src/DataType/macro.jl`). The bounds are used for optimization and validation, not just documentation:\n\n```julia\nLAI_max_o::FT = 4.5 | (0.1, 7.0)\n```\n\n**Default value sources.** Defaults flow from three places depending on the parameter type:\n\n| Source | File | Parameter Type |\n|--------|------|---------------|\n| Hardcoded constants | `src/DataType/Params/Params.jl` | Universal physical constants |\n| JSON lookup | `ParamVeg.json` | Vegetation-type-specific values |\n| Data table | `src/DataType/Params/GlobalData.jl` | Soil-type-specific values |\n\n**Core operations.** The three operations developers use most often:\n\n```julia\nmodel = ParamBEPS(VegType=6, SoilType=4; N=5)  # construct from lookup tables\ndf    = parameters(model)                        # extract as DataFrame\nupdate!(model, [:veg, :VCmax25], 100.0)         # update a nested field\n```\n\nThe `parameters()` / `update!()` interface comes from `ModelParams.jl`. The path `[:veg, :VCmax25]` navigates the nested struct hierarchy: `model.veg.VCmax25`.\n\n---\n\n## `docs/Structs.ipynb` — Interactive Struct Reference\n\nThis notebook demonstrates constructing and inspecting all major BEPS.jl structs. It serves as a runnable field reference — every field name, type, and default value is visible in the cell outputs.\n\n### What it demonstrates\n\n**1. `ParamBEPS` construction**\n\n```julia\nmodel = ParamBEPS(VegType, SoilType)\n```\n\n`ParamBEPS{Float64}` is the top-level parameter container. Its display shows three nested sub-structs:\n\n- `hydraulic::ParamSoilHydraulicLayers{Float64, N}` — per-layer Campbell retention parameters (`θ_vfc`, `θ_sat`, `K_sat`, `ψ_sat`, `b`)\n- `thermal::ParamSoilThermalLayers{Float64, N}` — per-layer thermal properties (`κ_dry`, `ρ_soil`, `V_SOM`)\n- `veg::ParamVeg{Float64}` — vegetation parameters (`LAI_max_o`, `Ω`, `VCmax25`, stomatal slope `g1_w`, etc.)\n\nAll layer arrays have length `N` (default 5), matching the soil discretization depth.\n\n**2. `StateBEPS` initialization via `setup`**\n\n```julia\nstate, model = setup(model; Ta, Tsoil, θ0, z_snow)\n```\n\n`StateBEPS` holds the mutable model state updated at every timestep. Key fields visible in the output:\n\n| Field | Description |\n|-------|-------------|\n| `θ[1:N]` | Volumetric soil moisture per layer |\n| `Tsoil_c[1:N]`, `Tsoil_p[1:N]` | Current/previous soil temperature |\n| `f_root[1:N]` | Root fraction per layer (exponential decay) |\n| `z_snow`, `z_water` | Snow and surface water depth |\n| `ψ[1:N]`, `Kmid[1:N]` | Matric potential, hydraulic conductivity (computed each step) |\n| `f_stress[1:N]` | Water stress factor per layer |\n\nArrays are allocated with length 10 (double the layer count) to allow future layer refinement; only indices `1:N` are active.\n\n**3. Per-timestep structs**\n\nThe notebook also constructs and displays the structs that are populated at each hourly timestep:\n\n- **`Met`** — meteorological forcing: `Rs`, `Rln_in`, `Tair`, `RH`, `Prcp`, `Uz`\n- **`Flux`** — carbon and energy outputs: `GPP`, `NPP`, `NEP`, `SH`, `LH`, `Trans`, `Evap`, snow/water depth\n- **`ETFlux`** — detailed ET partitioning: overstory/understory transpiration, interception evaporation (`Eil`, `EiS`), soil evaporation, sensible/latent heat by component\n- **`Radiation`** — shortwave partitioning into diffuse/direct, overstory/understory/ground components\n- **`LeafCache`** — per-leaf intermediate cache holding per-timestep conductances (`Gs_old/new`, `Gc`, `Gh`, `Gw`), temperatures, assimilation rates, and photosynthesis/aerodynamic constants\n\n---\n\n## Relationship to Source Code\n\n```\ndocs/README.md                    →  src/DataType/macro.jl          (@bounds syntax)\n                                  →  src/DataType/Params/Params.jl  (ParamVeg defaults)\n                                  →  src/DataType/Params/GlobalData.jl (soil tables)\n\ndocs/Structs.ipynb                →  src/DataType/BEPS_Param.jl     (ParamBEPS)\n                                  →  src/DataType/BEPS_State.jl     (StateBEPS)\n                                  →  src/DataType/setup.jl          (setup())\n                                  →  src/DataType/Met.jl            (Met)\n                                  →  src/DataType/OUTPUT.jl         (Flux, ETFlux)\n                                  →  src/DataType/CanopyLayer.jl    (Radiation)\n                                  →  src/DataType/LeafCache.jl      (LeafCache)\n```\n\nWhen `ParamBEPS` or `StateBEPS` field layouts change, the notebook output will need to be re-executed to stay current. The notebook is not tested in CI; treat its cell outputs as illustrative snapshots rather than guaranteed-current values.","other-manual":"# Other — manual\n\n# `docs/manual/` — BEPS.jl User Manual\n\nThe `docs/manual/` directory contains the user-facing documentation for the BEPS.jl model, authored in [Typst](https://typst.app/). It covers model inputs, physical algorithms, parameter tables, and a quick-start guide. A separate subdirectory (`软著/`) holds a software copyright registration document that auto-embeds source files from `src/`.\n\n## Document Map\n\n| File | Purpose |\n|---|---|\n| `manual.typ` | Primary user manual: overview, physics, API reference, quick-start |\n| `ch01_Input.typ` | Input data format and meteorological forcing variables |\n| `ch02_Param_SoilVeg_Types.typ` | Soil texture and vegetation type lookup tables |\n| `软著/ModernBEPS.typ` | Software copyright registration — reads and embeds `src/` files verbatim |\n\n## Toolchain\n\nAll documents use the `modern-cug-report` Typst package:\n\n```typst\n#import \"@local/modern-cug-report:0.1.3\": *\n#show: doc => template(doc, footer: \"BEPS.jl使用手册\", header: \"\")\n```\n\nThis is a **local** Typst package (not on the Typst Universe registry). It must be installed in the user's local Typst package cache before compilation. `codelst:2.0.2` is also used in `manual.typ` for styled code listings.\n\nBuild with:\n```bash\ntypst compile docs/manual/manual.typ\ntypst compile docs/manual/软著/ModernBEPS.typ\n```\n\n## `manual.typ` — Primary Manual\n\nThe main document is self-contained and covers six sections:\n\n1. **Model overview** — dual-layer canopy, sunlit/shaded leaf split, soil layers, module list\n2. **Data variables** — met forcing schema, state variables, canopy stratification\n3. **Model parameters** — vegetation (LAI, albedo, Vcmax25, g1, Ω) and soil (θ_sat, K_sat, b, κ)\n4. **Core physics** — full mathematical derivations for all sub-models (see below)\n5. **Runtime flow** — hourly time-step orchestration sequence\n6. **Usage guide** — `beps_main` call signature, key arguments, output structure\n\n### Physics Covered\n\nThe physics sections document the actual equations implemented in `src/`:\n\n| Section | Corresponding source file |\n|---|---|\n| Solar geometry, beam/diffuse split | `netRadiation.jl` |\n| Canopy radiation (Beer-Lambert + clumping) | `netRadiation.jl` |\n| Farquhar photosynthesis (Wc, Wj, cubic solve) | `photosynthesis.jl` |\n| Ball-Berry stomatal conductance | `photosynthesis.jl` |\n| Leaf energy balance iteration | `surface_temperature.jl` |\n| Penman-Monteith transpiration | `heat_H_and_LE.jl` |\n| Richards equation (Campbell model, implicit FD) | `SoilPhysics/UpdateSoilMoisture.jl` |\n| Soil heat diffusion (Johansen conductivity) | `SoilPhysics/UpdateHeatFlux.jl` |\n| Root water uptake + stress factor β_soil | `SoilPhysics/soil_water_factor_v2.jl` |\n| Snowpack (density evolution, melt energy) | `snowpack.jl` |\n\n### Usage API Documented\n\nThe manual documents `beps_main` (legacy API), not `simulate` (modern API). Key parameters described:\n\n```julia\nresults, et, t_soil, theta_soil = besp_main(\n  drivers, lai,\n  lon=120.0, lat=30.0,\n  VegType=4,\n  SoilType=8\n)\n```\n\n> **Note:** The manual uses `besp_main` (typo) and `df_out` / `Tsoil` / `θ` as output names. The actual modern API (`simulate`) returns `(df_flux, df_ET, states::StateSeries, caches::CacheSeries)`. When updating the manual, reconcile against `src/beps_main.jl` and `src/beps_modern.jl`.\n\n## `ch01_Input.typ` and `ch02_Param_SoilVeg_Types.typ`\n\nThese are chapter fragments, likely intended to be `#include`d into a larger document (they have standalone `#show` rules but are not referenced from `manual.typ`). They document:\n\n- **ch01**: The `par` named-tuple passed to the model (lon, lat, clumping, landcover, soil_type, initial Tsoil/soilwater/snowdepth) and the 7 forcing variables (doy, hour, rad, tem, hum, pre, wind)\n- **ch02**: Cross-reference tables mapping OpenLandMap soil texture IDs → BEPS `soil_type` integers, and IGBP land cover classes → BEPS `landcover` codes\n\nThe soil table in `ch02` is the human-readable version of the lookup data in `src/DataType/Params/GlobalData.jl`.\n\n## `软著/ModernBEPS.typ` — Source Embedding\n\nThis document serves a specific legal purpose (Chinese software copyright registration) and has an unusual structure: it uses a `disp()` helper that reads `src/` files at compile time and renders them as raw code blocks:\n\n```typst\n#let disp(f) = {\n  let dir-root = \"../../../src/\"\n  let file = dir-root + f\n  let evap = read(file)\n  raw(evap, lang: \"julia\")\n  v(1em)\n}\n```\n\nThe path `../../../src/` is relative to the `.typ` file's location (`docs/manual/软著/`), resolving to the repo root `src/`. This document must be compiled from within or with the repo root accessible. It currently embeds:\n\n- Main entry: `BEPS.jl`, `beps_main.jl`, `beps_modern.jl`\n- Physics modules: aerodynamics, evaporation, heat, radiation, photosynthesis, rainfall, snowpack, surface temperature\n- Data types: all files under `DataType/`\n- Parameters: all files under `DataType/Params/`\n\nWhen source files are added, renamed, or removed, update the `disp(...)` call list in this document to keep the registration current.\n\n## Maintenance Notes\n\n**When model physics change:** Update the equation in `manual.typ` section 4.x to match. The manual's equations are transcribed — they do not auto-generate from code.\n\n**When parameter tables change:** `src/DataType/Params/GlobalData.jl` is the source of truth. The tables in `ch02_Param_SoilVeg_Types.typ` and `manual.typ` §3 must be updated manually to match.\n\n**When the API changes:** Section 6 of `manual.typ` documents `beps_main` arguments and output field names. Cross-check against `src/beps_main.jl:beps_main` and `src/beps_modern.jl:simulate` after any signature change.\n\n**The manual carries a caveat:** `manual.typ` includes a prominent warning block that the content was auto-generated and should be used carefully (`*目前这个手册是：Claude帮忙自动生成的，需谨慎使用。*`). Remove or update this notice once the manual has been human-reviewed.","other-modules":"# Other — modules\n\n# Module-Level Test Suite (`test/modules/`)\n\nThis directory contains unit tests for individual BEPS physics modules. Each test file isolates one module and validates the Julia implementation against the C reference (`clang.*` wrappers from `libbeps`), establishing numerical equivalence at a tolerance of 1e-7 to 1e-12 depending on the module's sensitivity.\n\n## Purpose\n\nThe module tests serve two distinct goals:\n\n1. **C/Julia parity** — verify that the pure-Julia physics routines produce bit-equivalent (or near-equivalent) results to the original C library.\n2. **API compatibility** — verify that the newer `StateBEPS + ParamBEPS` struct pair produces the same results as the legacy `Soil` struct for every soil physics operation.\n\nTests in this suite are physics-level, not integration-level. They do not run `inter_prg_jl` or `simulate`; each test exercises a single function with controlled inputs.\n\n## Entry Point\n\n`modules.jl` loads `BEPS` and includes each sub-test file in order:\n\n```julia\ninclude(\"test-rainfall_stage1.jl\")\n# include(\"test-surface_temperature.jl\")  # disabled\ninclude(\"test-radiation.jl\")\ninclude(\"test-snowpack.jl\")\ninclude(\"test-photosynthesis.jl\")\ninclude(\"test-sensible_heat.jl\")\ninclude(\"test-aerodynamic_conductance.jl\")\ninclude(\"test-param.jl\")\ninclude(\"test-Soil.jl\")\ninclude(\"test-setup.jl\")\n```\n\n`test-surface_temperature.jl` is included but its `@testset` is commented out — the test logic exists but is not registered with the test harness. This is intentional; the file still documents the expected call signature for `surface_temperature_jl`.\n\n## Test Files\n\n### `test-param.jl` — Vegetation parameter loading\n\nTests `ReadParamVeg` for a representative set of land-cover codes. The C version is called via `clang.ReadParamVeg`. Known deviation: LAI_max_understory differs for `lc = -1` (Julia changed the default from 0.01 to 2.4), so that case is excluded. Tolerance is `atol=1e-3` rather than the usual tight tolerance, due to accumulated floating-point differences in `slope_Vc`.\n\n### `test-radiation.jl` — Net radiation partitioning\n\nTests `netRadiation_jl` against `clang.netRadiation_c`. Two physically distinct cases are covered: a mid-day high-radiation case and a lower-sun-angle case with non-zero snow fractions. The helper `test_net_radiation_case` constructs fresh `Leaf` and `Radiation` structs for both implementations and compares both the scalar return value and the four leaf-layer components (`o_sunlit`, `o_shaded`, `u_sunlit`, `u_shaded`). Tolerance is `atol=1e-10, rtol=1e-12` — the tightest in the suite.\n\n### `test-snowpack.jl` — Snowpack stages 1 and 3\n\n`snowpack_stage1` is tested across three temperature regimes (below zero, above zero, near-zero) because the near-zero path exercises floating-point edge cases noted in the source. Output from `call_jl()` and `call_c()` must be exactly equal (`==`), not just approximately equal.\n\n`snowpack_stage3` tests both melt (Tsnow goes positive) and refreezing (Tsnow goes negative) across multiple snow depths. A known bug in the C refreezing path is documented directly in the test:\n\n```julia\n#! C结冻存在错误\n@test r_jl == (0.06000080057281437, 0.019999453267346284)\n```\n\nThis is one of the few places where the Julia implementation is *not* expected to match C — the Julia result is the authoritative correct value.\n\n### `test-photosynthesis.jl` — Farquhar photosynthesis\n\nTests `photosynthesis_jl` vs `clang.photosynthesis_c` for both daytime (`Rs_leaf > 0`) and nighttime (`Rs_leaf = 0`) conditions. Uses fixed realistic inputs derived from an actual model run. Tolerance: `rtol=1e-7` (daytime), `rtol=1e-8` (nighttime).\n\n### `test-sensible_heat.jl` — Sensible heat fluxes\n\nExercises `sensible_heat_jl` vs `clang.sensible_heat_c` with a `Leaf` struct for multi-layer temperatures and conductances. Uses `all(r1 .≈ r2)` rather than `@test`, which means a failure silently returns `false` — this is a minor testing gap.\n\n### `test-aerodynamic_conductance.jl` — Aerodynamic resistance\n\nTests `aerodynamic_conductance_jl` vs `clang.aerodynamic_conductance_c`. Returns a 7-tuple `(ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)` and checks the maximum absolute difference across all outputs, tolerance `1e-8`.\n\n### `test-rainfall_stage1.jl` — Rainfall interception\n\nTests `rainfall_stage1_jl` vs `clang.rainfall_stage1` using `Layer2` inputs for canopy water state. Checks exact equality (`==`) on the ground precipitation scalar.\n\n### `test-Soil.jl` — Soil physics and API compatibility\n\nThe most comprehensive test file. It has two distinct sections:\n\n**C/Julia parity tests** (`Init_Soil_Parameters`, `UpdateHeatFlux`, `Init_Soil_var`): The helper `is_soil_equal` iterates over all `Soil` fields by position and applies unit conversions before comparing:\n\n| Julia unit | C unit | Conversion applied |\n|---|---|---|\n| cm/h | m/s | ÷ 360000 |\n| cm H₂O (negative) | m H₂O (positive) | ÷ -100 |\n| cm H₂O (positive) | m H₂O | ÷ 100 |\n\nSoil types 6 and 7 are skipped because they contain known divergences (`κ_dry` and `K_sat` respectively).\n\n**API compatibility tests** (`StateBEPS + ParamBEPS API 兼容性`, `Soil → StateBEPS 转换`, `Soil2Params! 逆函数测试`): These tests verify that the migration from the legacy `Soil` struct to the `StateBEPS + ParamBEPS` pair preserves numerical results. The full chain is tested:\n\n```\nParams2Soil!(soil, ps)           # params → legacy struct\nStateBEPS(soil)                  # legacy struct → new state\nState2Soil!(soil, st)            # write back: new state → legacy struct\nSoil2Params!(ps, soil)           # legacy struct → params (inverse of Params2Soil!)\n```\n\nFields verified include: `θ`, `Tsoil_c`, `Tsoil_p`, `ice_ratio`, `f_root`, `f_soilwater`, `f_stress`, `w_norm`, `κ`, `Cv`, `G`, `z_water`, `r_waterflow`, `ψ`, `Kmid`, `Kavg`, `ETi`. Post-moisture-update invariants are also checked: `ψ[1:5] .<= 0` and `Kmid[1:5] .>= 0`.\n\n### `test-setup.jl` — Model initialization\n\nTests all three `setup` overloads:\n\n| Signature | Returns | What it tests |\n|---|---|---|\n| `setup(VegType::Int, SoilType::Int; Ta, θ0)` | `(StateBEPS, ParamBEPS)` | Standard init from type indices |\n| `setup(ps::ParamBEPS; Ta, θ0)` | `(StateBEPS, ParamBEPS)` | Init from pre-built params; verifies `ps_ret === ps` (same object) |\n| `setup(soil::Soil)` | `(StateBEPS, ParamBEPS)` | Conversion from legacy struct |\n| `setup(; VegType, SoilType, ...)` | `(StateBEPS, ParamBEPS)` | Keyword-only form |\n\nAlso tests `setup_jl` (returns `(Soil, StateBEPS, ParamBEPS)`) and `setup_c` (returns `(Soil_c, Vector{Float64}[41], ParamBEPS)`). The `UpdateSoilMoisture compatibility` subtest verifies that a freshly set-up state can be immediately advanced with `r_rain_g` input.\n\n## How `Soil_c` and `Soil` field comparison works\n\n`is_soil_equal` accesses fields by index using a `Base.getindex` override on `Union{Soil, Soil_c}` — it calls `getfield(x, fieldnames(typeof(x))[i])`. This means field correspondence is positional, not by name. The `verbose=true` flag highlights mismatched names in red. Fields in `names_skip = [:θb, :ψb, :θ_vfc]` are excluded because they are Julia-only extensions with no C counterpart.\n\n## Adding a new module test\n\n1. Create `test/modules/test-<module>.jl`.\n2. Add an `include(...)` line in `modules.jl`.\n3. If comparing against C: call both `<func>_jl(...)` and `clang.<func>_c(...)` with identical inputs, then assert `maximum(abs.(r1 .- r2)) <= tol`.\n4. If there is a known C bug, test only the Julia result and add a `#!` comment explaining the divergence.\n5. If the function mutates a struct, construct fresh instances for both paths from the same initial state before calling.","other-path-mnt":"# Other — path_mnt\n\n# `path_mnt`\n\nThis file has no content and no code references — it appears to be an empty placeholder or artifact, likely a stale file from a path mapping or mount-point configuration experiment.\n\n**There is nothing to document.**\n\nIf this file was intentional, add its purpose as a comment or note inside the file itself. If it was created by accident, it can be safely deleted.","other-plan-md":"# Other — plan.md\n\n# ParamBEPS 2.0 Migration Plan\n\n## Overview\n\nThis document describes the migration of `ParamBEPS` from its current flat-field implementation to a layered parameter hierarchy that delegates to `ModelParams.jl` types (`HydraulicProfile`, `ThermalProfile`). The migration is entirely local to `BEPS.jl/src/`; `ModelParams.jl` is a read-only dependency.\n\n**Public API surface is preserved.** The type is still called `ParamBEPS`, existing constructors continue to work, and BEPS physics code is not algorithmically changed.\n\n---\n\n## Why This Migration\n\nThe current `BEPS_Param.jl` stores hydraulic and thermal parameters as flat vectors directly on `ParamBEPS`. `ModelParams.jl` provides composable soil-layer types (`HydraulicProfile`, `ThermalProfile`, `CampbellLayers`) with built-in `filter_params` and `update_params!` utilities. Migrating to these types:\n\n- Enables parameter optimization via `parameters(model)` / `update_params!` without custom boilerplate in BEPS.\n- Makes `ParamBEPS` composable with any solver that dispatches on `AbstractSoilModel{FT,N}`.\n- Eliminates duplicate retention-curve bookkeeping between BEPS and ModelParams.\n\n---\n\n## New Type Shape\n\n```julia\nstruct ParamBEPS{FT, N, H<:HydraulicProfile{FT,N}, T<:ThermalProfile{FT,N}} <: AbstractSoilModel{FT,N}\n    dz         :: Vector{FT}\n    r_drainage :: FT\n    ψ_min      :: FT\n    alpha      :: FT\n    hydraulic  :: H\n    thermal    :: T\n    veg        :: ParamVeg{FT}\nend\n```\n\nThe type parameter `N` encodes the number of soil layers at compile time. The four-parameter form `ParamBEPS{FT,N,H,T}` is the concrete type; call sites that dispatch on `ParamBEPS{FT}` will need broadened signatures or a compatibility alias.\n\n### Type hierarchy\n\n```mermaid\ngraph TD\n    A[AbstractSoilModel{FT,N}] --> B[ParamBEPS{FT,N,H,T}]\n    B --> C[hydraulic::HydraulicProfile{FT,N}]\n    B --> D[thermal::ThermalProfile{FT,N}]\n    B --> E[veg::ParamVeg{FT}]\n    C --> F[CampbellLayers / KvLayers]\n    D --> G[ThermalBaseLayers]\n```\n\n---\n\n## Scope of Changes\n\nAll edits are inside `BEPS.jl/src/`. File-level impact:\n\n| File | What changes |\n|---|---|\n| `DataType/BEPS_Param.jl` | Replace struct definition and constructors; add hydraulic/thermal access helpers; update `Params2Soil!`, `Soil2Params!` |\n| `DataType/Params/Params.jl` | Add `ModelParams` imports; re-export `filter_params`, `update_params!` |\n| `DataType/Params/Param_Init.jl` | `InitParam_Soil` returns `HydraulicProfile` + `ThermalProfile` instead of flat vectors |\n| `DataType/setup.jl` | Update `setup(soil::Soil)` and `_init_state(ps::ParamBEPS, ...)` |\n| `SoilPhysics/*.jl` | Replace direct `ps.K_sat` / `ps.θ_sat` access with helper calls |\n| `inter_prg.jl`, `Optim.jl`, `surface_temperature.jl` | Same field-access updates |\n\nTests in `test/` require assertion updates (field names, not numeric values).\n\n---\n\n## Compatibility Strategy\n\nTwo compatibility layers prevent a flag-day breakage:\n\n**1. BEPS-local access helpers** (preferred long-term)\n\n```julia\ntheta_vwp(ps::ParamBEPS)   # wilting point — not in ModelParams.Campbell, stored on ParamBEPS\ntheta_sat(ps::ParamBEPS)   # ps.hydraulic.profile.θ_sat\nk_sat(ps::ParamBEPS)       # ps.hydraulic.profile.Ksat  (note: Ksat not K_sat in ModelParams)\npsi_sat(ps::ParamBEPS)     # ps.hydraulic.profile.ψ_sat\ncampbell_b(ps::ParamBEPS)  # ps.hydraulic.profile.b\n```\n\nPhysics modules call these helpers instead of destructuring fields. No formula changes.\n\n**2. `BEPSHydraulicView` adapter** (short-term fallback)\n\nIf broad physics edits are deferred, `get_hydraulic(ps::ParamBEPS)` can return a lightweight view struct exposing `θ_vwp`, `θ_sat`, `K_sat`, `ψ_sat`, `b` under the old names. This costs one adapter layer but requires zero edits to physics modules initially.\n\n---\n\n## Unit Conventions\n\nThese must not silently change during migration:\n\n| Quantity | BEPS convention | ModelParams convention | Resolution |\n|---|---|---|---|\n| `K_sat` | cm h⁻¹ | (same in `KvLayers`) | No conversion needed |\n| `ψ_sat` | m, positive suction | cm (docs) | Derive `dz_cm = 100 .* dz`; verify table values |\n| `θ_vwp` | volumetric fraction | not present | Store as BEPS-only field on `ParamBEPS` |\n\n---\n\n## Implementation Order\n\nExecute in this sequence to keep tests green at each step:\n\n1. **Add helpers/adapters** in `BEPS_Param.jl` before touching the struct — existing code still compiles.\n2. **Replace the struct** and add outer constructors (`ParamBEPS(VegType, SoilType; N, FT, kw...)`).\n3. **Update `InitParam_Soil`** to return `HydraulicProfile` + `ThermalProfile`.\n4. **Update `Params2Soil!` / `Soil2Params!`** to copy arrays between old `Soil` layout and new profile types; call `UpdateRootFraction!` at the end.\n5. **Update field access** in `SoilPhysics/`, `inter_prg.jl`, `Optim.jl`, `surface_temperature.jl` to use helpers.\n6. **Update setup dispatches** in `setup.jl`.\n7. **Update tests** — field name assertions only; keep all numeric expectations unchanged. Add coverage for `filter_params(model, :hydraulic)` and `update_params!`.\n\nRun after each step:\n\n```bash\njulia --project=. test/test-BEPSmodel.jl\njulia --project=. test/modules/test-Soil.jl\njulia --project=. test/modules/test-setup.jl\n```\n\nFull suite at the end:\n\n```bash\njulia --project=. -e 'using Pkg; Pkg.test()'\n```\n\n---\n\n## Open Issues\n\nThese must be resolved before coding begins:\n\n1. **`θ_vwp` has no home in `ModelParams.Campbell`.** Options: store as a plain `Vector{FT}` field directly on `ParamBEPS`; subtype `CampbellLayers` to add it; or derive it from `ψ_sat`/`b` at runtime. The plain-field approach is lowest risk.\n\n2. **Four-parameter type breaks two-parameter dispatch.** Existing signatures like `foo(ps::ParamBEPS{FT})` stop matching the new `ParamBEPS{FT,N,H,T}`. Add `const ParamBEPS{FT} = ParamBEPS{FT, N, H, T} where {N,H,T}` UnionAll alias if Julia version allows, or widen dispatch signatures to `ps::ParamBEPS`.\n\n3. **`Ksat` vs `K_sat` naming.** `ModelParams` uses `Ksat`; BEPS code and tests use `K_sat`. The `k_sat(ps)` helper bridges this. No grep-and-replace.\n\n4. **Unit verification for `ψ_sat`.** Confirm BEPS soil lookup tables store `ψ_sat` in m before wiring `dz_cm = 100 .* dz`. A silent factor-of-100 error here would corrupt all retention-curve calculations.\n\n5. **`layers` cache synchronization.** `update_params!` updates `HydraulicProfile.profile` and rebuilds caches. Manual writes inside `Soil2Params!` must call the same cache rebuild path, or the `layers` field will be stale.\n\n---\n\n## Files Not Touched\n\n- `ModelParams.jl/` — read-only reference dependency.\n- `src/clang/` — C-path soil structs are independent; `Soil_c` is unaffected.\n- Numerical formulas in any physics module — this migration changes access patterns only.","other-project-toml":"# Other — Project.toml\n\n# Project.toml — Package Manifest\n\nBEPS.jl is declared as a Julia package via `Project.toml`. This file defines the package identity, all direct dependencies, and the compatibility constraints that the Julia package resolver enforces.\n\n## Package Identity\n\n```toml\nname = \"BEPS\"\nuuid = \"28258c4d-138d-4793-867a-264e32c782b1\"\nversion = \"0.1.11\"\n```\n\nThe registered name is `BEPS`. The UUID is the stable identifier used by the Julia registry and by any downstream package that declares `BEPS` as a dependency. The version follows semantic versioning; at `0.1.x` the public API is still considered unstable.\n\n## Dependency Inventory\n\nDependencies fall into a few functional groups:\n\n**Core scientific computing**\n- `StaticArrays` — stack-allocated arrays used throughout the physics kernels for fixed-size soil layer and canopy layer vectors; avoids heap allocation in the inner time-step loop.\n- `ComponentArrays` — named-field arrays that allow parameter vectors to carry axis metadata, used in optimization and calibration workflows.\n- `Statistics` — mean/std utilities used in output diagnostics.\n\n**Model parameterization**\n- `ModelParams` — the project's own companion library providing `@bounds`, `@units`, `@with_kw`, and related macros that define `ParamBEPS`, `ParamVeg`, `HydraulicProfile`, and `ThermalProfile`. Any change to the `ModelParams` compat bound directly affects what parameter-struct features are available.\n- `Parameters` — supplies `@unpack` / `@pack!` for destructuring structs in physics functions.\n- `UnPack` — thin complement to `Parameters`; provides `@unpack` in modules that avoid the heavier `Parameters` import.\n- `Functors` — enables recursive traversal of parameter structs (used for type promotion and device transfer).\n\n**I/O and data handling**\n- `DataFrames` — all `simulate` and `beps_main` outputs are returned as `DataFrame` objects (`df_flux`, `df_ET`).\n- `DelimitedFiles` — reading raw CSV/TSV forcing files.\n- `OrderedCollections` — `OrderedDict` used in output aggregation and lookup tables.\n- `Serialization` — saving/loading cached state or results to disk.\n- `JSON` — configuration and metadata serialization.\n- `Ipaper` — internal utility package (pretty-print, file I/O helpers).\n\n**Infrastructure**\n- `LazyArtifacts` — defers download of the `libbeps` shared library until it is first needed; keeps the package installable without a C toolchain.\n- `Libdl` — `dlopen`/`dlsym` calls that load `libbeps` at runtime for the C execution path.\n- `Reexport` — `@reexport` is used in submodule entry files (`BEPS_modules.jl`, `DataType/DataType.jl`) to surface symbols without explicit re-listing.\n- `DocStringExtensions` — `$(SIGNATURES)`, `$(TYPEDSIGNATURES)` interpolations in docstrings.\n- `Printf` — formatted diagnostic output.\n- `Dates` — `Date`/`DateTime` types used in the `simulate` API's `dates` argument.\n- `Test` — listed as a direct dep (not under `[extras]`) because test utilities are referenced from within `src/` for assertion helpers.\n- `ComponentArrays` — also used in the optimization module.\n\n## Compatibility Constraints\n\n```toml\njulia = \"1.8, 1.9, 1.10\"\n```\n\nThe package supports Julia 1.8 through 1.10. Julia 1.11+ is not yet in the compat list; if you are running 1.11, the resolver will warn but the code may still work.\n\nKey version pins to be aware of when upgrading dependencies:\n\n| Package | Pinned range | Why it matters |\n|---|---|---|\n| `ModelParams` | `0.2.0` | Exact lower-bound; `@bounds` macro API changed at 0.2 |\n| `StaticArrays` | `1.9.16` | MArray mutation API stabilised here |\n| `DataFrames` | `1.5, 1.6, 1.7` | Three-version range avoids 2.0 breaking changes |\n| `ComponentArrays` | `0.15.29` | Indexing semantics for named axes |\n| `Functors` | `0.5.2` | `fmap` signature used in parameter traversal |\n\n## Notes for Contributors\n\n- **Adding a new dependency:** run `julia --project -e \"using Pkg; Pkg.add(\\\"PkgName\\\")\"` — this updates both `Project.toml` and `Manifest.toml` atomically. Do not edit `Project.toml` by hand for new deps.\n- **Updating compat bounds:** use `Pkg.compat(\"PkgName\", \"X.Y\")` or edit the `[compat]` section and run `Pkg.resolve()`. A missing compat entry for a new dependency will cause a registry CI failure.\n- **`Test` in `[deps]` vs `[extras]`:** `Test` is listed under `[deps]` (not `[extras]`/`[targets]`) because test-assertion helpers are imported from within `src/`. This is intentional but unconventional — do not move it to `[extras]` without auditing all `src/` imports.\n- **`Manifest.toml` is not committed** (see `.gitignore`). Reproducible environments are pinned by `Project.toml` compat ranges, not a locked manifest.","other-readme-md":"# Other — README.md\n\n# BEPS.jl\n\nBoreal Ecosystem Productivity Simulator — a coupled soil-plant-atmosphere (SPAC) model that simulates carbon, water, and energy fluxes at hourly resolution. Ported entirely to Julia from the original C implementation, with a 4.3× speedup over the C version.\n\n## Two Execution Paths\n\nThe package exposes two runtime backends, selectable at call time:\n\n| Backend | Entry point | Notes |\n|---|---|---|\n| Pure Julia | `inter_prg_jl` | Default; 4.3× faster than C |\n| C binding | `inter_prg_c` | Requires `libbeps` shared library; cross-platform (Windows/Linux/Mac) |\n\nBoth backends share the same `beps_main` / `simulate` API surface — switching between them is a single keyword argument.\n\n## Quick Start\n\n```julia\nusing BEPS\nusing Serialization, DelimitedFiles\n\nd   = deserialize(\"data/p1_meteo\")\nlai = readdlm(\"examples/input/p1_lai.txt\")[:]\n\npar = (lon=120.5, lat=30.5, landcover=25, clumping=0.85,\n       soil_type=8, Tsoil=2.2,\n       soilwater=0.4115, snowdepth=0.0)\n\ndf_flux, df_ET, Tg = beps_main(d, lai, par; version=\"julia\")\n```\n\nFor the modern API with finer control over state and output variables, use `simulate` directly — see `beps_modern.jl` and the CLAUDE.md architecture section.\n\n## Parameter System\n\nParameters are organized into three layered structs, all annotated with `@bounds` for automatic optimization bounds and `@with_kw` for keyword construction:\n\n```\nParamBEPS{FT}\n├── hydraulic::ParamSoilHydraulicLayers  ← per-layer Campbell retention params\n├── thermal::ParamSoilThermalLayers      ← per-layer thermal conductivity params\n└── veg::ParamVeg                        ← vegetation + stomatal + Vcmax params\n```\n\nEach leaf struct (e.g., `ParamSoilHydraulic`, `ParamSoilThermal`, `ParamVeg`) uses the `field | (lo, hi)` syntax from `@bounds` to declare optimization search ranges alongside default values. This enables automatic parameter calibration without a separate configuration file.\n\nKey vegetation parameters in `ParamVeg`:\n\n| Field | Description |\n|---|---|\n| `VCmax25` | Max Rubisco capacity at 25°C [μmol m⁻² s⁻¹] |\n| `g1_w` / `g0_w` | Ball-Berry stomatal slope and intercept |\n| `LAI_max_o/u` | Peak LAI for overstory / understory |\n| `r_root_decay` | Exponential decay rate for root distribution |\n| `N_leaf`, `slope_Vc` | Leaf nitrogen → Vcmax scaling |\n\nKey soil hydraulic parameters in `ParamSoilHydraulic`:\n\n| Field | Description |\n|---|---|\n| `θ_vfc`, `θ_vwp`, `θ_sat` | Field capacity, wilting point, saturation |\n| `K_sat` | Saturated hydraulic conductivity [cm h⁻¹] |\n| `ψ_sat`, `b` | Campbell retention curve parameters |\n\n## Hourly Time-Step Architecture\n\n```mermaid\ngraph TD\n    A[MetSeries forcing] --> B[fill_met! timestep i]\n    B --> C[inter_prg_jl]\n    C --> D[Net radiation]\n    D --> E[Aerodynamics]\n    E --> F[Leaf temp / photosynthesis]\n    F --> G[ET: Trans + Eil + Evap]\n    G --> H[Rainfall / Snowpack stages]\n    H --> I[UpdateSoilMoisture!]\n    I --> J[UpdateHeatFlux!]\n    J --> K[Flux / ETFlux output]\n```\n\nEach call to `inter_prg_jl` advances one hour. The `simulate` wrapper manages the loop, state threading, and output collection into DataFrames.\n\n## Known Bugs Fixed\n\nSeveral physics errors have been corrected since the initial port. Notable ones:\n\n- **`LAMBDA` unit error** (2024-10): `lambda_ice` was `333 J/kg` instead of `333000 J/kg`.\n- **Snowpack accumulation** (2024-10): `snowpack_stage3` used `max` instead of `min` for frozen water mass; melt/freeze conditions were gated on sign-change transitions rather than current state, causing perpetual snow accumulation in summer.\n- **`surface_temperature_jl` formula** (2026-01): `T_weighted` was missing a `z_snow` factor in the numerator; `G_soil` used wrong conductivity formula (`G_snow * ΔT/z` instead of `κ_soil * ΔT/Δz`).\n- **`UpdateHeatFlux` argument** (2026-05): Was receiving `Tair` where `Ta_annual` was required, causing significant soil heat flux error.\n- **`V_SOM` domain** (2026-05): Organic matter fraction was initialized in [0, 100] instead of [0, 1].\n\n## Open Issues\n\n- Flux tower validation is the current top-priority task.\n- Leaf temperature passed into photosynthesis may have an incorrect argument (flagged 2025-10, not yet resolved).\n- Soil type parameter tables need review.\n- Clumping index data pipeline is incomplete.\n\n## References\n\n1. Hourly BEPS: Chen et al., [BEPS_hourly_site](https://github.com/JChen-UToronto/BEPS_hourly_site)\n2. Daily BEPS: Liu et al., [BEPS_D](https://github.com/JChen-UToronto/BEPS_D)\n3. Unified C/Julia interface: [CUG-hydro/BEPS.c](https://github.com/CUG-hydro/BEPS.c)","other-repomix-ps1":"# Other — repomix.ps1\n\n# `repomix.ps1` — Source Snapshot Script\n\nA one-shot PowerShell script that bundles the Julia source tree into a single Markdown file for sharing, review, or LLM context ingestion.\n\n## What It Does\n\nRuns [`repomix`](https://github.com/yamadashy/repomix) against the `src/` directory, concatenating all `.jl` files into a self-contained Markdown document at `BEPS_v0.1.9.md` in the repo root.\n\n## Usage\n\n```powershell\n# From the repo root\n./repomix.ps1\n```\n\nOutput: `BEPS_v0.1.9.md` — a flat Markdown file containing every included source file with fenced code blocks.\n\n## Scope\n\n| Setting | Value |\n|---|---|\n| Input root | `src/` |\n| Include pattern | `*.jl`, `*/*.jl`, `*/*/*.jl` (up to two levels deep) |\n| Excluded directories | `backup/`, `clang/` |\n| Output format | Markdown (`--style markdown`) |\n| Output file | `../BEPS_v0.1.9.md` (repo root) |\n\nThe `clang/` exclusion keeps C-binding boilerplate out of the snapshot — those files are auto-generated from the C headers and add noise without semantic value. The `backup/` exclusion avoids stale code copies.\n\n## Versioning\n\nThe output filename `BEPS_v0.1.9.md` is hardcoded. When the package version advances, update the filename to match (e.g. `BEPS_v0.2.0.md`). The script itself has no version-awareness.\n\n## Prerequisites\n\n```bash\nnode --version   # Node.js must be on PATH\nnpx repomix --version\n```\n\n`npx` fetches `repomix` on first run if it is not locally installed. Subsequent runs use the cached version.\n\n## Relationship to the Codebase\n\nThis script is a dev-tooling artifact — it has no runtime relationship to the Julia package, its tests, or any CI pipeline. It exists alongside `Project.toml` and `CLAUDE.md` at the repo root as a convenience for producing shareable source snapshots.","other-test":"# Other — test\n\n# Test Suite — BEPS.jl\n\nThe test suite validates BEPS.jl across four concerns: end-to-end simulation correctness (Julia vs C parity), physical plausibility of state variables, component-level unit behavior, and macro/type infrastructure. Tests are organized as independent files included by `runtests.jl`.\n\n## Running Tests\n\n```bash\n# Full suite\njulia --project test/runtests.jl\n\n# Single file\njulia --project test/test-beps_modern.jl\n```\n\nTests require the `p1_forcing` dataset and `p1_lai.txt` from `data/` and `examples/input/`, loaded via `path_proj(...)`.\n\n## File Map\n\n| File | What it tests |\n|---|---|\n| `test-beps_main.jl` | Legacy `beps_main` API; GPP/Evap regression against known annual totals; Julia vs C numerical parity |\n| `test-beps_modern.jl` | Modern `beps_modern` API; output shape, physical bounds, parameter sensitivity, state immutability |\n| `test-BEPSmodel.jl` | `ParamBEPS` construction across VegTypes and SoilTypes; float precision; soil hydraulic/thermal ranges |\n| `test-UpdateSoilMoisture_Q0.jl` | `UpdateSoilMoisture_Q0!` implicit Richards solver; physical bounds, rainfall/evaporation monotonicity, qualitative agreement with explicit solver |\n| `test-soil_sm.jl` | `soil_sm` standalone module; energy balance, diurnal cycles, snow, soil type variability |\n| `test-photosynthesis_standalone.jl` | Standalone Farquhar model; temperature/light/water stress responses; physical self-consistency |\n| `test-StateSeries.jl` | `StateSeries` construction, `save_state!`, scalar indexing to reconstruct state |\n| `test-macro.jl` | `@DefFluxSeries` / `@DefStateSeries` macros; field generation, `setindex!`, `getindex` slicing |\n| `test-ModelParams.jl` | `parameters()`, `update!()`, `filter_params()` on `ParamBEPS` |\n| `test-utilize.jl` | `Layer2`/`Layer3` structs, `VCmax`, `snow_density` |\n| `test-performance.jl` | Same regression as `test-beps_main.jl` plus `@profview` harness (not run in CI) |\n\n## Regression Anchors\n\nTwo hard numerical targets guard against silent physics regressions:\n\n```julia\n# test-beps_main.jl — full 2010 year, CRO/silty_clay_loam\n@test isapprox(r.GPP,  2146.110; atol=2)\n@test isapprox(r.Evap,  62.5378; atol=1)\n```\n\nThese are annual sums over 8760 hourly timesteps. The tolerances are intentionally loose enough to survive minor numerical drift while catching algorithmic changes.\n\n## Julia vs C Parity (`test-beps_main.jl`)\n\n`beps_main` dispatches to either `inter_prg_jl` (pure Julia) or `inter_prg_c` (C binding via `libbeps`). The test runs both and computes column-wise percent differences:\n\n```julia\ndf_diff_perc = abs.(df_jl .- df_c) ./ df_c .* 100\n```\n\nGPP, Evap, and Trans are compared; the current passing threshold is informally ~2.5% but the `@test true` guard means CI always passes — the output is printed for human review. Tighten the bound by uncommenting the `_nanmaximum(l) < 1.5` assertions when the C path is stable.\n\n## Physical Invariants\n\nTests across `test-beps_modern.jl`, `test-UpdateSoilMoisture_Q0.jl`, and `test-soil_sm.jl` share a common pattern: run the model for multiple timesteps, then assert bounds on state arrays.\n\nKey invariants checked:\n- `0 ≤ θ ≤ 1` for all soil layers at every timestep\n- `0 ≤ f_stress ≤ 1` (soil water stress factor)\n- `all(isfinite, ...)` on temperatures, conductances, GPP\n- `GPP ≥ 0`\n- Monotonicity: rainfall → `θ[1]` increases; evaporation-only → total `θ` decreases\n\n## Soil Moisture Solver Tests (`test-UpdateSoilMoisture_Q0.jl`)\n\nThe implicit Bonan-Q0 Richards solver (`UpdateSoilMoisture_Q0!`) is tested against both physical bounds and the existing explicit Darcy solver (`UpdateSoilMoisture`). The shared setup helper:\n\n```julia\nfunction make_state(; SoilType=4, VegType=\"DBF\", θ0=0.3, r_rain_g=0.0)\n    st, ps = setup(VegType, SoilType; θ0)\n    st.r_rain_g = r_rain_g\n    st, ps\nend\n```\n\nThe qualitative agreement test runs 24 hourly steps with both solvers and asserts the mean absolute difference in `θ` stays below 0.15 — a wide band that acknowledges the two discretization schemes are not expected to match numerically, only directionally.\n\n## Photosynthesis Standalone Tests (`test-photosynthesis_standalone.jl`)\n\nTests `BEPS.Photosynthesis.photosynthesis(T, RH, Srad, LAI, params)` independently of the canopy coupling. Key scenarios:\n\n- **Temperature response**: bell-shaped, with peak at ~15 °C (Bernacchi/Medlyn parameterization). Asserts `An(15°C) > An(-5°C)` and `An(15°C) > An(40°C)`.\n- **Light response**: `An(Srad=0) < 0` (dark respiration dominates); `An(Srad=500) > An(Srad=0)`.\n- **Water stress**: `An(β=1.0) > An(β=0.3)`.\n- **Cross-module comparison**: compares against the canopy-coupled `photosynthesis_jl` at five temperatures. Both must return `An ≥ 0` and `0 < Ci < 380 μmol/mol` in the 0–30 °C range. No tight numerical match is required because the two implementations use different Michaelis constants and temperature response functions (Harley & Baldocchi 1995 vs Bernacchi et al. 2001).\n\n## Macro Infrastructure Tests (`test-macro.jl`)\n\n`@DefFluxSeries` and `@DefStateSeries` generate time-series container structs from a source `Flux`/`State` type. The tests verify the full generated interface:\n\n```julia\n@DefFluxSeries MockSeries = MockFlux          # fields: ntime, a::Vector, b::Vector\n@DefFluxSeries ExtraSeries = MockFlux extra_c # adds extra_c::Vector\n@DefStateSeries MockStateSer = MockState\n```\n\nAssertions cover:\n- Field presence and `Vector{FT}` element types\n- `s[i] = flux` scatter (setindex!)\n- `s[r]` slice returns a new series with correct sub-vectors\n- `ss[i]` reconstruction returns the original State type\n- Subtype relationships: `MockSeries{Float64} <: AbstractFluxSeries{Float64}`\n\n## `ParamBEPS` Construction Tests (`test-BEPSmodel.jl`)\n\nCovers four PFTs (ENF, DBF, EBF, C4) × selected soil types, checking lookup-table values propagate correctly into nested structs. Notable checks:\n\n- DBF and EBF share `ψ_min = 1000.0`, `alpha = 1.5` (more isohydric); other PFTs get `ψ_min = 3300.0`.\n- C4 grass gets `g1_w = 4.0` vs C3 default of `8.0` (Ball-Berry slope).\n- `N` parameter controls `length(hydraulic.θ_sat)` — verified for `N=3` and `N=5`.\n- Float precision: `ParamBEPS{Float32}` propagates `Float32` through all sub-arrays.\n\n## `StateSeries` Tests (`test-StateSeries.jl`)\n\n```julia\nconst SF = (:z_water, :z_snow, :r_rain_g, :f_soilwater) |> Val  # scalar fields\nconst VF = (:θ, :Tsoil_c, :ETi, :G) |> Val                      # vector fields (N layers)\n```\n\n`StateSeries(SF, VF, NLAYER, ntime)` allocates scalar and vector time-series containers. `save_state!(out, state, t, SF, VF)` copies one timestep. `out[1]` reconstructs a named tuple with the saved values, verified against the initial `θ` profile from `setup(model)`.","other":"# Other — 结题\n\n# `docs/结题/` — Project Completion Reports\n\nThis directory contains the Typst-formatted final reports for the ModernBEPS.jl development project. The documents collectively demonstrate model accuracy, describe the data pipeline, and summarize the engineering and scientific work completed during the project.\n\n## Document Overview\n\n| File | Purpose |\n|---|---|\n| `ModernBEPS_精度证明.typ` | Model accuracy proof: site selection, calibration workflow, GOF metrics |\n| `数据方面.typ` | Data processing summary: ChinaFlux dataset compilation and QC |\n| `结题_PART2.typ` | Full completion report: development work, technical achievements, conclusions |\n\nAll three documents use the `@local/modern-cug-report:0.1.3` Typst template. `结题_PART2.typ` includes `数据方面.typ` directly via `#include \"数据方面.typ\"`, and also includes the module-level technical chapters from `../ch01_SoilMoisture/`, `../ch02_SoilTemperature/`, `../ch03_SurfaceFluxes/`, and `../ch04_Photosynthesis/`.\n\n## Validation Dataset\n\nThe project validates against **11 forest ecosystem flux stations** across China (12 were compiled; `MF_乔灌混交林_燕山` was excluded due to only half a year of data and magnitude errors in RE/GPP/ET). Sites span four vegetation types:\n\n- **DBF** (deciduous broadleaf): 宝天曼, 小浪底\n- **EBF** (evergreen broadleaf): 西双版纳热带雨林, 海南儋州橡胶林, 西双版纳橡胶林, 金佛山亚热带常绿阔叶林\n- **ENF** (evergreen needleleaf): 呼中, 千烟洲, 燕山\n- **MF** (mixed forest): 鼎湖山, 长白山\n\nEach station entry records `lon`, `lat`, `VegType`, `SoilType`, soil temperature depths (`z_TS`), soil moisture depths (`z_SM`), wind speed height (`z_Uz`), and canopy height (`z_overstory`). These fields map directly to `ParamBEPS` and `simulate` arguments.\n\n## Calibration Workflow\n\nThe full calibration pipeline is in `Project_ChinaFlux/case01_ChinaFlux.jl`. The `RunModel` function encapsulates the end-to-end flow:\n\n```\nLoad forcing (hourly) → Parse site metadata → Initialize ParamBEPS + StateBEPS\n  → simulate() → BEPS_GOF() [default params]\n  → optim() [SCE-UA, maxn=1000]\n  → goodness() [optimized params]\n  → jldsave() [persist gof_opt, gof, theta_opt, data_sim, data_obs]\n```\n\n**Key implementation details:**\n\n- Dates are converted from local CST (UTC+8) to UTC before passing to `simulate`: `dates_UTC = dates_local .- Hour(8)`\n- Hourly forcing is truncated to match the daily validation data length when there is a mismatch\n- `sanitize_forcing!` performs quality control on the forcing DataFrame before use\n- Missing values are replaced with `replace_missing!` after reading CSVs\n\n**Optimized parameters** (5 selected for sensitivity to GPP/ET):\n\n| Parameter | Path | Controls |\n|---|---|---|\n| `r_drainage` | `[:r_drainage]` | Lateral drainage rate |\n| `Ω` | `[:veg, :Ω]` | Clumping index |\n| `g1_w` | `[:veg, :g1_w]` | Stomatal slope (Ball-Berry) |\n| `g0_w` | `[:veg, :g0_w]` | Minimum stomatal conductance |\n| `VCmax25` | `[:veg, :VCmax25]` | Maximum carboxylation rate at 25°C |\n\nThe `optim` function accepts `paths` as a vector of symbol-path vectors, and `goal` switches the objective function between `:NSE`, `:KGE`, and `:R2` (with `goal_multiplier=-1` for minimization).\n\n## Performance Results\n\nPost-calibration median metrics across 11 sites (daily GPP and ET):\n\n- **NSE** improvement: +0.247 (GPP), +0.086 (ET)\n- **R² before optimization**: GPP = 0.69, ET = 0.53\n- **R² after optimization**: GPP = 0.65, ET = 0.66\n\nThe R² target of ≥ 0.6 is met for both GPP and ET post-calibration. NSE and KGE show the largest gains, indicating the optimization primarily corrects systematic bias rather than correlation structure. The slight R² drop in GPP after calibration is expected when the optimizer targets NSE (which penalizes bias and variance jointly) rather than R² directly.\n\n## Model Engineering Achievements (`结题_PART2.typ`)\n\nThis document is the most comprehensive and records the technical progress made during the project lifecycle. Key restructuring results:\n\n**Signature reduction (C → Julia):**\n\n- `sensible_heat`: 19 arguments → 7 (via `Leaf`, `Gheat` structs)\n- `transpiration`: 16 arguments → 5\n\n**Seven development tracks completed:**\n\n1. Julia reimplementation and modular reorganization of the BEPS main model\n2. State/parameter/cache variable system (`StateBEPS`, `ParamBEPS`, `LeafCache`)\n3. Soil moisture physics (Richards equation, Campbell hydraulics, numerical stability)\n4. Soil temperature and surface energy boundary corrections\n5. Snow, rainfall, and ET process integration\n6. Radiation, aerodynamic resistance, and photosynthesis module cleanup\n7. Parameter optimization and multi-site evaluation framework\n\n**Computational performance:** The Julia implementation runs at ~4.3× the speed of the C version after cache mechanism optimization. A 13-year hourly simulation completes in seconds.\n\n## Data Pipeline (`数据方面.typ`)\n\nThe upstream data processing uses the R package `ChinaFlux2026.R`, which provides:\n\n- Automated variable name normalization, unit unification, and timestamp standardization\n- Quality control, gap filling, and LAI augmentation\n- Batch plotting for data verification\n\nThe full compiled dataset covers **39 flux stations** across China (CRO, DBF, EBF, ENF, GRA, MF, SAV, WET, WSA), with DOIs recorded for each station in `数据方面.typ`. The forest subset (12 stations) feeds directly into the ModernBEPS validation pipeline.\n\n## Rendering\n\nThese documents are compiled with the Typst toolchain. The template `@local/modern-cug-report:0.1.3` must be installed in the local Typst package registry. Figures referenced (e.g., `../Figures/Figure1_GOF_NSE.png`, `../Figures/数据核验/*.pdf`) are generated by the R analysis scripts and must exist before compilation. The `#include` directives in `结题_PART2.typ` pull in sibling `.typ` files at compile time, so the full document tree must be present.","overview":"# BEPS.jl — Wiki\n\n# BEPS.jl\n\n**Boreal Ecosystem Productivity Simulator** — a coupled soil-plant-atmosphere (SPAC) model that simulates carbon, water, and energy fluxes at hourly resolution. The pure-Julia implementation runs **4.3× faster than the original C version**, while a C-binding path (`libbeps`) remains available for comparison and legacy workflows.\n\n## What it does\n\nGiven hourly meteorological forcing (air temperature, radiation, humidity, wind, precipitation) and a leaf area index time series, BEPS.jl simulates the full surface energy balance: gross primary productivity, net ecosystem exchange, transpiration, interception, snowpack dynamics, and soil heat and moisture — hour by hour over a multi-year run. The [ChinaFlux case studies](chinaflux-case-studies.md) show the model calibrated and validated against ~12 Chinese flux tower sites.\n\n## Architecture at a glance\n\n```mermaid\ngraph TD\n    F[MetSeries<br/>Forcing] --> S[simulate / beps_main]\n    P[setup<br/>ParamBEPS + StateBEPS] --> S\n    S --> I[inter_prg_jl<br/>hourly engine]\n    I --> SE[Surface Energy Fluxes]\n    I --> PS[Photosynthesis]\n    I --> SP[Soil Physics]\n    I --> HY[Hydrology & Snowpack]\n    SE --> SPAC[SPAC & Canopy]\n    SP --> O[Flux / ETFlux output]\n```\n\nThe diagram captures the entire runtime call chain. Each box below is a self-contained module with its own wiki page.\n\n## Getting started\n\n```julia\n] add https://github.com/jl-pkgs/ModernBEPS.jl\n```\n\n```julia\nusing BEPS\n\n# Initialize parameters and state for a site\nstate, ps = setup(\"evergreen_needleleaf\", \"silty_clay_loam\")\n\n# Run simulation\ndf_flux, df_ET, states, caches = simulate(forcing, lai, dates;\n    ps, state,\n    lon = 120.0, lat = 40.0)\n```\n\n`setup` builds a [`ParamBEPS` and `StateBEPS`](data-types-and-state.md) pair from built-in vegetation and soil lookup tables. `simulate` returns tidy DataFrames of hourly fluxes alongside optional time-series of internal state variables.\n\n## The two execution paths\n\nBEPS.jl ships two implementations of the hourly physics kernel:\n\n- **`inter_prg_jl`** — the native Julia path. This is the default and the one being actively developed.\n- **`inter_prg_c`** — thin Julia wrappers around `libbeps`, the original compiled C library. Documented in the [C Library Bindings](c-library-bindings.md) module. Useful for cross-validation and Windows/Linux/Mac compatibility is ensured.\n\nBoth paths share the same [data types](data-types-and-state.md) and are dispatched through the same `beps_main` / `simulate` API.\n\n## Inside the hourly engine\n\nEach call to `inter_prg_jl` runs one hour of physics in this order:\n\n1. **Radiation & aerodynamics** — net radiation partitioned across overstory, understory, and ground layers; aerodynamic conductances computed from wind and LAI. Handled by the [SPAC & Canopy](spac-and-canopy.md) module.\n2. **Leaf temperature & conductance** — iterative solution of the leaf energy balance coupling stomatal conductance to vapour pressure deficit. Part of [Surface Energy Fluxes](surface-energy-fluxes.md).\n3. **Photosynthesis** — Farquhar-Ball-Berry model solved for each of four leaf classes (overstory/understory × sunlit/shaded). See the [Photosynthesis](photosynthesis.md) module.\n4. **Evapotranspiration** — sensible and latent heat separated per canopy layer; interception and throughfall tracked. [Hydrology & Snowpack](hydrology-and-snowpack.md) handles the rain/snow cascade.\n5. **Soil physics** — heat diffusion with freeze-thaw, soil moisture transport (Richards equation or bucket), and root uptake. Documented in [Soil Physics](soil-physics.md).\n\n## Parameters and state\n\nAll physical constants, vegetation traits, and soil hydraulic/thermal parameters live in [Model Parameters](model-parameters.md) (`src/DataType/Params/`). Every vegetation parameter carries explicit bounds metadata via `@bounds`, enabling the [Optimization & Calibration](optimization-and-calibration.md) module to run SCE-UA parameter estimation against observed fluxes without any manual bookkeeping.\n\n## Running the tests\n\n```bash\n# Full suite\njulia --project test/runtests.jl\n\n# Single focused test\njulia --project test/test-beps_modern.jl\njulia --project test/test-soil_sm.jl\n```\n\n## Where to go next\n\n- New to the physics? Start with [Surface Energy Fluxes](surface-energy-fluxes.md) and [Photosynthesis](photosynthesis.md).\n- Setting up a new site run? See [Data Types & State](data-types-and-state.md) for `setup()` and forcing conventions, then [Utilities](utilities.md) for forcing QC helpers.\n- Calibrating parameters? Go to [Optimization & Calibration](optimization-and-calibration.md), then read the [ChinaFlux Case Studies](chinaflux-case-studies.md) for a worked example.\n- Working with the C path? See [C Library Bindings](c-library-bindings.md).","photosynthesis":"# Photosynthesis\n\n# Photosynthesis Module\n\nThe photosynthesis module implements leaf-level carbon assimilation using the Farquhar-von Caemmerer-Berry (FvCB) model coupled with Ball-Berry stomatal conductance. It exists in two forms: a **coupled implementation** tightly integrated with the BEPS canopy model, and a **standalone implementation** for isolated use or testing.\n\n---\n\n## Two Implementations\n\n### Coupled (`src/photosynthesis.jl` + `src/photosynthesis_helper.jl`)\n\nThe primary path used during model simulation. `photosynthesis_jl` is the single-leaf solver; `photosynthesis` (in `photosynthesis_helper.jl`) fans it out across the four leaf classes in the canopy:\n\n| Leaf class | Description |\n|---|---|\n| `o_sunlit` | Overstory, sun-exposed |\n| `o_shaded` | Overstory, shaded |\n| `u_sunlit` | Understory, sun-exposed |\n| `u_shaded` | Understory, shaded |\n\nEach class receives its own radiation (`R.o_sunlit`, etc.), temperature (`Tc_old`), and `Vcmax`. Boundary layer conductances differ between overstory (`Gb_o`) and understory (`Gb_u`). The dispatch via `version=\"julia\"` or `version=\"c\"` selects between `photosynthesis_jl` and `photosynthesis_c` (C library binding).\n\n### Standalone (`src/standalone/Photosynthesis/`)\n\nA self-contained module for use outside the full SPAC loop. It exposes a simpler `photosynthesis(Tair, RH, Srad, LAI, params)` entry point that internally iterates to solve the coupled An–gs–ci system. Results are returned as a `PhotoResult{FT}` struct containing `LeafPhoto` (sunlit/shaded) fields for `An`, `Gs`, `Ci`, `Gc`, and `Rd`.\n\n---\n\n## Core Algorithm\n\nBoth implementations follow the same FvCB logic:\n\n**1. Enzyme kinetics (Rubisco-limited rate)**\n```\nWc = Vcmax × (ci − Γ) / (ci + K)\n```\nwhere `K` is the effective Michaelis-Menten constant for CO₂ in the presence of O₂, and `Γ` is the CO₂ compensation point.\n\n**2. Electron transport (light-limited rate)**\n```\nJₓ = Jmax × PPFD / (PPFD + 2.1 × Jmax)      # Chen 1999, Eq. 6\nWj = Jₓ × (ci − Γ) / (4 ci + 8 Γ)\n```\n\n**3. Rate selection** — whichever of `Wj` or `Wc` is smaller determines the coefficients `(a, b, e)` used in the cubic/quadratic solver.\n\n**4. Coupled An–gs–ci solution**\n```\ngs = g0 + g1 × RH_leaf × β_soil × Ag / cs     # Ball-Berry\n```\nSubstituting the Ball-Berry model into the diffusion equations yields a cubic in `An` (when `An > 0`) solved by `solve_cubic`, or a quadratic solved by `solve_quad` when `An ≤ 0` or the light/Rubisco rates fall below `Rd`.\n\n**5. Sucrose limitation** caps the result: `An ≤ Vcmax/2 − Rd`.\n\n---\n\n## Temperature Response\n\n`init_photo_consts(T_leaf_K)` (in `DataType/PhotoConsts.jl`) pre-computes the factors `Γ`, `K`, `Rd_factor`, `Jmax_factor`, `Vcmax_factor` for a given leaf temperature. These can be passed in as a `PhotoConsts{T}` to avoid recomputation when the solver is called in a loop.\n\nThe standalone module uses explicit Medlyn-form peaked Arrhenius functions:\n\n- `fTv(T, Vcmax25, evc, toptvc)` — Vcmax, normalized to 25 °C with deactivation terms (Hd = 200 kJ mol⁻¹, S = 640 J mol⁻¹ K⁻¹)\n- `fTj(T, Jmax25, ejm, toptjm)` — same form for Jmax\n- `fTd(T, Rd25, erd)` — simple Arrhenius for dark respiration\n- `TBOLTZ(T, E)` — bare Arrhenius factor, used for Kc, Ko, τ\n\nThe coupled module uses a different formulation (`TBOLTZ` with explicit `topt` in `photosynthesis_helper.jl`) and the older linear relationship `jmopt = 2.39 × Vcmax25 − 14.2`.\n\n---\n\n## Solver Details\n\n```mermaid\ngraph TD\n    A[photosynthesis_jl] --> B{Wj or Wc < Rd?}\n    B -- No --> C[solve_cubic: cubic An]\n    B -- Yes --> D[solve_quad: quadratic An]\n    C --> E{An <= 0?}\n    E -- Yes --> D\n    E -- No --> F[apply sucrose cap]\n    D --> F\n    F --> G[compute gs_w, ci]\n```\n\n**`solve_cubic`** implements the Vieta / Numerical Recipes trigonometric method for three real roots. After sorting, `findroot` picks the physically meaningful root: the minimum positive root when all roots are positive, or the sign-change-spanning root otherwise.\n\n**`solve_quad`** assumes `gs = g0` (cuticular conductance only), reducing the system to a standard quadratic.\n\n---\n\n## Unit Conventions\n\n| Quantity | Unit | Notes |\n|---|---|---|\n| Conductances (`gb_w`, `gs_w`) | s m⁻¹ | Used internally; converted via `umol_m` / `m_umol` |\n| Conductances (`gb_c_mol`, `gs_w_mol`) | µmol m⁻² s⁻¹ | Molar form for FvCB equations |\n| CO₂ concentrations | ppm (µmol mol⁻¹) | `ca`, `ci`, `cs` |\n| Assimilation | µmol m⁻² s⁻¹ | `An`, `Rd` |\n| Radiation (`Rsn_leaf`) | W m⁻² | Converted: `PPFD = 4.55 × 0.5 × Rsn_leaf` |\n\nThe 1.6 factor (`g_H2O / g_CO2`) appears throughout when converting between H₂O and CO₂ conductances.\n\n---\n\n## Integration with BEPS\n\n`photosynthesis` in `src/photosynthesis_helper.jl` is called from `inter_prg_jl` (the hourly time-step orchestrator) after leaf temperature and aerodynamic resistances are resolved. It writes results into mutable `Leaf` structs:\n\n```julia\nphotosynthesis(Tc_old, R, Ci_old, leleaf, Ta, ea, f_soilwater,\n               g0_h2o, g1_hio, Gb_o, Gb_u, Vcmax_sunlit, Vcmax_shaded,\n               Gs_new, Ac, Ci_new;          # output structs, modified in-place\n               version=\"julia\", pc=nothing)\n```\n\n`β_soil` (`f_soilwater`) scales stomatal conductance to represent water stress; it is computed upstream in `soil_water_factor_v2.jl`. The `pc::PhotoConsts` kwarg lets the caller pass pre-computed temperature constants to avoid redundant work across the four leaf classes.\n\n---\n\n## Key Differences Between the Two Implementations\n\n| Aspect | Coupled (`photosynthesis_jl`) | Standalone (`farquhar_model`) |\n|---|---|---|\n| Stomatal coupling | Cubic/quadratic solver, fully coupled | Iterative ci loop (up to 15 iterations) |\n| Jmax relationship | `jmopt = 2.39 Vcmax25 − 14.2` | Explicit `Jmax25` parameter |\n| Temperature response | `init_photo_consts` / `TBOLTZ` with topt | Medlyn peaked Arrhenius (`fTv`, `fTj`) |\n| TPU limitation | Sucrose cap: `Vcmax/2 − Rd` | Explicit `Wp = 0.5 Vcmax` |\n| Entry point | Four leaf classes via `photosynthesis()` | Single leaf via `farquhar_model()` |\n| RH at leaf surface | Computed from LE via `SFC_VPD` | Passed directly as forcing |","soil-physics":"# Soil Physics\n\n# Soil Physics Module\n\n`src/SoilPhysics/`\n\nThis module implements the coupled soil heat and moisture physics for BEPS. It handles heat diffusion with freeze-thaw, two alternative soil moisture solvers, root water uptake, and soil water stress. It is called once per timestep from `inter_prg_jl`, after canopy fluxes are resolved but before output collection.\n\n## State and Parameter Access\n\nFunctions in this module accept either `StateBEPS`/`ParamBEPS` (the modern JAX-style API) or the legacy `Soil` struct, which carries both state and parameters in one object. Adapter functions at the bottom of `SoilPhysics.jl` handle the dispatch:\n\n```julia\nget_hydraulic(ps::ParamBEPS) = ps.hydraulic   # → HydraulicProfile\nget_hydraulic(ps::Soil)      = ps             # Soil is its own hydraulic source\nget_thermal(ps::ParamBEPS)   = ps.thermal\nget_thermal(ps::Soil)        = ps\nget_root_decay(ps::ParamBEPS) = ps.veg.r_root_decay\nget_root_decay(ps::Soil)      = ps.r_root_decay\n```\n\nA separate conversion normalizes the sign convention for matric potential: `ParamBEPS` stores `ψ_sat` in negative cm (ModelParams convention), while legacy code used positive m. `_get_ψ_sat_m` converts both to positive m for callers that need it.\n\n---\n\n## Heat Subsystem\n\n### Thermal conductivity — `UpdateThermal_κ`\n\nCalled from `surface_temperature.jl` before the heat flux update. Uses the LHE equation (eq. 8) blending dry-soil, ice, and water conductivities:\n\n```\nκ[i] = (κ_dry^(1-θ_sat) · ki^(1.2·θ·ice_ratio) · kw^(θ·(1-ice_ratio)) − 0.15) · Sr + 0.15\n```\n\nwhere `Sr = θ/θ_sat` is the relative saturation and the 0.15 floor is the Juweimin (2005) lower bound. Constants: `ki = 2.1` (ice), `kw = 0.61` (water).\n\n### Volumetric heat capacity — `UpdateThermal_Cv`\n\nAlso called from `surface_temperature.jl`. Follows Chen Baozhang (2007, *Ecol. Modelling* 209) eq. 18:\n\n```\nCv[i] = 2.0e6 · ρ_soil/2650   # mineral solid\n      + 1.0e6 · θ · (4.2·(1−ice_ratio) + 2.09·ice_ratio)  # liquid + ice\n      + 2.5e6 · V_SOM          # organic matter\n```\n\nUnits: MJ m⁻³ K⁻¹.\n\n### Heat flux and temperature update — `UpdateHeatFlux`\n\nFinite-difference heat diffusion across `n` layers plus one bottom boundary. The flux between layers uses harmonic-mean resistance:\n\n```\nG[i] = 2(Tsoil_p[i-1] − Tsoil_p[i]) / (dz[i-1]/κ[i-1] + dz[i]/κ[i])\n```\n\nThe bottom boundary layer drives toward `Tair_annual_mean` over a fixed depth `DEPTH_F`. Temperatures are updated explicitly:\n\n```\nTsoil_c[i] = Tsoil_p[i] + (G[i] − G[i+1]) / (Cv[i] · dz[i]) · Δt\n```\n\nPass `fix_Tsoil=true` to skip the temperature update (observation-prescribed mode) while still running `Update_ice_ratio`.\n\n### Freeze-thaw — `Update_ice_ratio`\n\nRuns after every heat flux step. When a layer crosses 0 °C, latent heat of fusion (`Lf0 = 3.34×10⁵ J kg⁻¹`) determines how much of the temperature excursion goes into phase change instead of sensible heat. The layer is pinned to 0 °C until all water has changed phase. After the phase update, `ice_ratio` is rescaled by `θ_prev/θ` to account for concurrent moisture changes.\n\nThe simplified scalar version `get_ice_ratio(Tsoil)` maps temperature linearly: fully frozen below −1 °C, fully liquid above 0 °C, linear between.\n\n---\n\n## Moisture Subsystem\n\n### Surface water and infiltration\n\nBoth moisture solvers share the same surface entry point:\n\n**`update_surface_water!(st, ps, kstep) → inf`**\n\n1. Calls `update_SoilWaterFrac!` to compute `f_water[i]` — the liquid fraction of pore water — as a piecewise linear function of `Tsoil_c`: 1.0 above 0 °C, 0.1 below −1 °C, linearly blended between.\n2. Computes maximum infiltration capacity (`inf_max`) from the Green-Ampt-like expression using layer-1 Campbell hydraulics.\n3. Clamps actual infiltration to available rainfall + ponded water.\n4. Updates `st.z_water` (ponded depth) after applying `r_drainage` runoff fraction.\n\nReturns infiltration in cm h⁻¹.\n\n### Explicit adaptive solver — `UpdateSoilMoisture` / `solve_SM_beps`\n\nThe original solver. Uses a variable-timestep explicit loop: each sub-step length is chosen by `guess_step(max_Fb)` based on the maximum inter-layer flux (Darcy velocity in cm h⁻¹):\n\n| `max_Fb` threshold | Sub-step `Δt` |\n|--------------------|---------------|\n| > 3.6 cm h⁻¹ | 1 s |\n| > 0.36 cm h⁻¹ | 30 s |\n| else | 360 s |\n\nWithin each sub-step, hydraulic conductivity follows Campbell (1974):\n```julia\nK[i]  = K_sat[i] · (θ/θ_sat)^(2b+3)      # cal_K\nψ[i]  = ψ_sat[i] · (θ/θ_sat)^(-b)         # cal_ψ  [cm, negative]\n```\n\nInter-layer flux uses a weighted-harmonic average conductivity and a Darcy gradient that includes the gravitational head:\n```\ngrad_ψ = 1 − (ψ[i+1] − ψ[i]) / Δz_cm\nQ[i]   = Kavg[i] · grad_ψ                  # [cm h⁻¹]\n```\n\n`Q` is capped at `Q_max` to prevent over-draining the receiving layer. The bottom layer has zero outflow (zero-flux lower boundary, following CLM 3.5).\n\nθ is updated as:\n```\nθ[i] += (inflow − outflow − ETi[i]) · Δt / dz[i]\n```\nthen clamped to `[θ_res, θ_sat]`.\n\n### Implicit Crank-Nicolson solver — `UpdateSoilMoisture_Q0!` / `solve_SM_beps_Q0!`\n\nA drop-in replacement for the explicit solver (same signature, same `update_surface_water!` call). Delegates to `soil_moisture_Q0!` from `ModelParams.jl`, which solves one implicit tridiagonal system per timestep using the `TriSolver` workspace pre-allocated in `StateBEPS`.\n\n```julia\nfunction UpdateSoilMoisture_Q0!(st, ps, kstep; fix_sm=false)\n    θ_prev .= θ\n    inf = update_surface_water!(st, ps, kstep)\n    fix_sm && return\n    solve_SM_beps_Q0!(st, ps, inf, kstep)\n    # rescale ice_ratio for moisture change\nend\n```\n\nThe top boundary condition is `Q0 = -inf` (negative = downward). Sink terms (`ETi`, in m s⁻¹) are converted to cm h⁻¹ before passing to the solver. After the solve, `ψ` may be positive under saturated conditions — this is physically valid and must not be forced negative.\n\n### Choosing a solver\n\n| | `UpdateSoilMoisture` | `UpdateSoilMoisture_Q0!` |\n|---|---|---|\n| Method | Explicit adaptive | Implicit Crank-Nicolson |\n| Sub-steps | Many (1–360 s each) | One per timestep |\n| Stability | CFL-limited | Unconditionally stable |\n| Backend | Self-contained | Requires `ModelParams.TriSolver` workspace |\n\nFor typical hourly BEPS timesteps, the implicit solver is preferred for stiff conditions (high K_sat, near-saturated layers). Both are exercised in `test/test-UpdateSoilMoisture_Q0.jl`.\n\n---\n\n## Root Water Uptake and Stress\n\n### `UpdateRootFraction!`\n\nDistributes roots across layers using an exponential decay profile:\n\n```\nf_root[i] = β^z[i-1] − β^z[i]     # cumulative depth z in cm\nf_root[n] = β^z[n-1]               # remainder to deepest layer\n```\n\n`β = r_root_decay` is a per-vegetation decay constant (from `ParamVeg`). This is called once at setup by `_init_state`, not per timestep.\n\n### `Root_Water_Uptake`\n\nSplits total canopy transpiration (`Trans_o + Trans_u`) and soil evaporation across layers. The surface layer receives both transpiration (weighted by `w_norm[1]`) and soil evaporation:\n\n```julia\nETi[1] = Trans/ρ_w * w_norm[1] + Evap_soil/ρ_w\nETi[i] = Trans/ρ_w * w_norm[i]     # deeper layers\n```\n\n`w_norm` is the normalized root weight from `soil_water_factor_v2`. Units: m s⁻¹.\n\n### `soil_water_factor_v2`\n\nComputes the scalar soil water stress factor `f_soilwater` and per-layer weights `w_norm`. Two stress terms are multiplied:\n\n**Water stress** (He 2017, *JGR-B* eq. 4):\n```\nf_stress[i] = 1 / (1 + ((|ψ[i]| − ψ_min) / ψ_min)^α)    if |ψ[i]| > ψ_min\n            = 1.0                                            otherwise\n```\n\n**Temperature stress**:\n```\nf_temp[i] = 1 − exp(−0.02 · Tsoil_p[i]²)    if Tsoil_p[i] > 0\n           = 0                                 otherwise\n```\n\nRoot weight per layer: `w_root[i] = f_root[i] · f_stress[i]` (when `FW_VERSION == 1`). Normalized: `w_norm[i] = w_root[i] / Σw_root`. The scalar `f_soilwater` is the weighted mean of `f_stress`, floored at 0.1.\n\n---\n\n## Dataflow Within a Timestep\n\n```mermaid\ngraph TD\n    A[inter_prg_jl] --> B[Root_Water_Uptake<br/>fill ETi]\n    A --> C[UpdateThermal_κ<br/>update κ]\n    A --> D[UpdateThermal_Cv<br/>update Cv]\n    C --> E[UpdateHeatFlux<br/>solve T, ice_ratio]\n    D --> E\n    B --> F[UpdateSoilMoisture<br/>or _Q0!]\n    F --> G[soil_water_factor_v2<br/>update f_soilwater, w_norm]\n```\n\n`UpdateThermal_κ` and `UpdateThermal_Cv` are called from `surface_temperature.jl` (before `UpdateHeatFlux`). `Root_Water_Uptake` must run before the moisture solver so `ETi` is populated.\n\n---\n\n## Notes for Contributors\n\n- **Unit conventions**: `ψ` is stored in cm and is negative (ModelParams convention). `cal_ψ` enforces `ψ ≤ ψ_sat`. `K` in `solve_SM_beps` uses cm h⁻¹; in the implicit solver it follows ModelParams units internally.\n- **`fix_sm=true`**: Both moisture solvers accept this flag. When set, only `update_surface_water!` runs — soil θ is unchanged. Used when prescribing observed soil moisture.\n- **`fix_Tsoil=true`** in `UpdateHeatFlux`: skips the temperature update but still calls `Update_ice_ratio`. Always use this when prescribing observed temperatures.\n- **`ice_ratio` rescaling**: After any step that changes θ, `ice_ratio` is scaled by `θ_prev/θ` to conserve ice mass. Both solvers do this at the end of their outer function.\n- **Bottom boundary**: Zero-flux (no drainage out the bottom layer). This follows CLM 3.5 and is not configurable.","spac-canopy":"# SPAC & Canopy\n\n# SPAC & Canopy Module (`src/SPAC/`)\n\nThis module provides the canopy-physics layer for BEPS: solar geometry, LAI partitioning, conductance assembly, Vcmax-nitrogen coupling, and thermodynamic helpers. It supplies the building blocks consumed each timestep by `inter_prg_jl`.\n\n---\n\n## The `Leaf` Struct\n\n`Leaf <: FieldVector{4, Cdouble}` is the universal container for any per-leaf-class scalar. Four named fields cover the full two-layer × two-light-class decomposition used throughout the canopy:\n\n```\nLeaf\n├── o_sunlit   overstory, sunlit fraction\n├── o_shaded   overstory, shaded fraction\n├── u_sunlit   understory, sunlit fraction\n└── u_shaded   understory, shaded fraction\n```\n\nBecause `Leaf` is a `FieldVector`, it supports arithmetic operators element-wise without additional dispatch. Practical constructors:\n\n```julia\nLeaf()                    # zeros\nLeaf(x::Cdouble)          # broadcast scalar to all four slots\nLeaf(o::Cdouble, u)       # overstory value, understory value (sunlit = shaded per layer)\ninit_leaf_dbl2(x, ov, un) # in-place version of the two-argument constructor\n```\n\n`Leaf` is reused for LAI, PAI, conductances (Gw, Gc, Gs), absorbed radiation, and anything else that must vary by layer and light class. `multiply!(Z, X, Y)` does element-wise `Z = X .* Y` in-place.\n\n---\n\n## Solar Geometry\n\n```julia\nCosZs = s_coszs(jday, hour, lat, lon)\n```\n\nReturns `cos(solar zenith angle)` — negative at night, positive during the day. The calculation follows the standard declination approximation (Fourier series in day-of-year) plus an hour-angle derived from UTC time (`hour + lon/15`). Several downstream computations gate on `CosZs > 0` to suppress nighttime photosynthesis and LAI partitioning.\n\n---\n\n## LAI Partitioning (`lai2.jl`)\n\nSeparates total LAI and PAI into sunlit/shaded fractions for each layer using Beer's law with clumping index Ω:\n\n```julia\nLAI, PAI = lai2(Ω, CosZs, stem_o, stem_u, lai_o, lai_u)\n# or in-place:\nlai2!(Ω, CosZs, stem_o, stem_u, lai_o, lai_u, LAI::Leaf, PAI::Leaf)\n```\n\nThe sunlit fraction of overstory LAI:\n\n```\nLAI.o_sunlit = 2 · CosZs · (1 − exp(−0.5 · Ω · lai_o / CosZs))\n```\n\nUnderstory sunlit accounts for overstory shading by cumulating through the overstory first. Shaded fractions are the remainder (`lai − sunlit`), floored at zero. `partition_lai` is a simpler single-layer version used for quick checks.\n\n**Key invariant:** all sunlit fractions return 0 when `CosZs ≤ 0`.\n\n---\n\n## Vcmax from Leaf Nitrogen (`VCmax.jl`)\n\n```julia\nVCmax_sunlit, VCmax_shaded = VCmax(lai, Ω, CosZs, VCmax25, N_leaf, χ)\n```\n\nScales `VCmax25` by a nitrogen-weighted, light-weighted integral over canopy depth, following Chen et al. (2012, GBC). The extinction coefficient for nitrogen `Kn = 0.3` and the light extinction `K = 0.5·Ω/CosZs` are used to integrate analytically. Returns `0, 0` immediately at night (`CosZs ≤ 0`).\n\n---\n\n## Conductance Assembly (`BEPS_helper.jl`)\n\nTwo functions in `BEPS_helper.jl` combine aerodynamic, boundary-layer, and stomatal resistances into total leaf conductances for water vapour and CO₂:\n\n```julia\nupdate_Gw!(Gw, Gs_new, Ga_o, Ga_u, Gb_o, Gb_u)\nupdate_Gc!(Gc, Gs_new, Ga_o, Ga_u, Gb_o, Gb_u)\n```\n\nBoth write into a pre-allocated `Leaf`. The resistance formulas differ by diffusivity ratio:\n\n| Quantity | Boundary-layer factor | Stomatal factor |\n|---|---|---|\n| Water vapour (Gw) | `1/Gb` | `1/Gs` |\n| CO₂ (Gc) | `1.4/Gb` | `1.6/Gs` |\n\nAerodynamic resistance `1/Ga` is identical for both. These are called directly from `inter_prg_jl` after each stomatal conductance update.\n\n---\n\n## Thermodynamic Helpers (`helper.jl`)\n\nScalar functions exported for use across the model:\n\n| Function | Returns | Units |\n|---|---|---|\n| `cal_es(Ta)` | Saturation vapour pressure | kPa |\n| `cal_ea(Ta, RH)` | Actual vapour pressure | kPa |\n| `cal_slope(Ta)` | `d(es)/dT` | kPa K⁻¹ |\n| `cal_lambda(Ta)` | Latent heat of vaporisation | J kg⁻¹ |\n| `cal_cp(q)` | Specific heat of moist air | J kg⁻¹ K⁻¹ |\n| `ea2q(ea)` | Specific humidity from `ea` | kg kg⁻¹ |\n| `RH2q(Ta, RH)` | Specific humidity from RH | kg kg⁻¹ |\n| `q2RH(q, Ta)` | Relative humidity from `q` | % |\n| `blackbody(T)` | Blackbody emission | W m⁻² |\n| `cal_Rln(emiss, T)` | Longwave emission | W m⁻² |\n\n`meteo_pack_jl(Ta, RH)` bundles these into a named tuple `(ρₐ, cp, VPD, λ, Δ, γ, es, ea, q)` — a convenient way to compute all atmospheric scalars in one call.\n\n---\n\n## Snow Density (`snow_density.jl`)\n\n```julia\nρ_snow = snow_density(Ta, U10=NaN; tfrz=0.0, method=\"LoTmpDnsSlater2017\")\n```\n\nReturns fresh-snow density [kg m⁻³] as a function of air temperature and 10 m wind speed. Three temperature regimes:\n\n- `Ta > tfrz + 2` → wet/warm snow, fixed 50 + 1.7·17^1.5 ≈ 191 kg m⁻³  \n- `tfrz − 15 < Ta ≤ tfrz + 2` → polynomial in `(Ta − tfrz + 15)^1.5`  \n- `Ta ≤ tfrz − 15` → either truncated at 50 kg m⁻³ (Anderson 1976) or quadratic (Slater 2017, default)\n\nWind compaction adds up to ~267 kg m⁻³ via a `tanh`-based saturating function. Reference: van Kampenhout et al. 2017 / CLM5.\n\n---\n\n## How It Connects to `inter_prg_jl`\n\n```mermaid\ngraph TD\n    A[inter_prg_jl] -->|\"s_coszs()\"| B[Solar geometry]\n    A -->|\"lai2!()\"| C[LAI partitioning]\n    C --> D[Leaf structs: LAI, PAI]\n    A -->|\"VCmax()\"| E[Vcmax per leaf class]\n    A -->|\"update_Gw! / update_Gc!\"| F[Conductance assembly]\n    D --> F\n    F --> G[photosynthesis / ET]\n```\n\nAt each timestep in `inter_prg_jl`:\n\n1. `s_coszs` sets `CosZs` for all light-dependent calculations.\n2. `lai2!` fills `LAI` and `PAI` `Leaf` structs.\n3. `VCmax` derives per-class Vcmax from leaf nitrogen and the partitioned canopy.\n4. After stomatal conductance iteration, `update_Gw!` and `update_Gc!` assemble total leaf-to-atmosphere conductances used by photosynthesis and ET routines.\n\nThe thermodynamic helpers (`cal_es`, `cal_slope`, `cal_lambda`, etc.) are called both here and in `netRadiation.jl` independently.","surface-energy-fluxes":"# Surface Energy Fluxes\n\n# Surface Energy Fluxes\n\nThis module implements the coupled energy balance at the land surface: radiation partitioning across three canopy layers, aerodynamic transport, leaf temperature equilibration, and evaporation/sensible heat fluxes from leaves, soil, and snow. Together these files form the heart of the SPAC physics that runs inside `inter_prg_jl` every hour.\n\n---\n\n## Canopy structure and notation\n\nAll flux computations use a three-layer geometry: **overstory** (o), **understory** (u), and **ground** (g). Leaf-level quantities further split each vegetated layer into sunlit and shaded fractions, giving four leaf classes stored in the `Leaf` struct: `o_sunlit`, `o_shaded`, `u_sunlit`, `u_shaded`. The `Layer2{FT}` and `Layer3{FT}` types hold two or three scalar layer values, and `Radiation` accumulates intermediate SW/LW terms computed during `netRadiation_SW!`.\n\n---\n\n## Hourly time-step structure\n\n```mermaid\ngraph TD\n    A[Forcing: Rs, Rln, Ta, wind, RH] --> B[aerodynamic_conductance_jl]\n    A --> C[netRadiation_SW!]\n    B --> D[Ga, Gb for o/u/g]\n    C --> E[Rns_o, Rns_u, Rns_g + per-leaf Rns]\n    D & E --> F{Iteration over Tc}\n    F --> G[netRadiation_LW! → Rn per leaf]\n    G --> H[Leaf_Temperatures_jl → Tc converged]\n    H --> F\n    H --> I[sensible_heat_jl / transpiration_jl]\n    H --> J[evaporation_canopy_jl]\n    H --> K[surface_temperature! → T_surf, G]\n    I & J & K --> L[Flux output]\n```\n\nSW radiation (`netRadiation_SW!`) does not depend on temperature, so it runs once. LW radiation (`netRadiation_LW!`) depends on leaf temperatures and is recomputed each iteration until `Tc` converges.\n\n---\n\n## Aerodynamic conductance\n\n**Files:** `aerodynamic_conductance.jl` (V1), `aerodynamic_conductance_V2.jl` (V2)\n\nBoth versions expose the same signature and return the same seven outputs:\n\n```julia\nra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =\n    aerodynamic_conductance_jl(h, h_u, z_wind, clumping, Tair, u, H, lai_o, lai_u)\n```\n\n| Output | Meaning |\n|--------|---------|\n| `ra_o` | Aerodynamic resistance, reference height → overstory [s m⁻¹] |\n| `ra_u` | Resistance, overstory canopy base → understory top |\n| `ra_g` | Cumulative resistance, reference height → ground (includes `ra_o + ra_u`) |\n| `Ga_o`, `Gb_o` | Aerodynamic and leaf boundary conductances, overstory [m s⁻¹] |\n| `Ga_u`, `Gb_u` | Same for understory |\n\n**V1** uses a simplified Monin-Obukhov stability term: `ra_o = 1/(k·u*) · (ln((z-d)/z0) + n·(z-d)·L)`, where `L` is the Obukhov length derived from sensible heat flux `SH_o_p`.\n\n**V2** is physically more complete. It applies the Bonan (2019) stability correction functions `cal_Ψh` and `cal_ϕh` and separates momentum (`z0m`) and heat (`z0h = 0.1·z0m`) roughness lengths. The eddy diffusivity `kh_o` is handed off layer by layer — it decays from overstory top through the trunk space to the understory top using separate extinction coefficients (`gamma_o_m`, `gamma_o_h`, `gamma_u_m`), so each segment has a physically distinct turbulence regime.\n\n**Leaf boundary layer resistance** is computed by `ra_leaf_boundary(Tair, u)` in V2 (inline in V1 via `cal_Nu`): Reynolds number from leaf-width scale (0.1 m), yielding `rb = 0.5 · 0.1 / (α_w · Nu)` clamped to 40 s m⁻¹.\n\nFallback values (`Ga = 1/200`, `ra_g = 300`) are returned when `u == 0`.\n\n---\n\n## Net radiation\n\n**File:** `netRadiation.jl`\n\nEntry point: `netRadiation_jl` calls `netRadiation_SW!` then `netRadiation_LW!`.\n\n### Shortwave (`netRadiation_SW!`)\n\n1. Computes effective snow-weighted albedo for each layer from the bare and snow albedos and snow fraction.\n2. Separates global radiation into direct and diffuse using the cloud index `ratio_cloud = Rs / (1367·CosZs)` (Luo 2018 Eq. A4).\n3. Exponential Beer-Lambert transmittances for direct (`τ_dir = exp(-0.5·Ω·LAI/CosZs)`) and diffuse (`τ_df = exp(-0.5·Ω·LAI/cosQ)`) separately for each layer.\n4. Layer net SW = incoming that is absorbed at that layer after subtracting what passes through to the layer below.\n5. Per-leaf-class irradiance: sunlit leaves receive direct + diffuse; shaded leaves receive diffuse only.\n\n### Longwave (`netRadiation_LW!`, `cal_Rln_Longwave`)\n\nStefan–Boltzmann emission `cal_Rln(ϵ, T) = ϵ·σ·(T+273.15)⁴`. Atmospheric emissivity follows `ϵ_air = 1 - exp(-(ea·10)^((T+273)/1200))`. The three-layer exchange accounts for multiple reflections between layers using layer transmittances `τ_o_df` and `τ_u_df`. Canopy emissivities: `ϵ_o = ϵ_u = 0.98`, `ϵ_g = 0.96`.\n\nPer-leaf-class longwave is distributed by `Rnl_o / lai_os` (area-weighted from the layer total).\n\n---\n\n## Leaf temperature and heat fluxes\n\n**File:** `heat_H_and_LE.jl`\n\n### Leaf temperature (`Leaf_Temperature_jl`)\n\nSingle-leaf energy balance solved analytically:\n\n```\nTc = Ta + (Rn - VPD·ρa·cp·p*) / (ρa·cp·(Gh + Δ·p*))\n```\n\nwhere `p* = (Gw + Gw_wet·Xc) / γ` combines the dry-leaf and wet-leaf water conductances. The result is clamped to `[Ta - 3, Ta + 5]` by default.\n\n`Leaf_Temperatures_jl` calls this for all four leaf classes in sequence, mutating the `Tc::Leaf` struct in place.\n\n### Sensible heat\n\nTwo methods: scalar form `sensible_heat(T_w, T_a, ρa, cp, gH)` and leaf-struct form that returns a `Leaf`. `sensible_heat_jl` integrates over the leaf-area distribution to produce `SH_o`, `SH_u`, `SH_g`.\n\n### Latent heat / transpiration\n\n`latent_heat(Ta, Ts, gw, met)` uses the Penman linearisation:\n\n```\nLE = ρa·cp·gw·(VPD + Δ·(Ts - Ta)) / γ\n```\n\n`transpiration_jl` integrates this over all four leaf classes and divides by `λ` (latent heat of vaporisation) to return `trans_o` and `trans_u` in kg m⁻² s⁻¹.\n\n`latent_heat!` mutates a `Leaf` struct in place for use inside the inter-program iteration loop (`inter_prg_jl`).\n\n---\n\n## Ground surface temperature\n\n**File:** `surface_temperature.jl`\n\n`surface_temperature!` is the public entry point. It first updates soil thermal conductivity and heat capacity (`UpdateThermal_κ`, `UpdateThermal_Cv`), then dispatches to `surface_temperature_jl!` which handles three snow regimes:\n\n| Case | Snow depth | Approach |\n|------|-----------|----------|\n| 1 | ≤ 2 cm | Bare/thin snow — single implicit solve for `T_surf` |\n| 2 | 2–5 cm | Mixed snow-soil — separate bare and snow surface temperatures, area-weighted `T_mix0` |\n| 3 | > 5 cm | Deep snow — three-layer snow model (top 2 cm / next 2 cm / remainder) |\n\nThe core solver is `solve_imp`, an implicit Euler step for the 1D thermal diffusion equation:\n\n```\nT_new = (T_old·I + G·ra·z_rad + ρCp·T_bnd·z + η_c·ra·κ_bot·T_bot) /\n        (ρCp·z + η_c·ra·κ_bot + I)\n```\n\nwhere `I = ΔM·ra·z` is the thermal inertia term. The explicit helper `step_exp` is used for interior snow layers in Case 3.\n\nPhase-change checks (`check_phase`) clamp temperature to 0 °C when crossing the freezing point and liquid water is present.\n\nSnow thermal conductivity follows Jordan (1991): `κ_snow = 0.021 + 4.2e-4·ρ + 2.2e-9·ρ³`.\n\n---\n\n## Evaporation\n\n### Canopy interception (`evaporation_canopy.jl`)\n\n`evaporation_canopy_jl` computes evaporation of intercepted liquid water and sublimation of intercepted snow for each of the four leaf classes. Both use `latent_heat` from `heat_H_and_LE.jl`; the only difference is the denominator: liquid uses `λ` (vaporisation), snow uses `λ_snow` (sublimation). Area integration over `lai` produces `E_water_o`, `E_water_u`, `E_snow_o`, `E_snow_u`.\n\n### Soil/ground (`evaporation_soil.jl`)\n\n`evaporation_soil_jl` handles three mutually exclusive ground states in priority order:\n\n1. **Standing water, no snow:** Penman-Monteith with soil resistance `Gwater_g = 0` (open water, rs → 0).\n2. **Snow present:** Sublimation proportional to snow coverage fraction `perc_snow.g`.\n3. **Bare moist soil:** PM with soil surface resistance `Gwater_g = 1/(4·exp(8.2 - 4.2·θ/n))` where `θ/n` is relative saturation.\n\nBoth water depth (`z_water`) and snow mass (`mass_snow.g`) are updated in-place to conserve mass within the time step. Evaporation is bounded above by available water/snow storage.\n\n---\n\n## Connection to `inter_prg_jl`\n\n`inter_prg_jl` calls into this module in the following order each hour:\n\n1. `aerodynamic_conductance_jl` — compute resistances from canopy geometry and stability.\n2. `netRadiation_SW!` — partition solar radiation (once, before the iteration).\n3. Iterative loop (~5 iterations):\n   - `netRadiation_LW!` — update LW given current `Tc`.\n   - `Leaf_Temperatures_jl` — solve new `Tc`.\n   - `latent_heat!` — update LE for convergence check.\n4. `sensible_heat_jl`, `transpiration_jl` — final flux outputs.\n5. `evaporation_canopy_jl` — intercepted water/snow evaporation.\n6. `evaporation_soil_jl` — ground evaporation.\n7. `surface_temperature!` — advance `T_surf` and compute soil heat flux `G`.","utilities":"# Utilities\n\n# Utilities Module\n\nTwo source files — `src/ultilize.jl` and `src/check_forcing.jl` — provide the data preparation and aggregation infrastructure that sits between raw inputs and the model core. Neither file contains physics; together they handle the unglamorous work of cleaning forcing data, interpolating gaps, and collapsing hourly outputs to daily summaries.\n\n---\n\n## Temporal Aggregation (`src/ultilize.jl`)\n\n### `parse_time`\n\n```julia\nparse_time(x::AbstractString) -> DateTime\n```\n\nParses ISO 8601 strings with a trailing `Z` (`\"2010-06-15T14:00:00Z\"`) into `DateTime`. Used wherever forcing files carry UTC timestamps as strings.\n\n---\n\n### `agg` — generic grouped reduction\n\n```julia\nagg(A::AbstractArray, by; dims::Integer=3, fun=mean) -> Array{Float64}\n```\n\nReduces an N-dimensional array along one axis by grouping indices according to `by`. The output array has the same rank as `A`; the target dimension shrinks to `length(unique(by))`.\n\n**How it works:**\n\n1. Sorts the unique group labels to fix output ordering.\n2. Builds an index map `Dict{group → [row indices]}` in a single pass over `by`.\n3. For each group, calls `fun(selectdim(A, dims, inds[g]); dims=dims)` and writes the result into the corresponding slice of `out`.\n\n`@views` and `@inbounds` are applied throughout. The `selectdim` + `dropdims` pattern avoids materializing intermediate arrays.\n\n`fun` defaults to `mean` but can be any function with the signature `f(x; dims)` — pass `sum` for precipitation totals, for example.\n\n---\n\n### `agg_daily` — hourly → daily\n\nTwo methods cover the two common call sites:\n\n```julia\n# 1. For a raw matrix (rows = timesteps)\nagg_daily(mat::AbstractMatrix, dates) -> Matrix\n\n# 2. For a flux DataFrame — returns a NamedTuple of daily vectors\nagg_daily(df_fluxes::AbstractDataFrame, dates) -> NamedTuple\n```\n\nThe `DataFrame` overload hard-codes the BEPS flux variables:\n\n| Output field | Source columns | Reduction |\n|---|---|---|\n| `GPP_sim` | `df_fluxes.GPP` | `sum` (gC m⁻² d⁻¹) |\n| `ET_sim` | `Trans + Evap` | `sum` (mm d⁻¹) |\n| `Hs_sim` | `SH` | `mean` (W m⁻²) |\n\n`dates_day` is the sorted unique set of `Date.(dates)`, returned alongside the aggregated vectors so the caller always has aligned date and value arrays.\n\n---\n\n## Forcing Quality Control (`src/check_forcing.jl`)\n\n### Pipeline overview\n\n```mermaid\ngraph TD\n    A[Raw forcing DataFrame] --> B[sanitize_forcing!]\n    B --> C[_sanitize_column! ×5]\n    C --> D{bad values?}\n    D -- yes --> E[set to NaN or fill]\n    E --> F[_interp_invalid!]\n    F --> G[clamp to physical range]\n    D -- no --> G\n    G --> H[Clean forcing]\n```\n\n---\n\n### `_interp_invalid!` (internal)\n\n```julia\n_interp_invalid!(x::AbstractVector{<:AbstractFloat}) -> x\n```\n\nIn-place linear interpolation over non-finite gaps. Edge behaviour:\n\n- Leading gap → backfilled with the first finite value.\n- Trailing gap → forward-filled with the last finite value.\n- Interior gaps → linearly interpolated between the nearest finite neighbours.\n\nA gap of length one (`j == i + 1`) is skipped since there is nothing to interpolate. The function modifies `x` in place and returns it.\n\n---\n\n### `_sanitize_column!` (internal)\n\n```julia\n_sanitize_column!(d, col; lo, hi, fill, interp, clamp_lo, clamp_hi) -> Int\n```\n\nOperates on a single named column of a DataFrame. Returns the count of bad values replaced. Steps:\n\n1. Mark bad: `!isfinite(x) || x < lo || x > hi`.\n2. Replace bad values with `fill` (default `NaN`).\n3. Call `_interp_invalid!` unless `interp=false`.\n4. Optionally clamp to `[clamp_lo, clamp_hi]`.\n\nSilently no-ops if `col` is absent from the DataFrame (returns `0`).\n\n---\n\n### `sanitize_forcing!`\n\n```julia\nsanitize_forcing!(d::AbstractDataFrame) -> NamedTuple\n```\n\nApplies `_sanitize_column!` to all five standard meteorological columns:\n\n| Column | Physical range | Fill strategy | Clamp |\n|---|---|---|---|\n| `Tair` | −50 to 50 °C | interpolate | — |\n| `RH` | 0 to 100 % | interpolate | [0, 100] |\n| `Uz` | ≥ 0 m s⁻¹ | interpolate | ≥ 0.01 (avoids division by zero) |\n| `Rs` | ≥ 0 W m⁻² | interpolate | ≥ 0 |\n| `Prcp` | ≥ 0 mm | `0.0`, **no interpolation** | ≥ 0 |\n\nPrecipitation uses `fill=0.0` and `interp=false` — a missing precipitation record should be treated as no rain, not interpolated from neighbours. Wind speed is clamped to a minimum of 0.01 m s⁻¹ to prevent downstream aerodynamic resistance from blowing up.\n\nReturns a `NamedTuple` `(; bad_Tair, bad_RH, bad_Uz, bad_Rs, bad_Prcp)` with the count of bad values per column — useful for logging data quality.\n\n---\n\n### `normalize_flux_obs!`\n\n```julia\nnormalize_flux_obs!(d::AbstractDataFrame) -> d\n```\n\nSome eddy-covariance datasets store GPP as a negative flux (net ecosystem exchange sign convention). If the column `GPP_obs` exists and its mean is negative, this function negates it in place. Called in `LoadData` before any calibration or comparison against simulated GPP.\n\n---\n\n### `_align_daily_data` and `_daily_obs_dates` (internal)\n\nThese two helpers serve `BEPS_GOF` in `src/Optim.jl` during calibration:\n\n```julia\n_daily_obs_dates(d) -> Vector{Date} | nothing   # extracts Date.(d.date); nil-safe\n_align_daily_data(data_sim, data_obs) -> (sim′, obs′)  # inner join on :date\n```\n\n`_align_daily_data` computes the intersection of simulation and observation dates, then reindexes both DataFrames to the common sorted date vector. It errors if the intersection is empty, which catches mismatched time periods early rather than silently producing an empty goodness-of-fit.\n\n---\n\n## Integration Points\n\n`sanitize_forcing!` and `normalize_flux_obs!` are called once per site run inside `LoadData` (`Project_ChinaFlux/case01_ChinaFlux.jl`), before the forcing array is ever passed to `simulate`. `_align_daily_data` and `_daily_obs_dates` are used exclusively inside `BEPS_GOF` during parameter optimisation. `agg_daily` is the standard post-processing step for converting hourly `simulate` output to daily totals for GOF computation or plotting."};
 var TREE = [{"name":"Model Core & API","slug":"model-core-api","files":["src/BEPS.jl","src/BEPS_modules.jl","src/beps_main.jl","src/beps_modern.jl","src/inter_prg.jl"]},{"name":"Data Types & State","slug":"data-types-state","files":["src/DataType/DataType.jl","src/DataType/BEPS_Param.jl","src/DataType/BEPS_State.jl","src/DataType/CanopyLayer.jl","src/DataType/Constant.jl","src/DataType/LeafCache.jl","src/DataType/Met.jl","src/DataType/OUTPUT.jl","src/DataType/PhotoConsts.jl","src/DataType/AeroConsts.jl","src/DataType/StateSeries.jl","src/DataType/macro.jl","src/DataType/setup.jl"]},{"name":"Model Parameters","slug":"model-parameters","files":["src/DataType/Params/Params.jl","src/DataType/Params/GlobalData.jl","src/DataType/Params/ParamPhoto.jl","src/DataType/Params/Param_Init.jl","src/DataType/Params/deprecated/Init_Soil_Parameters.jl","src/DataType/Params/deprecated/ReadParamVeg.jl","src/DataType/Params/deprecated/VegHelper.jl"]},{"name":"Soil Physics","slug":"soil-physics","files":["src/SoilPhysics/SoilPhysics.jl","src/SoilPhysics/Solve_SM.jl","src/SoilPhysics/UpdateHeatFlux.jl","src/SoilPhysics/UpdateSoilMoisture.jl","src/SoilPhysics/soil_water_factor_v2.jl","src/standalone/UpdateSoilMoisture.jl"]},{"name":"SPAC & Canopy","slug":"spac-canopy","files":["src/SPAC/SPAC.jl","src/SPAC/BEPS_helper.jl","src/SPAC/Leaf.jl","src/SPAC/VCmax.jl","src/SPAC/helper.jl","src/SPAC/lai2.jl","src/SPAC/snow_density.jl","src/SPAC/ultilize.jl"]},{"name":"Surface Energy Fluxes","slug":"surface-energy-fluxes","files":["src/aerodynamic_conductance.jl","src/aerodynamic_conductance_V2.jl","src/heat_H_and_LE.jl","src/netRadiation.jl","src/surface_temperature.jl","src/evaporation_canopy.jl","src/evaporation_soil.jl"]},{"name":"Photosynthesis","slug":"photosynthesis","files":["src/photosynthesis.jl","src/photosynthesis_helper.jl","src/standalone/Photosynthesis/core.jl","src/standalone/Photosynthesis/helper.jl","src/standalone/Photosynthesis/photosynthesis.jl","src/standalone/Photosynthesis/radiation.jl","src/standalone/Photosynthesis/stomatal.jl","src/standalone/Photosynthesis/temperature.jl","src/standalone/Photosynthesis/types.jl"]},{"name":"Hydrology & Snowpack","slug":"hydrology-snowpack","files":["src/rainfall_stage.jl","src/snowpack.jl"]},{"name":"C Library Bindings","slug":"c-library-bindings","files":["src/clang/BEPS_c.jl","src/clang/SOIL_c.jl","src/clang/module.jl","src/clang/snowpack_stage.jl","src/clang/struct_SOIL.jl"]},{"name":"Optimization & Calibration","slug":"optimization-calibration","files":["src/Optim.jl"]},{"name":"Utilities","slug":"utilities","files":["src/ultilize.jl","src/check_forcing.jl"]},{"name":"ChinaFlux Case Studies","slug":"chinaflux-case-studies","files":["Project_ChinaFlux/case01_ChinaFlux.jl","Project_ChinaFlux/main_vis.jl","Project_ChinaFlux/Figure1_GOF.Rmd","Project_ChinaFlux/Table1_GOF.R","Project_ChinaFlux/check_Input.Rmd","Project_ChinaFlux/ex_model.ipynb","Project_ChinaFlux/template.jl"]},{"name":"Other","slug":"other","files":[],"children":[{"name":"Other — AGENTS.md","slug":"other-agents-md","files":["AGENTS.md"]},{"name":"Other — Artifacts.toml","slug":"other-artifacts-toml","files":["Artifacts.toml"]},{"name":"Other — CLAUDE.md","slug":"other-claude-md","files":["CLAUDE.md"]},{"name":"Other — Project.toml","slug":"other-project-toml","files":["Project.toml"]},{"name":"Other — README.md","slug":"other-readme-md","files":["README.md"]},{"name":"Other — case01.ipynb","slug":"other-case01-ipynb","files":["case01.ipynb"]},{"name":"Other — codecov.yml","slug":"other-codecov-yml","files":["codecov.yml"]},{"name":"Other — deps","slug":"other-deps","files":["deps/build.ps1","deps/gen_artifacts.jl"]},{"name":"Other — 01_参数化方案","slug":"other-01","files":["docs/01_参数化方案/ch01_params.typ"]},{"name":"Other — docs","slug":"other-docs","files":["docs/README.md","docs/Structs.ipynb"]},{"name":"Other — ch01_SoilMoisture","slug":"other-ch01-soilmoisture","files":["docs/ch01_SoilMoisture/Soil_struct.md","docs/ch01_SoilMoisture/ch01_SoilMoisture.typ","docs/ch01_SoilMoisture/soil_hydrology.md","docs/ch01_SoilMoisture/wilting_point.md","docs/ch01_SoilMoisture/推导_界面Kmid加权.typ"]},{"name":"Other — ch02_SoilTemperature","slug":"other-ch02-soiltemperature","files":["docs/ch02_SoilTemperature/TS_Case02.md","docs/ch02_SoilTemperature/TS_Case03.md","docs/ch02_SoilTemperature/ch02_SoilTemperature.typ","docs/ch02_SoilTemperature/surface_temperature.md"]},{"name":"Other — ch03_SurfaceFluxes","slug":"other-ch03-surfacefluxes","files":["docs/ch03_SurfaceFluxes/V2.typ","docs/ch03_SurfaceFluxes/aerodynamic_conductance.typ","docs/ch03_SurfaceFluxes/forcing_reference_height.md"]},{"name":"Other — ch04_Photosynthesis","slug":"other-ch04-photosynthesis","files":["docs/ch04_Photosynthesis/photosynthesis.typ","docs/ch04_Photosynthesis/photosynthesis_comparison.typ"]},{"name":"Other — manual","slug":"other-manual","files":["docs/manual/ch01_Input.typ","docs/manual/ch02_Param_SoilVeg_Types.typ","docs/manual/manual.typ","docs/manual/软著/ModernBEPS.typ"]},{"name":"Other — 结题","slug":"other","files":["docs/结题/ModernBEPS_精度证明.typ","docs/结题/数据方面.typ","docs/结题/结题_PART2.typ"]},{"name":"Other — path_mnt","slug":"other-path-mnt","files":["path_mnt"]},{"name":"Other — plan.md","slug":"other-plan-md","files":["plan.md"]},{"name":"Other — repomix.ps1","slug":"other-repomix-ps1","files":["repomix.ps1"]},{"name":"Other — debug","slug":"other-debug","files":["test/debug/debug-Rn.qmd","test/debug/debug.md","test/debug/debug_Rln.jl","test/debug/debug_Rln2.jl","test/debug/debug_snow.jl"]},{"name":"Other — modules","slug":"other-modules","files":["test/modules/modules.jl","test/modules/test-Soil.jl","test/modules/test-aerodynamic_conductance.jl","test/modules/test-param.jl","test/modules/test-photosynthesis.jl","test/modules/test-radiation.jl","test/modules/test-rainfall_stage1.jl","test/modules/test-sensible_heat.jl","test/modules/test-setup.jl","test/modules/test-snowpack.jl","test/modules/test-surface_temperature.jl"]},{"name":"Other — test","slug":"other-test","files":["test/runtests.jl","test/test-BEPSmodel.jl","test/test-ModelParams.jl","test/test-StateSeries.jl","test/test-UpdateSoilMoisture_Q0.jl","test/test-beps_main.jl","test/test-beps_modern.jl","test/test-macro.jl","test/test-performance.jl","test/test-photosynthesis_standalone.jl","test/test-soil_sm.jl","test/test-utilize.jl"]}]}];
 var META = {"fromCommit":"a8b7fbee4f05eb2400852618aba6ed02e5f344ef","generatedAt":"2026-05-29T07:46:52.466Z","model":"","lang":"","moduleFiles":{"Model Core & API":["src/BEPS.jl","src/BEPS_modules.jl","src/beps_main.jl","src/beps_modern.jl","src/inter_prg.jl"],"Data Types & State":["src/DataType/DataType.jl","src/DataType/BEPS_Param.jl","src/DataType/BEPS_State.jl","src/DataType/CanopyLayer.jl","src/DataType/Constant.jl","src/DataType/LeafCache.jl","src/DataType/Met.jl","src/DataType/OUTPUT.jl","src/DataType/PhotoConsts.jl","src/DataType/AeroConsts.jl","src/DataType/StateSeries.jl","src/DataType/macro.jl","src/DataType/setup.jl"],"Model Parameters":["src/DataType/Params/Params.jl","src/DataType/Params/GlobalData.jl","src/DataType/Params/ParamPhoto.jl","src/DataType/Params/Param_Init.jl","src/DataType/Params/deprecated/Init_Soil_Parameters.jl","src/DataType/Params/deprecated/ReadParamVeg.jl","src/DataType/Params/deprecated/VegHelper.jl"],"Soil Physics":["src/SoilPhysics/SoilPhysics.jl","src/SoilPhysics/Solve_SM.jl","src/SoilPhysics/UpdateHeatFlux.jl","src/SoilPhysics/UpdateSoilMoisture.jl","src/SoilPhysics/soil_water_factor_v2.jl","src/standalone/UpdateSoilMoisture.jl"],"SPAC & Canopy":["src/SPAC/SPAC.jl","src/SPAC/BEPS_helper.jl","src/SPAC/Leaf.jl","src/SPAC/VCmax.jl","src/SPAC/helper.jl","src/SPAC/lai2.jl","src/SPAC/snow_density.jl","src/SPAC/ultilize.jl"],"Surface Energy Fluxes":["src/aerodynamic_conductance.jl","src/aerodynamic_conductance_V2.jl","src/heat_H_and_LE.jl","src/netRadiation.jl","src/surface_temperature.jl","src/evaporation_canopy.jl","src/evaporation_soil.jl"],"Photosynthesis":["src/photosynthesis.jl","src/photosynthesis_helper.jl","src/standalone/Photosynthesis/core.jl","src/standalone/Photosynthesis/helper.jl","src/standalone/Photosynthesis/photosynthesis.jl","src/standalone/Photosynthesis/radiation.jl","src/standalone/Photosynthesis/stomatal.jl","src/standalone/Photosynthesis/temperature.jl","src/standalone/Photosynthesis/types.jl"],"Hydrology & Snowpack":["src/rainfall_stage.jl","src/snowpack.jl"],"C Library Bindings":["src/clang/BEPS_c.jl","src/clang/SOIL_c.jl","src/clang/module.jl","src/clang/snowpack_stage.jl","src/clang/struct_SOIL.jl"],"Optimization & Calibration":["src/Optim.jl"],"Utilities":["src/ultilize.jl","src/check_forcing.jl"],"ChinaFlux Case Studies":["Project_ChinaFlux/case01_ChinaFlux.jl","Project_ChinaFlux/main_vis.jl","Project_ChinaFlux/Figure1_GOF.Rmd","Project_ChinaFlux/Table1_GOF.R","Project_ChinaFlux/check_Input.Rmd","Project_ChinaFlux/ex_model.ipynb","Project_ChinaFlux/template.jl"],"Other":["AGENTS.md","Artifacts.toml","CLAUDE.md","Project.toml","README.md","case01.ipynb","codecov.yml","deps/build.ps1","deps/gen_artifacts.jl","docs/01_参数化方案/ch01_params.typ","docs/README.md","docs/Structs.ipynb","docs/ch01_SoilMoisture/Soil_struct.md","docs/ch01_SoilMoisture/ch01_SoilMoisture.typ","docs/ch01_SoilMoisture/soil_hydrology.md","docs/ch01_SoilMoisture/wilting_point.md","docs/ch01_SoilMoisture/推导_界面Kmid加权.typ","docs/ch02_SoilTemperature/TS_Case02.md","docs/ch02_SoilTemperature/TS_Case03.md","docs/ch02_SoilTemperature/ch02_SoilTemperature.typ","docs/ch02_SoilTemperature/surface_temperature.md","docs/ch03_SurfaceFluxes/V2.typ","docs/ch03_SurfaceFluxes/aerodynamic_conductance.typ","docs/ch03_SurfaceFluxes/forcing_reference_height.md","docs/ch04_Photosynthesis/photosynthesis.typ","docs/ch04_Photosynthesis/photosynthesis_comparison.typ","docs/manual/ch01_Input.typ","docs/manual/ch02_Param_SoilVeg_Types.typ","docs/manual/manual.typ","docs/manual/软著/ModernBEPS.typ","docs/结题/ModernBEPS_精度证明.typ","docs/结题/数据方面.typ","docs/结题/结题_PART2.typ","path_mnt","plan.md","repomix.ps1","test/debug/debug-Rn.qmd","test/debug/debug.md","test/debug/debug_Rln.jl","test/debug/debug_Rln2.jl","test/debug/debug_snow.jl","test/modules/modules.jl","test/modules/test-Soil.jl","test/modules/test-aerodynamic_conductance.jl","test/modules/test-param.jl","test/modules/test-photosynthesis.jl","test/modules/test-radiation.jl","test/modules/test-rainfall_stage1.jl","test/modules/test-sensible_heat.jl","test/modules/test-setup.jl","test/modules/test-snowpack.jl","test/modules/test-surface_temperature.jl","test/runtests.jl","test/test-BEPSmodel.jl","test/test-ModelParams.jl","test/test-StateSeries.jl","test/test-UpdateSoilMoisture_Q0.jl","test/test-beps_main.jl","test/test-beps_modern.jl","test/test-macro.jl","test/test-performance.jl","test/test-photosynthesis_standalone.jl","test/test-soil_sm.jl","test/test-utilize.jl"],"Other — AGENTS.md":["AGENTS.md"],"Other — Artifacts.toml":["Artifacts.toml"],"Other — CLAUDE.md":["CLAUDE.md"],"Other — Project.toml":["Project.toml"],"Other — README.md":["README.md"],"Other — case01.ipynb":["case01.ipynb"],"Other — codecov.yml":["codecov.yml"],"Other — deps":["deps/build.ps1","deps/gen_artifacts.jl"],"Other — 01_参数化方案":["docs/01_参数化方案/ch01_params.typ"],"Other — docs":["docs/README.md","docs/Structs.ipynb"],"Other — ch01_SoilMoisture":["docs/ch01_SoilMoisture/Soil_struct.md","docs/ch01_SoilMoisture/ch01_SoilMoisture.typ","docs/ch01_SoilMoisture/soil_hydrology.md","docs/ch01_SoilMoisture/wilting_point.md","docs/ch01_SoilMoisture/推导_界面Kmid加权.typ"],"Other — ch02_SoilTemperature":["docs/ch02_SoilTemperature/TS_Case02.md","docs/ch02_SoilTemperature/TS_Case03.md","docs/ch02_SoilTemperature/ch02_SoilTemperature.typ","docs/ch02_SoilTemperature/surface_temperature.md"],"Other — ch03_SurfaceFluxes":["docs/ch03_SurfaceFluxes/V2.typ","docs/ch03_SurfaceFluxes/aerodynamic_conductance.typ","docs/ch03_SurfaceFluxes/forcing_reference_height.md"],"Other — ch04_Photosynthesis":["docs/ch04_Photosynthesis/photosynthesis.typ","docs/ch04_Photosynthesis/photosynthesis_comparison.typ"],"Other — manual":["docs/manual/ch01_Input.typ","docs/manual/ch02_Param_SoilVeg_Types.typ","docs/manual/manual.typ","docs/manual/软著/ModernBEPS.typ"],"Other — 结题":["docs/结题/ModernBEPS_精度证明.typ","docs/结题/数据方面.typ","docs/结题/结题_PART2.typ"],"Other — path_mnt":["path_mnt"],"Other — plan.md":["plan.md"],"Other — repomix.ps1":["repomix.ps1"],"Other — debug":["test/debug/debug-Rn.qmd","test/debug/debug.md","test/debug/debug_Rln.jl","test/debug/debug_Rln2.jl","test/debug/debug_snow.jl"],"Other — modules":["test/modules/modules.jl","test/modules/test-Soil.jl","test/modules/test-aerodynamic_conductance.jl","test/modules/test-param.jl","test/modules/test-photosynthesis.jl","test/modules/test-radiation.jl","test/modules/test-rainfall_stage1.jl","test/modules/test-sensible_heat.jl","test/modules/test-setup.jl","test/modules/test-snowpack.jl","test/modules/test-surface_temperature.jl"],"Other — test":["test/runtests.jl","test/test-BEPSmodel.jl","test/test-ModelParams.jl","test/test-StateSeries.jl","test/test-UpdateSoilMoisture_Q0.jl","test/test-beps_main.jl","test/test-beps_modern.jl","test/test-macro.jl","test/test-performance.jl","test/test-photosynthesis_standalone.jl","test/test-soil_sm.jl","test/test-utilize.jl"]},"moduleTree":[{"name":"Model Core & API","slug":"model-core-api","files":["src/BEPS.jl","src/BEPS_modules.jl","src/beps_main.jl","src/beps_modern.jl","src/inter_prg.jl"]},{"name":"Data Types & State","slug":"data-types-state","files":["src/DataType/DataType.jl","src/DataType/BEPS_Param.jl","src/DataType/BEPS_State.jl","src/DataType/CanopyLayer.jl","src/DataType/Constant.jl","src/DataType/LeafCache.jl","src/DataType/Met.jl","src/DataType/OUTPUT.jl","src/DataType/PhotoConsts.jl","src/DataType/AeroConsts.jl","src/DataType/StateSeries.jl","src/DataType/macro.jl","src/DataType/setup.jl"]},{"name":"Model Parameters","slug":"model-parameters","files":["src/DataType/Params/Params.jl","src/DataType/Params/GlobalData.jl","src/DataType/Params/ParamPhoto.jl","src/DataType/Params/Param_Init.jl","src/DataType/Params/deprecated/Init_Soil_Parameters.jl","src/DataType/Params/deprecated/ReadParamVeg.jl","src/DataType/Params/deprecated/VegHelper.jl"]},{"name":"Soil Physics","slug":"soil-physics","files":["src/SoilPhysics/SoilPhysics.jl","src/SoilPhysics/Solve_SM.jl","src/SoilPhysics/UpdateHeatFlux.jl","src/SoilPhysics/UpdateSoilMoisture.jl","src/SoilPhysics/soil_water_factor_v2.jl","src/standalone/UpdateSoilMoisture.jl"]},{"name":"SPAC & Canopy","slug":"spac-canopy","files":["src/SPAC/SPAC.jl","src/SPAC/BEPS_helper.jl","src/SPAC/Leaf.jl","src/SPAC/VCmax.jl","src/SPAC/helper.jl","src/SPAC/lai2.jl","src/SPAC/snow_density.jl","src/SPAC/ultilize.jl"]},{"name":"Surface Energy Fluxes","slug":"surface-energy-fluxes","files":["src/aerodynamic_conductance.jl","src/aerodynamic_conductance_V2.jl","src/heat_H_and_LE.jl","src/netRadiation.jl","src/surface_temperature.jl","src/evaporation_canopy.jl","src/evaporation_soil.jl"]},{"name":"Photosynthesis","slug":"photosynthesis","files":["src/photosynthesis.jl","src/photosynthesis_helper.jl","src/standalone/Photosynthesis/core.jl","src/standalone/Photosynthesis/helper.jl","src/standalone/Photosynthesis/photosynthesis.jl","src/standalone/Photosynthesis/radiation.jl","src/standalone/Photosynthesis/stomatal.jl","src/standalone/Photosynthesis/temperature.jl","src/standalone/Photosynthesis/types.jl"]},{"name":"Hydrology & Snowpack","slug":"hydrology-snowpack","files":["src/rainfall_stage.jl","src/snowpack.jl"]},{"name":"C Library Bindings","slug":"c-library-bindings","files":["src/clang/BEPS_c.jl","src/clang/SOIL_c.jl","src/clang/module.jl","src/clang/snowpack_stage.jl","src/clang/struct_SOIL.jl"]},{"name":"Optimization & Calibration","slug":"optimization-calibration","files":["src/Optim.jl"]},{"name":"Utilities","slug":"utilities","files":["src/ultilize.jl","src/check_forcing.jl"]},{"name":"ChinaFlux Case Studies","slug":"chinaflux-case-studies","files":["Project_ChinaFlux/case01_ChinaFlux.jl","Project_ChinaFlux/main_vis.jl","Project_ChinaFlux/Figure1_GOF.Rmd","Project_ChinaFlux/Table1_GOF.R","Project_ChinaFlux/check_Input.Rmd","Project_ChinaFlux/ex_model.ipynb","Project_ChinaFlux/template.jl"]},{"name":"Other","slug":"other","files":[],"children":[{"name":"Other — AGENTS.md","slug":"other-agents-md","files":["AGENTS.md"]},{"name":"Other — Artifacts.toml","slug":"other-artifacts-toml","files":["Artifacts.toml"]},{"name":"Other — CLAUDE.md","slug":"other-claude-md","files":["CLAUDE.md"]},{"name":"Other — Project.toml","slug":"other-project-toml","files":["Project.toml"]},{"name":"Other — README.md","slug":"other-readme-md","files":["README.md"]},{"name":"Other — case01.ipynb","slug":"other-case01-ipynb","files":["case01.ipynb"]},{"name":"Other — codecov.yml","slug":"other-codecov-yml","files":["codecov.yml"]},{"name":"Other — deps","slug":"other-deps","files":["deps/build.ps1","deps/gen_artifacts.jl"]},{"name":"Other — 01_参数化方案","slug":"other-01","files":["docs/01_参数化方案/ch01_params.typ"]},{"name":"Other — docs","slug":"other-docs","files":["docs/README.md","docs/Structs.ipynb"]},{"name":"Other — ch01_SoilMoisture","slug":"other-ch01-soilmoisture","files":["docs/ch01_SoilMoisture/Soil_struct.md","docs/ch01_SoilMoisture/ch01_SoilMoisture.typ","docs/ch01_SoilMoisture/soil_hydrology.md","docs/ch01_SoilMoisture/wilting_point.md","docs/ch01_SoilMoisture/推导_界面Kmid加权.typ"]},{"name":"Other — ch02_SoilTemperature","slug":"other-ch02-soiltemperature","files":["docs/ch02_SoilTemperature/TS_Case02.md","docs/ch02_SoilTemperature/TS_Case03.md","docs/ch02_SoilTemperature/ch02_SoilTemperature.typ","docs/ch02_SoilTemperature/surface_temperature.md"]},{"name":"Other — ch03_SurfaceFluxes","slug":"other-ch03-surfacefluxes","files":["docs/ch03_SurfaceFluxes/V2.typ","docs/ch03_SurfaceFluxes/aerodynamic_conductance.typ","docs/ch03_SurfaceFluxes/forcing_reference_height.md"]},{"name":"Other — ch04_Photosynthesis","slug":"other-ch04-photosynthesis","files":["docs/ch04_Photosynthesis/photosynthesis.typ","docs/ch04_Photosynthesis/photosynthesis_comparison.typ"]},{"name":"Other — manual","slug":"other-manual","files":["docs/manual/ch01_Input.typ","docs/manual/ch02_Param_SoilVeg_Types.typ","docs/manual/manual.typ","docs/manual/软著/ModernBEPS.typ"]},{"name":"Other — 结题","slug":"other","files":["docs/结题/ModernBEPS_精度证明.typ","docs/结题/数据方面.typ","docs/结题/结题_PART2.typ"]},{"name":"Other — path_mnt","slug":"other-path-mnt","files":["path_mnt"]},{"name":"Other — plan.md","slug":"other-plan-md","files":["plan.md"]},{"name":"Other — repomix.ps1","slug":"other-repomix-ps1","files":["repomix.ps1"]},{"name":"Other — debug","slug":"other-debug","files":["test/debug/debug-Rn.qmd","test/debug/debug.md","test/debug/debug_Rln.jl","test/debug/debug_Rln2.jl","test/debug/debug_snow.jl"]},{"name":"Other — modules","slug":"other-modules","files":["test/modules/modules.jl","test/modules/test-Soil.jl","test/modules/test-aerodynamic_conductance.jl","test/modules/test-param.jl","test/modules/test-photosynthesis.jl","test/modules/test-radiation.jl","test/modules/test-rainfall_stage1.jl","test/modules/test-sensible_heat.jl","test/modules/test-setup.jl","test/modules/test-snowpack.jl","test/modules/test-surface_temperature.jl"]},{"name":"Other — test","slug":"other-test","files":["test/runtests.jl","test/test-BEPSmodel.jl","test/test-ModelParams.jl","test/test-StateSeries.jl","test/test-UpdateSoilMoisture_Q0.jl","test/test-beps_main.jl","test/test-beps_modern.jl","test/test-macro.jl","test/test-performance.jl","test/test-photosynthesis_standalone.jl","test/test-soil_sm.jl","test/test-utilize.jl"]}]}]};

 (function() {
  var activePage = 'overview';

  document.addEventListener('DOMContentLoaded', function() {
    mermaid.initialize({ startOnLoad: false, theme: 'neutral', securityLevel: 'loose' });
    renderMeta();
    renderNav();
    document.getElementById('menu-toggle').addEventListener('click', function() {
      document.getElementById('sidebar').classList.toggle('open');
    });
    if (location.hash && location.hash.length > 1) {
      activePage = decodeURIComponent(location.hash.slice(1));
    }
    navigateTo(activePage);
  });

  function renderMeta() {
    if (!META) return;
    var el = document.getElementById('meta-info');
    var parts = [];
    if (META.generatedAt) {
      parts.push(new Date(META.generatedAt).toLocaleDateString());
    }
    if (META.model) parts.push(META.model);
    if (META.fromCommit) parts.push(META.fromCommit.slice(0, 8));
    el.textContent = parts.join(' \u00b7 ');
  }

  function renderNav() {
    var container = document.getElementById('nav-tree');
    var html = '<div class="nav-section">';
    html += '<a class="nav-item overview" data-page="overview" href="#overview">Overview</a>';
    html += '</div>';
    if (TREE.length > 0) {
      html += '<div class="nav-group-label">Modules</div>';
      html += buildNavTree(TREE);
    }
    container.innerHTML = html;
    container.addEventListener('click', function(e) {
      var target = e.target;
      while (target && !target.dataset.page) { target = target.parentElement; }
      if (target && target.dataset.page) {
        e.preventDefault();
        navigateTo(target.dataset.page);
      }
    });
  }

  function buildNavTree(nodes) {
    var html = '';
    for (var i = 0; i < nodes.length; i++) {
      var node = nodes[i];
      html += '<div class="nav-section">';
      html += '<a class="nav-item" data-page="' + escH(node.slug) + '" href="#' + encodeURIComponent(node.slug) + '">' + escH(node.name) + '</a>';
      if (node.children && node.children.length > 0) {
        html += '<div class="nav-children">' + buildNavTree(node.children) + '</div>';
      }
      html += '</div>';
    }
    return html;
  }

  function escH(s) {
    var d = document.createElement('div');
    d.textContent = s;
    return d.innerHTML;
  }

  function navigateTo(page) {
    activePage = page;
    location.hash = encodeURIComponent(page);

    var items = document.querySelectorAll('.nav-item');
    for (var i = 0; i < items.length; i++) {
      if (items[i].dataset.page === page) {
        items[i].classList.add('active');
      } else {
        items[i].classList.remove('active');
      }
    }

    var contentEl = document.getElementById('content');
    var md = PAGES[page];

    if (!md) {
      contentEl.innerHTML = '<div class="empty-state"><h2>Page not found</h2><p>' + escH(page) + '.md does not exist.</p></div>';
      return;
    }

    contentEl.innerHTML = marked.parse(md);

    // Rewrite .md links to hash navigation
    var links = contentEl.querySelectorAll('a[href]');
    for (var i = 0; i < links.length; i++) {
      var href = links[i].getAttribute('href');
      if (href && href.endsWith('.md') && href.indexOf('://') === -1) {
        var slug = href.replace(/\.md$/, '');
        links[i].setAttribute('href', '#' + encodeURIComponent(slug));
        (function(s) {
          links[i].addEventListener('click', function(e) {
            e.preventDefault();
            navigateTo(s);
          });
        })(slug);
      }
    }

    // Convert mermaid code blocks into mermaid divs
    var mermaidBlocks = contentEl.querySelectorAll('pre code.language-mermaid');
    for (var i = 0; i < mermaidBlocks.length; i++) {
      var pre = mermaidBlocks[i].parentElement;
      var div = document.createElement('div');
      div.className = 'mermaid';
      div.textContent = mermaidBlocks[i].textContent;
      pre.parentNode.replaceChild(div, pre);
    }
    try { mermaid.run({ querySelector: '.mermaid' }); } catch(e) {}

    window.scrollTo(0, 0);
    document.getElementById('sidebar').classList.remove('open');
  }
 })();

 </script>
 </body>
 </html>
No results found