nltecool.F

路径

LMDZ.MARS\libf\phymars\nltecool.F

所属目录 / 模块

libf/phymars

文件定位

nltecool.F 提供旧版 nltemodel=0/1 的 CO2 15 微米 NLTE cool-to-space 冷却参数化。它把文件内置的 68 点压强表、两组 escape function 表和三组成分 VMR 表插值到 GCM 层中，然后在每个水平格点和垂直层上求解两个激发能级、两条辐射通道的简化非 LTE 平衡，输出 dtnlte(ngrid,nlayer)，单位在源码注释中标为 K/s。

这个文件与 nlte_tcool_mod 是互斥的 physiq_mod.F 分支：physiq_mod.F 在 nltemodel.eq.0.or.nltemodel.eq.1 时调用 nltecool，在 nltemodel.eq.2 时调用 nlte_tcool 并把其 K/day 输出除以 86400。

定义的符号

符号	类型	行号	作用
`nltecool_mod`	module	1	封装旧版 NLTE 冷却入口、内置逃逸函数/成分表和两个本地插值 helper。
`np`	integer parameter	9	内置压强、escape function 和 VMR 表的点数，固定为 68。
`pnb`	saved real array	10	初始保存为压强指数；首次调用时转换为 Pa，并声明为 OpenMP `THREADPRIVATE`。
`ef1`, `ef2`	real parameter arrays	36, 61	两条 CO2 15 微米通道的预计算 escape function 表。
`co2vmr`, `o3pvmr`, `n2covmr`	real parameter arrays	86, 111, 136	`nltemodel=0` 使用的 CO2、O 和 N2+CO 体积分数表。
`nltecool`	subroutine	164	主入口；从 GCM 压强、温度和 tracer 混合比计算 NLTE 温度倾向。
`interp1`	subroutine	422	把单个压强表剖面线性插值到 GCM 层压强。
`interp3`	subroutine	451	同时把三组成分 VMR 表线性插值到 GCM 层压强。

依赖的模块

use 模块	only 列表	用途	待确认
`tracer_mod`	`igcm_co2`, `igcm_co`, `igcm_o`, `igcm_n2`, `mmol`	`nltemodel=1` 时从 `pq` 中按 tracer 索引取 CO2、CO、O、N2 质量混合比，并用摩尔质量换算为体积分数。	否
`conc_mod`	`mmean`	`nltemodel=1` 时把 tracer 质量混合比转换为体积分数。	否
`callkeys_mod`	`nltemodel`	在表成分方案 `0` 与动态成分方案 `1` 之间切换。	否

调用的关键例程

被调用例程	所在模块 / 文件	调用位置	作用
`interp1`	本文件	314-315	将 `ef2` 和 `ef1` 从内置 `pnb` 压强表插值到当前水平格点的 `pyy` 层压强。
`interp3`	本文件	316-318	`nltemodel=0` 时将内置 CO2、O、N2+CO VMR 表插值到当前层压强。

输入

输入	来源	类型 / 维度	单位	含义
`ngrid`	`physiq_mod.F`	integer	-	水平格点数。
`nlayer`	`physiq_mod.F`	integer	-	GCM 垂直层数。
`nq`	`physiq_mod.F`	integer	-	tracer 数量，用于第三维 `pq`。
`pplay(ngrid,nlayer)`	`physiq_mod.F::zplay`	real array	Pa	层中心压强；`nltecool` 不再做 Pa/atm 转换。
`pt(ngrid,nlayer)`	`physiq_mod.F::pt`	real array	K	层中心温度。
`pq(ngrid,nlayer,nq)`	`physiq_mod.F::pq`	real array	kg/kg_air	tracer 质量混合比；仅 `nltemodel=1` 分支读取 CO2、O、CO、N2。
`mmean(ngrid,nlayer)`	`conc_mod`	real array	与 `mmol` 一致	动态成分体积分数换算使用的平均分子量。
`mmol(*)`	`tracer_mod`	real array	与 `mmean` 一致	tracer 摩尔质量。

输出

输出	去向	类型 / 维度	单位	含义
`dtnlte(ngrid,nlayer)`	`physiq_mod.F::zdtnlte`	real array	K/s	CO2 15 微米 NLTE 温度倾向；压强范围外置零。

共享状态与副作用

pnb 是 save 数组，首次调用时从指数值转换为 Pa：pnb(i)=1.0e-4*exp(pnb(i))。它同时是 OpenMP THREADPRIVATE，每个线程需要依靠自己的 firstcall 副本完成一次转换。
firstcall 也是 save + THREADPRIVATE 逻辑量，用于防止同一线程重复转换 pnb。如果复现时手动重入或复制本模块，必须保留这对状态的一致性。
文件不读外部 NLTEDAT 表，也不调用 nlte_setup。其表值直接内嵌在源码中，来源注释为历史 nlte_escape.dat / nltedata.h。
本模块没有显式日志输出；旧的 write(7,25) 调试输出已注释。

核心逻辑

