尘埃循环主题

主题目标

尘埃循环是 LMDZ.MARS 物理包中影响最广的过程之一。火星大气尘埃决定了大气的辐射加热/冷却，进而驱动温度、环流和其他循环。本主题串联尘埃从"注入大气—被风/对流输运—改变粒径分布—产生光学厚度—被辐射模块使用—沉降回地表"的完整流程，目标是说明：

1. 尘埃示踪剂（质量、数量、暴尘、地形尘）如何在大气中被抬升、注入与沉降。
2. 尘埃粒径分布如何由二阶矩（质量+数量）方案计算。
3. 尘埃光学厚度如何计算，并如何被缩放以匹配观测尘暴场景（dust scenario）。
4. 尘埃如何与辐射主题耦合。

源码依据：dust_param_mod.F90 第 1 行注释 "This module contains flags and saved variables for the dust cycle"。

入口

尘埃循环没有单一入口，而是分散在 physiq() 主物理例程（physiq_mod.F）的多个阶段中按顺序调用。

源码依据：physiq_mod.F 第 21-99、1379-1530、2107-2135 行的 use/CALL 语句。

参与文件

rocketduststorm、topmons、dust_coagulation 的内部数值方案已由对应文件页覆盖；本主题只保留跨模块调用链、配置门控和复现检查点。

尘埃示踪剂

尘埃以一组示踪剂（tracer）形式存在，具体哪些被激活取决于配置开关。

源码依据：tracer_mod 中定义的 igcm_* 索引；updatereffrad_mod.F 第 13-23 行、aeropacity_mod.F 第 20-23 行的 use tracer_mod。

配置开关

dustscaling_mode 三种取值（dust_scaling_mod.F90 与 conf_phys.F 注释）：

• 0：freedust，tauscaling=1，不重缩放。
• 1：GCM v5.3 风格，计算 tauscaling 使 GCM 光学厚度匹配 tau_pref_scenario。
• 2：GCM v6 风格，tauscaling=1，仅对辐射传输用 dust_rad_adjust 缩放。

源码依据：conf_phys.F 第 559-578 行；dust_scaling_mod.F90 第 42-107 行。

复现风险：dustscaling_mode 的默认值依赖 freedust，复现时必须先确定 freedust、doubleq 的组合，否则会触发 abort_physic（conf_phys.F L553-578 有多处互斥校验）。

尘埃场景（iaervar）

read_dust_scenario 根据 iaervar 决定读取哪个 NetCDF 场景文件：

源码依据：read_dust_scenario_mod.F90 第 72-95 行。

场景变量类型：dustop（可见光消光）或 cdod/cdodext（红外吸收/消光）。若为红外（IRscenario=.true.），最终乘 IRtoVIScoef(ig) 转换为可见光消光。旧 TES 文件额外乘 1.3*(odpref/700.) 校正；swrtype==1（Fouquart）时再乘 0.825/1.3。

源码依据：read_dust_scenario_mod.F90 第 139-174、387-411 行。

状态变量

参考压力常数 odpref = 610. Pa（dust_param_mod.F90 第 19 行）。

主调用链

physiq()  (physiq_mod.F)
│
├── [辐射前] callradite()
│   ├── updatereffrad()          更新尘埃粒径 reffrad
│   │   └── updaterdust()        doubleq: rdust = (r3n_q*qdust/ndust)^(1/3)
│   ├── aeropacity()             计算尘埃光学厚度
│   │   ├── read_dust_scenario() 读取并插值场景 tau_pref_scenario
│   │   │   └── (IR场景) × IRtoVIScoef
│   │   ├── [累加] aerosol(:,:,iaerdust)  及 dsodust, tau_pref_gcm
│   │   └── compute_dustscaling()
│   │       ├── mode 1: tauscaling = tau_pref_scenario*ps/odpref/taudust
│   │       │           aerosol *= tauscaling
│   │       └── mode 2: compute_dust_rad_adjust()  (跨 sol 插值)
│   │                   aerosol *= dust_rad_adjust
│   └── (lwmain/swmain 使用 aerosol)   → 辐射主题
│
├── [3.x 输运] rocketduststorm() / topmons()   暴尘/地形尘垂直输运
├── [3.4 注入] compute_dtau()        → dustliftday (注入率)
│
├── [扩散 vdifc] dust_windstress_lift()  风应力抬升 (lifting=T)
│                 使用 dustliftday 注入 (dustinjection>0)
├── [尘卷风] dustdevil()              callddevil=T
└── [沉降] callsedim()                尘埃/水冰重力沉降

