moldiff.F

路径

LMDZ.MARS\libf\aeronomars\moldiff.F

概述

moldiff_mod 是 LMDZ.MARS 的原版（legacy）分子扩散模块，在 65 km 以上热层区域求解 14 种中性物种的多组分分子扩散。它采用 Dickinson & Ridley (1972) 的 α 矩阵表述，经 LU 分解/求逆后转为隐式三对角系统求解混合比倾向，并在顶部用 Jeans 逃逸方程处理 H/H₂ 逃逸通量。

当前状态：无直接调用方。 thermosphere_mod.F 通过 moldiff_scheme 开关（默认 2）选择 moldiff_red（legacy scheme, moldiff_scheme==1）或 moldiff_MPF（MPF scheme, moldiff_scheme==2），两者均不引用 moldiff_mod。moldiff.F 为遗留死代码，changelog 记录了从活跃使用到逐步被 moldiff_red/moldiff_MPF 取代的完整历史。

定义符号

物种编号（moldiff 内部）

共享状态（SAVE + OMP THREADPRIVATE）

use 依赖

call 调用

moldiff 参数表

核心逻辑

firstcall 初始化（第 117–213 行）

1. 调用 moldiffcoeff(dij) 获取 14×14 二元扩散系数矩阵。
2. 逐一检查 14 个必需 tracer 的 igcm_* 索引；任一为 0 即 stop。
3. 建立 gcmind(1:14) 映射（本地物种索引 → GCM tracer 索引）。

逐列求解（第 222–459 行，do ig=1,ngrid）

1. 前向推算温度和混合比（第 224–249 行）：tt = pt + pdt·dt + pdteuv·dt + pdtconduc·dt；qq = pq + pdq·dt，钳位 ≥ 1e-30。
2. 标高和平均分子量梯度（第 233–255 行）：hp = -log(p(l+1)/p(l-1))；dmmeandz = (mmean(l+1)-mmean(l-1))/hp，边界层用 3 点差分。
3. 构造 α 矩阵（第 259–286 行，Dickinson & Ridley 1972）：
   • 对 13×13 子矩阵（排除 O），每对物种 (nn, n) 按 dij(nn,n)·ptfac（ptfac = (1e5/p)·(T/273)^1.75）计算扩散系数。
   • alf(nn,n) = qq(nn)/ntot/1.66e-27 · (1/(mmol(n)·d(nn)) - 1/(mmol(O)·d_O))。
   • 对角线修正 alfdiag。
   • b(l,n) = -(dmmeandz/mmean + mmol(n)/mmean - 1) 重力分离项。
4. LU 求逆 α 矩阵（第 291–310 行）：ludcmp_sp + lubksb_sp，得 alfinv(l,nn,n) = y(nn,n)/hh（hh = rnew·T/g 标高）。
   • 若 LU 失败（奇异矩阵），全列 pdqdiff = 0 并 return。
5. 三对角系数（第 316–334 行）：
   • del1 = hp·p/(g·dt)，del2 = (hp/2)²·p/(g·dt)。
   • aa/bb/cc 分别为下/主/上对角线系数，包含 alfinv 垂直梯度和 b 项。
6. Jeans 逃逸通量（第 339–376 行，Hunten 1973 eq. 5）：
   • H：wi(H) = √(2kT/(m·1.6605e-27)) / (2√π) · (1+pote)·exp(-pote)；pote = (R_mars+z)/(kT/(m·g))。
   • H₂：同公式。
   • flux(H/H₂) = qq(nz)·dens/m_H · wi，abfac(H/H₂) = 0（逃逸边界），其余物种 abfac = 1（扩散平衡）。
7. 三对角求解（第 385–435 行，逐物种）：
   • 组装 atri/btri/ctri/rtri（层 2 到 nz-1），含交叉物种耦合项。
   • 顶部边界（第 408–414 行）：逃逸物种用 flux；非逃逸物种用 abfac·ctri/(3-2·hp·b) 扩散平衡条件。
   • 底部边界（第 417 行）：完全混合 btri(1) += atri(1)。
   • tridag_sp 求解得 qnew(l,nn)，钳位 ≥ 1e-30。
   • 顶层 qnew(nz,nn) 由逃逸通量或扩散平衡外推。
   • 底层 qnew(1,nn) = qnew(2,nn)。
