extract.F90

路径

LMDZ.MARS\util\extract.F90

所属目录/模块

util

文件定位

extract.F90 是 LMD Mars GCM 输出的 ASCII 点值/廓线抽取工具。它读取经过 zrecast 处理的 diagfi 或 concat 类 NetCDF 文件，要求目标变量具有 (longitude, latitude, altitude, Time) 四维结构，然后按用户给定的经度、纬度、垂直坐标、太阳经度 Ls 和当地真太阳时 LT 做四维线性插值，输出 infile_var.dat。

本文件是独立命令行 program extract，通过 stdin 或 extract.points.def / extract.profile.def 重定向输入运行。源码注释和 util/README 都说明它服务于两种模式：逐点抽取和固定经纬时刻的垂直廓线抽取。

定义的符号

依赖的模块

无 LMDZ 运行时模块依赖。

调用的关键例程

输入

示例 extract.points.def：

diagfi_P.nc
0.0
1
temp
-20 25 100 35 2

示例 extract.profile.def：

diagfi_S.nc
0.0
2
temp
100 25 10 0
0
500
1000

输出

共享状态与副作用

• stdin 交互：所有运行参数从标准输入读取，.def 文件必须按源码读取顺序排列。
• ASCII 文件副作用：使用 Fortran unit 42 打开输出文件；同名文件会被覆盖或重写，具体行为取决于编译器默认 open 状态。
• 整场读入内存：目标变量被完整读入 field(lonlen,latlen,altlen,timelen)，大文件会占用较多内存。
• 硬停止：缺维度、缺变量、变量不是四维、维度顺序不对、缺 missing_value 属性、非法坐标都会 stop。
• 插值例程的 SAVE 变量：extraction 声明 prev_lon/lat/alt/sol 和索引为 SAVE，但源码未在例程末尾更新 prev_*；因此当前行为近似每次重新定位索引，而不是有效缓存。

核心逻辑

1. 读取输入文件和模式：打开 NetCDF；读 starttimeoffset、extract_mode、varname。模式只能是 1 或 2。
2. 读取坐标：按固定名称读取 latitude、longitude、Time、altitude；随后执行 time(:)=time(:)+starttimeoffset。
3. 判断垂直坐标类型：读取 altitude 的 units 属性；Pa 视为压力坐标 alttype="p"，m 视为几何高度/本地高度坐标 alttype="z"，其他单位直接停止。
4. 读取目标变量：目标变量必须存在，必须是四维，且维度顺序严格等于 longitude、latitude、altitude、Time；读取完整四维场和 missing_value 属性。
5. 创建输出文件：输出名由输入文件基本名和变量名拼接为 *_var.dat。
6. 读取查询坐标：
   • 模式 1：每行读 lon lat alt Ls LT。
   • 模式 2：先读固定 lon lat Ls LT，循环读多个 alt。
   • 经度允许 [-360,360]，随后规整到 [-180,180]；纬度必须 [-90,90]；Ls 必须 [0,360]；LT 必须 [0,24]。
7. Ls/LT 到 sol：ls2sol(Ls,sol) 得到年内 sol；再用 sol=floor(sol)+(LT-lon/15)/24 把当地时和经度折算到模型时间轴；超过 [0,669] 时加/减 669。
8. 四维插值：
   • 定位经纬度、垂直坐标、时间的下/上包围索引。压力坐标按大到小寻找，几何高度按小到大寻找。
   • 任一需要的 16 个角点等于 missing_value 时返回缺测。
   • 先沿时间线性插值，再沿垂直线性插值；当变量名为 rho 或 pressure 且高度坐标为 z 时，垂直方向对变量值取 log 插值。
   • 最后沿纬度和经度线性插值，返回 value。
9. 写输出：模式 1 写 6 列，模式 2 写 2 列；空行或非法输入文本结束循环。

伪代码

program extract:
  read infile, starttimeoffset, extract_mode, varname
  open infile
  read latitude, longitude, Time, altitude
  time = time + starttimeoffset
  alttype = "p" if altitude units == Pa else "z" if units == m
  require target variable dims == (longitude, latitude, altitude, Time)
  read field and missing_value
  open output infile_var.dat

  if mode 2:
    read fixed lon, lat, Ls, LT
  loop over input lines:
    if mode 1: read lon, lat, alt, Ls, LT
    if mode 2: read alt
    validate ranges; normalize lon to [-180,180]
    ls2sol(Ls, sol)
    sol = floor(sol) + (LT - lon/15) / 24
    wrap sol into current Martian year
    extraction(...) -> value
    write output row

subroutine extraction:
  value = missing_value
  find surrounding lon/lat/alt/time indexes
  if outside alt/time range: return missing_value
  if any of 16 corner values == missing_value: return missing_value
  interpolate time -> altitude -> latitude -> longitude

参与的主题流程

写法特点

• nf90 接口：不同于 concatnc 的 NF_* 旧接口，本文件使用 use netcdf 和 nf90_*。
• 严格变量契约：只接受 (lon,lat,alt,time) 四维变量；二维/三维变量不会被处理。
• 时间输入仍是 sol：用户输入 Ls/LT，但文件 Time 必须是 sol；源码注释明确这一要求。
• 垂直单位决定分支：仅识别 Pa 和 m，压力坐标和高度坐标的排序假设不同。
• 缺测传播保守：16 个角点任一缺测则整个查询点缺测，不做部分角点插值。

复现要点

• 输入文件应先用 zrecast 转成压力或高度坐标，否则 altitude:units 和垂直插值语义可能不对。
• starttimeoffset 要让文件 Time=0 对齐到 Ls=0 后的 sol；示例文件用 0.0。
• 查询经度最终被规整到 [-180,180]；输入 NetCDF 的 longitude 坐标也应采用兼容范围，否则索引定位可能失败。
• 对 rho 和 pressure 在高度坐标下抽取时，垂直方向按 log(value) 插值；其他变量线性插值。
• 输出没有表头；后处理脚本需要按模式自行解释列。

待确认

• extraction 中 prev_lon/prev_lat/prev_alt/prev_sol 未被更新，注释暗示的缓存优化实际没有生效；这不改变结果，但影响性能预期。
• 若 longitude 坐标使用 0..360 而查询经度被规整到负值，ilon_inf 可能保持非法初值；实际输入文件经度约定需运行前确认。
• missing_value 属性是硬要求；若变量只提供 _FillValue 或其他缺测属性，程序会停止。
• sol 包裹使用 669，而 ls2sol 内火星年是 668.6；边界 Ls/LT 查询可能存在小差异。

复现风险

• 目标变量完整读入内存，大型四维输出可能内存占用明显。
• open(42,file=outfile,form="formatted") 没有显式 status，重复运行的覆盖行为依赖编译器默认。
• missing_value 逐值等号比较，若输入变量缺测约定经过数值转换，可能无法识别。
• 对极点、边界经度/时间端点的查询可能触发上/下索引越界或返回缺测，应避免贴边抽取。

相关页面

• util/index.md - util 后处理工具总览。
• zrecast.md - 本工具通常读取其输出的压力/高度坐标文件。
• concatnc.md - 上游时间拼接工具。
• simu_MCS.md - 更复杂的观测分箱抽样工具。
• extractcolumnfor1D.F90 - 1D 列抽取工具。