simu_MCS.F90

路径

LMDZ.MARS\util\simu_MCS.F90

所属目录/模块

util

文件定位

MRO/MCS 观测模拟器（observer simulator）离线后处理工具。把 GCM 输出（diagfi.nc/concat.nc 或 stats.nc）按 MRO/MCS 观测数据（Luca Montabone 分箱版）的时空采样和分箱方式重采样，生成一个与 MCS 文件同格式、同维度声明顺序的 netcdf 文件，用于 GCM 与 MCS 观测的逐 bin 对比。

核心思路：MCS 数据按 Ls 区间分箱。程序把每个 Ls 区间换算为 GCM 的 sol 区间，在该区间内对每个 MCS 空间 bin（经/纬/高）的若干 sol×当地时（LT）采样点做 GCM 场的四维线性插值，再平均成一个输出 bin 值。昼/夜（dayside/nightside）分两遍处理。还能计算柱积分尘埃光学厚度比 tau_ratio = tau_GCM/tau_MCS，用于归一化对比。

作者 Antoine Bierjon（2019-2020）。本文件是一个独立的 program（非 GCM 运行时模块），编译为命令行工具，通过 stdin 交互式读入文件名和变量列表。

定义的符号

依赖的模块

无 LMDZ 物理模块依赖——本工具自包含，仅依赖 netcdf 库。

调用的关键例程

输入

输出

共享状态与副作用

• stdin 交互：通过 READ(*,*) 读取 obsfile、gcmfile、变量列表、outfile——是交互式命令行工具，非批处理库。
• 大量 stdout 提示：欢迎语、变量清单、进度、跳过提示等。
• stop 硬终止：任何 netcdf 错误（经 status_check）、非预期坐标值、无变量可处理时直接 stop。
• extraction 的 SAVE 索引缓存：prev_lon/lat/alt/sol 与 ilon_inf/sup 等 SAVE 变量缓存上次的包围索引；坐标未变时跳过二分查找——性能优化，但使该子程序非纯函数、非线程安全。
• netcdf 文件句柄：obsfid/gcmfid（只读 nf90_nowrite）、outfid（读写）全程保持打开，section 6 关闭。

核心逻辑

1. 打开文件 + 读维度（section 1.1-1.2）：交互读 obsfile/gcmfile；分别读 lon/lat/alt/time 维度。MCS 时间轴是 Ls，GCM 时间轴是 sol。由 altitude 的 units 属性判定坐标类型 'p'（Pa，压力）或 'z'（m，高度），两文件须一致。OBSdeltaLs = OBSLs(2)-OBSLs(1) 为 Ls bin 宽。

2. 创建输出文件（section 1.3）：用 inidim 按 MCS 维度顺序建 outfile 坐标。

3. 变量管理（section 2）：

   • 读 GCM 变量清单，跳过 15 个非处理变量（坐标/混合系数等）；用户选变量或 all。
   • 检测能否算不透明度：dso→dustok1、dsodust→dustok2、dustq→dustok3、h2o_ice→wiceok；若可，加载 rho（stats 文件只有部分 sol，循环填满所有 sol），并把 dust/wice 追加进待处理变量。
   • 昼/夜大循环开始：dayside 用 d* 变量、nightside 用 n* 变量。从 obsfile 读 LT 的均值/max/min（OBSLT/OBSLTmax/OBSLTmin）和各 bin 样本数 numbin（temp/dust/wice）。

4. 变量循环（section 2.4-2.6）：逐个待处理变量，通用读取 GCM 4D 场（stats 文件同样循环填充）；dust/wice 不透明度由 dso/dsodust/dustq(×rho) 和 h2o_ice×rho 算出；读对应 MCS 参考变量与 numbin；在 outfile 定义该变量。

5. 提取与分箱（section 3，四重坐标循环 Ls→lat→lon→alt）：

   • 3.1：把 Ls bin [Ls-δ/2, Ls+δ/2] 经 ls2sol 换算成 sol 区间，定位 GCM 时间索引 m_minsol/m_maxsol（±0.5 sol 余量覆盖全经度）；取区间内的整数 sol 列表 int_sol_list（共 solcount 个）。无 GCM 数据则该 Ls 填缺测。
   • lon 规整到 [-180,180]；取该 bin 的观测 LT 均值 LT_val；缺测则跳过。
   • 3.2：LTcount = floor(numbin) 为该 bin 的样本数；gen_sol_list 生成 solcount*LTcount 个 sol 采样点（每个整数 sol 配 LTcount 个围绕 LT_val 对称分布的 LT，换算到 lon=0° 的 sol：sol = int_sol + (LT - lon/15)/24）。
   • 3.3：对每个采样点调 extraction 插值；非缺测值累加，最后除以有效样本数得 solbinned_value 写入 outvar。
   • 3.4：整变量写入 outfile，释放数组。

6. CDOD ratio（section 4，仅 OBSalttype='p'）：在 outfile 找尘埃不透明度（d/ndust→d/nopa_dust→d/nopadust 依次回退）、obsfile 找温度；用 OBS 温度+压力算 OBS 密度，按层积分得 GCM 与 MCS 柱光学厚度，算比值 tau_ratio 写入 outfile。

7. 昼夜循环结束 + 关闭文件（section 5-6）。

伪代码

