dust_coagulation_main

所在文件

LMDZ.MARS\libf\phymars\dust_coagulation_mod.F90

SUBROUTINE dust_coagulation_main（源码第 93-400 行，共 308 行）。

例程类型

subroutine

职责

源码依据（第 4-7 行注释 "Dust coagulation scheme, Based on Bertrand et al., 2022"）：在每个物理时间步中，对所有格点的所有层计算**尘埃凝并（coagulation）**倾向。

凝并的物理含义：两个尘埃粒子碰撞后粘连形成更大的粒子，改变尘埃粒径分布。本例程：

1. 将 mass/number 双矩量展开为 nres_coag=10 个 bin 的对数正态分布；
2. 用离散化的凝并方程（Smoluchowski 方程）计算每个 bin 的数浓度变化；
3. 通过体积守恒归一化修正质量守恒；
4. 将 bin 分布重新合并回 mass/number 的双矩量；
5. 输出质量/数浓度的倾向 pdqcoag。

涉及三种尘埃示踪物：dust、stormdust、topdust（各自含 mass/number 双矩量）。

参数

使用的 module 变量

调用的其他例程

控制流程

1. 初始化（第 142-213 行）：无条件执行
2. 凝并分支（第 220 行 if (full_coag_equations)）：
   • 完整方程模式（第 220-266 行）：逐格点逐层逐 bin 直接计算 kernel 函数
   • 查表模式（第 268-357 行）：先通过 findval 定位 (T, mfp) 或 (T, P) 角点，再用 boxinterp 双线性插值 kernel
3. 凝并激活阈值（第 223、273 行）：if (ntot(i,j) > 1000.) — 总粒子浓度 > 1000 nb/m³ 时才触发
4. 质量守恒归一化（第 261-262、350-351 行）：在每个格点层执行
5. 新矩量计算（第 362-398 行）：将 bin 分布重新合并为 mass/number 双矩量
6. 非负保护（第 391-396 行 where 语句）：防止混合比变负

核心计算步骤

阶段 1：初始化（第 142-213 行）

1. 分布常数（第 145-147 行）：

   dens = rho_dust         ! 2500 kg/m³
   sig0 = dev_dt           ! 0.63676
   cst  = 0.75 / (π*dens) * exp(-4.5*sig0²)

2. 示踪物场推进（第 153-154 行）：

   zq = pq + pdqfi * ptimestep
   zt = pt + pdtfi * ptimestep

3. 大气密度（第 157-159 行）：rho = pplay / (r * pt)

4. 初始对数正态分布（第 173-193 行）：对三种尘埃示踪物（mass, mass-1=number）：

   mass_ini = max(zq(mass), 1e-15)
   numb_ini = max(zq(number), 1e-15)
   mtot += mass_ini                           ! 累积总质量用于守恒
   r0 = (0.75 * mass_ini / (π * numb_ini * dens))^(1/3) * exp(-1.5*sig0²)
   
   对每个 bin i:
     ndis_tab(i) = numb_ini * rho * deltar_coag(i) / (2*rad(i)*√(2π)*sig0)
                   * exp(-0.5*(ln(rad(i)/r0))² / sig0²)

   这是标准的对数正态分布 PDF，归一化为数浓度。

5. 各模式数浓度比（第 196-208 行）：

   rat_tab(iq) = ndis_tab(iq) / sum(ndis_tab(all_iq))
   rat_tab(iq) /= sum(rat_tab(all_iq))         ! 二次归一化确保 Σ=1

   此比值在凝并后用于将总分布重新拆分回各尘埃模式。

6. 总初始分布（第 211-213 行）：

   ndis = sum(ndis_tab(all_iq))     ! 所有模式合并
   ntot = sum(ndis(all_bins)) / rho ! nb/kg

阶段 2：凝并（第 215-357 行）

对每个格点层，若 ntot > 1000：

1. 环境参数（第 224-226 行）：pt0=pt, rho0=rho, mfp0=mfp(pt0,rho0)

