moldiff_MPF.F90

路径

LMDZ.MARS\libf\aeronomars\moldiff_MPF.F90

所属目录/模块

libf\aeronomars

文件定位

moldiff_MPF_mod 是 LMDZ.MARS 的当前默认（moldiff_scheme==2）分子扩散模块，在 Pdiff=15 Pa 以上热层区域求解动态数量（最多 16 种）中性物种的分子扩散。它基于 Modified Pass Flow (MPF) 方法（Parshev et al. 1987，2023 年 6 月 JYC 更新），逐物种独立求解含分子扩散、涡旋扩散和热扩散的隐式一维垂直输运方程，并在顶部用 Jeans 逃逸处理 H/H₂/D 逃逸通量。与 moldiff_red（legacy scheme）的主要区别是：可选柱质量修正器（call_mass_fixer_moldiff_MPF）和 MPF 算法框架。由 thermosphere_mod.F 在 callmoldiff 且 moldiff_scheme==2 时调用。

定义的符号

扩散物种列表（ListeDiffNb=16）

模块级参数

共享状态（SAVE + OMP THREADPRIVATE）

use 依赖

call 调用

moldiff_MPF 参数表

共享状态与副作用

• call_mass_fixer_moldiff_MPF 为 module 级 SAVE 变量，由 conf_phys.F（第 1258–1264 行）通过 getin_p("call_mass_fixer_moldiff_MPF", ...) 读取。
• firstcall 翻转后 dij、gcmind 和可选 masscorrfac 被 SAVE 缓存，后续调用跳过初始化。
• 所有 SAVE 变量均有 !$OMP THREADPRIVATE 声明，支持 OpenMP 并行。
• 每次 firstcall 打印 'specie X diffused in moldiff_MPF' 和 'number of diffused species: N'。
• 非 MESOSCALE 编译时，调用 planetwide_sumval 对逃逸通量做全行星求和。
• MESOSCALE 编译时，打印 "no planetwide in MESOSCALE" 并跳过求和。
• Zmax > 4000 km 或 NaN 时 stop。
• GCMGRID_P2 中的外推循环设 compteur >= 100000 上限防无限循环。

核心逻辑

firstcall 初始化（第 150–267 行）

1. 构造 16 物种 ListeDiff 列表，逐 tracer 按 noms 名匹配，统计实际存在的 ncompdiff 种。
2. 记录 i_h/i_h2/i_d/i_hd 在扩散列表中的位置（默认 1000 表示不存在）。
3. master 上找 il0（pplay > 15 Pa 之上第一层）和可选 il1（pplay > 1 Pa 之上第一层），再 bcast。
4. 调用 moldiffcoeff_red(dij, indic_diff, gcmind, ncompdiff) 获取二元扩散系数矩阵。
5. 分配所有 SAVE 工作数组（qq/qnew/qint/FacMass/rhoK/rhokinit 和可选质量修正数组 Mtot1-3/Mhi1-2）。

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

1. 前向推算（第 304–323 行）：tt = pt + pdt·dt（NaN 检查）；qq = pq + pdq·dt（max(1d-30) 钳位）。
2. 大气结构（第 327–344 行）：HSCALE 气压标高、MMOY 平均分子质量、DMMOY 梯度、ZVERT 静力高度、RHOTOT 密度。
3. 质量守恒前测（第 352–369 行，仅 call_mass_fixer）：柱积分质量 Mtot1 和 1 Pa 以上质量 Mhi1。
4. 安全检查（第 377–384 行）：Zmax > 4000 km 则 stop。
5. 高度网格构造（第 390–415 行）：nlraf = (Zmax-Zmin)/dzres + 1（最少 40 层），UPPER_RESOL 将 GCM 气压层插值到等间距高度网格（Rrafk 用幂律插值保持正性）。
6. 涡旋扩散（第 422–426 行）：Krasnopolsky 2002 方案 K = K0·sqrt(T·T_TOA·kB/P)，K0=1.2e11。
7. 质量修正因子（第 430–435 行）：GCMGRID_P 回插→CORRMASS 得 FacMass = qq/qint。
8. Jeans 逃逸速度（第 487–501 行）：对 H/H₂/D，wi = Uthermal/(2√π)·exp(-Λ)·(Λ+1)；热扩散系数 alphaT=-0.25（H/H₂/D/HD）。
9. 扩散系数和时间步（第 513–541 行）：逐物种 DCOEFF（Wilke 多组分公式）、HSCALEREAL、VELVERT（含热扩散项）、TIMEDIFF。最终 tdiff = ptimestep/5，ntime = ptimestep/tdiff = 5。
10. 时间循环（第 552–678 行，do istep=1,ntime）：
    • 逐物种非无量纲化：H0 = kT₀/(m·g)，D0 = Draf(TOA)，Time0 = H0²/D0。
    • DIFFPARAM：构造 delta/ksi/eps/zeta/prod/loss 含温度梯度和热扩散。
    • SEQUENCY：从顶到底递推 alpha/beta。
    • 从底到顶正向求解 Xtri(l) 得新密度分布。
    • NaN/负值检查（CO₂ 负值 stop）。
    • 更新 Rrafk、总密度、混合比、气压。
