newsedim_mod.F

路径

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)"）：

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

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

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 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）。

核心逻辑

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

1. firstcall 初始化（第 82-93 行）：b = (2/9)·g/visc，即 Stokes 速度公式 Vstokes = b·ρ·r² 的系数。
2. 平均自由程系数（第 97-100 行，每次调用）： a = 0.707·8.31 / (4π·molrad²·N_A)，再乘形状因子：a = a·beta。 气体平均自由程为 λ = a·T/P（第 96 行注释）。
3. 逐层逐列算下落速度（第 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 行）。
4. 算每个界面被输运的大气质量 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。
5. 调用输运格式（第 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

参与的主题流程

复现要点

• 下落速度公式必须同时含 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 文献指引，未展开。

相关页面

• callsedim_mod.md：上层调度器，准备粒径/密度并分类调用本例程。
• vlz_fi.md：本例程调用的垂直输运内核。
• dust-cycle.md：尘埃沉降主题流程。