newsedim_mod.F

快速理解

它做什么： 单示踪物重力沉降核心计算单元，被 callsedim（callsedim_mod.F:569/593/654）和 co2cloud 调用；用 Stokes 定律 + Cunningham 滑移修正计算下落速度，再交给 vlz_fi 做垂直输运。

基本过程： 输入粒子半径/密度/形状因子 → 计算 Stokes 速度 + Cunningham 修正 → 换算界面输运质量 w → 调用 vlz_fi 更新混合比 → 返回界面通量 wq。

关键结果： 更新的示踪物混合比 pqi 和界面通量 wq，供 callsedim 汇总为 pdqsed/pdqs_sed。

路径

LMDZ.MARS\libf\phymars\newsedim_mod.F

所属目录

libf/phymars

文件定位

newsedim_mod.F（全文 224 行，固定格式 Fortran）只定义一个例程 newsedim，是 LMDZ.MARS 中单示踪物重力沉降的核心计算单元。

它做两件事（源码依据：文件头注释第 13-18 行 "Compute sedimentation of 1 tracer of radius rd (m) and density rho (kg.m-3)"）：

用 Stokes 定律 + Cunningham 滑移修正计算每层的粒子下落速度 vstokes，再换算出本时间步内穿过每个层间界面的"被输运大气质量" w（kg·m⁻²）。
把 w 交给 vlz_fi（半拉格朗日垂直输运），由后者实际更新示踪物混合比 pqi 并返回界面通量 wq。

也就是说：下落速度公式在本文件，输运数值格式在 vlz_fi.F。上层调度（按示踪物类别准备粒径/密度）在 callsedim_mod.F 和 co2cloud_mod.F90。

定义的符号

符号	类型	行号	作用
`newsedim_mod`	module	1	容器模块，无模块级变量
`newsedim`	subroutine	7	单示踪物沉降：算下落速度 → 算输运质量 → 调 `vlz_fi`

依赖的模块

use 模块	only 列表	用途	待确认
`vlz_fi_mod`	`vlz_fi`	垂直输运数值格式（第 9 行）	无
`comcstfi_h`	`r`, `g`	比气体常数 r 与重力加速度 g（第 10 行）；`comcstfi_h.F90` 中 r、g 为运行时由物理初始化赋值的 SAVE 变量	无

调用的关键例程

被调用例程	所在模块/文件	调用位置	作用
`vlz_fi`	`vlz_fi_mod`（`vlz_fi.F`）	第 216 行	用 `w` 实施"pseudo upstream"垂直平流，更新 `pqi`，输出 `wq`；坡度限制参数固定传 `2.`

输入

输入	来源	类型/维度	单位	含义
`ngrid`	调用者	integer	-	水平列数
`nlay`	调用者	integer	-	垂直层数
`naersize`	调用者	integer	-	`rd` 数组长度：1 = 全场单一半径；`ngrid*nlay` = 每格点每层各自半径
`nrhosize`	调用者	integer	-	`rho` 数组长度，语义同上
`ptimestep`	调用者	real	s	物理时间步（或微物理子步）
`pplev`	调用者	real(ngrid,nlay+1)	Pa	层间界面气压
`masse`	调用者	real(ngrid,nlay)	kg（注释原文；推断：实为 kg·m⁻²，见"待确认"）	每层大气质量
`epaisseur`	调用者	real(ngrid,nlay)	m	每层厚度
`pt`	调用者	real(ngrid,nlay)	K	层中心温度
`rd`	调用者	real(naersize)	m	粒子半径
`rho`	调用者	real(nrhosize)	kg·m⁻³	粒子密度
`beta`	调用者	real	-	粒子形状修正因子（Murphy et al. JGR 1990）：1 球形、0.85 不规则、0.5 盘状（第 41-45 行注释）

输出

