fit_Iceinertia_MONSicedepth.F

路径

LMDZ.MARS\util\xvik\fit_Iceinertia_MONSicedepth.F

所属目录/模块

util/xvik

文件定位

fit_Iceinertia_MONSicedepth.F 定义 Fortran 程序 SUPERFIT，是 util/xvik 下的 Viking Lander 1 压力拟合辅助程序。它读取四组 xvik.F 输出，分别代表冰深系数和冰热惯量的最小/最大组合；再用本目录的 VL1 观测曲线作为约束，暴力扫描北/南半球冰热惯量、北/南半球冰深系数和总 CO2 压力，寻找 VL1 压力均方差最小的组合。

本文件是历史拟合脚本，不是通用 NetCDF 工具。源码顶部明确写有“THE FOLLOWING LINES MUST BE ADAPTED”，且数组长度、搜索步长、压力范围和文件命名都硬编码在源码或 .def 输入中。

定义的符号

符号	类型	行号	作用
`SUPERFIT`	program	1	主程序：读取四组 xvik 输出和 `VL1`，扫描 CO2 cycle 参数并写出拟合结果。
`RUNAVE`	subroutine	698	对周期坐标上的数据做滑动平均。
`ls2sol`	subroutine	753	把 Ls 转换为 sol；用于读取 Ls-only xvik 输出时重建 sol 坐标。
`sol2ls`	subroutine	798	把 sol 转换为 Ls；用于最终 `pvlfit` 和 `xprestotfit` 输出。

配套文件

文件	作用
`README`	说明本目录 5 个文件、编译脚本、输入样例和输出文件。
`compile_fit`	用 `gfortran fit_Iceinertia_MONSicedepth.F -o fit_Iceinertia_MONSicedepth.e` 编译。
`fit_Iceinertia_MONSicedepth.def`	标准输入样例，包含冰深/热惯量范围、四组 xvik 输出路径、年份、时间格式和 plot 模式。
`VL1`	669 行 Viking Lander 1 压力观测/谐波曲线，程序要求与可执行文件在同一目录。

输入

输入	来源	类型/维度	单位	含义
`icedmin` / `icedmax`	stdin / `.def`	scalar	coefficient	冰深系数搜索下界/上界；示例为 `3.2e-4` 到 `1.3e-3`。
`iceimin` / `iceimax`	stdin / `.def`	scalar	thermal inertia	冰热惯量搜索下界/上界；示例为 `500` 到 `2000`。
`dset_DminImin`	stdin / `.def`	path	-	最小冰深 + 最小冰热惯量的 xvik 输出目录。
`dset_DminImax`	stdin / `.def`	path	-	最小冰深 + 最大冰热惯量的 xvik 输出目录。
`dset_DmaxImin`	stdin / `.def`	path	-	最大冰深 + 最小冰热惯量的 xvik 输出目录。
`dset_DmaxImax`	stdin / `.def`	path	-	最大冰深 + 最大冰热惯量的 xvik 输出目录。
`year_xvik`	stdin / `.def`	integer	-	追加到 `xpsol1N` 和 `xprestotN` 文件名的年份编号。
`time_unit`	stdin / `.def`	1/2/3	-	xvik 文件列格式：1=sol，2=Ls，3=sol+Ls。
`plot_test`	stdin / `.def`	1/2	-	`minimization.txt` 输出切片：1 写 `COST(DN,DS)`，2 写 `COST(IN,IS)`。
`xpsol1N`	四组 xvik 输出目录	text table	Pa	VL1 站点模拟压力；列格式由 `time_unit` 决定。
`xprestotN`	四组 xvik 输出目录	text table	Pa	`patm, pcapn, pcaps`；列格式由 `time_unit` 决定。
`VL1`	当前工作目录	669 行 text table	sol, Ls, Pa	VL1 观测压力曲线。

硬编码参数

参数	源码值	影响
`ngcm`	8028	每组 xvik 原始记录数；注释保留 10704 和 2672 的替代值。
`nvl1` / `nsol`	669 / 669	VL1 观测长度和每日平滑输出长度。
`refrun`	1	以第一组 `Dmin/Imin` 作为线性化参考。
`maxiceidiff`	3000	允许南北冰热惯量差异的上限。
`deltaicei`	20	冰热惯量扫描步长。
`maxiceddiff`	`30.e-4`	允许南北冰深系数差异的上限。
`deltaiced`	`0.5e-4`	冰深系数扫描步长。
`pmin/pmax/deltap`	`690/720/1.0`	总 CO2 压力扫描范围和步长。
平滑窗口	20 sol、60 sol	原始压力用 20 sol 运行平均；敏感度导数再用 60 sol 平滑。

输出

输出	去向	类型/维度	含义
`minimization.txt`	当前工作目录	text table	每个切片上的最小 RMS、对应另一组参数和总压力；格式由 `plot_test` 决定。
stdout	终端	text	最优 `Ptot/In/IS/Dn/DS` 和 RMS。
`pvlfit`	当前工作目录	669 行 text table	最优参数下的合成 VL1 压力：`sol, Ls, pvl1_fit`。
`xprestotfit`	当前工作目录	669 行 text table	最优参数下的 CO2 库存分解：`sol, Ls, patm, pcapn, pcaps, pfit`。

核心逻辑

