check_co2cycle.F90

路径

LMDZ.MARS\util\check_co2cycle.F90

所属目录/模块

util

文件定位

check_co2cycle.F90 是一个离线 CO2 cycle 诊断工具。它读取一个或多个 diagfi、concat 或 Xhistins 类 NetCDF 输出文件，把地表压力 ps 双线性插值到 Viking Lander 1、Viking Lander 2 和 InSight 三个硬编码站点，按站点真实高度做静力外推，同时计算全球平均地表压力、南北半球地表 CO2 冰等效压力和总 CO2 库存。

程序还会在第一轮输出之后重新读取站点压力 ASCII 文件，生成逐日平均文件，并在完整火星年上做 6 阶 Fourier 谐波重建。

定义的符号

依赖

调用的关键例程

输入

硬编码站点

程序把这些高度乘以 g 转成地形势，再与四邻点 phisinit 加权得到的模式地形势比较。

输出

副作用

• 输出文件用固定文件名直接打开，已有同名 ASCII 文件会被覆盖。
• 程序不创建输出目录；pathsor 必须已经存在。
• 会二次读取自己刚写出的 ps_*_yearN 文件来计算日平均，再读取 _diurnal 文件做谐波重建。
• day0 在读取 controle(3) 后被强制重置为 0，跨文件时间连续性由每个文件末尾的 day0 = day0 + day 维护。

核心逻辑

1. 读取交互输入：输入目录、输出目录、首个文件名和输出时间轴类型；后续文件名在每个文件处理完后继续从 stdin 读取。
2. 读取首个 NetCDF 文件建立网格：
   • 若存在 time_counter 维度，认为是 DYNAMICO/XIOS 插值输出；
   • 否则认为是 native lon-lat 输出，并把最后一个冗余经度排除在 iend=lonlength-1 外。
3. 读取静态场和常数：
   • phisinit 提供模式地形；
   • controle(3/5/7/8) 分别用于初始日、半径、重力和平均摩尔质量；
   • constR=1000*8.314511/controle(8)；
   • area 或 aire 用于面积权重。
4. 定位三站点四邻点：
   • 在经纬度数组中寻找包围硬编码站点坐标的格盒；
   • 用相对位置 zalpha/zbeta 构造四个双线性权重；
   • 用四邻点 phisinit 得到模式站点地形势 phisim。
5. 逐输入文件和逐时间步读取变量：
   • 读取 Time 或 time_instant；
   • 每个时间步读取 ps、co2ice，并尝试读取 temp；temp 缺失时尝试 temp7，再缺失则 t7=0。
6. 计算全球 CO2 监控量：
   • pstot=sum(area*ps)；
   • captotN=sum(area*co2ice) for j=1..jend/2+1；
   • captotS=sum(area*co2ice) for j=jend/2+1..jend；
   • 写出时乘 airtot1=1/sum(area)，CO2 冰再乘 g 转 Pa。
7. 计算站点压力：
   • 四邻点上用 sum(weight*log(ps)) 插值，随后 zp1=exp(sum)；
   • 第 7 层温度加权为 zt，gh=constR*zt；
   • 真实高度压力 zp2=zp1*exp(-(phisite-phisim)/gh)；
   • 若没有温度导致 gh=0，使用 10 km 标尺高度退化公式。
8. 跨年输出管理：当 sol >= 669 时递增年号；下一时间步关闭上一年文件并打开 yearN 新文件。
9. 日平均：对每个站点、每年，把同一 floor(sol) 的连续记录平均，输出 _diurnal。
10. 谐波重建：若日平均最后一个 Ls >= 360-1.e-3，计算 0 到 6 阶 Fourier 系数，并用 1 到 669 sol 的 Ls 重建压力序列。

伪代码

read input directory, output directory, first file stem, time_unit
open first NetCDF
read lon, lat, alt length, phisinit, controle, area
compute constants and station interpolation weights
open year1 prestot and ps_SITE files
while file stem is not blank:
  open file.nc
  read time vector
  for each time index:
    if year changed: rotate output files
    read ps, co2ice, temp or temp7
    sol = time + day0, fold sol into 0..669
    write global pressure and CO2 ice inventory
    for each site:
      bilinear interpolate log(ps) and temp7
      extrapolate pressure from model terrain to site terrain
      write ps_SITE_yearN
    if sol reaches year end: increment year
  day0 += last time
  read next file stem
for each completed year and site:
  read ps_SITE_yearN and average by integer sol
  write ps_SITE_yearN_diurnal
  if full year is present:
    compute Fourier coefficients through n=6
    reconstruct one value per sol
    write ps_SITE_yearN_harmonics

参与的主题流程

复现要点

• 当前源码读取顺序是：输入目录、输出目录、首个文件名、time_unit，然后在每个文件处理完后继续读取下一个文件名，空行结束。
• 仓库中的 check_co2cycle.def 文件把 time_unit 放在文件列表末尾；按当前源码字面读取会把第二个文件名读入 time_unit，自动化运行前应重排或交互输入。
• 输入文件名不带 .nc；程序自行拼接 pathchmp/nomfich.nc。
• 输入 NetCDF 必须至少有 ps、co2ice、phisinit、controle、area 或 aire。
• 经度、纬度排序被源码写死为经度递增、纬度递减；不符合此约定会导致站点定位失败或权重错误。
• 若使用 native lon-lat 输出，最后一个冗余经度不参与面积积分和站点定位。
• 完整年判定依赖最后的 Ls >= 360-1.e-3；不完整年只会写日平均，不写有效谐波重建。

待确认

• sol2ls 使用 year_day=669.，ls2sol 使用 year_day=668.6d0；两者常数不完全一致，页面按源码记录。
• 南北半球 CO2 冰积分都包含 jend/2+1 纬圈，赤道或中间纬圈是否有意双计需开发者确认。
• temp 读取成功时直接取 t(:,:,7)；若垂直层少于 7 层，源码没有保护。
• time_unit 非 1/2 时都进入 “both” 分支，没有拒绝非法值。
• check_co2cycle.def 的说明把 time_unit 写在文件列表之后，但主程序先读取首个文件名、再读取 time_unit，后续文件名在文件循环末尾继续读取；页面按源码记录，样例文件是否仍能直接用于当前源码需实机确认。

复现风险

• 输出目录缺失、输入文件缺少任一必需变量、站点坐标不落入网格都会 stop。
• 程序会覆盖同名输出文件。
• nbmax=999999 是静态数组上限；超长数据集可能越界，源码没有动态扩容。
• 终端输出非常详细，大批量处理时日志会很长。

相关页面

• util/index.md - util 后处理工具总览。
• co2-cycle.md - CO2 cycle 主题页。
• simu_MCS.md - 另一个观测比较后处理工具。