jthermcalc_util.F

路径

LMDZ.MARS\libf\aeronomars\jthermcalc_util.F

所属目录/模块

libf/aeronomars

文件定位

热层光吸收计算的共享工具模块：提供柱密度积分（column，含球面大气光学路径几何）、快速查表线性插值（interfast）、几何光程计算（espesor_optico_A）、网格索引查找（grid_R8）以及一个遗留的太阳通量太阳周期修正例程（flujo，无调用方）。被 jthermcalc.F 和 jthermcalc_e107_mod.F 共同依赖。

定义的符号

依赖的模块

column

espesor_optico_A

flujo

调用的关键例程

column 子例程

输入

输出

核心逻辑

1. 重力计算（第 134–137 行）：grav(i) = gg * masa / (radio + iz(i))²，随高度递减。
2. 标高计算（第 140–158 行）：在顶层用 kT/(mg) 计算各物种标高 H*。O₃ 需 chemthermod ≥ 1；N₂/N/NO/NO₂ 需 chemthermod ≥ 2。
3. 数密度提取（第 160–194 行）：从 rm 矩阵按热层化学槽位索引（i_co2=1, i_co=2, i_o=3, ...）提取各物种密度，同样受 chemthermod 门控。
4. 柱密度积分（第 196–329 行）：从顶向下逐层调用 espesor_optico_A 获取光学路径长度：
   • 若返回 ilayesp(nlayesp) == -1（无日照），所有柱密度置 1e25（光学厚极限）。
   • 顶层（jj == nlayer）使用标高 × 数密度 × 路径长度积分；SZA ≤ 60° 用 1/cos(SZA) 平面近似，SZA > 60° 用球面弦长。
   • 其余层用梯形积分：esp(j) × (density(jj) + density(jj+1)) / 2。

interfast 子例程

输入

输出

核心逻辑

1. 可选健全性检查（extra_sanity_checks=.false.，编译时常量关闭）。
2. 遍历每层：若 p(n1) 超出 [limdown, limup] 则权重置零。
3. 否则在 pin 表中从上次位置 nini 起顺序搜索包含 p(n1) 的区间，计算线性权重 wm = |pin(n) - p| / (pin(n+1) - pin(n))。

espesor_optico_A 子例程

核心逻辑

计算从高度 z 出发、太阳天顶角 szadeg 方向的光线穿过的各层几何路径长度（球面大气）。三个区间：

• SZA < 60°：简单 Δz / cos(SZA) 平面平行近似。
• 60° ≤ SZA ≤ 90°：球面弦长 √((R+z_top)² - (R+z_min)²) - √((R+z_bot)² - (R+z_min)²)。
• SZA > 90°：光线先向下到达最低点 z_min，再向上穿过；分 5 个区域处理。若 z_min < iz(1)（地表以下），标记 ilayesp = -1 表示无日照。

grid_R8 函数

在单调递增网格 zgrid(nz) 中查找 z 所在区间索引。越界时钳位并打印警告（不中止），网格非递增或维度 < 2 时 stop。搜索从中间点开始线性扫描。

flujo 子例程

遗留太阳通量 11 年周期修正：用 date（年份，钳位到 1985–2001）和正弦函数 sin(2π/11 × (date - 1985 - π)) 对前 24 个光谱区间的 fluxtop 做修正，末尾加 (1.52/dist_sol)² 日火距离修正。当前源码中未发现任何调用方。

共享状态与副作用

• 读：param_v4_h::radio（火星半径）、gg/masa（重力和质量）、kboltzman/n_avog（物理常数）；tracer_mod::igcm_*（tracer 索引）、mmol（分子量）。
• 写：column 无模块级副作用，只写调用方提供的输出参数；interfast 同理。
• 标准输出：espesor_optico_A 在 SZA > 90° 且 z_min 接近 z 时打印舍入警告（第 549–551 行）；grid_R8 越界时打印警告（第 659–670 行）。

参与的主题流程

写法特点

• 固定格式 Fortran 77 风格（.F 扩展名但使用 module/end module 结构）。
• 热层化学物种槽位索引硬编码（i_co2=1, i_co=2, i_o=3, ..., i_n2=17），注释要求必须与 chemthermos.F90 保持一致。
• 注释掉的旧索引方案（第 116–128 行）显示历史重构痕迹。
• column 的 iz 维度为 nlayer+1（界面高度），与其他热层例程的 iz(nlayer) 不同。
• interfast 的 extra_sanity_checks 为编译时常量（.false.），开启后做单调性和索引范围检查。
• flujo 中 date 钳位范围 1985–2001 对应一个太阳周期，参考文献 Gonzalez-Galindo et al. 2005。

复现要点

• column 依赖 tracer_mod 中 igcm_* 索引正确映射到化学 tracer 列表，若 traceur.def 改变需同步更新。
• espesor_optico_A 的 iz 需要 nlayer+1 个界面高度值（从底到顶单调递增），单位为 km。
• interfast 要求输入 pin 单调递增，输入 p 也需单调递增（extra_sanity_checks 可校验）。
• column 返回 1e25 作为无日照标志（espesor_optico_A 返回 ilayesp=-1），调用方需正确处理该极端值。

待确认

• iz 在 column 中的物理含义与单位（km？几何高度还是位势高度？）。column 的 iz 维度为 nlayer+1，而 jthermcalc/jthermcalc_e107 传入的 iz 维度为 nlayer——需确认调用方如何处理维度差异。
• flujo 是否确实已废弃（源码中未发现任何调用方），是否可从模块中移除。
• espesor_optico_A 中 SZA > 90° 的 5 区域划分是否在所有高度和 SZA 组合下都能正确覆盖（注释提及 SZA 接近 90° 时存在舍入误差风险）。
• grid_R8 越界时钳位到端点而非中止——是否可能导致下游计算的静默误差。

相关页面

• jthermcalc.F（基础版调用方）
• jthermcalc_e107_mod.F（E107 增强版调用方）
• hrtherm.F（间接上游）
• chemthermos.F90（热层化学主例程，物种槽位索引一致性要求）
• param_v4_h.F90（物理常数和火星参数）
• tracer_mod.F90（tracer 索引和分子量）