读搜索范围和四组输出目录：从 stdin 读取冰深系数、冰热惯量范围、四个 xvik 输出目录、年份、时间列格式和 plot_test。
构造四个实验角点：
- run 1: Dmin/Imin
- run 2: Dmin/Imax
- run 3: Dmax/Imin
- run 4: Dmax/Imax
打开 xvik 输出文件：按 write(filename1,'(a6,i1)') 'xpsol1',year_xvik 和 write(filename2,'(a8,i1)') 'xprestot',year_xvik 构造文件名。
读观测和模拟数据：
- VL1 提供 solvl, lsvl, pvl_obs。
- xpsol1N 读入四组 pvl_gcm。
- xprestotN 读入四组 patm, pcapn, pcaps。
- 若 time_unit=2，调用 ls2sol 把 Ls 转为 sol。
总 CO2 压力一致性检查：每组 run 在首个时间点定义 ptot，后续若 patm+pcapn+pcaps 偏离超过 3 Pa 则打印警告。
平滑到每日 sol 网格：RUNAVE 用周期 669 sol 和 20 sol 窗口把四组 VL1 压力、patm、pcapn、pcaps 平滑到 sol=1..669。
估计敏感度导数：
- dpdicein/dpdiceis 由 run 1-2 与 3-4 差分估计北/南 CO2 cap 对冰热惯量的敏感度。
- dpden/dpdes 由 run 1-3 与 2-4 差分估计北/南 CO2 cap 对冰深系数的敏感度。
- alphavl1 是四组 pvl_sm/patm_sm 的平均。
- 四条导数再用 60 sol 窗口平滑。
暴力扫描最优参数：
- plot_test=1 时外层扫描 DN,DS，内层扫描 IN,IS,Ptot，输出适合画 COST(DN,DS) 的表。
- plot_test=2 时外层扫描 IN,IS，内层扫描 DN,DS,Ptot，输出适合画 COST(IN,IS) 的表。
- 每个候选点用线性化公式构造 pcapn_new/pcaps_new，并限制为非负。
- VL1 压力预测为 alphavl1*(ptry-pcapn_new-pcaps_new)。
- cost 是 669 个 sol 上观测和预测压力差平方和。
写最优合成曲线：用最优参数重新生成 pvlfit 和 xprestotfit，并调用 sol2ls 输出 Ls。

伪代码

read icedmin, icedmax, iceimin, iceimax
read four xvik output directories, year_xvik, time_unit, plot_test
assign run1..run4 to D/I corner combinations
open xpsol1N and xprestotN in each directory
read VL1 observed pressure curve
for n in 1..ngcm:
  read four simulated VL1 pressure series
  read four patm/pcapn/pcaps series
  convert Ls to sol if needed
  check total CO2 pressure consistency
smooth pvl, patm, pcapn, pcaps to 669-sol grid
compute sensitivity derivatives to ice inertia and ice depth
smooth derivatives with 60-sol window
for each allowed DN, DS, IN, IS, Ptot:
  construct pcapn_new and pcaps_new from reference run and sensitivities
  pvl1 = alphavl1 * (Ptot - pcapn_new - pcaps_new)
  cost = sum((VL1_obs - pvl1)^2)
  update global best fit and plot-slice best fit
write minimization.txt
write pvlfit and xprestotfit for best fit

与 xvik.F 的关系

项目	源码事实
预期上游	本程序的提示、README 和 `.def` 都说输入来自 `xvik.F` 输出。
当前文件名	本程序硬编码读取 `xpsol1N` 和 `xprestotN`。
已有 xvik 页面记录	xvik.F 页面记录当前源码输出 `ps_VL1_yearN`、`ps_VL2_yearN`、`ps_INS_yearN` 和 `prestot_yearN`。
兼容性风险	`xpsol1N/xprestotN` 与当前 `xvik.F` 页面记录的命名不一致，可能是历史分支或需重命名/预处理。

复现要点

编译用 util/xvik/compile_fit，不会链接 NetCDF；本程序只读写 ASCII 文本。
运行目录必须能读取 VL1，否则 open(30,file='VL1') 没有显式 iostat 检查，会在后续 read 路径失败。
四个 xvik 输出目录必须分别代表 Dmin/Imin、Dmin/Imax、Dmax/Imin、Dmax/Imax。
year_xvik 只按一位整数拼入文件名；多位年份编号会被格式截断或不匹配。
ngcm=8028 必须与输入 xvik 文件行数一致，不一致会读错或提前 EOF。
plot_test 只接受 1 或 2；其他值会停止。

待确认

xpsol1N/xprestotN 与当前 xvik.F 输出 ps_VL*_yearN/prestot_yearN 的历史对应关系。
fit 逻辑里仍有“cap albedo”旧注释，但当前实际扫描的是冰热惯量、冰深系数和总压力；是否保留 albedo 拟合路径需查历史版本。
fonc、fonc2n、fonc2s 当前都是恒等函数，旧注释中的指数/log 变换是否仍有科学用途需确认。
VL1 文件第三列是压力曲线；README 称其为 VL1 surface pressure harmonics，但实际文件内容需结合生成流程确认。

复现风险

硬编码数组长度和搜索步长需要针对每组数据手工修改。
RUNAVE 内部 xx(nbig) 使用 99999999 个 real，内存占用约数百 MB。
minimization.txt、pvlfit、xprestotfit 会覆盖同名文件。
搜索是五重循环，步长变小会显著增加运行时间。
输入文件打开大多有 iostat 检查，但 VL1 没有显式打开失败检查。