输出	去向	类型/维度	单位	含义
`pqi`	inout，调用者	real(ngrid,nlay)	如 kg/kg	示踪物混合比，被 `vlz_fi` 原地更新为沉降后的值
`wq`	out，调用者	real(ngrid,nlay+1)	注释写 ?/m-2（推断：kg·m⁻²，即本时间步穿过各界面的示踪物量）	界面示踪物通量；`wq(ig,1)` 即落到地表的量，调用者用 `wq(ig,1)/ptimestep` 得地表通量倾向

共享状态与副作用

firstcall（SAVE，THREADPRIVATE，第 53-55 行）：首调时计算 Stokes 常数 b = 2/9 · g/visc。
a, b（SAVE，THREADPRIVATE，第 73-75 行）：b 仅 firstcall 算一次；a 每次调用都重算（第 97-100 行在 firstcall 块之外），因为 a 要乘以本次调用传入的 beta。
物理常数硬编码在本文件：CO2 气体分子粘度 visc = 1.e-5 N·s·m⁻²、CO2 有效分子半径 molrad = 2.2e-10 m（第 68-70 行）。
无文件读写；无诊断输出（仅注释掉的 write）。

核心逻辑

按执行顺序（行号为源码行号）：

firstcall 初始化（第 82-93 行）：b = (2/9)·g/visc，即 Stokes 速度公式 Vstokes = b·ρ·r² 的系数。
平均自由程系数（第 97-100 行，每次调用）： a = 0.707·8.31 / (4π·molrad²·N_A)，再乘形状因子：a = a·beta。气体平均自由程为 λ = a·T/P（第 96 行注释）。
逐层逐列算下落速度（第 113-136 行）：
- 按 naersize/nrhosize 取标量或场值：rfall、rhofall。场值索引方式为 i = ngrid*(l-1)+ig，即把 (ngrid,nlay) 数组按列优先展平后取第 i 个元素。
- Stokes + Cunningham 滑移修正（第 130-131 行）： vstokes = b·ρ·r²·(1 + 1.333·(a·T/P)/r) 注释注明该修正按 Rossow (Icarus 36, 1-50, 1978)。
- 穿层时间 traversee = epaisseur/vstokes（第 134 行）。
算每个界面被输运的大气质量 w（第 145-215 行）：
- 下落距离 dztop = vstokes·ptimestep（第 148 行）。
- 简单法（第 161-171 行，总是先算）： w = (1 − exp(−dztop·g/(r·T)))·pplev/g 即静力学下 dztop 厚度内的气柱质量。数值保护（JF+AS 2005/11 注释）：当指数太小使 1−exp(−x) 圆整为 0 时，改用一阶泰勒展开 w = (dztop·g/(r·T))·pplev/g。
- 复杂法（第 176-201 行，仅当 dztop > epaisseur，即一个时间步穿越多层）：用 while 循环逐层下推（第 182-187 行）：累计已穿越厚度 Ep 和已耗时间 Stra，剩余时间用下一层的 vstokes(ig,l+k) 继续下落，直到 dztop ≤ Ep 或到达 l+k+1 > nlay。最终落点气压 ptop = pplev(ig,l+k)·exp(−(dztop−Ep)·g/(r·T(l+k))) （同样有指数过小时的泰勒保护，第 193-198 行），则 w(ig,l) = (pplev(ig,l) − ptop)/g。
调用输运格式（第 216 行）：call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)，原地更新 pqi，输出 wq。

伪代码

