extract.F90

快速理解

它做什么： 独立命令行程序，从 NetCDF 中按四维线性插值抽取点值或垂直廓线。

基本过程： stdin/def 读模式（逐点或廓线）→ 指定 lon/lat/alt/Ls/LT → 四维线性插值。

关键结果： ASCII infile_var.dat，服务于逐点抽取和廓线抽取。

路径

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 都说明它服务于两种模式：逐点抽取和固定经纬时刻的垂直廓线抽取。

定义的符号

符号	类型	行号	作用
`extract`	program	1	主程序：读 NetCDF 输入、变量名、模式和查询坐标，输出 ASCII 抽取结果。
`extraction`	subroutine	406	四维插值核心：定位 lon/lat/alt/time 包围索引，检查缺测，再按 time→alt→lat→lon 插值。
`ls2sol`	subroutine	621	将用户输入的太阳经度 Ls 换算为 sol，用于和 LT/longitude 合成 NetCDF `Time` 查询点。

依赖的模块

use 模块	用途
`netcdf`	使用 Fortran 90 `nf90_*` 接口打开文件、查维度/变量、读取属性和四维场。

无 LMDZ 运行时模块依赖。

调用的关键例程

被调用例程	所在位置	作用
`nf90_open`	NetCDF 库	只读打开输入文件。
`nf90_inq_dimid` / `nf90_inquire_dimension`	NetCDF 库	查 latitude、longitude、altitude、Time 维度和长度。
`nf90_inq_varid` / `nf90_inquire_variable`	NetCDF 库	查目标变量 ID、维数和维度顺序。
`nf90_get_var`	NetCDF 库	读取坐标数组和完整四维目标变量。
`nf90_get_att`	NetCDF 库	读取 `altitude:units` 和目标变量 `missing_value`。
`ls2sol`	本文件	把 Ls 转换为年内 sol。
`extraction`	本文件	对一个查询点返回插值值或 `missing_value`。

输入

输入	来源	类型/维度	单位	含义
输入文件名	stdin / `.def`	NetCDF 文件	-	经过 `zrecast` 的 `diagfi`/`concat` 类文件。
`starttimeoffset`	stdin / `.def`	real	sol	输入文件 `Time=0` 相对 Ls=0 的 sol 偏移；读入后加到整个 `time(:)`。
`extract_mode`	stdin / `.def`	1 或 2	-	1=逐点抽取；2=垂直廓线抽取。
`varname`	stdin / `.def`	字符串	-	要抽取的四维变量名。
查询点	stdin / `.def`	模式 1: `lon lat alt Ls LT`；模式 2: 先 `lon lat Ls LT`，再逐行 `alt`	deg, deg, m/Pa, deg, hour	插值目标坐标。
`longitude` / `latitude` / `altitude` / `Time`	输入 NetCDF	1D 坐标	随文件	插值网格。
目标变量	输入 NetCDF	REAL(lon,lat,alt,time)	随变量	被抽取的四维场。
`missing_value`	目标变量属性	real	随变量	缺测输出值；任何包围角点缺测则返回它。

示例 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

输出

输出	去向	类型/维度	单位	含义
`infile_var.dat`	磁盘 ASCII	文本	随变量	`outfile=输入文件去掉最后 .nc + "_" + varname + ".dat"`。
模式 1 行	输出文件	`lon lat alt Ls LT value`	混合	每个输入查询点一个插值结果。
模式 2 行	输出文件	`alt value`	垂直坐标 + 变量单位	固定 lon/lat/Ls/LT 下的垂直廓线。

共享状态与副作用

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

核心逻辑

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

参与的主题流程

主题	参与方式
util 后处理链	常接在 zrecast 后，对压力或高度坐标上的 GCM 输出抽点/抽廓线。
观测或现场点比较	把模型四维场按给定经纬度、Ls、当地时和高度/压力抽成 ASCII 表，便于和观测点位或剖面比较。
一维输入准备辅助	可抽取垂直廓线；更专用的 1D 列抽取见 extractcolumnfor1D.F90。

写法特点

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 逐值等号比较，若输入变量缺测约定经过数值转换，可能无法识别。
对极点、边界经度/时间端点的查询可能触发上/下索引越界或返回缺测，应避免贴边抽取。