program simu_MCS:
  读 obsfile, gcmfile (stdin); 读 lon/lat/alt/time 维度
  由 altitude units 定 alttype ('p'/'z'); OBSdeltaLs = OBSLs(2)-OBSLs(1)
  inidim(outfile)                          // 按 MCS 维度顺序建输出
  读 GCM 变量清单 (跳过 15 个非处理变量); 用户选变量或 all
  检测 dust/wice 不透明度可行性; 若可: 加载 rho (stats 循环填充)

  for dayornight in [dayside, nightside]:   // 昼/夜大循环
    从 obsfile 读 OBSLT(ave/max/min), numbin(temp/dust/wice)
    for var in 待处理变量:                   // 变量循环
      读 GCM 4D 场 (stats 循环填充); 算 dust/wice 不透明度 (×rho)
      读对应 MCS 参考变量 + numbin; 在 outfile 定义该变量
      for l in Ls:                          // 坐标循环
        [Ls±δ/2] --ls2sol--> sol 区间 -> m_minsol/m_maxsol
        int_sol_list = 区间内整数 sol (solcount 个)
        for j in lat, i in lon:
          lon 规整 [-180,180]; LT_val = OBSLT(i,j,l)
          for k in alt:
            LTcount = floor(numbin); 若 0 -> 填缺测
            sol_list = gen_sol_list(int_sol_list, LTcount, LT_val, ...)  // solcount*LTcount 点
            for sol in sol_list:
              extraction(lon,lat,alt,sol, GCM_var) -> extr_value
              累加非缺测值
            outvar(i,j,k,l) = 累加 / 有效样本数
        outfile <- outvar
      if alttype=='p' 且有尘埃+温度:          // CDOD ratio
        tau_ratio = 柱积分(tau_GCM) / 柱积分(tau_MCS) -> outfile
  关闭文件

subroutine extraction(lon,lat,alt,sol, ..., field, ..., value):
  若坐标变化: 二分查找包围索引 (SAVE 缓存); 越界 -> value=missing 返回
  16 个角点任一为 missing -> 返回
  时间线性插值 -> 垂直插值 (rho/pressure 在 z 坐标下取 log) -> 纬度 -> 经度

参与的主题流程

写法特点

• 独立 program：本文件首行即 program simu_MCS，无内嵌 module，编译为独立命令行工具。
• 交互式 stdin：通过 READ(*,*)/READ(*,'(a50)') 读文件名和变量列表，适合人工或脚本管道驱动。
• 逐 netcdf 调用 + status_check：每个 netcdf 操作后立即查错并 stop，无优雅恢复。
• extraction 的 SAVE 索引缓存：坐标连续时复用上次二分查找结果，是性能优化但非线程安全。
• stats 文件特殊处理：stats.nc 只覆盖一个 sol 周期，按 modulo 循环填充到整个观测期（l 为 GCMstatstimelen 整数倍时特判）。
• 昼/夜双遍：同一套逻辑跑两次，变量名前缀 d/n 切换。
• 垂直 log 插值：rho/pressure 在高度坐标下按对数插值（更符合指数廓线）。
• 硬编码火星常数：r_atm=191、g=3.72（CDOD 积分）；ls2sol 内 year_day=668.6、peri_day=485.35、e_elips=0.0934 等轨道参数。

复现要点

• 工具交互输入顺序：obsfile → gcmfile → 变量列表（逐行，空行或 all 结束）→ outfile。
• obsfile 与 gcmfile 的 altitude 坐标类型（'p'/'z'）必须一致，否则插值无意义；类型由 altitude 的 units 属性前两字符判定（Pa/m）。
• GCM 文件必须完整覆盖观测期的 sol 区间，否则对应 Ls bin 填缺测。
• 不透明度计算依赖 GCM 文件含 rho 及 dso/dsodust/dustq（尘埃）或 h2o_ice（水冰）之一；缺 rho 则关闭不透明度计算。
• ls2sol 用开普勒方程反解（火星轨道：year_day=668.6 sol、peri_day=485.35、timeperi=1.90258...、e_elips=0.0934），与 LMDZ 主模式的 Ls↔︎sol 换算需一致。
• sol 采样点公式 sol = int_sol + (LT - lon/15)/24：把当地时 LT 在给定经度换算为 lon=0° 参考的 sol（经度每 15° 差 1 小时）。
• CDOD ratio 仅在压力坐标（OBSalttype='p'）且 outfile 有尘埃不透明度、obsfile 有温度时计算；尘埃变量名按 d/ndust→d/nopa_dust→d/nopadust 回退查找（兼容旧版命名）。
• 夜间 LT 用 LTmod 映射到 [-12,12) 以保证跨午夜的均值/距离计算连续。

待确认

• extraction 的越界/缺测分支大量 write 被注释掉（仅返回 missing）——调试输出默认关闭，排障时需手动开启。
• inidim 的 londimid/latdimid 等为 intent(inout)，进入时携带"哪个维度是 altitude"的标记（do d=1,4 中 altdimid.eq.d 判断）——调用侧如何预置这些 ID 需结合 section 1.3 完整确认（本次未逐行读 1.3 定义段）。
• gen_sol_list 中 N=floor(LTcount/2)，当 LTcount=1 时 N=0、循环不执行，仅靠末尾 if mod(LTcount,2)==1 补入 LTave——单样本 bin 退化为仅取均值 LT，符合预期但边界值得注意。

复现风险

• 工具为交互式 stdin 驱动，自动化复现需用管道/重定向喂入固定顺序的输入行。
• stop 式错误处理：输入文件缺变量或维度不匹配会直接终止，无部分结果保留。
• extraction 的 SAVE 缓存使其在并行/多次独立调用场景下行为依赖调用历史——单线程顺序调用安全。
• stats 文件循环填充假设观测期为 sol 周期的整数延拓，长期外推的物理意义有限。

相关页面

• zrecast.md - 同为 util 下的 GCM 输出后处理工具，负责垂直坐标重投影。
• util-toolchain - util 离线后处理和观测比较工具链。