chemthermos.F90

路径

LMDZ.MARS\libf\aeronomars\chemthermos.F90

所属目录 / 模块

libf/aeronomars

文件定位

chemthermos.F90 定义 chemthermos_mod，包装热层化学主例程 chemthermos。该例程接收单个物理网格列的化学示踪物 mole fraction、热层温度、数密度和太阳几何，先在 firstcall 阶段校验热层化学需要的 tracer 组合并缓存 GCM tracer 索引，再把 zycol 转成 rm(nlayer,nesptherm) 数密度矩阵，调用 jthermcalc_e107 计算光吸收/光解相关系数，调用 paramfoto_compact 执行 C/O/H/N/ion 化学源汇，最后把 lswitch:nlayer 的结果写回 zycol。

调用方是 calchim_mod 的 .not.unichim 分支。calchim_mod 源码默认 unichim=.true.，因此常规统一化学路径不会进入本例程；只有切到分离的低层/高层化学方案时才会先读取 chemthermos_reactionrates.def，再逐列调用 chemthermos。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• nesptherm、firstcall 和所有 g_* tracer 索引为 save 且 THREADPRIVATE，每个 OpenMP 线程首次进入时独立初始化。
• firstcall 会打印 tracer/方案信息，并在 tracer 缺失或物种数不匹配时 write(*,*) 后直接 stop；这里没有走 abort_physic。
• chemthermod 没有 intent，且 firstcall 中当调用方传入 <1 但存在 igcm_o3 时会把 chemthermod=1，因此该修改可能对调用方可见。
• 每次调用都会 allocate(rm(nlayer,nesptherm))，结束时 deallocate(rm)。
• paramfoto_compact 会写 param_v4_h 中的共享 P*、L*、tau* 等数组；本例程随后读取 Pno、Po2。
• zpress 在本文件中未被读取。

核心逻辑

1. 首次进入时读取 tracer_mod::igcm_*，校验基础 C/O/H 物种必备：CO2, CO, O, O1D, O2, H, H2, OH, HO2, H2O2, H2O。
2. 按 chemthermod 扩展物种集合：>=1 需要 O3；>=2 需要 N2,N,NO,NO2,N2D；==3 需要 CO2+,O+,O2+,CO+,C+,N+,NO+,N2+,H+,HCO2+,e-。
3. 校验 nesptherm 是否等于 11、12、17 或 28；否则打印错误并 stop。
4. 分配 rm(nlayer,nesptherm)，把 zycol(l,g_*)*zdens(l) 转成数密度，所有值用 max(...,1.e-30) 加下限。
5. 调 jthermcalc_e107 计算光吸收/光解相关系数。
6. 调 paramfoto_compact 执行热层化学积分，输入输出都在 rm。
7. 仅对 l=lswitch,nlayer 把 rm(l,i_*)/zdens(l) 写回对应 zycol(l,g_*)，并设置 em_no(l)=Pno(l,45)、em_o2(l)=Po2(l,10)*0.75。
8. 释放 rm。

伪代码

chemthermos(args):
  if firstcall:
    load igcm_* tracer indexes into g_*
    require base C/O/H species
    if chemthermod >= 1 require O3
    else if O3 exists, include it and set chemthermod = 1
    if chemthermod >= 2 require N, N2, NO, NO2, N2D
    if chemthermod == 3 require ion/electron species
    check nesptherm is 11/12/17/28
    firstcall = false

  allocate rm(nlayer, nesptherm)
  for each layer:
    rm(layer, species_slot) = max(zycol(layer, tracer_index) * zdens(layer), 1.e-30)

  call jthermcalc_e107(...)
  call paramfoto_compact(...)

  for layer = lswitch..nlayer:
    zycol(layer, tracer_index) = max(rm(layer, species_slot) / zdens(layer), 1.e-30)
    em_no(layer) = Pno(layer,45)
    em_o2(layer) = Po2(layer,10) * 0.75

  deallocate rm

参与的主题流程

写法特点

• chemthermod 同时控制物种清单、rm 宽度、jthermcalc_e107 和 paramfoto_compact 的下游分支；复现实验需要记录它的实际值。
• 内部物种槽位硬编码，源码注释要求与 jthermcal.F、column.F、paramfoto_compact.F、euvheat.F、hrtherm.F 同步；本仓库实际还可见 jthermcalc_e107_mod.F 和 jthermcalc_util.F 使用相同要求。
• rm(l,i_co) 在 chemthermod>=2 的 N 化学块中被再次赋值一次，与前面的基础赋值相同，属于重复写入。
• em_no 和 em_o2 只在 lswitch:nlayer 被赋值；低层值需要由调用方初始化或保持原值。

复现要点

• 默认 calchim_mod 中 unichim=.true.，不会进入 chemthermos；要复现本路径需使用 .not.unichim 的配置/源码路径。
• 运行目录需要有 chemthermos_reactionrates.def，因为 calchim_mod 会在 .not.unichim firstcall 中调用 chemthermos_readini。
• traceur.def 必须包含与 chemthermod 对应的 tracer 集合；缺任一必需 tracer 会直接 stop。
• chemthermod 若在首次调用中被 O3 自动提升为 1，后续调用沿用同线程缓存的 nesptherm 和 g_*；不要在一次运行中假设可安全切换不同物种集合。
• zdens 不能为 0，因为写回阶段除以 zdens(l)。
• 复现数值时要保留 1.e-30 下限，否则低浓度层会与源码偏离。

待确认

• zpress 作为参数保留但未使用，可能是历史接口遗留。
• chemthermod 由 chemthermos firstcall 修改是否为预期接口行为；源码未声明 intent。
• em_no/em_o2 低于 lswitch 的值是否始终由 calchim_mod 或其他调用方预先清零。

相关页面

• calchim_mod：唯一已确认调用方，说明 .not.unichim、chemthermos_readini 和 chemthermos 调用点。
• chemthermos_readini：读取 chemthermos_reactionrates.def 到 param_v4_h::rcoef。
• reaction-rates：61 个热层化学反应速率常数。
• euvheat：共享热层物种槽位的 EUV 加热入口。
• hrtherm：共享热层物种槽位的 EUV 加热率计算。
• paramfoto_compact：展开本例程调用的热层光化学/离子化学积分细节。