thermcell_main_mars.F90

路径

LMDZ.MARS\libf\phymars\thermcell_main_mars.F90

所属目录 / 模块

libf\phymars

文件定位

thermcell_main_mars.F90 定义 Martian thermal plume model 的主体子步长求解器。它由 calltherm_interface_mod 在 nsplit_thermals 子步长循环内调用，负责从当前温度、风、压力、层几何和 CO2 tracer 状态出发，估计热羽流源层、保守羽流高度、卷入/退卷、上升质量通量、下沉气流、温度 tendency、CO2 tracer tendency、覆盖率和热通量诊断。

本文件不处理非 CO2 tracer 或风的被动输送；那些变量在接口层用 thermcell_dqup 后处理。CO2 tracer 是例外：它会改变平均分子量和修正位温，因此在本主求解器内部用同一个输送算子更新。

定义的符号

依赖的模块

调用的关键例程

被调用位置

输入

输出

共享状态与副作用

• 本例程读取 comtherm_h 中的编译期参数和少量模块变量；不写这些共享参数。
• 非 MESOSCALE 构建会调用 planetwide_maxval，并在并行物理网格上触发 gather/bcast 式全行星规约。
• 例程不读写文件、不分配模块级内存、不直接写 NetCDF。
• 当 limz >= nlayer 时会打印警告并把 limz 限制为 nlayer。随后若存在 CO2 tracer，第 1193-1195 行仍会用该 limz 调用 thermcell_dqup；而 thermcell_dqup 的环境更新会读取 q(:,limz+1)，所以触顶截断路径需要按顶层边界/数组维度风险处理。当质量通量修正次数过多且 thermverbose 为真时会打印修正统计。
• entr 会在 CO2 tracer 输送前随 detrmod<0 修正块被再次调整，因此下游使用的是清理后的质量通量。

核心逻辑

1. 清零 entr/detr/fm，用 zhc=pt/zpopsk 得到位温，并令内部 CO2 输送子步数 ndt=1。
2. 若 igcm_co2.ne.0，用 CO2 与非凝结气体分子量计算 A/B，把位温改写为依赖 CO2 混合比的修正位温；否则直接用普通位温。
3. 识别源层、保守羽流速度和羽流高度，估计无量纲夹卷/退卷率，并迭代更新上升流温度、垂直速度和源层夹卷。
4. 根据 zw2 加权积分得到 zmax；非 MESOSCALE 下用 planetwide_maxval(lmax,limz) 后令 limz=limz+2，MESOSCALE 下直接设为 nlayer-5。
5. 用源层闭合得到质量通量尺度 f，再把无量纲 entr_star/detr_star/alim_star 转为真实 entr/detr。
6. 重建 fm，随后逐层执行负质量通量、过大退卷、过强夹卷和 alphamax 覆盖率限制等修正，保证 fm(lmax+1)=0。
7. 按上升流温度差计算 zdthladj，再按 fdfu 和高度比例构造下沉质量通量 fm_down 与下沉温度廓线 ztvd，计算 zdthladj_down。
8. 用清理后的 fm/(rhobarz*zw2) 得到界面覆盖率 fraca。
9. 若存在 CO2 tracer，先由质量守恒重算 detrmod=fm(k)-fm(k+1)+entr(k)，再调用 thermcell_dqup 输送 pq(:,:,igcm_co2)，并用输送后的 CO2 重新计算 ratiom。如果第 786-790 行已因热羽流触顶把 limz 截断为 nlayer，这里传入 thermcell_dqup 的仍是 nlayer；下游会读取 q(:,limz+1)，复现时必须检查 pq(:,:,igcm_co2) 的顶层边界是否安全。
10. 用 (zdthladj+zdthladj_down)*zpopsk*ratiom 得到 pdtadj，最后在界面插值环境、上升流和下沉流位温，输出 heatFlux、heatFlux_down、buoyancyOut 和 buoyancyEst。