2. 逐 bin 凝并方程（第 229-258 行）：

   Term 1（增益项，第 231-243 行）：所有 ii < k, jj ≤ k 的碰撞中，产物落入 bin k 的部分

   term1 = Σ_{jj=1..k} Σ_{ii=1..k-1} frac(ii,jj,k) * vol(ii) * kernel(ii,jj) * ndis_new(ii) * ndis(jj)
   term1 *= ptimestep / vol(k)

   Term 2（损耗项，第 246-255 行）：bin k 与所有 bin 碰撞而消失

   term2 = Σ_{jj=1..nres_coag} (1 - frac(k,jj,k)) * kernel(k,jj) * ndis(jj)
   term2 *= ptimestep

   半隐式更新（第 257 行）：

   ndis_new(k) = (ndis(k) + coal_fac * term1) / (1 + coal_fac * term2)

   这是经典的半隐式格式：增益项用部分更新值（ndis_new for ii<k），损耗项隐式处理。

3. 体积守恒归一化（第 261-262 行）：

   norm = Σ(ndis * vol) / Σ(ndis_new * vol)
   ndis_new *= norm

完整方程 vs 查表：

• 完整方程（第 220 分支）：直接调用 betab/betag/betade/betati 函数计算 kernel。
• 查表模式（第 268 分支）：通过 findval 在 table_pt/table_mfp/table_pres 中定位角点（第 279-290 行），再用 boxinterp 做 2D 双线性插值（第 300-340 行）。Brownian/Gravitational 按 (T, mfp) 插值；Diffusion Enhancement 按 (T, P) 插值。Turbulent kernel 标注 "not implemented yet"。

阶段 3：新矩量计算（第 359-398 行）

1. 按模式拆分（第 364-380 行）：用初始比值分配

   ndis_tab_new(iq) = ndis_new * rat_tab(iq)
   numb_new(iq) = Σ(ndis_tab_new(iq)) / rho
   mass_new(iq) = Σ(ndis_tab_new(iq) * vol) * dens / rho

2. 质量守恒（第 382-398 行）：

   mass_new(iq) *= mtot / mtotnew     ! 按总量比缩放确保总质量守恒
   pdqcoag(mass)  = (mass_new - mass_ini) / ptimestep
   pdqcoag(number)= (numb_new - numb_ini_dis) / ptimestep

3. 非负保护（第 391-396 行）：

   where (zq(mass) + pdqcoag(mass)*ptimestep ≤ 0):
     pdqcoag(mass) = -zq(mass)/ptimestep + 1e-15
   where (zq(number) + pdqcoag(number)*ptimestep ≤ 0):
     pdqcoag(number) = -zq(number)/ptimestep + 1e-14

关键公式或算法

凝并方程（Smoluchowski 方程的离散化）

dn_k/dt = ½ Σ_{i+j→k} β(i,j) n_i n_j - Σ_{j} β(k,j) n_k n_j

本代码的半隐式离散形式：

n'_k = (n_k + Σ_{ii<k,jj≤k} frac(ii,jj,k)·vol(ii)·β(ii,jj)·n'_ii·n_jj · Δt/V_k)
       / (1 + Σ_j (1-frac(k,j,k))·β(k,j)·n_jj · Δt)

对数正态分布初始化

n(r) = N·ρ·ΔD / (2·r·√(2π)·σ₀) · exp(-0.5·(ln(r/r₀))²/σ₀²)

其中 r₀ = (0.75·Q/(π·N·ρ_d))^(1/3) · exp(-1.5·σ₀²) 是几何中值半径。

frac 函数（体积分数分配）

当 bin i 和 bin j 的粒子凝并，产物体积 v = vol(i) + vol(j)。frac(i,j,k) 表示产物落入 bin k 的比例：

• 若 v ∈ [vol(k), vol(k+1))：frac = (vol(k+1)-v)/(vol(k+1)-vol(k)) · vol(k)/v
• 若 k = nres_coag 且 v ≥ vol(k)：frac = 1（溢出到最大 bin）
• 若 v ∈ (vol(k-1), vol(k))：frac = (v-vol(k-1))/(vol(k)-vol(k-1)) · vol(k)/v
• 否则 frac = 0

Kernel 函数

边界条件与保护逻辑

