moldiff_red.F90

路径

LMDZ.MARS\libf\aeronomars\moldiff_red.F90

所属目录/模块

libf/aeronomars

文件定位

简化版（"red"）多组分分子扩散求解器，是 moldiff.F（legacy，死代码）的现代化精简版本。moldiff_red_mod 对 16 种扩散物种逐物种独立做隐式三对角扩散求解，含 Jeans H/H₂/D 逃逸和行星球积分通量输出。与 moldiff_MPF.F90 共享同一 thermosphere_mod 调用入口，通过 moldiff_scheme 开关切换（moldiff_scheme==1 选本模块，==2 选 MPF，默认 2）。July 2014 JYC 加入弹道输运耦合以计算 H/H₂ 的上边界逃逸速度。

定义的符号

模块参数

SAVE / THREADPRIVATE 状态

全部标记 !$OMP THREADPRIVATE。

16 种扩散物种（ListeDiff）

co2, o, n2, ar, co, h2, h, d2, hd, d, o2, h2o_vap, o3, n, he, hdo_vap

与 moldiff_MPF 的 ListeDiff 完全相同。

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 读 tracer_mod::noms/mmol：firstcall 物种匹配和持续分子量计算。
• 读 geometry_mod::cell_area：逃逸通量面积加权。
• 写 planetwide_sumval 全局和：输出 PhiEscH/H2/D。
• SAVE/THREADPRIVATE 缓存数组在 firstcall 分配后持续存在。
• 无数值文件 I/O。

核心逻辑

1. firstcall：匹配 16 种 ListeDiff 中存在于 noms 的物种，记 ncompdiff、gcmind、i_h/i_h2/i_d；找 il0（pplay>Pdiff 的最高层+1）并广播；调 moldiffcoeff_red 得 dij；分配所有缓存数组。

2. 逐列循环（ig=1,ngrid）：

   • 温度前向推进：tt=pt+pdt*ptimestep（内联，不调 TMNEW）。
   • 混合比前向推进：qq=pq+pdq*ptimestep，下限 1d-30（内联，不调 QMNEW）。
   • 大气结构：HSCALE 压力标高 → MMOY 平均分子量 → DMMOY 梯度 → ZVERT 高度 → RHOTOT 密度。
   • 记录扩散层以上的初始柱质量 Mtot1。
   • 安全检查：Zmax>4000 km 则 stop。
   • 精细网格：nlraf=max((Zmax-Zmin)/dzres+1, 40)，调 UPPER_RESOL 插值到均匀 dz 网格。
   • GCMGRID_P 回插检查质量，CORRMASS 计算 FacMass=qq/qint。
   • 记录精细网格初始质量 Mraf1。
   • qnew=qq（il0 以下无扩散效果）。

3. Jeans 逃逸速度：

   • Uthermal=sqrt(2·kB·T_top/m)，Lambdaexo=m·G·M_mars/(R_mars+Z_top)/kB/T_top。
   • 对 H/H₂/D：wi=Uthermal/2/sqrt(pi)*exp(-Lambda)*(Lambda+1)。
   • 其他物种：wi=0。
   • 注意：无热扩散 αT 修正（MPF 版加了 αT=−0.25）。

4. 逐物种扩散参数：

   • DCOEFF：有效扩散系数 D=(1-Nk/N)/Σ(xj/Dij)，用 dij 二元系数和 (P0/P)·(T/T0)^1.75 缩放。
   • HSCALEREAL：密度标高 H=-1/(dlnNk/dz)。
   • VELVERT：扩散速度 W=-D·(1/H-1/Hmol-dT/T)。
   • TIMEDIFF：扩散时间 |H/W|。
   • 取全局 Tdiff=min(Tdiffrafmol)，下限 tdiffmin·Mraf(nlraf)。

5. 时间分裂：ntime=round(ptimestep/tdiff)。

   • 逐物种、逐时间步：
     • 无量纲化：H0=kB·T0/m/g，D0=Draf(nlraf)，Time0=H0²/D0。
     • DIFFPARAM：差分系数 alpha/beta/gama/delta/eps。
     • MATCOEFF：三对角矩阵 A/B/C/D。底边界 Dirichlet，顶边界扩散廓线含 FacEsc=exp(-w·dz/D)。
     • tridag：Thomas 算法求解 → Xtri → 更新 Rrafk(l,nn)=rho0·Xtri(l)。
     • NaN/负值检查：若 Rrafk 异常则 stop（对 nn=16 回退到初始值）。

6. 更新总密度和混合比：Rraf=ΣRrafk，Qraf/Rraf/Nrafk/Mraf/Nraf/Praf 依次更新。

7. 质量诊断：记录 Mraf2（扩散后精细网格柱质量）。

8. 回插 GCM 网格：GCMGRID_P2 用 FacMass 修正，顶层以上用标高外推。RHOTOT 更新密度。

9. 逃逸通量：PhiauxH/H2/D(ig)=wi·Nrafk(nlraf)·cell_area(ig)，ig 循环后调 planetwide_sumval 全球求和。

10. 扩散 tendency：pdqdiff(ig,l,gcmind(nn))=(qnew(l,nn)-qq(l,nn))/ptimestep。

11. 释放精细网格数组（每列释放，每列重新分配）。

伪代码

if firstcall:
  match 16 species against tracer noms → ncompdiff, gcmind, i_h/i_h2/i_d
  find il0: highest layer where pplay > Pdiff, +1, broadcast
  call moldiffcoeff_red → dij(ncompdiff,ncompdiff)
  allocate all SAVE arrays
  firstcall = .false.

