zrecast.F90

快速理解

它做什么： 独立命令行程序，把 GCM 输出从 sigma/hybrid 垂直坐标重投影到物理垂直坐标（压力/高度/半径）。

基本过程： 静力方程积分建立几何高度 → 按用户给定新垂直层插值 → 写新 NetCDF。

关键结果： 物理垂直坐标下的 NetCDF 场，是制图前的必要步骤。

路径

LMDZ.MARS\util\zrecast.F90

所属目录/模块

util

文件定位

zrecast.F90 是 LMD Mars GCM 输出的垂直坐标重投影工具。它读取 diagfi.nc、concat.nc 或 stats.nc 一类 lon-lat-alt-time 四维 NetCDF 场，把原始 GCM sigma 或 hybrid sigma 垂直坐标转换到物理垂直坐标：压力、areoid 以上高度、本地地表以上高度，或距行星中心半径。

util/README 明确说明 GCM 原始 hybrid 坐标不对应任何物理垂直坐标，zrecast 是制作可发表科学图像前的必要步骤。程序通过静力方程积分建立 GCM 层的几何高度/压力关系，再按用户给定的新垂直层进行插值，并把结果写入一个新的 NetCDF 文件。

本文件是独立命令行 program zrecast，不是运行时物理模块；它通过 stdin 交互或 zrecast.*.def 重定向输入。

定义的符号

依赖的模块

主程序主体还包含 include 'netcdf.inc'，大量使用旧式 NetCDF Fortran NF_* 接口读写文件。

调用的关键例程

输入

示例输入文件：

zrecast.auto.def:
../stats.nc
temp

2
yes
-2000 2000
3

zrecast.manual.def:
../stats.nc
temp

2
no
3
-2000
0
1000

输出

共享状态与副作用

• 文件副作用：读取输入 NetCDF 和可能的回退文件 diagfi.nc、diagfi1.nc、phisinit.nc；创建新的 *_P.nc、*_A.nc、*_S.nc 或 *_R.nc。
• stdout 交互提示：程序通过 write(*,*) 输出大量选择菜单、查找提示和错误信息。
• 硬停止：多数 NetCDF 错误、缺少必要维度/变量、非法输入值都会直接 stop。
• 全局常数：planet_const 中的 a0、g0、Rmean 由 controle 初始化，之后被高度积分和 areoid/半径计算使用。
• 内置 areoid 数据：半径坐标分支依赖文件末尾内嵌的 Lemoine/MOLA 系数，而不是外部数据文件。

核心逻辑

1. 读输入文件和维度：打开用户指定的 GCM 输出，读取 longitude、latitude、altitude、Time 维度，收集变量列表、ps、温度、sigma 或 aps/bps、phisinit、可选 rho 和 controle。
2. 初始化行星常数：init_planet_const 从 controle 读行星半径、参考重力和平均气体常数；若当前输入文件没有 controle，继续查找 diagfi.nc 和 diagfi1.nc。
3. 选择目标垂直坐标：用户选择压力、areoid 高度、本地地表高度或半径坐标；层数可由 min/max/nlevels 自动生成，也可逐层手动输入。本地地表高度还有 mcd 特殊构建模式。
4. 构建原始 GCM 层压力：若有 sigma，按 P=sigma*ps；否则按 P=aps+bps*ps。
5. 构建几何高度关系：需要高度坐标时，程序积分静力方程。若有 rho，用 R=P/(rho*T)；若没有 rho，使用 Rmean 的近似分支。
6. 定义输出 NetCDF：复制维度、坐标、phisinit、原始 sigma/hybrid 层和可得的 controle；根据输出模式定义 pressure 或 zareoid 辅助变量。
7. 逐变量垂直插值：对用户选择的变量读取四维场；rho* 和 num_* 前缀变量设置为 log-space 插值；其他变量按线性或压力对数坐标规则插值。
8. 写结果并关闭文件：每个变量写入输出文件，最后写辅助坐标变量、关闭输入/输出 NetCDF。

伪代码