• 凝并激活阈值（第 223、273 行）：ntot(i,j) > 1000 nb/m³——低于此阈值凝并无效，倾向为零。
• 混合比下限（第 176-177 行）：mass_ini/numb_ini = max(zq, 1e-15)，避免除零。
• 非负保护（第 391-396 行）：where 语句防止 zq + pdqcoag*ptimestep ≤ 0，但 mass 和 number 的下限常数不同（1e-15 vs 1e-14）。
• 查表插值角点保护（第 280-290 行）：findval 返回 0（超出范围）时截断到最远端点；返回 1 时下推到 2，保证插值区间有效。
• 体积守恒归一化（第 261-262、350-351 行）：确保凝并前后总粒子体积不变。
• coal_fac 因子（第 257、346 行）：控制整体凝并强度，默认 1.0，可通过 callphys.def 修改。

输入数据来源

• 上游倾向：pdqfi/pdtfi 来自物理时间步中凝并之前的所有过程（辐射、湍流等）。
• 模块状态：rads_coag/vols_coag/deltar_coag 由 dust_coagulation_init 在首调时填充。
• 查表数据：table_b/g/de/ti 由 make_tables() 在 dust_coagulation_init 中预计算（仅查表模式）。
• 运行配置：full_coag_equations、dust_coag_kernel_b/g/de/ti、coal_coeff_mode、coal_fac 均来自 callphys.def（通过 conf_phys.F 读取）。

输出与副作用

• pdqcoag（唯一输出）：real(ngrid,nlayer,nq)，凝并倾向。仅尘埃类（mass 和 number）分量有非零值，其他分量保持初始化时的 0。
• 无文件 I/O、无诊断输出。
• 无模块变量写入：本例程只读模块变量。
• firstcall（第 132 行）：声明但未在 dust_coagulation_main 中使用——推断为 dust_coagulation_init 的 firstcall 已处理了初始化逻辑，本例程中的声明是残留代码。

伪代码

subroutine dust_coagulation_main(ngrid, nlayer, nq, ptimestep, pq, pdqfi, pt, pdtfi, pplay, pplev, pdqcoag):

  #--- 阶段 1: 初始化 ---
  sig0 = dev_dt; dens = rho_dust; cst = 0.75/(π*dens)*exp(-4.5*sig0²)
  zq = pq + pdqfi*ptimestep
  zt = pt + pdtfi*ptimestep
  rho = pplay / (r * pt)
  pdqcoag = 0

  for each dust mode iq in {dust_mass, stormdust_mass, topdust_mass}:
    mass_ini(iq) = max(zq(mass), 1e-15)
    numb_ini(iq) = max(zq(number), 1e-15)
    mtot += mass_ini(iq)
    r0 = (0.75*mass/(π*numb*dens))^(1/3) * exp(-1.5*sig0²)
    for each bin i:
      ndis_tab(i,iq) = numb*rho*deltar/(2*rad*√(2π)*sig0) * exp(-0.5*(ln(rad/r0))²/sig0²)
    numb_ini_dis(iq) = sum(ndis_tab(iq)) / rho

  rat_tab(iq) = normalize(ndis_tab(iq) / sum(ndis_tab))
  ndis = sum(ndis_tab(all_iq))
  ntot = sum(ndis) / rho

  #--- 阶段 2: 凝并（逐格点逐层）---
  for each (i,j):
    if ntot(i,j) > 1000:
      pt0 = pt(i,j); rho0 = rho(i,j); mfp0 = mfp(pt0, rho0)

      # 查表模式: findval + boxinterp 确定 kernel
      # 完整方程模式: 直接调 betab/betag/betade/betati

      for each bin k = 1..nres_coag:
        term1 = Σ_{ii<k, jj≤k} frac(ii,jj,k)*vol(ii)*K(ii,jj)*ndis_new(ii)*ndis(jj) * Δt/V_k
        term2 = Σ_j (1-frac(k,j,k))*K(k,j)*ndis(jj) * Δt
        ndis_new(k) = (ndis(k) + coal_fac*term1) / (1 + coal_fac*term2)

      # 体积守恒
      norm = Σ(ndis*vol) / Σ(ndis_new*vol)
      ndis_new *= norm

  #--- 阶段 3: 新矩量 ---
  for each dust mode iq:
    ndis_tab_new(iq) = ndis_new * rat_tab(iq)     # 按初始比值拆分
    numb_new(iq) = sum(ndis_tab_new(iq)) / rho
    mass_new(iq) = sum(ndis_tab_new(iq)*vol)*dens / rho
    mtotnew += mass_new(iq)

  for each dust mode iq:
    mass_new(iq) *= mtot / mtotnew                  # 质量守恒
    pdqcoag(mass, iq)   = (mass_new - mass_ini) / ptimestep
    pdqcoag(number, iq) = (numb_new - numb_ini_dis) / ptimestep
    where (zq + pdqcoag*ptimestep ≤ 0):             # 非负保护
      pdqcoag = -zq/ptimestep + small_value