源码依据：调用顺序见 physiq_mod.F；updatereffrad → updaterdust 见 updatereffrad_mod.F 第 99-106 行；aeropacity → compute_dustscaling 见 aeropacity_mod.F 第 39 行 use dust_scaling_mod。

推断：rocketduststorm、topmons、compute_dtau 在 physiq 中的精确先后次序按 physiq_mod.F 第 1379-1530 行的出现顺序判断；其与扩散/沉降段的相对位置以源码行号为准，但完整控制流（含各 IF 分支条件）未逐一展开。

计算阶段

阶段 1：粒径更新（updatereffrad → updaterdust）

doubleq.AND.active 时，对每层每点由质量与数量混合比计算几何平均半径：

rdust = ( r3n_q * qdust / ndust )^(1/3)

并对极小/极大值用阈值 r3dustmin/r3dustmax（对应 rdustmin=1e-8 m、rdustmax=500e-6 m）截断；否则（非 doubleq）固定 rdust=0.8e-6 m、nueffdust=0.3。随后转为辐射有效半径：dust_doubleq 用 reffrad = rdust*ref_r0，dust_conrath 用 reffrad = rdust*(1+nueffdust)^2.5。

源码依据：updaterad.F90 第 203-238 行；updatereffrad_mod.F 第 99-113、228-234 行。

阶段 2：光学厚度计算（aeropacity）

每层各类尘埃光学厚度按下式累加（doubleq）：

aerosol(ig,l,iaer) = 0.75 * QREFvis3d / (rho_dust * reffrad) * qmass * dp / g

同时计算红外密度缩放光学厚度 dsodust（用 QREFir3d，供与 TES/MCS 观测比较）和 GCM 柱光学厚度 tau_pref_gcm（换算到 610 Pa 参考压力）。

源码依据：aeropacity_mod.F 第 423-455、888-910、726-736 行。

阶段 3：尘埃缩放（compute_dustscaling）

• mode 1：先对所有尘埃类累加柱光学厚度 taudust，再令 tauscaling = tau_pref_scenario * pplev(:,1) / odpref / taudust， 然后 aerosol *= tauscaling（下限 1e-20）。这样 GCM 光学厚度被强制匹配观测场景。
• mode 2：tauscaling=1，调用 compute_dust_rad_adjust 算 dust_rad_adjust，仅对辐射传输用的 aerosol 缩放。
• mode 0（freedust）：tauscaling=1，不缩放。

源码依据：dust_scaling_mod.F90 第 40-108 行。

compute_dust_rad_adjust（mode 2）逻辑：在每个格点本地时间跨越 t_scenario_sol=14/24 sol（即 14:00）时更新目标，目标 = tau_pref_scenario_next*ps/odpref/taudust（上限 5），再在前后两天的 14:00 端点间按本地时间线性插值。

源码依据：dust_rad_adjust_mod.F90 第 66-119 行；dust_param_mod.F90 第 28 行 t_scenario_sol。

阶段 4：地表抬升与注入

• 风应力抬升（dust_windstress_lift）：在 GCM 非 doubleq 分支中由 vdifc 调用。当摩擦速度 us = pcdh/sqrt(pcdh_true) 超过阈值 ust = sqrt(stress_seuil/rho)（stress_seuil=0.0225 N/m²）时，水平通量 fhoriz = 2.61*(rho/g)*(us-ust)*(us+ust)^2（White 1979），垂直注入通量 dqslift(ig,iq) = -alpha_lift(iq)*fhoriz。有 CO₂ 冰覆盖（co2ice>0）的点不抬升。

  源码依据：dust_windstress_lift.F90 第 102-162 行；vdifc_mod.F 第 813-883 行。