11. 回插到 GCM 网格（第 693–694 行）：GCMGRID_P2 含 FacMass 修正。
12. 逃逸通量（第 702–706 行）：Phiaux(ig) = wi·Nrafk(nlraf)·cell_area(ig)。
13. 质量修正器（第 724–806 行，仅 call_mass_fixer）：
    • 柱积分 Mtot2/Mhi2（扩散后）。
    • 非 H/N/O 物种：masscorrfac = (Mtot1 - Mhi2)/(Mtot2 - Mhi2)，缩放 il1 以下混合比。
    • O 损失 → 补偿 O₂（masscorrfac1）；H → H₂；N → N₂。
14. 写倾向（第 816–848 行）：
    • call_mass_fixer：非 CO₂ 物种 pdqdiff = (qnew-qq)/ptimestep，CO₂ 倾向 = -Σ(其他)。
    • 传统：所有物种直接 pdqdiff = (qnew-qq)/ptimestep。

全行星逃逸通量汇总（第 867–873 行）

非 MESOSCALE：planetwide_sumval(PhiauxH, PhiEscH) 等。MESOSCALE 跳过。

伪代码

moldiff_MPF(ngrid, nlayer, nq, ..., pdqdiff, PhiEscH, PhiEscH2, PhiEscD):
    if firstcall:
        build ListeDiff(16), match present species → ncompdiff
        find i_h, i_h2, i_d, i_hd
        find il0 (P > 15 Pa boundary), il1 (P > 1 Pa boundary)
        moldiffcoeff_red(dij)
        allocate SAVE arrays
        firstcall = .false.

    PhiEscH = PhiEscH2 = PhiEscD = 0

    for ig = 1, ngrid:
        tt = pt + pdt * dt                        # forward temperature
        qq = pq + pdq * dt                        # forward mixing ratios, clip ≥ 1e-30
        HSCALE, MMOY, DMMOY, ZVERT, RHOTOT        # atmospheric structure
        if mass_fixer: Mtot1, Mhi1                 # pre-diffusion column mass
        nlraf = max(40, (Zmax-Zmin)/dzres+1)       # refined altitude grid
        UPPER_RESOL → Praf/Traf/Qraf/Zraf/...      # interpolate to equal-dz grid
        Kraf = K0 * sqrt(T*T_TOA*kB/P)             # eddy diffusion (Krasnopolsky 2002)
        GCMGRID_P + CORRMASS → FacMass             # mass correction factors

        for nn = 1, ncompdiff:
            DCOEFF (Wilke), HSCALEREAL, VELVERT, TIMEDIFF
        tdiff = ptimestep / 5; ntime = 5

        for istep = 1, ntime:
            for nn = 1, ncompdiff:
                non-dimensionalize (H0, D0, Time0)
                DIFFPARAM → delta/ksi/eps/zeta/prod/loss
                SEQUENCY → alpha/beta (top-down recurrence)
                Xtri = forward solve (bottom-up)     # new density profile
                update Rrafk(l,nn) = rho0 * Xtri(l)
            update Rraf, Qraf, Mraf, Praf

        GCMGRID_P2 → qnew (back to GCM pressure grid, with FacMass)
        PhiauxH/H2/D(ig) = wi * Nrafk(nlraf) * cell_area(ig)

        if mass_fixer:
            Mtot2, Mhi2                             # post-diffusion column mass
            scale qnew(1:il1) for non-H/N/O species  # column mass conservation
            compensate O2/H2/N2 below 1 Pa for O/H/N losses
            pdqdiff(non-CO2) = (qnew-qq)/dt
            pdqdiff(CO2) = -sum(other tendencies)    # ensure total = 0
        else:
            pdqdiff = (qnew-qq)/dt for all species

    if not MESOSCALE:
        planetwide_sumval(PhiauxH/H2/D → PhiEscH/H2/D)