首次进入 nltecool 时，把 pnb 从压强指数转换为 Pa，之后保持转换后的表压强。
对每个水平格点 j，把 pplay(j,:) 复制到一维 pyy，再把 ef2、ef1 插值为 escf2(:)、escf1(:)。
若 nltemodel=0，调用 interp3 从内置 co2vmr/o3pvmr/n2covmr 表得到各层 VMR；若 nltemodel=1，后续逐层从 tracer pq、mmean 和 mmol 动态换算成 VMR。
逐层检查压强：pyy(i)>1.0 Pa 或 pyy(i)<4.0e-6 Pa 时把 hr(i) 和 dtnlte(j,i) 置零。
对有效压强层，用 nt = p/(1.381e-17*T) 计算数密度，再把 CO2、O、N2+CO VMR 转成 cm-3 数密度。
用固定同位素丰度 imr1=0.987、imr2=0.00408+0.0112 得到两个基态族群 n1/n2，并根据温度计算 k19、k20、k21 碰撞速率族。
把碰撞损失/产生项与 Einstein 自发辐射项 ae1/ae2 合并，求解两个能级耦合的代数系统，得到 el1/el2 和 pl1/pl2。
计算两条通道冷却 hr1/hr2，求和后通过 0.1*hr*T/(4.4*p) 转成 dtnlte。

伪代码

if firstcall for this thread:
  for each table point:
    pnb = 1e-4 * exp(pnb)
  firstcall = false

for each horizontal grid point j:
  pyy = pplay(j,:)
  escf2 = interp1(ef2, pnb -> pyy)
  escf1 = interp1(ef1, pnb -> pyy)

  if nltemodel == 0:
    co2, o3p, n2co = interp3(tabulated VMRs, pnb -> pyy)

  for each layer i:
    if pyy(i) > 1 Pa or pyy(i) < 4e-6 Pa:
      dtnlte(j,i) = 0
    else:
      nt = pyy(i) / (k_B_cgs * pt(j,i))
      if nltemodel == 1:
        derive CO2, O, CO+N2 VMRs from pq, mmean and mmol
      convert VMRs to number densities
      build two-level collisional production/loss terms
      add radiative escape terms
      solve the two coupled excitation equations
      hr = hr1 + hr2
      dtnlte(j,i) = 0.1 * hr * T / (4.4 * p)

参与的主题流程

主题	参与方式
CO2 15 微米 NLTE 冷却	作为 `nltemodel=0/1` 的旧版 cool-to-space 参数化，直接输出 K/s 温度倾向。
辐射倾向合成	`physiq_mod.F` 把 `zdtnlte` 传给 `blendrad`，后者在 LTE/NLTE 过渡层附近把长波 LTE、NIR CO2 与 NLTE 倾向合成。
高层大气成分敏感性	`nltemodel=0` 使用内置标准成分表，`nltemodel=1` 使用当前 tracer 组成，因此 O、CO、N2 和 CO2 剖面对冷却率有直接影响。

写法特点

固定格式 Fortran 文件内含大段初始化数组；这些表是计算逻辑的一部分，不是运行时外部数据。
pnb 的含义会在首次调用后改变：源码顶部注释称其为 pressure exponent，首次调用后它变成 Pa。文档或调试时不能把转换后的 pnb 仍解释为指数。
interp1 / interp3 对压强范围使用 3.5 Pa 与 4.0e-6 Pa，而主计算层过滤使用 1.0 Pa 与 4.0e-6 Pa；有效计算区间由主循环的 1.0 Pa 上限决定。
nltemodel=1 的成分换算使用 pq*mmean/mmol，这是从质量混合比到体积分数的源代码字面行为。

复现要点

调用 nltecool 时 pplay 必须是 Pa；不要像 nlte_tcool 分支那样先转为 atm。
dtnlte 在本模块内已经是 K/s；不要再像 nltemodel=2 输出那样除以 86400。
callnlte 打开且 nltemodel 为 0 或 1 时才会走本模块；nltemodel=2 使用 nlte_tcool_mod。
若要复现 nltemodel=1，igcm_co2/igcm_o/igcm_co/igcm_n2 对应 tracer 必须存在，且 mmean/mmol 必须与 pq 的质量混合比一致。
复现风险：主循环允许 pyy(i) >= 4.0e-6 Pa，但转换后 pnb 表最后一点约为 1e-4*exp(-3) Pa，略高于 4.0e-6 Pa。若输入压强落在这两个下限之间，插值循环可能找不到包围区间。

待确认

内置 ef1/ef2/co2vmr/o3pvmr/n2covmr 表的外部物理来源仅由源码注释追溯到 nlte_escape.dat / Miguel Lopez-Valverde 方案；本页未核验原始数据文件。
0.1*hr*T/(4.4*p) 中 0.1 和 4.4 的单位换算来源未在本文件内解释；本页按源码公式记录。
interp1 / interp3 在有效压强但未命中任何表区间时没有显式初始化 nm/np/wm/wp，该边界行为需要运行或外部方案文档确认。