• GCM doubleq 分支边界：vdifc_mod.F 中 lifting 先进入 doubleq 判断。doubleq=.true. 且 dustinjection=0 时使用 old injection scheme，直接给 dust_mass/dust_number 设置 -alpha_lift 源项；dustinjection=1 时用 compute_dtau 产生的 dustliftday 注入。只有落到非 doubleq 的 else 分支时才调用 dust_windstress_lift()。因此在当前 GCM two-moment dust 配置下，lifting=.true. 不能等同于风应力阈值起尘。

  源码依据：vdifc_mod.F 第 813-883 行；conf_phys.F 第 553-604 行。

• 尘卷风（dustdevil，callddevil=T）：基于 Renno et al. 1998，尘卷风活动 devila = max(0,(1-b)*Fs - seuil)，其中 Fs 为地表感热通量（用固定风速 15 m/s），抬升通量 pdqs_dev = -alpha_devil*devila，并把尘埃注入到边界层顶。

  源码依据：dustdevil.F90 第 137-221 行。

• 场景注入（compute_dtau，dustinjection>0）：在本地 14:00 时计算 GCM 与场景目标光学厚度差 dtau，当 dtau>0 时给出注入率 dustliftday = coeff_injection*(dtau*ps/odpref)/(daysec*(tf-ti))，该率随后在 vdifc 中用于注入。

  源码依据：compute_dtau_mod.F90 第 87-102 行；vdifc_mod.F 第 813-857 行。

阶段 5：沉降（callsedim）

尘埃（与水冰）按重力沉降回地表。

沉降调度见 callsedim_mod，单粒径下落速度公式和 Stokes-Cunningham 修正见 newsedim_mod，守恒输运格式见 vlz_fi。

数据依赖

源码依据：read_dust_scenario_mod.F90 第 100-110 行；dust_windstress_lift.F90 第 60-84 行（#ifdef MESOSCALE）。

与其他主题的耦合

源码依据：vdifc_mod.F 第 31、43 行 use dust_param_mod/use dust_windstress_lift_mod；updaterad.F90 第 255-292 行 updaterccn/updaterccnCO2。

输出结果

源码依据：compute_dtau_mod.F90 第 83-108 行的 write_output；aeropacity_mod.F 输出参数列表第 95-117 行。

伪代码总览

! ---- 每个物理步，按 physiq 顺序 ----

! A. 辐射前：粒径与光学厚度
CALL updatereffrad(...)       ! rdust = (r3n_q*qdust/ndust)^(1/3) ; reffrad
CALL aeropacity(...)
  IF (iaervar/=1) CALL read_dust_scenario(... tau_pref_scenario)
  DO 各层各尘类:
    aerosol = 0.75*QREFvis3d/(rho_dust*reffrad)*qmass*dp/g
    累加 dsodust, tau_pref_gcm
  CALL compute_dustscaling(...)
    SELECT CASE (dustscaling_mode)
      CASE(1): tauscaling = tau_pref_scenario*ps/odpref/taudust
               aerosol = max(1e-20, aerosol*tauscaling)
      CASE(2): CALL compute_dust_rad_adjust(...)   ! 跨 sol 14:00 插值
               aerosol = max(1e-20, aerosol*dust_rad_adjust)
      CASE(0): tauscaling = 1                       ! freedust
! (随后 lwmain/swmain 使用 aerosol —— 见辐射主题)

! B. 输运与注入
IF (rdstorm)  CALL rocketduststorm(...)
IF (topflows) CALL topmons(...)
IF (dustinjection>0) CALL compute_dtau(... dustliftday)

! C. 抬升（在 vdifc 内）
IF (lifting) THEN
  IF (doubleq) THEN
    IF (dustinjection==0) 使用 old constant lifting 源项
    IF (dustinjection==1) 使用 dustliftday 注入
  ELSE
    CALL dust_windstress_lift(... dqslift)  ! 风应力阈值路径
  ENDIF