for each column ig:
  tt = pt + pdt * ptimestep  (inline)
  qq = pq + pdq * ptimestep  (inline, floor 1d-30)

  atmospheric structure: HSCALE → MMOY → DMMOY → ZVERT → RHOTOT
  record Mtot1 (column mass above il0)
  stop if Zmax > 4000 km

  nlraf = max((Zmax-Zmin)/dzres + 1, 40)
  allocate fine-grid arrays
  UPPER_RESOL: interpolate to uniform dz grid
  GCMGRID_P: back-interpolate for mass check
  CORRMASS: FacMass = qq / qint
  record Mraf1 (fine-grid mass)

  qnew = qq  (no diffusion below il0)

  Jeans escape: Uthermal, Lambdaexo → wi for H/H2/D (no alphaT correction)

  for each species nn:
    DCOEFF → Draf
    HSCALEREAL → Hraf
    VELVERT → Wraf
    TIMEDIFF → Tdiffraf
  Tdiff = min(Tdiffrafmol), floor tdiffmin * Mraf(nlraf)
  ntime = round(ptimestep / Tdiff)

  for istep = 1 to ntime:
    for each species nn:
      non-dimensionalize: H0, D0, Time0
      DIFFPARAM → alpha/beta/gama/delta/eps
      MATCOEFF → A/B/C/D  (Dirichlet bottom, diffusion-profile top with FacEsc)
      tridag → Xtri
      Rrafk = rho0 * Xtri  (NaN/negative check)
    update total: Rraf, Qraf, Nrafk, Mraf, Nraf, Praf

  record Mraf2 (post-diffusion mass)
  GCMGRID_P2: back-interpolate to GCM grid with FacMass
  RHOTOT: update densities

  PhiauxH/H2/D(ig) = wi * Nrafk_top * cell_area(ig)

  pdqdiff = (qnew - qq) / ptimestep

  deallocate fine-grid arrays

planetwide_sumval(PhiauxH → PhiEscH)
planetwide_sumval(PhiauxH2 → PhiEscH2)
planetwide_sumval(PhiauxD → PhiEscD)

参与的主题流程

写法特点

• JulY 2014 JYC 注释：标题注释提到 "ADD BALISTIC Transport coupling to compute wup for H and H2"，但实际源码中 Jeans 逃逸公式是标准形式，未见弹道输运的显式耦合。
• TMNEW/QMNEW 遗留：这两个子例程定义在模块内但主例程用内联代码替代，保留了源码注释中的 CALL TMNEW/QMNEW 行。
• tridagbloc 遗留：块三对角求解器保留在模块中但未被主例程调用（pre-refactor 残留）。
• DCOEFF 简化：用 (1-Nk/N)/interm 而非 Wilke 完整多组分公式（MPF 版用 Wilke），注释中有 "Temporary: eliminate modification to include Wilke's formulation"。
• 无涡流扩散：与 MPF 版不同，本模块不包含 Krasnopolsky 涡流扩散廓线。
• 无热扩散：Jeans 逃逸不含 αT=−0.25 热扩散修正（MPF 版有）。
• 无质量守恒修正器：MPF 版有可选 call_mass_fixer_moldiff_MPF，本模块无对应功能。
• 硬编码常数：g=3.72d0、Rmars=3390000d0、Mmars=6.4d23、pi=3.141592653、ij0=6000（测试列号）均在子例程内部重复定义。
• MESOSCALE 条件编译：planetwide_mod 的 use 和 planetwide_sumval 调用在 #ifndef MESOSCALE 下，mesoscale 模式打印警告。
• step 变量：firstcall 中设为 1 后未被使用。

与 moldiff.F / moldiff_MPF.F90 对比

调用方

复现要点

• moldiff_scheme 在 thermosphere_mod.F 中默认为 2（选 MPF），需设为 1 才能执行本模块。
• 依赖 moldiffcoeff_red 提供 dij，需确保 tracer 列表中包含 ListeDiff 中的物种。
• 需要 MCD5/MCD6 配置开启 callmoldiff（conf_phys.F）。
• Pdiff=15 Pa 决定了扩散起始高度，典型火星条件下约 65–80 km。
• 逃逸通量输出依赖 planetwide_sumval，mesoscale 模式下不可用。

待确认

• 标题注释提到 "BALISTIC Transport coupling" 但实际 Jeans 逃逸公式是标准形式，未见弹道输运的显式实现，可能指 Uthermal/Lambdaexo 的推导方式。
• DCOEFF 注释 "Temporary: eliminate modification to include Wilke's formulation" 暗示曾经尝试加入 Wilke 公式后又回退，具体原因未说明。
• step 变量 firstcall 设 1 后未使用，可能为遗留。
• ij0=6000 硬编码测试列号，实际 ngrid 可能远小于此值，诊断输出永远不会触发。
• Checkmass/Checkmass2 被注释掉，质量守恒诊断当前不活跃。

相关页面

• moldiff：原版 legacy 14 物种分子扩散模块（死代码）。
• moldiff_MPF：MPF 增强版，默认 moldiff_scheme==2。
• moldiffcoeff_red：简化扩散系数矩阵，由本模块 firstcall 调用。
• moldiffcoeff：原版扩散系数矩阵，供 legacy moldiff 使用。
• aeronomars/index：aeronomars 目录总览，本文件在其"分子扩散"分类下。
• tracer_mod：提供 noms/mmol/nqmx。