dust_coagulation_mod.F90

路径

LMDZ.MARS\libf\phymars\dust_coagulation_mod.F90

所属目录/模块

libf/phymars

文件定位

dust_coagulation_mod.F90 定义尘埃粒子凝并方案。源码文件头注释写明该方案基于 Bertrand et al., 2022，作者为 T. Bertrand。它接收当前尘埃质量/数量混合比、温度、压力和已有物理倾向，把背景尘埃、rocket storm dust、mountain top dust 的连续对数正态粒径谱离散到 10 个半径 bin，在每个格点/层上用凝并核更新粒子数分布，再把结果折回到每个尘埃模式的质量和数量 tendency pdqcoag。

源码依据：dust_coagulation_mod.F90:1-16 定义模块与 public 初始化例程；dust_coagulation_mod.F90:93-400 定义主例程 dust_coagulation_main；physiq_mod.F:1477-1520 在 dust_coagulation 开关打开时调用本模块并把输出 tendency 加回 pdq。

定义的符号

符号	类型	行号	作用
`dust_coagulation_mod`	module	1	尘埃凝并模块容器。
`dust_coagulation_init`	subroutine	63	分配并初始化凝并半径/体积 bin；查表模式下构造凝并核查找表。
`dust_coagulation_main`	subroutine	93	主计算：由质量/数量示踪剂构造粒径分布，执行凝并，输出 mass/number tendency。
`dynvis`	function	403	用 Sutherland 公式计算 CO2 动力粘度。
`kinvis`	function	413	由动力粘度和密度计算运动粘度。
`thervel`	function	419	计算气体分子或尘埃粒子的热速度。
`mfp`	function	427	由温度和密度估计平均自由程。
`cun`	function	435	Cunningham slip-flow 修正。
`diff`	function	442	粒子扩散系数。
`lambda`	function	449	粒子平均自由程。
`delta`	function	459	Brownian transition regime 中的平均距离修正。
`fallv`	function	469	带 Cunningham 修正的粒子下落速度。
`reyn`	function	478	Reynolds 数。
`stokes`	function	485	以两粒子相对沉降速度定义的 Stokes 数。
`schmidt`	function	495	Schmidt 数。
`coal`	function	502	coalescence efficiency 公式；当前主流程未直接调用。
`betab`	function	520	Brownian coagulation kernel。
`betag`	function	540	gravitational coagulation kernel。
`betade`	function	561	Brownian diffusion enhancement kernel。
`betade_tabfunc`	function	575	构造查表时使用的 diffusion enhancement kernel。
`betati`	function	591	turbulent inertial kernel；源码注释标为 not implemented yet。
`betati_tabfunc`	function	601	turbulent inertial 查表函数；当前查表构造未填充 `table_ti`。
`frac`	function	613	把两个 bin 凝并后的中间体积分配到目标 bin。
`make_tables`	subroutine	644	构造温度、压力/平均自由程网格上的 kernel 查找表。
`boxinterp`	function	686	对查找表做二维线性插值。
`findval`	function	718	找到数组中第一个大于目标值的位置。

依赖的模块

use 模块	only 列表	用途	待确认
`comcstfi_h`	`pi`, `g`, `r`, `mugaz`	`pi` 用于体积/分布公式，`g` 用于落速，`r` 用于空气密度和 pressure kernel；`mugaz` 在本文件中未使用。	`mugaz` 是否为历史遗留依赖待确认。
`tracer_mod`	`igcm_dust_mass`, `igcm_dust_number`, `igcm_stormdust_mass`, `igcm_stormdust_number`, `igcm_topdust_mass`, `igcm_topdust_number`	标识三类尘埃的质量/数量示踪剂索引。	无。
`microphys_h`	`kbz`, `m0co2`	Boltzmann 常数和 CO2 分子质量，用于热速度、扩散和平均自由程。	无。

源码依据：dust_coagulation_mod.F90:9-13。

调用的关键例程

被调用例程	所在模块/文件	调用位置	作用
`make_tables`	本文件	`dust_coagulation_mod.F90:83`	当 `full_coag_equations=.false.` 时预计算 kernel 查找表。
`mfp`	本文件	`dust_coagulation_mod.F90:225,278,585,606`	由温度/密度或 pressure 派生平均自由程。
`betab`	本文件	`236,249,520-538` 等	Brownian kernel，直接模式和查表构造均使用。
`betag`	本文件	`237,250,540-559` 等	gravitational kernel，受 `coal_coeff_mode` 影响。
`betade` / `betade_tabfunc`	本文件	`238,251,575-589,675-677`	diffusion enhancement kernel；查表模式用 `betade_tabfunc` 预计算。
`betati`	本文件	`239,252,591-599`	turbulent kernel；源码标注未实现，且查表路径中对应调用被注释。
`boxinterp`	本文件	`300-341`	查表模式下对 kernel 表做二维插值。
`findval`	本文件	`281,285,288`	查表模式下定位温度、平均自由程、压力相邻格点。
`frac`	本文件	`240,254,316,343`	计算凝并后体积落入目标 bin 的比例。

