streamfunction.F90

路径

LMDZ.MARS\util\streamfunction.F90

所属目录/模块

util

文件定位

streamfunction.F90 是一个离线 NetCDF 诊断程序，用于从 diagfi.nc、concat.nc 或 stats.nc 类文件计算纬向平均的经向质量流函数、总角动量、纬向平均风、温度、密度、地表压力和初始地形势。

它不是逐时刻输出工具。源码先对输入的 u/v/temp/rho/ps 做全时段平均，再做纬向平均或垂直积分，输出文件固定为单经度、单时间维度的 *_stream.nc。

定义的符号

依赖

输入

输出

副作用和错误处理

• 使用 NF_CLOBBER 创建输出，同名 *_stream.nc 会被覆盖。
• 程序整块读取 u/v/temp/rho/ps，内存随 lon*lat*alt*time 增长。
• 必需变量或维度缺失时直接 stop。
• aire 和 phisinit 可从同目录 diagfi.nc 或 diagfi1.nc 回退读取；cv 缺失时会额外交互读取一个文件名，再回退 start.nc。
• 源码行 852 使用 ierr=nf_close(nid) 关闭输入，但主路径打开的是 infid；nid 未在主路径赋值。待确认：运行时是否因此出现 NetCDF close 错误或仅被忽略。

核心逻辑

1. 读坐标和必需字段：读取 latitude/longitude/Time/altitude、aps/bps、ps、u/v、aire、phisinit、cv、temp 和可选 rho。
2. 密度回退：若无 rho，用 hybrid 压力和理想气体近似：

rho(lon,lat,alt,time) = (aps(alt) + bps(alt) * ps(lon,lat,time)) / (Rgas * temp(lon,lat,alt,time))

3. 时间平均：对每个 (lon,lat,alt) 计算 u/v/temp/rho 的时间平均，对每个 (lon,lat) 计算 ps 的时间平均。
4. 经向风转换：对 ilat=1..latlength-1，使用相邻纬度平均的经向风除以 cv：

vcont(lon,lat,alt) = 0.5 * (vcum(lon,lat,alt) + vcum(lon,lat+1,alt)) / cv(lon,lat)
vcont(lon,latlength,alt) = 0

5. hybrid 压力和层质量：

Plev = aps + bps * pscum
mass = aire * Plev / g0
dmass(alt) = mass(alt) - mass(alt+1)
dmass(top) = 0

6. 流函数积分：从顶层向下积分，每层先继承上一层 psi，再减去所有经度上的 vcont*(dmass(lat)+dmass(lat+1))。源码还计算 vcontcum，但后续未写出。
7. 角动量：先对每个经度计算：

mom = a0*cos(lat) * (omega*a0*cos(lat) + u_time_mean)

随后又用纬向平均风 uzm 重算 momave，因此最终输出的 momave 对应 zonal mean wind 路径。 8. 纬向平均诊断：对 u/v/temp/rho/ps/phisinit 沿经度平均。 9. 写输出：创建单经度、单时间输出文件，写 psi/momave/u/v/rho/temp/ps/phisinit。

伪代码

read infile from stdin
open infile
read lon, lat, alt, Time, aps, bps, ps, u, v
read aire and phisinit, falling back to diagfi.nc or diagfi1.nc
read cv, asking for a start-like file if missing
read temp
if rho exists:
  read rho
else:
  rho = hybrid_pressure / (Rgas * temp)
average u, v, temp, rho, ps over all input times
vcont = adjacent-lat average of v / cv
Plev = aps + bps * mean_ps
mass = aire * Plev / g0
dmass = vertical difference of mass
integrate psi downward from top using vcont and neighboring-lat dmass
compute zonal means of u, v, temp, rho, ps, phisinit
momave = a0*cos(lat) * (omega*a0*cos(lat) + zonal_mean_u)
create *_stream.nc with lon=0 and Time=0
write psi, momave, u, v, rho, temp, ps, phisinit

参与的主题流程

复现要点

• 输入必须包含 u/v/temp/ps/aps/bps，否则程序停止；rho 可以缺失但会用常数 Rgas 估算。
• aire、phisinit 和 cv 对结果不是可选物理量；若主文件没有，需要同目录有 diagfi.nc、diagfi1.nc 或用户提供的 start 类文件。
• 输出纬度维和高度维保留输入坐标，但经度和时间均压缩为单点 0。
• psi 的单位属性为空；复现时应从 mass/g0 与 vcont 的源码公式推导实际量纲。
• dmass(top)=0、vcont(latlength)=0 和 psi(latlength,:)=0 是源码边界条件。

待确认

• cv 的精确定义和单位需结合生成该字段的 start/archive 代码确认。
• psi 是否缺少对经度宽度、纬向平均或传统流函数归一化因子的显式处理，需要和目标诊断定义比对。
• momave 输出是否应乘层质量或层厚；源码当前输出的是每单位质量的绝对角动量形式。
• 结尾 nf_close(nid) 是否是 infid 拼写错误。
• 输出变量定义为 4D，但写入数组为 2D 或 1D 纬向平均数组；需用实际 NetCDF 库接口确认写入布局是否总是按预期。

复现风险

• 同名 _stream.nc 会被覆盖。
• 大文件全量读入内存。
• 缺少 cv 时会阻塞等待交互输入。
• 主路径没有复制输入文件的时间范围信息，输出 Time=0 只表示平均诊断。
• 行星常数由编译宏控制；跨 Mars/Jupiter 编译时结果会不同。

相关页面

• util/index.md - util 后处理工具总览。
• concatnc.md - 常见上游时间拼接工具。
• localtime.md - 另一个时间相关后处理工具。
• lslin.md - Ls 线性化工具。