xvik.F

路径

LMDZ.MARS\libf\dynphy_lonlat\phymars\xvik.F

所属目录/模块

libf/dynphy_lonlat/phymars

文件定位

xvik.F 定义交互式后处理程序 xvik。它逐个读取用户给出的 diagfi NetCDF 文件，提取 Viking Lander 1、Viking Lander 2 和 InSight 站点附近的地表压力，按真实站点高度做静力订正，同时输出全球平均地表压、南北极 CO2 冰等效压力、每日平均压力以及最多 6 阶的季节谐波拟合。它依赖 readhead_NC 从 diagfi 文件读取网格头信息和地表位势，再用 ps/temp/co2ice/Time 做诊断。

定义的符号

符号	类型	行号	作用
`xvik`	program	1	交互式读取 diagfi 文件并生成 Viking/InSight 压力、全球压力/CO2 库存、日平均和谐波输出。
`sol2ls`	subroutine	827	将 sol 日期转换为太阳经度 `Ls`，可累计跨年 360 度。
`ls2sol`	function	906	将太阳经度 `Ls` 反算到火星年内 sol，用于按 Ls 输入的日平均阶段。
`DiscreetFourierHn`	subroutine	974	计算离散周期函数的第 `n` 阶 Fourier 系数 `a/b`。

依赖的模块

use 模块	only 列表	用途	待确认
`filtreg_mod`	`inifilr`	被导入但本文件未直接调用。	是否为历史残留。
`comconst_mod`	`dtvr,g,r,pi`	使用/打印动力时间步、重力、气体常数和圆周率；`pi` 在程序开始重设。	无。
`comvert_mod`	`pa,preff`	初始化 `pa=20`、`preff=610.`，配合读头/动力垂直坐标状态。	`pa` 具体意义需结合垂直坐标模块确认。

include 依赖

include	用途
`dimensions.h`	提供 `iim/jjm/llm`。
`paramet.h`	提供 `iip1/jjp1/ip1jmp1` 等派生维度。
`comdissip.h`	动力耗散公共块；本文件未直接读写。
`comgeom2.h`	提供 `rlonu/rlatv/rlonv/rlatu/aire` 等二维几何数组。
`netcdf.inc`	提供 NetCDF API。

调用的关键例程

被调用例程	所在模块/文件	调用位置	作用
`readhead_NC`	`readhead_NC`	行 229	读取 diagfi 头信息、几何、地表位势、`day0` 和 `constR`。
`conf_gcm`	外部配置例程	行 233	读取/设置 GCM 配置。
`iniconst`、`inigeom`	动力初始化链	行 234-235	初始化常数和几何状态。
`NF_OPEN/NF_INQ_DIMID/NF_INQ_DIMLEN/NF_INQ_VARID/NF_GET_VAR*`	NetCDF API	多处	打开 diagfi，读取 `Time`、`ps`、`temp`/`temp7`、`co2ice`。
`SSUM`	外部求和函数	行 322	计算网格总面积倒数 `airtot1`。
`sol2ls`	本文件	多处	将 sol 转为 Ls 输出或谐波重建。
`ls2sol`	本文件	日平均阶段	当输出时间轴为 Ls 时反算 sol 用于按日分组。
`DiscreetFourierHn`	本文件	谐波阶段	对日平均压力序列计算 0-6 阶 Fourier 系数。

输入

输入	来源	类型/维度	单位	含义
diagfi 目录	标准输入	path	-	用户输入的 diagfi 文件目录。
输出目录	标准输入	path	-	用户输入 `pathsor`，但源码后续未用于 `OPEN` 文件路径。
`Time_unit`	标准输入	integer	-	`1` 输出 sol，`2` 输出 Ls，其他值同时输出 sol 和 Ls。
diagfi 文件名	标准输入循环	name without `.nc`	-	每次循环读取一个 `nomfich`，空行结束。
NetCDF `Time`	diagfi	1D	sol	时间坐标。
NetCDF `ps`	diagfi	`iip1,jjp1,Time`	Pa	地表压力。
NetCDF `temp` 或 `temp7`	diagfi	3D 或 2D	K	用第 7 层温度做高度订正；无 `temp` 时尝试 `temp7`。
NetCDF `co2ice`	diagfi	`iip1,jjp1,Time`	kg m-2（推断）	CO2 冰面密度，用于南北极等效压力。

输出

输出	去向	类型	含义
`ps_VL1_yearN`、`ps_VL2_yearN`、`ps_INS_yearN`	当前工作目录	文本	每个站点逐时间输出真实高度订正压力 `zp2` 与 GCM 高度压力 `zp1`。
`prestot_yearN`	当前工作目录	文本	全球平均地表压、南北极 CO2 冰等效压力和总 CO2 等效压力。
`ps_*_yearN_diurnal`	当前工作目录	文本	按整数 sol 聚合的站点日平均压力。
`ps_*_yearN_harmonics`	当前工作目录	文本	若年内覆盖到 `Ls>350`，输出 0-6 阶 Fourier 系数和 669 sol 重建序列。