输入

输入	来源	类型/维度	单位	含义
`ngrid`	`physiq`	integer	-	水平格点数。
`nlayer`	`physiq`	integer	-	垂直层数。
`nq`	`physiq`	integer	-	示踪剂数。
`ptime`	`physiq`	real	sol 或模型时间（待确认）	传入但本文件未使用。
`ptimestep`	`physiq`	real	s	物理时间步，用于把 tendency 转成当前状态并推进凝并。
`pq(ngrid,nlayer,nq)`	`physiq`	real	kg/kg 或 nb/kg	当前 advected tracer field。
`pdqfi(ngrid,nlayer,nq)`	`physiq`	real	tracer unit/s	上游物理过程已累积的 tracer tendency。
`pt(ngrid,nlayer)`	`physiq`	real	K	层中温度。
`pdtfi(ngrid,nlayer)`	`physiq`	real	K/s	上游温度 tendency；用于构造 `zt`，但后续未使用。
`pplay(ngrid,nlayer)`	`physiq`	real	Pa	层中压力，用于空气密度、查表 pressure 坐标。
`pplev(ngrid,nlayer+1)`	`physiq`	real	Pa	层界面压力；传入但本文件未使用。
`full_coag_equations`	`conf_phys`/模块 SAVE	logical	-	真：每步直接计算 kernel；假：用初始化时的查找表。
`dust_coag_kernel_b/g/de/ti`	`conf_phys`/模块 SAVE	logical	-	是否启用 Brownian、gravitational、diffusion enhancement、turbulent inertial kernel。
`coal_coeff_mode`	`conf_phys`/模块 SAVE	integer	-	gravitational kernel 中 coalescence efficiency 的选择模式。

配置依据：conf_phys.F:392-408 设置 dust_coagulation 默认 .false.、full_coag_equations 默认 .false.、Brownian kernel 默认 .true.、其他 kernel 默认 .false.、coal_coeff_mode=0。同一段源码没有看到 getin_p 读取 kernel 子开关，因此本文件页只记录源码默认值。

输出

输出	去向	类型/维度	单位	含义
`pdqcoag(ngrid,nlayer,nq)`	`physiq`	real	tracer unit/s	尘埃凝并产生的 tracer tendency。
`pdqcoag(:,:,igcm_dust_mass)`	`physiq` → `pdq`	real	kg/kg/s	背景尘埃质量 tendency。
`pdqcoag(:,:,igcm_dust_number)`	`physiq` → `pdq`	real	nb/kg/s	背景尘埃数量 tendency。
`pdqcoag(:,:,igcm_stormdust_*)`	`physiq` → `pdq` when `rdstorm`	real	kg/kg/s, nb/kg/s	rocket storm dust 的质量/数量 tendency。
`pdqcoag(:,:,igcm_topdust_*)`	`physiq` → `pdq` when `topflows`	real	kg/kg/s, nb/kg/s	mountain top dust 的质量/数量 tendency。

源码依据：dust_coagulation_mod.F90:116,387-395 生成 pdqcoag；physiq_mod.F:1484-1518 把背景尘埃、topdust、stormdust 的对应 tendency 加入 pdq。

共享状态与副作用

full_coag_equations、dust_coag_kernel_b/g/de/ti、coal_coeff_mode 是 SAVE 且 THREADPRIVATE 的模块变量，由 conf_phys.F 设置默认值并打印。
rads_coag、vols_coag、deltar_coag 在 dust_coagulation_init 中分配，后续主例程依赖这些数组。
查表模式下 table_pt、table_pres、table_mfp、table_b、table_g、table_de 在 make_tables 中填充；table_ti 声明但未填充。
dev_dt 和 rho_dust 在本模块中为 public，默认分别为 0.63676 和 2500.；dust_coagulation_main 用它们确定粒径谱宽度和尘埃密度。
本模块不直接做文件 I/O；副作用主要是分配模块数组、填充 lookup table、输出 pdqcoag。