8. 写倾向（第 437–447 行）：仅在 zlocal > 65000 m 时写入：
   • pdqdiff(ig,l,n) = (qnew(l,n) - qq(l,n)) / ptimestep。
   • O 倾向由质量守恒反算：pdqdiff(ig,l,g_o) = -Σ_{n≠O} pdqdiff(ig,l,n)。

辅助例程

tridag_sp（第 467–492 行）

Thomas 算法三对角求解器。a/b/c 为下/主/上对角线，r 为右端向量，输出 u。主元 b(j)-a(j)·gam(j) == 0 时 stop。

LUBKSB_SP（第 498–531 行）

Numerical Recipes 风格 LU 回代。输入 A(NP,NP) 和 INDX 置换向量，原地修改 B。

LUDCMP_SP（第 537–616 行）

Numerical Recipes 风格 LU 分解 + 部分主元选取。奇异矩阵（行最大元素为 0）时设 ierr=1 并 return（不再 stop），由调用方决定后续处理。

伪代码

moldiff(ngrid, nlayer, nq, pplay, ..., pdqdiff):
    if firstcall:
        moldiffcoeff(dij)
        for each of 14 tracers: check igcm_* != 0
        build gcmind(1:14)
        firstcall = .false.

    for ig = 1, ngrid:
        tt = pt + (pdt + pdteuv + pdtconduc) * dt    # forward temperature
        qq = pq + pdq * dt                            # forward mixing ratios
        hp = -log(p(l+1)/p(l-1))                      # log-pressure thickness
        dmmeandz = d(mmean)/hp                         # mean molecular weight gradient

        for l = 1, nz:
            build alf(13×13) from dij, qq, mmol, ntot  # Dickinson-Ridley alpha matrix
            ludcmp_sp(alf) + lubksb_sp → alfinv        # invert alpha
            build b(l,n) gravity separation term

        build aa/bb/cc tridiagonal coefficients from alfinv gradients
        compute Jeans escape flux for H and H2 at TOA (Hunten 1973)

        for nn = 1, 13:                                # solve each species
            assemble tridiagonal system with cross-coupling
            apply escape BC (H/H2) or diffusive equilibrium BC (others) at top
            apply perfect mixing BC at bottom
            tridag_sp → qnew(:,nn)

        for l where z > 65 km:
            pdqdiff(l, n) = (qnew - qq) / dt          # tendency for each species
            pdqdiff(l, O) = -sum(pdqdiff(l, n≠O))     # O by mass conservation

调用关系

取代关系

thermosphere_mod.F（第 18–19 行）引入 moldiff_red_mod::moldiff_red（legacy scheme）和 moldiff_MPF_mod::moldiff_MPF（MPF scheme），通过 moldiff_scheme 开关（默认 2，第 67–68 行）选择。callmoldiff 在 conf_phys.F 读取（默认 .false.，MCD5/MCD6 配置开启 .true.，MCD6 设 moldiff_scheme=1），要求 callthermos=.true.。

与 moldiff_red/moldiff_MPF 的区别

待确认

1. moldiff.F 是否在某个未编译的测试路径中仍被引用（当前 libf 全目录 grep 无 use moldiff_mod 或 call moldiff）。
2. 与 moldiff_red 的具体算法差异——两者均标注 legacy，但 moldiff.F 的 α 矩阵维数固定为 13×13（排除 O），而 moldiff_red 使用动态维度。
3. 第 370–376 行注释块 "TEMPORAIRE: no escape for h and h2" 是否曾被启用（当前被注释，逃逸生效）。
4. dij(h,o) = 0.000114 与 moldiffcoeff.F 第 195 行的覆盖旧值 0.0000144（~10× 差异）的物理依据。

相关页面

• aeronomars/index：aeronomars 目录总览，本文件在其"分子扩散"分类下。
• moldiffcoeff：提供 dij(14,14) 二元扩散系数矩阵，由本文件 firstcall 调用。
• moldiffcoeff_red：简化变体扩散系数例程，供 moldiff_red/moldiff_MPF 使用。
• tracer_mod：提供 igcm_* 索引与 mmol 摩尔质量。
• conc_mod：提供 rnew/mmean 混合气体属性。
• comcstfi_h：提供重力常数 g。
• conduction：热层分子热传导，与分子扩散共同影响热层温度/成分剖面。
• molvis：热层分子黏性速度扩散，与分子扩散并列为热层三大输运过程。