hrecast.F90

快速理解

它做什么： 水平重网格工具，把 GCM 输出按用户给定新经纬网格做面积权重映射。

基本过程： 读输入和新网格定义 → interp_horiz 计算格盒球面面积交叠 → 面积权重重分配 → 写 *_h.nc。

关键结果： 保持单位面积强度量总量一致的重网格文件（非点值双线性插值）。

路径

LMDZ.MARS\util\hrecast.F90

所属目录/模块

util

文件定位

hrecast.F90 是 LMD Mars GCM 输出的水平重网格工具。它读取 diagfi.nc、stats.nc、concat.nc 等 NetCDF 输出，按用户给定的新经纬度规则网格，把 3D (lon,lat,Time) 和 4D (lon,lat,alt,Time) 变量水平映射到新网格，输出文件名自动生成为 input_h.nc。

本工具不是点值双线性插值。核心例程 interp_horiz / iniinterp_h 计算旧网格格盒与新网格格盒的球面面积交叠，用 intersec/airen 面积权重把旧格点标量场重分配到新格点，源码注释说明目标是对单位面积强度量保持总量一致。

定义的符号

依赖的模块

调用的关键例程

输入

示例 hrecast.def 的顺序是：输入文件、变量名列表、输出经度数量和值、输出纬度数量和值。

输出

共享状态与副作用

• 固定输出命名：outfile=infile(1:len_trim(infile)-3)//"_h.nc"；NF_CLOBBER 会覆盖同名输出。
• stdin 交互：变量列表和输出网格完全由 stdin 提供，hrecast.def 必须严格按源码读取顺序。
• 整变量读入内存：3D 和 4D 变量分别整块读入，长时间序列或大网格可能占用大量内存。
• 缓存插值权重：interp_horiz 用 SAVE 保存 ktotal/iik/jjk/ik/jk/intersec/airen，网格尺寸未变时复用权重。
• 硬停止/exit：NetCDF 错误多为 stop；插值尺寸超过静态上限时调用 exit(1)。

核心逻辑

1. 读输入和变量列表：打开输入文件，列出 3D/4D 变量。用户可指定变量或输入 all；all 分支只收集 4D 变量。
2. 读取输入坐标：读取 latitude、longitude、altitude、Time。输入纬度必须首尾接近 90 和 -90，经度必须首尾接近 -180 和 180。
3. 读取垂直辅助量：优先读取 sigma；若无，则读取 aps 和 bps。可选读取 phisinit。
4. 构造输入格盒边界：输入经度边界用相邻经度中点，最后一个边界为第一个边界加 2π；输入纬度边界用相邻纬度中点。
5. 读输出网格并构造工作网格：
   • 若输出经度首尾已相差约 360，则视为周期网格；否则追加一个 lon(1)+360 工作点。
   • 若输出纬度不含 90/-90，则在工作网格中补极点，并把近极边界放到接近极点的位置。
6. 创建输出文件：定义 longitude、latitude、altitude、Time 维度和坐标变量；复制 altitude 属性、Time units、sigma 或 aps/bps、变量 long_name/title、units、missing_value。
7. 水平映射：
   • 对 phisinit 调 interp_horiz，得到工作网格结果，再按是否补极点裁剪回用户网格。
   • 对 3D 变量逐 time 调 interp_horiz。
   • 对 4D 变量逐 time 调 interp_horiz，一次传入全部 altitude 层作为 lm。
8. 面积交叠权重：iniinterp_h 为每个新格盒找所有相交旧格盒，计算交集面积 (bb-aa)*(sin(dd)-sin(cc))；interp_horiz 累加 old_value*intersec/new_area，并对新网格南北极点做面积加权平均。

伪代码

program hrecast:
  read infile; open NetCDF
  list candidate variables with ndims 3 or 4
  read selected variables
  read input lon/lat/alt/Time, sigma or aps/bps, optional phisinit
  require input lon starts -180 and ends 180; lat starts 90 and ends -90
  build input lon/lat cell boundaries
  read output lon/lat coordinate lists
  build periodic/pole-complete work grid and boundaries
  create infile_h.nc and define coords/variables
  write coordinates and vertical metadata
  if phisinit exists:
    interp_horiz(phisinit) -> output phisinit
  for each selected 3D variable:
    for each time: interp_horiz(surface field)
    write output variable
  for each selected 4D variable:
    for each time: interp_horiz(all altitude layers)
    write output variable

subroutine interp_horiz:
  if first call or dimensions changed: iniinterp_h(...)
  zero output
  for each old/new grid-box intersection:
    new_cell += old_cell * intersection_area / new_cell_area
  area-average north and south pole rows

参与的主题流程

写法特点

• 面积交叠加权：intersec/airen 权重来自格盒边界交集，不是点值插值。
• 输入网格假设强：必须是包含 -180/180 和 90/-90 的 LMDZ GCM 标量网格。
• 输出网格可不含端点：程序会为内部计算补周期经度或极点，最后裁剪回用户指定网格。
• 静态数组上限：interp_horiz/iniinterp_h 对新旧网格维度和交集数有固定上限。
• 属性复制有限：只显式复制或转写 long_name/title、units、missing_value；其他属性会提示未转移。

复现要点

• hrecast.def 中变量列表后必须用空行结束；若输入 all，源码只会处理 4D 变量，不会自动包含 3D surface+time 变量。
• 输出经度必须递增；输出纬度必须从北到南递减。
• 输入文件需包含 sigma 或同时包含 aps/bps；否则创建输出垂直坐标会停止。
• 若目标变量没有 missing_value，输出会写默认 -9.99e+33，但插值本身不显式跳过缺测角点。
• 大文件一次读完整变量，内存需求随 lonlatalt*time 增长。

待确认

• all 分支只收集 4D 变量，是否应包含 3D 变量需与使用者意图确认。
• miss_val 只写属性，不参与 interp_horiz 缺测跳过；含缺测输入时面积平均可能把缺测值当普通数值传播。
• 代码中多个错误信息把 longitude 写成 latitude，属诊断文本问题，不影响计算。
• NF_DEF_DIM Time 失败时错误信息写成 latitude dimension，属诊断文本问题。

复现风险

• 输出 _h.nc 会被 NF_CLOBBER 覆盖。
• 超过 imnmx2=360/jmnmx2=190、imomx=361/jmomx=180 或 kmax=360*179+200000 会 exit(1)。
• 输入经纬度若不是源码假设的 GCM 端点网格会直接停止。
• 对非单位面积强度量使用面积交叠平均时，物理意义需要按变量性质确认。

相关页面

• util/index.md - util 后处理工具总览。
• concatnc.md - 常见上游时间拼接工具。
• zrecast.md - 常见上游垂直坐标转换工具。
• extract.md - 下游点值/廓线抽取工具。