复现风险：dust_coagulation_init 无 allocated() 保护，若被重复调用会再次 allocate(rads_coag/vols_coag/deltar_coag)；源码调用点 physiq_mod.F:806-808 表明它在初始化阶段由 dust_coagulation 开关保护调用。

核心逻辑

初始化粒径 bin
dust_coagulation_init 分配 10 个 bin，最小半径 r1_coag=0.01e-6 m，最大半径 rn_coag=40.e-6 m，体积比 vrat_coag=(rn/r1)^(3/(nres-1))。每个 bin 的半径和体积按几何级数生成，deltar_coag 存放直径宽度。查表模式下调用 make_tables。
把上游 tendency 并入当前状态
主例程中 zq = pq + pdqfi*ptimestep，zt = pt + pdtfi*ptimestep。源码后续计算空气密度与 kernel 时使用 pt 而不是 zt；待确认这是有意使用未加热后的温度，还是保留变量未使用。
构造各尘埃模式的初始对数正态数分布
对 igcm_dust_mass、igcm_stormdust_mass、igcm_topdust_mass 三类 mass tracer，源码假设对应 number tracer 为 iq-1。质量和数量都设下限 1e-15。几何均值半径：
```
r0 = (3/4 * mass / (pi * number * rho_dust))^(1/3) * exp(-1.5*dev_dt^2)
```
每个 bin 的数浓度分布：
```
ndis_bin = number * rho_air * deltar
           / (2*r_bin*sqrt(2*pi)*dev_dt)
           * exp(-0.5*(log(r_bin/r0))^2/dev_dt^2)
```
ndis_tab 单位按源码变量关系推断为 nb/m3/bin。
合并多个尘埃模式并保存初始比例
ndis=sum(ndis_tab, mode) 得到所有尘埃模式合并后的总 bin 分布。rat_tab 保存每个模式在每个 bin 中的初始占比，并归一化，使凝并后的总分布能再按初始比例分回背景尘埃、stormdust、topdust。
执行凝并
仅当 ntot=sum(ndis)/rho_air > 1000 时启用凝并。两条路径：
- full_coag_equations=.true.：逐格点/层/目标 bin 直接调用 betab/betag/betade/betati 计算 kernel。
- full_coag_equations=.false.：用 findval 定位查表角点，再用 boxinterp 插值 table_b/table_g/table_de。
每个目标 bin k 计算两个项：
```
term1 = sum_{jj<=k, ii<k} frac(ii,jj,k) * vols(ii) * kernel(ii,jj)
        * ndis_new(ii) * ndis(jj) * dt / vols(k)

term2 = sum_j (1-frac(k,j,k)) * kernel(k,j) * ndis(j) * dt

ndis_new(k) = (ndis(k) + coal_fac*term1) / (1 + coal_fac*term2)
```
随后按 sum(ndis*vols)/sum(ndis_new*vols) 归一化，保持总体积（质量）守恒。
把凝并后的总分布折回各尘埃模式
ndis_tab_new = ndis_new * rat_tab。每个模式的新 number mixing ratio 为 sum(ndis_tab_new)/rho_air，新 mass mixing ratio 为 sum(ndis_tab_new*vols)*rho_dust/rho_air。然后用 mass_new/mtotnew*mtot 再次校正各模式质量，保持总质量。
输出 tendency 并做非负保护
mass tendency 为 (mass_new-mass_ini)/dt；number tendency 为 (numb_new-numb_ini_dis)/dt。若加上 tendency 后质量或数量非正，则用 -zq/dt + epsilon 保护到小正值。

伪代码

init:
  allocate 10 coagulation bins between 0.01 um and 40 um
  compute bin radii, volumes, diameter widths
  if using lookup tables:
    build temperature, mean-free-path/pressure grids
    fill Brownian, gravitational, diffusion-enhancement kernels

