callsedim

气溶胶/冰晶重力沉降统一入口例程。

所在文件

LMDZ.MARS\libf\phymars\callsedim_mod.F

MODULE callsedim_mod → SUBROUTINE callsedim（源码第 1-754 行）。

例程类型

subroutine（模块内唯一 CONTAINS 例程）。

职责

源码依据（第 35-50 行注释 + 整体结构）：在每个物理时间步内，对所有“有半径”的气溶胶/冰晶示踪物施加重力沉降（gravitational sedimentation），计算：

1. 各层因沉降产生的混合比倾向 pdqsed；
2. 沉降到地表的通量 pdqs_sed；
3. 沉降后更新的粒子几何平均半径 rdust / rstormdust / rtopdust / rice。

该例程是尘埃循环、水循环、（部分）CO2 云循环共享的底层沉降调度器：它本身不实现下落速度公式，而是按示踪物类别准备粒径分布与密度，然后调用 newsedim 完成单示踪物的垂直输运。

源码依据（第 39-48 行 history 注释）：

• F. Forget 1999 原始版本；
• J.-B. Madeleine 2010 引入 doubleq（双矩）技术，使 physiq 只需一次调用 callsedim；
• J. Audouard 2016/09 加入 co2clouds 分支——CO2 冰与 CCN 的沉降不在本例程，而在 co2cloud.F 的微时间步内完成。

参数

实参对齐来自 physiq_mod.F:2143（见下方“调用的其他例程/被调用关系”），方向取自源码第 60-85 行 intent 声明。

源码依据：第 7-12 行参数列表；第 60-85 行 intent 声明。

复现风险：rstormdust / rtopdust 仅在 rdstorm / topflows 开关打开时被赋值（第 704-722 行）。关闭时其值未在本例程初始化，调用方需注意未定义内存（源码未对其做无条件初始化）。

使用的 module 变量

来自 USE 列表（源码第 14-32 行）。

本例程内的 SAVE 状态

源码依据（第 105-186 行）：以下为 SAVE/THREADPRIVATE 变量，仅在 firstcall 时初始化：

• beta：冰粒形状修正因子（默认 0.75，可由 ice_shape 覆盖）；
• rd(nr)、rdi(nr+1)、rr(ninter,nr)：doubleq 离散粒径分箱（nr=12，ninter=4）；
• 一系列示踪物索引 idust_mass, idust_number, iccn_mass, iccn_number, istormdust_*, itopdust_*, iccnco2_*, ico2_ice；
• firstcall：首调标志。

调用的其他例程

被调用关系（源码依据 physiq_mod.F:2143）：

physiq  →  callsedim  →  newsedim  →  vlz_fi
                      →  updaterdust / updaterice_micro / updaterice_typ
                      →  vlz_fi (HDO 同位素)

控制流程

源码依据（结构性阅读第 193-751 行）：

1. firstcall 初始化块（第 193-402 行）：仅首次执行
   • doubleq：构造 12 个对数等比分箱半径 rd、边界 rdi、子积分点 rr；按 noms 定位 dust_mass/dust_number 索引，缺失则 abort_physic。
   • microphys：定位 ccn_mass/ccn_number。
   • co2clouds：定位 CO2 相关 CCN/冰索引（含 co2useh2o、meteo_flux 子情形）。
   • water：读取 ice_shape → beta。
   • rdstorm：定位 stormdust_*。
   • topflows：定位 topdust_*。
2. 初始化（第 404-428 行）：把上游倾向并入局部 zqi/zq0/zt；计算层质量 masse 和厚度 epaisseur。
3. 预计算 doubleq 半径（第 433-467 行）：对 doubleq/rdstorm/topflows 各自调用 updaterdust 得到 r0dust/r0stormdust/r0topdust。
4. 示踪物主循环 do iq=1,nq（第 469-686 行），仅对 radius(iq)>1e-9 且非 CO2 冰/CCN/HDO 的示踪物处理，按三类分支：
   • DOUBLEQ 分支（第 482-584 行）：尘埃/stormdust/topdust 质量与数浓度；
   • WATER CYCLE 分支（第 588-649 行）：CCN 与 H2O 冰，含 HDO 同位素子情形；
   • GENERAL 分支（第 653-661 行）：其余单一示踪物。
   • 末尾统一计算 pdqsed（第 667-682 行）。
5. 半径更新（第 690-751 行）：对 doubleq/rdstorm/topflows 调 updaterdust；对 water 按 microphys 与否调 updaterice_micro/updaterice_typ。

提前跳过条件（第 470-478 行）：CO2 冰、各类 CO2 CCN、HDO 冰被显式排除——CO2 在微时间步内沉降，HDO 由 H2O 携带。

核心计算步骤

1. 层质量（第 425 行，源码依据）：

   masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
   epaisseur(ig,l) = zlev(ig,l+1) - zlev(ig,l)