program zrecast:
  read infile and variable list from stdin
  open NetCDF input; read lon/lat/alt/time and ps
  read temp from preferred or fallback names
  read sigma, else aps/bps
  read phisinit from input, else diagfi.nc, diagfi1.nc, phisinit.nc
  optionally read rho and controle
  init_planet_const(controle)

  ask output vertical coordinate:
    1 pressure -> build target pressures
    2 above areoid -> build target zareoid levels
    3 above local surface -> build target zsurface levels
    4 radius -> build target radius levels

  compute original pressure on GCM levels:
    if sigma: press = sigma * ps
    else:     press = aps + bps * ps

  if target needs height:
    if rho exists: use hydrostatic integration with R=P/(rho*T)
    else:          use crude integration with Rmean
    if radius target: add areoid radius from geoid()

  create output file and metadata
  for each requested variable:
    read 4D data
    flag = log interpolation if variable starts with "rho" or "num_"
    if pressure target:
      p_coord_interp(...)
    else if areoid or radius target:
      z_coord_interp(...)
    else if local-surface target:
      zs_coord_interp(...)
    write output variable
  write pressure or zareoid auxiliary field; close files

参与的主题流程

写法特点

• 混合 NetCDF API：主程序使用 netcdf.inc 里的 NF_* 旧接口，init_planet_const 使用 use netcdf 的 nf90 接口。
• 交互式 stdin 驱动：与 util 目录其他工具一致，既可回答屏幕问题，也可用 .def 文件重定向。
• 变量名前缀规则：rho* 和 num_* 变量被强制按对数空间插值，用来处理密度或数密度廓线。
• 多层回退查找：phisinit 和 controle 不只从当前输入文件读取，还会在常见伴随文件中查找。
• 常数和经验层构建：build_zs 使用 P_ref=610、低层/高层尺度高度和 tanh 过渡，为 mcd 模式生成本地高度层。
• 半径坐标自带 areoid 模型：build_gcm_ra 通过 areoid_ini/geoid 使用文件内嵌球谐系数计算 areoid 半径。

复现要点

• 输入文件必须包含 ps、温度变量、sigma 或 aps/bps，并能提供 phisinit；否则程序会在回退文件中查找，失败则停止。
• 温度变量名最好显式给出；源码有候选名回退，但不同历史输出的命名差异会影响自动识别。
• 若缺少 rho，高度积分会退化为常数 Rmean 近似；同一输出坐标下，含 rho 和不含 rho 的结果并不完全等价。
• 对密度、数密度一类变量，应确保变量名以 rho 或 num_ 开头，否则不会触发 log-space 垂直插值。
• 压力坐标输出会额外写 zareoid；高度/本地高度/半径坐标输出会额外写 pressure，后续工具应读取正确的辅助变量。
• zrecast.auto.def 的自动层模式输入顺序是：输入文件、变量列表空行、坐标类型、yes、min/max、层数；手动层模式则是：坐标类型、no、层数、逐层值。

待确认

• 源码注释说明压力坐标输出时缺少温度可被容忍；本页尚未逐行验证所有缺温分支的行为。
• temp 的候选变量名回退已按源码搜索记录，具体优先级和所有提示文本未在本文完整展开。
• 半径坐标的 areoid 系数来自文件内嵌常数，未与外部 MOLA/Lemoine 数据版本做数值核对。

复现风险

• 程序大量使用 stop，批处理链中任一缺失变量或 NetCDF 错误都会终止而不生成可用的部分结果。
• 回退读取 diagfi.nc、diagfi1.nc、phisinit.nc 依赖当前工作目录，移动 .def 或在其他目录运行时容易找不到伴随文件。
• stats.nc、diagfi.nc 和 concat.nc 的时间覆盖、变量维度和命名可能不同；重投影前应先用 ncdump 或同类工具确认。
• 本地地表高度的 zs_coord_interp 在最低层以下有外推逻辑，而最高层以上可能填缺测或做指数外推，边界层外结果需要谨慎解释。

相关页面

• util/index.md - util 后处理工具总览。
• simu_MCS.md - 同属 util 的观测模拟和分箱比较工具。
• extract.md：常接在 zrecast 后抽取点值或廓线。