共享状态与副作用

交互式从标准输入读取路径、输出时间轴和多个 diagfi 文件名。
打开并关闭多个文本文件；pathsor 变量被读取但未拼接到 OPEN(file=...)，实际输出在当前工作目录。
修改 pi、pa、preff 等模块变量。
使用固定站点经纬度和高度：Viking 1、Viking 2、InSight；站点高度先按 3.73 转为位势。
第一条时间记录不输出站点压力，因为 firstcal 初始为 .true.，站点循环在 .NOT.firstcal 后执行。

核心逻辑

设置 pi、pa=20、preff=610.，定义火星年长度 unanj=669.。
初始化三站点经纬度和真实高度，打印站点信息并把经纬度转为弧度。
从标准输入读取 diagfi 目录、输出目录、时间轴选项和第一个文件名。
对每个非空文件名，打开 pathchmp/nomfich.nc，调用 readhead_NC 读取头信息和 phis，再初始化 GCM 常数/几何。
读取 Time 维长度和时间数组。
对三个站点，利用 rlonu/rlatv 找到包围站点的 2x2 GCM 网格，计算双线性权重 zw(0:1,0:1,iv)。
用权重插值得到每个站点的模型地表位势 phisim。
对每个时间记录读取 ps、temp 或 temp7、co2ice。
计算当前 sol 和火星年编号；跨年时打开新的站点输出文件和全球压力输出文件，按 Time_unit 写表头。
用面积加权得到全球平均地表压；按南北半球累加 co2ice 并乘 g 得等效压力。
对每个站点，在 2x2 网格上对 log(ps) 做双线性插值，得到 GCM 高度压力 zp1；用第 7 层温度估计 gh=constR*t7，将压力订正到真实站点高度 zp2。无温度时使用 10 km 标高的备用尺度高度。
写逐时站点文件和 prestot_yearN。
所有 diagfi 文件处理完后，读取逐时站点文件，按 floor(sol) 分组写日平均文件。
读取日平均文件；若年内 Ls 覆盖到 350 度以上，计算 0-6 阶 Fourier 系数，并重建 1-669 sol 的谐波压力序列。

伪代码

read input directory, output directory, time axis option
read first diagfi basename
while basename is non-empty:
    open basename.nc
    readhead_NC -> day0, phis, constR, geometry
    read Time
    locate 2x2 GCM cells around Viking1, Viking2, Insight
    interpolate model topography at sites
    for each Time record:
        read ps, temp/temp7, co2ice
        open new yearly output files if needed
        compute planetary mean pressure and CO2 ice equivalents
        if not first record:
            interpolate log(ps) around each site
            correct pressure from model altitude to true site altitude
            write station pressure files
        increment year at sol ~= 669
    close NetCDF
    read next basename

for each produced yearly station file:
    group rows by integer sol and write diurnal means

for each diurnal file:
    if year is complete enough:
        compute Fourier coefficients n=0..6
        write coefficients and reconstructed annual curve

参与的主题流程

主题	参与方式
Viking/InSight 压力后处理	从 diagfi 网格输出生成站点压力、日平均和谐波拟合。
CO2 库存诊断	输出全球平均地表压、南北极 CO2 冰等效压力和总 CO2 等效压力。
`util/xvik` 拟合链	`util/xvik/fit_Iceinertia_MONSicedepth.F` 的配置和 README 说明它消费 xvik 输出。

写法特点

固定格式 Fortran，主程序后直接附带三个辅助例程。
输出路径变量 pathsor 被读取但没有实际用于文件名。
sol2ls 使用 year_day=669.，ls2sol 使用 year_day=668.6d0，两者常数略不一致。
temp 缺失时不终止，而是尝试 temp7；若也不存在，则温度置零并用 10 km 标高备用公式。
谐波拟合只在最后一个日平均点 Ls>350 时执行，否则写 “not computed because year is not complete”。

复现要点

输入 diagfi 文件必须至少有 Time、ps、co2ice，并应有 temp 或 temp7。
readhead_NC 要求同一文件中有 controle、经纬度、面积、phisinit、aps/bps。
需要从标准输入依次提供目录、输出时间轴和一个或多个文件名。
输出文件实际写到当前工作目录；运行前应切到希望保存结果的目录，或在源码/脚本层处理路径。

待确认

pathsor 未使用是否为 bug 或历史残留。
sol2ls 和 ls2sol 火星年常数不一致是否会影响 Time_unit==2 的日平均分组。
util/xvik/fit_Iceinertia_MONSicedepth.F 文档提到 xpsol/xprestot，而当前 xvik.F 写 ps_VL*_yearN/prestot_yearN，两者文件名约定是否已分叉。
firstcal 跳过首个站点压力输出的原因需结合诊断时间轴约定确认。