参与的主题流程

写法特点

• 自由格式 F90（& 续行），但部分内部辅助例程使用固定格式风格。
• 大量 SAVE + OMP THREADPRIVATE 工作数组，避免每次调用重新分配。
• 扩散物种数动态（按 noms 名匹配），而非像 moldiff.F 那样固定 14 种。
• 使用 moldiffcoeff_red（简化变体）而非 moldiffcoeff（原版），gcmind 由本文件构造后传入。
• tdiff = ptimestep/5 硬编码（第 541 行），前面计算的物理扩散时间 Tdiff 被覆盖。注释保留多种 K0 备选值（6e10 到 1.2e12）。
• 涡旋扩散 K0=1.2e11（第 282 行），注释中保留 Krasnopolsky 2002 原始值和其他测试值。
• TMNEW 和 QMNEW 定义但未调用——主例程内联了等价逻辑。
• tridagbloc 和 tridag 定义但主例程未调用——MPF 方案使用 SEQUENCY 递推 + 显式前向求解替代了三对角直接求解。
• GCMGRID_P 和 GCMGRID_P2 功能相似，但 P2 版本额外应用 FacMass 并在 call_mass_fixer 路径中按 sum(qq(il,:)) 归一化。
• 外推区使用气压标高指数衰减逐层推进（do while (pnew2 >= pp(il))），设 100000 次上限防无限循环。

复现要点

• conf_phys.F 读取 call_mass_fixer_moldiff_MPF（默认 .false.）。设为 .true. 启用质量修正器。
• thermosphere_mod.F 通过 moldiff_scheme（默认 2）选择本文件。callmoldiff 须为 .true.（由 deftank 配置设定）。
• 需要 tracer_mod 中存在 ListeDiff 中的物种名（如 'co2','h2','h','d','hd' 等）才能正确设置 i_h/i_h2/i_d/i_hd。缺失不影响扩散计算，仅跳过逃逸通量。
• moldiffcoeff_red 要求 noms 中存在 'h2'/'h'/'o'（由 moldiffcoeff_red 内部查找）；缺失会导致 mmol 越界。
• cell_area 来自 geometry_mod，用于逃逸通量面积加权。
• 非 MESOSCALE 编译时依赖 planetwide_mod::planetwide_sumval。

调用关系

deftank 配置

与 moldiff.F/moldiff_red 的区别

待确认

1. tdiff = ptimestep/5 硬编码（第 541 行）覆盖了前面计算的物理扩散时间 Tdiff，是否有物理依据或仅为数值稳定性。
2. K0=1.2e11 值相对 Krasnopolsky 2002 原文的适用范围——注释中列出 6e10 到 1.2e12 多个备选值，最终选择依据。
3. tridagbloc（块三对角求解器）和 tridag（标准三对角求解器）已定义但主例程未调用——是否为 MPF 重构前的遗留代码。
4. TMNEW/QMNEW 已定义但未调用——同上，主例程内联了等价逻辑。
5. GCMGRID_P2 中 call_mass_fixer 分支用 sum(qq(il,:)) 归一化（第 1720 行），与传统分支的 rhoknew/rhonew 不同——其物理含义是否为补偿非扩散 tracer 的质量份额。

相关页面

• aeronomars/index：aeronomars 目录总览，本文件在其"分子扩散"分类下。
• moldiff：原版（legacy）14 物种分子扩散模块（死代码），含 α 矩阵 + LU 求逆方法对比。
• moldiffcoeff_red：简化扩散系数例程，由本文件 firstcall 调用。
• moldiffcoeff：原版扩散系数例程，供 moldiff.F 使用。
• tracer_mod：提供 noms/mmol。
• geometry_mod：提供 cell_area（网格面积）。
• planetwide_mod：提供 planetwide_sumval（全行星求和）。
• conduction：热层分子热传导，与分子扩散共同影响热层温度/成分剖面。
• molvis：热层分子黏性速度扩散，与分子扩散并列为热层三大输运过程。
• conf_phys：读取 call_mass_fixer_moldiff_MPF 标志。
• moldiff_red.F90：legacy scheme 变体。
• thermosphere_mod：热层物理过程总调度，本文件调用方。