main:
  zq = pq + pdqfi * dt
  rho_air = pplay / (r * pt)
  initialize pdqcoag = 0

  for each dust mass tracer in {dust, stormdust, topdust}:
    mass_ini = max(zq(mass), 1e-15)
    number_ini = max(zq(number index = mass index - 1), 1e-15)
    mtot += mass_ini
    r0 = moment radius from mass/number/rho_dust/dev_dt
    for each coag bin:
      ndis_tab(bin, mode) = lognormal number distribution
    number_ini_dis = sum(ndis_tab over bins) / rho_air

  rat_tab = ndis_tab / sum(ndis_tab over dust modes)
  ndis = sum(ndis_tab over dust modes)
  ntot = sum(ndis over bins) / rho_air

  for each column and layer:
    if ntot <= 1000: keep ndis_new at initialized value
    else:
      for each target bin:
        kernel = selected Brownian + gravitational + diffusion + turbulent terms
        compute gain term from smaller-bin collisions
        compute loss term from target-bin collisions
        ndis_new = (old + coal_fac*gain) / (1 + coal_fac*loss)
      normalize ndis_new volume to initial ndis volume

  for each dust mode:
    ndis_tab_new = ndis_new * initial mode/bin fraction
    number_new = sum(ndis_tab_new) / rho_air
    mass_new = sum(ndis_tab_new * volume) * rho_dust / rho_air
  rescale mass_new so total mass equals initial total mass

  pdqcoag(mass) = (mass_new - mass_ini) / dt
  pdqcoag(number) = (number_new - number_ini_dis) / dt
  protect against negative post-tendency values

参与的主题流程

主题	参与方式
尘埃循环	可选微物理阶段；减少/重分布 dust number，保持质量守恒，并影响后续由 mass/number 反算的粒径。
辐射	推断：凝并改变 number tendency，从而影响 doubleq 粒径和后续辐射有效半径；本文件不直接调用辐射模块。
rocket storm dust	当 `rdstorm` 打开时，`physiq_mod.F:1507-1518` 把 `igcm_stormdust_*` 的 `pdqcoag` 加回 `pdq`。
mountain top dust	当 `topflows` 打开时，`physiq_mod.F:1494-1505` 把 `igcm_topdust_*` 的 `pdqcoag` 加回 `pdq`。

写法特点

模块变量大量 SAVE + THREADPRIVATE，说明该方案按 OpenMP 线程保存独立状态。
public 声明只显式列出 dust_coagulation_init，但模块变量和 dust_coagulation_main 未设 private，因此默认仍可被 use 访问；physiq_mod.F:27-28 只导入 dust_coagulation_main 和 dust_coagulation_init。
dust_coag_kernel_ti 注释为 Not implemented yet；直接计算路径可调用 betati，但 betati 中 eps 未赋值，查表路径的 turbulent 调用被注释，table_ti 也未填充。
ptime、pplev、局部 zt、rn、若干局部变量声明后未在有效计算中使用。待确认：这些是否为预留接口或历史遗留。
查表路径中 p2=findval(table_pres,pres0) 后若 p2==0 则设为 table_numt，而不是 table_nump。待确认：由于两者当前都等于 25，数值上不表现差异，但复现时不应把它泛化为不同维度。

复现要点

必须先调用 dust_coagulation_init，否则 rads_coag/vols_coag/deltar_coag 和查表数组没有初始化。
该方案依赖 mass tracer 与 number tracer 的索引相邻关系：代码用 iq-1 读取 number tracer。复现或重排 tracer 表时必须保持 *_number 紧邻并位于对应 *_mass 前一位。
默认配置只打开 Brownian kernel，且 dust_coagulation 本身默认关闭；本仓库 deftank/datadir/doc 中未找到显式 dust_coagulation 配置样例。
full_coag_equations=.false. 时，kernel 使用初始化阶段的查表和双线性插值；=.true. 时每步直接计算。
ntot <= 1000 nb/kg 的格点/层不会进入凝并更新。复现风险：源码未显式把这些层的 ndis_new 设为初始 ndis，后续仍参与折回 moments；该行为需用运行测试确认是否符合预期。
质量守恒通过两层归一化实现：先对总 bin 体积分布归一化，再对各模式 mass_new 按初始总质量 mtot 重缩放。
pdqcoag 的 number tendency 用离散化后的 numb_ini_dis 而非原始 numb_ini 作基准；这是避免 bin 离散误差立即进入 tendency 的关键。

待确认

dust_coag_kernel_*、full_coag_equations、coal_coeff_mode 在当前源码中只有默认赋值和打印，未见 getin_p 读取；是否在其他分支或运行框架中覆盖，待确认。
zt=pt+pdtfi*dt 已计算但未用于 kernel；是否应使用更新后温度待确认。
betati/betati_tabfunc 中 eps 未赋值，且源码注释 turbulent kernel 未实现；不要在复现中启用 dust_coag_kernel_ti，除非先补全物理定义。
coal 函数当前未被 betag 直接调用；betag 使用 coal_coeff_mode 的三种简化效率。完整 coalescence efficiency 公式的设计用途待确认。
ntot<=1000 时 ndis_new 保持初始化零值，后续 mass/number 折回的结果需数值测试验证。