2. doubleq 粒径分布积分（第 521-560 行）：对每个分箱 ir，在 ninter=4 个子点上用梯形法积分对数正态分布，得到该分箱权重 qr，再归一化到该示踪物总量 zqi：
   • radpower = 2（质量）或 -1（数浓度）；
   • qr(ig,l,ir) = zqi(ig,l,iq) * qr / Sq(ig,l)。
3. 逐分箱沉降（第 568-584 行）：对 doubleq 的每个固定半径 rd(ir) 调 newsedim，密度恒为 rho_dust，beta=0.5；累加地表通量与各层新混合比。
4. 水冰沉降（第 593-602 行）：微物理时用 rsedcloud+rhocloud、beta；非微物理时用 rsedcloud+rho_q(iq)、beta。
5. 同位素输运（第 611-649 行）：用母体输运质量 w 通过 vlz_fi 输运 HDO/H2O 比值，再换算 HDO 冰混合比与地表通量。
6. 通用示踪物（第 654-656 行）：用 radius(iq)+rho_q(iq)、beta=1.0 调一次 newsedim。
7. 最终倾向（第 669-670 行，源码依据）：

   pdqsed(ig,l,iq) = (zqi(ig,l,iq) - (pq + pdqfi*ptimestep)) / ptimestep

关键公式或算法

源码依据（第 198-210 行）—— doubleq 分箱构造（对数等比）：

rd(ir)  = rdmin * (rdmax/rdmin)^((ir-1)/(nr-1))           ! nr=12, rdmin=1e-8, rdmax=30e-6
rdi(1)  = rdimin                                          ! 1e-8
rdi(ir) = sqrt(rd(ir-1)*rd(ir))   (ir=2..nr)
rdi(nr+1) = rdimax                                        ! 1e-4
rr(iint,ir) = rdi(ir) * (rdi(ir+1)/rdi(ir))^((iint-1)/(ninter-1))   ! ninter=4

对数正态分布积分核（第 528-542 行）：

weight(r) = r^radpower * exp( -(ln(r/r0))^2 / (2*sigma0^2) )

其中 sigma0 = varian（tracer_mod），radpower=2（质量矩）或 -1（数浓度矩）。

newsedim 内部的下落速度公式（Stokes 定律 + Cunningham 滑移修正等）见 newsedim_mod.md。本页只确认 callsedim 把粒径 rd、密度 rho、形状因子 beta 传入 newsedim。

待确认：beta（形状因子）在 doubleq 分支取 0.5、water 分支取 beta（默认 0.75）、general 分支取 1.0；其物理含义（对下落速度的修正方式）需在 newsedim 专页确认。

边界条件与保护逻辑

源码依据：

• 示踪物缺失保护：firstcall 时若开启某方案却找不到对应示踪物索引，调用 abort_physic 中止（第 236, 258, 315, 372, 397 行）。
• 半径阈值：仅 radius(iq)>1e-9 的示踪物参与沉降（第 470 行），过滤气体示踪物。
• CO2 显式跳过：CO2 冰与各类 CO2 CCN 在主循环条件中被排除（第 470-477 行）。
• 同位素母体阈值：zq0(iq)>qparentmin 才计算同位素比，否则比值置 0（第 620-624 行）；masseq 下限为 masseqmin（第 626, 635 行）。
• 同位素非法保护：若有 child 的示踪物不是 igcm_h2o_ice，则 abort_physic（第 615, 676 行）。

输入数据来源

• 配置开关：callkeys_mod（water/microphys/doubleq/rdstorm/topflows/co2clouds/...），由 callphys.def 经初始化填充（推断：根据 callkeys_mod 的常规填充路径；本文件未直接读取 def）。
• 运行参数 ice_shape：getin_p（第 342 行），默认 0.75。
• 示踪物属性：tracer_mod（radius/rho_q/rho_dust/varian/noms/...）。
• 上游倾向：由 physiq 传入的 pdqfi（沉降前累计倾向）与 pdt。

复现风险：ice_shape 的实际运行值由运行目录的 def 文件决定；若复现实验未提供该项，则采用默认 0.75。沉降通量对 beta 敏感。

输出与副作用

• 写出参数：pdqsed（各层沉降倾向）、pdqs_sed（地表通量）、rdust/rstormdust/rtopdust/rice、rhocloud（inout，微物理时更新）。
• SAVE 状态：firstcall 时设置 beta、rd/rdi/rr、各示踪物索引、firstcall=.false.。
• 标准输出：firstcall 时打印各示踪物索引与水循环参数（第 223, 227, 246, 250, 340-348 行等）；可能产生 abort_physic 终止。
• 无文件 I/O：本例程不直接读写数据文件。

复现风险：CO2 冰与 CO2 CCN 的沉降不在此处完成；复现 CO2 云循环时必须在 co2cloud.F 微时间步内单独处理，不能假设 callsedim 已覆盖它们。