复现注意事项

• nres_coag = 10（第 35 行，parameter）：bin 数量硬编码为 10。改变此值需同步修改 rads_coag/vols_coag/deltar_coag 的初始化逻辑（dust_coagulation_init 中）。
• bin 半径等比分布（dust_coagulation_init 第 72-77 行）：vrat_coag = (rn/r1)^(3/(nres-1))，rads(i) = r1 * vrat^((i-1)/3)。范围 r1=0.01e-6 到 rn=40e-6 m。
• dev_dt = 0.63676（第 53 行）：对数正态分布标准差，影响 cst 和 r0 计算。重新实现时必须使用此精确值。
• frac 函数的体积分数分配（第 613-638 行）：决定了凝并产物如何在 bin 间分配。公式含 vol(k)/vint 校正因子，容易实现错误。
• 半隐式格式（第 257/346 行）：增益项用 ndis_new(ii)（ii < k，已更新），损耗项分母隐式处理。这种格式比纯显式稳定但比纯隐式简单。
• coal_fac 因子（第 257/346 行）：全局缩放凝并效率。若设为 0 则凝并完全关闭。
• 凝并阈值 ntot > 1000（第 223/273 行）：低于此浓度不触发凝并。这是性能优化而非物理要求。
• 质量守恒双重保护：体积守恒归一化（阶段 2）+ 总质量缩放（阶段 3 第 385 行）。两个机制共同保证质量守恒。
• rat_tab 二次归一化（第 204-208 行）：先算比值再做归一化确保 Σ=1。若省略第二步，浮点误差会导致模式间质量泄漏。
• 非负保护下限不一致（第 392 vs 395 行）：mass 用 1e-15，number 用 1e-14。重新实现时需保留此差异（推断：number 分量需要略高的下限以防止数值不稳定）。
• 复现风险：firstcall 变量在 dust_coagulation_main 中声明但未使用（第 132-133 行）。若误用为初始化开关，可能导致首步不计算凝并。
• 复现风险：betati（湍流 kernel，第 591-599 行）标注 "not implemented yet"。若 dust_coag_kernel_ti = .true.，完整方程模式会调用但结果可能不可靠。查表模式中湍流 kernel 表 table_ti 从未被填充。
• 复现风险：查表模式的插值——Brownian/Gravitational 按 (T, mfp) 对数插值（log10(table_mfp)），Diffusion Enhancement 按 (T, P) 对数插值（log10(table_pres)）。插值轴的对数变换必须正确实现。

待确认

• firstcall（第 132 行）是否曾有用途——当前代码中未被引用。
• coal_fac 是否通过 callphys.def 可调——模块声明为 real :: coal_fac = 1.，未在 conf_phys.F 中找到对应 getin_p 调用（待确认）。
• 非负保护的 mass/number 下限差异（1e-15 vs 1e-14）是否有物理依据，还是数值调试产物。
• table_pres 的初始化公式（第 655 行 1.e-8 * 10^(j/2)）：从约 0.03 Pa 到约 3e4 Pa 的对数等距分布——具体覆盖范围和物理意图待确认。
• betati 中 eps 变量（第 594、598 行）——被声明但从未赋值，推断为湍流耗散率，但当前未实现。