ENDIF
IF (callddevil) CALL dustdevil(... pdqdev,pdqs_dev)

! D. 沉降
CALL callsedim(...)            ! 尘埃 + 水冰

复现检查点

必须保留的逻辑

1. 二阶矩半径公式：rdust = (r3n_q*qdust/ndust)^(1/3)，并保留 r3dustmin/max 阈值截断。
2. 光学厚度公式：0.75*QREF/(rho_dust*reffrad)*qmass*dp/g，可见光用 QREFvis3d，红外（dsodust）用 QREFir3d。
3. 参考压力 610 Pa：所有柱光学厚度换算都用 odpref=610（旧 TES 文件需 700→610 校正）。
4. 缩放模式选择：dustscaling_mode 与 freedust/doubleq/iaervar 的互斥约束必须保留，否则运行会中止。
5. 抬升阈值与公式：White 1979 水平通量公式与 stress_seuil=0.0225；CO₂ 冰覆盖点不抬升。注意该路径在当前 GCM doubleq=.true. 分支中不会被 dustinjection=0 调到。
6. dust_rad_adjust 跨 sol 插值：以本地 14:00（t_scenario_sol=14/24）为锚点，上限 5。

关键公式

1. 几何平均半径：rdust = (r3n_q*qdust/ndust)^(1/3)
2. 层光学厚度：aerosol = 0.75*QREFvis3d/(rho_dust*reffrad)*qmass*(p_l - p_{l+1})/g
3. mode1 缩放：tauscaling = tau_pref_scenario*ps/odpref/taudust
4. 风应力水平通量：fhoriz = 2.61*(rho/g)*(us-ust)*(us+ust)^2，ust = sqrt(stress_seuil/rho)
5. 注入率：dustliftday = coeff_injection*(dtau*ps/odpref)/(daysec*(tf-ti))

验证检查点

1. 光学厚度合理性：背景尘埃可见光柱光学厚度典型 0.1-1（晴期），尘暴期可达数倍。dustscaling_mode=0 的 10 天短试验已见 dust-free-evolution-tests：dust devil 候选 E 的 tau_pref_gcm 10 天均值约 0.526，末值约 0.482；仍需长积分验证。
2. 缩放一致性：mode1 下 tau_pref_gcm 经缩放后应≈tau_pref_scenario。
3. 质量守恒：抬升/注入为负趋势（离开地表），沉降为正向地表通量，长期收支应闭合。待确认：未核验收支闭合实现。
4. 自由尘埃配置边界：当前 def-only 候选为 freedust=.true., dustscaling_mode=0, lifting=.false., dustinjection=0, callddevil=.true.；若要 two-moment 风应力阈值起尘，需要源码改造或新增分支。

相关模块

• 高分辨率稳定性：高分辨率下的沙尘数值稳定性
• 辐射：callradite_mod.F、lwmain_mod.F、swmain_mod.F、aeroptproperties.F、suaer.F90（见 radiation.md）
• 粒径：updaterad.F90、updatereffrad_mod.F
• 抬升/注入：dust_windstress_lift.F90、dustdevil.F90、compute_dtau_mod.F90
• 缩放：dust_scaling_mod.F90、dust_rad_adjust_mod.F90、read_dust_scenario_mod.F90
• 暴尘/地形尘：rocketduststorm_mod.F90、topmons_mod.F90
• 微物理：dust_coagulation_mod.F90、watercloud_mod.F、improvedclouds_mod.F90
• 沉降：callsedim_mod.F

待确认

• submicron（亚微米尘）分布在抬升/输运中的完整处理（vdifc 第 815 行起 doubleq.AND.submicron 分支）只读到片段。
• 各尘埃光学厚度的典型量级已有 10 天 free-dust 试验补充，但长积分稳定性、季节尺度偏差与收支闭合仍未完成数值验证。
• physiq 中 rocketduststorm/topmons/compute_dtau 与沉降段的完整控制流分支条件未全部展开（仅按行号顺序判断）。
• GCM doubleq=.true. 分支下接入风应力阈值 lifting 的设计尚未实现；现有 dustinjection=0 是 old constant lifting。