伪代码

procedure callsedim(ngrid, nlay, ptimestep, ...):
  if firstcall:
    if doubleq:   build rd[nr], rdi[nr+1], rr[ninter,nr]; locate dust_mass/number
    if microphys: locate ccn_mass/number
    if co2clouds: locate co2 ccn/ice indices (+ h2o/meteor sub-cases)
    if water:     beta = getin("ice_shape", 0.75)
    if rdstorm:   locate stormdust_mass/number
    if topflows:  locate topdust_mass/number
    firstcall = false

  zqi = pq + pdqfi*ptimestep          # apply upstream tendencies
  zt  = pt + pdt*ptimestep
  masse = (pplev[l] - pplev[l+1]) / g
  epaisseur = zlev[l+1] - zlev[l]

  if doubleq:   r0dust     = updaterdust(dust_mass, dust_number, tauscaling)
  if rdstorm:   r0stormdust = updaterdust(stormdust_*)
  if topflows:  r0topdust   = updaterdust(topdust_*)

  for iq in 1..nq:
    if radius[iq] <= 1e-9: continue
    if iq is CO2 ice / CO2 ccn / HDO ice: continue

    if doubleq and iq in {dust,storm,top}_{mass,number}:
      r0 = matching r0*
      radpower = 2 if mass else -1
      for each bin ir:                # log-normal trapezoid integration
        qr[ir] = integrate(rr[:,ir], r0, sigma0, radpower)
      normalize qr to zqi[iq]
      zqi[iq] = 0; pdqs_sed[iq] = 0
      for each bin ir:
        newsedim(radius=rd[ir], rho=rho_dust, beta=0.5) -> wq
        pdqs_sed[iq] += wq[1]/ptimestep
        zqi[:,iq]    += qr[:,ir]

    else if iq in {ccn_mass, ccn_number, h2o_ice}:
      if microphys: newsedim(rsedcloud, rhocloud, beta)
      else:         newsedim(rsedcloud, rho_q[iq], beta)
      pdqs_sed[iq] = wq[1]/ptimestep
      if iq == h2o_ice and has child (HDO):
        Ratio = zq0[hdo]/zq0[h2o]   (guarded by qparentmin)
        vlz_fi(Ratio, masseq, w=wq) -> transport HDO with H2O mass
        update zqi[hdo], pdqs_sed[hdo]

    else:   # general tracer
      newsedim(radius[iq], rho_q[iq], beta=1.0)
      pdqs_sed[iq] = wq[1]/ptimestep

    pdqsed[:,iq] = (zqi[:,iq] - (pq + pdqfi*ptimestep)) / ptimestep
    if HDO child: pdqsed[:,hdo] likewise

  if doubleq:  rdust      = updaterdust(zqi dust_*)
  if rdstorm:  rstormdust = updaterdust(zqi storm_*)
  if topflows: rtopdust   = updaterdust(zqi top_*)
  if water:
    if microphys: rice = updaterice_micro(zqi h2o_ice, ccn_*, -> rhocloud)
    else:         rice = updaterice_typ(zqi h2o_ice, tau, zlay)
end procedure

复现注意事项

• 沉降是“先并入上游倾向 → 沉降 → 反算倾向”的算子分裂模式；复现时 pdqsed 必须按 (zqi_final - (pq+pdqfi*dt))/dt 计算，而非直接覆盖混合比。
• doubleq 用 12 个固定半径分箱近似对数正态谱的沉降，每箱单独调用 newsedim；不能简化为单一有效半径，否则下落速度的尺寸依赖丢失。
• 三类分支的 beta 取值不同（doubleq=0.5，water=beta，general=1.0），是复现下落速度修正的关键。
• HDO 不独立沉降，而是用 H2O 冰的输运质量 w 经 vlz_fi 携带，保证同位素比守恒。
• 半径更新（rdust/rice 等）发生在沉降之后，使用沉降后的 zqi；这些输出半径供下一物理过程/辐射使用。

待确认

• callkeys_mod 各开关如何由 def 文件填充（本文件未直接读取，标推断）。
• beta 在平均自由程系数中生效的物理文献依据仍以 newsedim_mod.md 的待确认项为准。

复现风险

• 复现风险：rstormdust/rtopdust 在对应开关关闭时不被赋值。
• 复现风险：CO2 冰/CCN 沉降不在此例程，需在 co2cloud.F 微时间步内完成。
• 复现风险：ice_shape（beta）默认 0.75，但实际由运行 def 决定，影响水冰沉降通量。

相关页面

• phymars 模块总览
• 尘埃循环主题
• newsedim_mod.md：核心沉降数值方案。
• vlz_fi.md：限坡半拉格朗日垂直输运。
• updaterad.md：尘埃/水冰有效半径反演模块