subroutine newsedim(..., rd, rho, pqi, wq, beta):
  if firstcall: b = (2/9)*g/visc
  a = 0.707*8.31/(4*pi*molrad^2*Avogadro) * beta

  for each layer l, column ig:
    rfall  = rd(1)  or rd(ngrid*(l-1)+ig)    # 按 naersize
    rhofall= rho(1) or rho(ngrid*(l-1)+ig)   # 按 nrhosize
    vstokes = b*rhofall*rfall^2 * (1 + 1.333*(a*T/P)/rfall)
    traversee = epaisseur/vstokes

  for each layer l, column ig:
    dztop = vstokes*ptimestep
    w = (1 - exp(-dztop*g/(r*T))) * pplev/g     # 数值下溢时用泰勒展开
    if dztop > epaisseur(l):                    # 一步穿多层
      沿 l+1, l+2, ... 累计厚度 Ep 与时间 Stra，
      用各层各自的 vstokes 推进剩余时间，
      ptop = pplev(l+k)*exp(-(dztop-Ep)*g/(r*T(l+k)))
      w = (pplev(l) - ptop)/g

  call vlz_fi(ngrid, nlay, pqi, 2., masse, w, wq)  # 更新 pqi，输出 wq

参与的主题流程

主题	参与方式
尘埃循环	`callsedim_mod.F` 第 569 行（doubleq 分箱，`beta=0.5`，密度 `rho_dust`）
水循环	`callsedim_mod.F` 第 593/599 行（水冰/CCN，`beta` 来自 run 参数 `ice_shape`，默认 0.75，源码依据 callsedim 第 341-342 行）
通用示踪物沉降	`callsedim_mod.F` 第 654 行（`beta=1.0`）
CO2 云	`co2cloud_mod.F90` 共 8 处调用（第 819-884 行区间），均为 `naersize=nrhosize=ngrid*nlay` 的场模式，`beta=0.85`（co2cloud 第 178 行参数常量），对 CO2 冰与各类 CCN 示踪物在微时间步内沉降

复现要点

下落速度公式必须同时含 Stokes 项与 Cunningham 滑移项：v = (2/9)(g/visc)·ρ·r²·(1 + 1.333·λ/r)，其中 λ = a·T/P，且 a 已乘形状因子 beta。注意 beta 进的是滑移修正项而不是 Stokes 项。
粘度 visc=1e-5、分子半径 molrad=2.2e-10 是 CO2 的硬编码值；重新实现时不能换成地球空气值。
naersize/nrhosize 的 1 vs ngrid*nlay 双模式必须保留，且场模式索引是 ngrid*(l-1)+ig 的展平约定——调用方传 rsedcloud/rhocloud 这类 (ngrid,nlay) 数组时依赖该内存布局。
两处 1-exp(-x) / exp(-x) 的小指数泰勒保护（第 167-171、193-198 行）影响极小粒子/极短时间步的结果，需保留。
跨层 while 循环上界条件 l+k+1.le.nlay：粒子最多推到模式顶下一层，超出部分按最后到达层的温度做静力外推。
w 的含义是"该界面在本时间步被示踪物穿越的大气质量（kg·m⁻²）"，正方向向下；真正的守恒输运与坡度限制（pente_max=2）在 vlz_fi 内完成。
复现风险：复杂法只在 dztop > epaisseur 时覆盖简单法的 w；两种方法在临界处不连续是否影响结果未验证。
复现风险：firstcall 只算 b 不算 a，若误把 a 也放进 firstcall（如早期版本结构），不同 beta 的多次调用会互相污染。

待确认

待确认：masse 的注释单位写 "kg"，但 vlz_fi.F 第 25 行同名参数注释为 "mass of atmospheric layer delta(P)/g"，推断实际单位为 kg·m⁻²（柱密度）。
待确认：wq 注释单位写 "?/m-2"，推断对 kg/kg 混合比示踪物而言是 kg·m⁻²。
待确认：0.707 即 1/√2（来自分子运动论平均自由程公式 λ = kT/(√2·π·d²·P) 的变形），此为推断；8.31 为通用气体常数，6.023e23 为 Avogadro 常数（源码字面值，精度较低）。
待确认：形状因子 beta 乘在平均自由程系数 a 上（而非直接乘速度）的物理依据，注释只给出 Murphy et al. JGR 1990 文献指引，未展开。