伪代码

thermcell_main_mars(...):
    entr, detr, fm = 0
    theta = pt / zpopsk

    if igcm_co2 != 0:
        theta_model = theta * (A * pq[:,:,igcm_co2] + B)
    else:
        theta_model = theta

    find source layers and conservative plume velocity
    estimate entrainment/detrainment and updraft profile
    compute lmax and zmax

    if MESOSCALE:
        limz = nlayer - 5
    else:
        limz = planetwide_max(lmax) + 2
    limz = min(limz, nlayer)

    f = closure_mass_flux(...)
    entr = f * (entr_star + alim_star)
    detr = f * detr_star
    reconstruct fm and apply correction blocks

    zdthladj = upward_temperature_transport(fm, theta_model)
    fm_down, ztvd = downdraft_parameterization(fm, zmax)
    zdthladj_down = downdraft_temperature_transport(fm_down, ztvd)
    fraca = fm / (rhobarz * zw2)

    if igcm_co2 != 0:
        detrmod = fm - fm_next + entr
        thermcell_dqup(..., pq[:,:,igcm_co2], pdqadj[:,:,igcm_co2], ndt=1, limz)
        ratiom = 1 / (A * transported_CO2 + B)

    pdtadj = (zdthladj + zdthladj_down) * zpopsk * ratiom
    heatFlux, heatFlux_down, buoyancyOut, buoyancyEst = diagnostics(...)

参与的主题流程

写法特点

• 文件头明确说明本例程在 calltherm_interface 的子步长循环中运行，质量通量、温度和 CO2 MMR 都在这里处理。
• MESOSCALE 与 GCM 路径的 limz 不同：mesoscale 固定用 nlayer-5，GCM 用全行星最大 lmax+2。
• CO2 只在本例程内部作为主动耦合 tracer 输送；非 CO2 tracer 和风在接口层按平均质量通量被动输送。
• 质量通量重建后有多段数值保护和修正计数，结果复现时不能只看闭合公式，还要保留这些修正分支。

复现要点

• igcm_co2 会改变热力变量定义和温度 tendency 换算；没有 CO2 tracer 的运行不能与有 CO2 tracer 的 polar mixing 结果直接等价比较。
• limz 是全域后续输送范围。并行 GCM 中若把 planetwide_maxval 替换成本地最大值，会造成不同进程覆盖层数不一致。
• comtherm_h 中 ad/bd/ae/be/fdfu/omega/a1/b1/a1inv/b1inv 是主公式调参入口；修改这些常数会同时影响高度、质量通量、下沉气流和热通量。
• thermverbose=.false. 时质量通量修正通常静默执行；调试复现实验若要查看修正频率，需要改编译期参数。
• limz>=nlayer 的警告表示热羽流触及模型顶层，源码会强制截断为 nlayer，这类案例不应视作普通边界层高度结果。
• 触顶截断不是孤立诊断：limz=nlayer 会继续传入第 1193 行 thermcell_dqup，而该算子读取 q(:,limz+1)。移植、数组边界检查或重写 CO2 输送时必须显式处理 q(:,nlayer+1) 访问风险。

待确认

• 源码中 pdqadj 注释写作局部变量，但它是子程序形参之一；本页按接口层实际使用解释为 CO2 输出 tendency。
• 部分质量通量修正计数变量如 ncorecfm4/5/7/8 的触发语义需要结合完整修正分支和调试输出进一步命名化；本页只按源码行为描述。

相关页面

• calltherm_interface_mod：热羽流与 physiq 的接口层，负责子步长调度和非 CO2 被动输送。
• thermcell_dqup：质量通量已知后的被动变量输送算子，本例程用于 CO2 tracer。
• comtherm_h：热羽流编译期参数和调参常数。
• planetwide_mod：非 MESOSCALE 下 limz/llmax 的全行星最大值规约工具。
• comcstfi_h：提供 g 和 r 等物理常数。