dyn1d/disvert_noterre.F

路径

LMDZ.MARS\libf\phymars\dyn1d\disvert_noterre.F

文件定位

dyn1d/disvert_noterre.F 在当前源码树中不是独立实现，而是一行相对路径：

../../../../LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F

它把 1D testphys1d 编译路径中的 disvert_noterre 和内部辅助 sig_hybrid 解析到公共非地球垂直离散例程。该例程读取 esasig.def 或 z2sig.def，根据 hybrid 开关生成 comvert_mod 中的 ap/bp/aps/bps/presnivs/pseudoalt，供 dyn1d 后续构造 plev/play。

定义的符号

符号	类型	行号	作用
`disvert_noterre.F`	path stub	dyn1d line 1	指向 `LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F`。
`disvert_noterre`	subroutine	common target line 2	读取非地球垂直坐标输入并写入 `comvert_mod` 的 hybrid/sigma 系数。
`sig_hybrid`	subroutine	common target line 277	从目标 sigma 层位置反解 hybrid 坐标中的 `newsig`。

依赖模块和输入状态

依赖	使用符号	作用
`IOIPSL` 或 `ioipsl_getincom`	`getin`	读取 `hybrid` 配置键。
`comvert_mod`	`ap`, `bp`, `aps`, `bps`, `presnivs`, `pseudoalt`, `nivsig`, `nivsigs`, `pa`, `preff`, `scaleheight`	读写垂直坐标共享状态。
`comconst_mod`	`kappa`	`esasig.def` 分支中构造能量守恒相关的 `s` 权重。
`logic_mod`	`hybrid`	控制使用 hybrid 坐标还是纯 sigma 坐标。
include	`dimensions.h`, `paramet.h`, `iniprint.h`	提供 `llm/llmp1`、维度参数和 `lunout` 输出单元。

输入文件选择

例程按固定顺序打开输入文件：

先尝试 esasig.def。
若失败，关闭单元 99 后尝试 z2sig.def。
若两个文件都不存在，向 lunout 写出缺文件提示并 stop。

文件	读取内容	生成方式
`esasig.def`	`scaleheight`, `dz0`, `dz1`, `nhaut`	先把 `dz0/dz1` 除以 `scaleheight`，通过指数/tanh 组合生成 `sig(1:llm+1)`，并用 `kappa` 构造归一化权重 `s`。
`z2sig.def`	`scaleheight` 后跟 `llm` 行 `zsig(l)`	`sig(1)=1`，`sig(l)=0.5*(exp(-zsig(l)/scaleheight)+exp(-zsig(l-1)/scaleheight))`，`sig(llm+1)=0`。

testphys1d.F90 的文件头注释也提示 1D 运行需要一个描述 sigma layers 的文件，例如 z2sig.def。

主算法流程

设置 hybrid=.true.，再通过 getin('hybrid', hybrid) 允许配置覆盖，并打印结果。
读取 esasig.def 或 z2sig.def，生成界面 sigma 数组 sig(1:llm+1)。
写入 nivsigs(l)=real(l) 和 nivsig(l)=real(l) 作为层号坐标。
若 hybrid=.true.：
- 对每个 l=1,llm 调用 sig_hybrid(sig(l), pa, preff, newsig)。
- 写入 bp(l)=exp(1 - 1/newsig**2)。
- 写入 ap(l)=pa*(newsig-bp(l))。
- 顶界 ap(llmp1)=0、bp(llmp1)=0。
若 hybrid=.false.：
- 写入 ap(l)=0、bp(l)=sig(l)。
- 顶界 ap(llmp1)=0、bp(llmp1)=0。
对 l=1,llm-1，层中系数取相邻界面平均：
- aps(l)=0.5*(ap(l)+ap(l+1))
- bps(l)=0.5*(bp(l)+bp(l+1))
顶层层中值单独外推：
- hybrid: aps(llm)=aps(llm-1)**2/aps(llm-2)，bps(llm)=0.5*(bp(llm)+bp(llm+1))
- sigma: bps(llm)=bps(llm-1)**2/bps(llm-2)，aps(llm)=0
生成参考层压和伪高度：
- presnivs(l)=aps(l)+bps(l)*preff
- pseudoalt(l)=-scaleheight*log(presnivs(l)/preff)

`sig_hybrid` 求解

sig_hybrid(sig, pa, preff, newsig) 要求解：

(1 - pa/preff) * exp(1 - 1/newsig**2) + (pa/preff) * newsig = sig

源码中的分支：

条件	行为
`sig >= 1`	直接 `newsig=sig`。
`sig*preff/pa >= 0.25`	在 `[0,1]` 上用最多 9999 次二分迭代，根据 `F` 相对 1 的大小移动上下界；当 `abs(10*log(F)) < 1.E-5` 时退出。
其他	使用近似 `newsig=sig*preff/pa`。

这个辅助例程没有显式错误返回；若参数导致收敛慢，循环结束后使用最后一次 newsig。

dyn1d 调用关系

调用方	行号	调用	作用
`init_testphys1d_mod.F90`	337	`call disvert_noterre`	在读取 `psurf/pa/preff/hybrid` 后生成垂直坐标数组。

调用之后，init_testphys1d_mod.F90 第 339 行把 ap/bp/aps/bps/presnivs/pseudoalt 复制到 vertical_layers_mod，第 341-342 行构造 plev=ap+psurf*bp 和 play=aps+psurf*bps。

输出和副作用

项目	说明
写入模块状态	`comvert_mod::ap/bp/aps/bps/presnivs/pseudoalt/nivsig/nivsigs/scaleheight`。
读取模块状态	`comvert_mod::pa/preff`、`comconst_mod::kappa`、`logic_mod::hybrid`。
文件 I/O	读取 `esasig.def` 或 `z2sig.def`；向 `lunout` 打印分支、系数和层压/伪高度。
程序控制	两个垂直坐标文件都缺失时执行 `stop`。
诊断文件	`testhybrid.tab` 生成块全部被注释，不会实际写出。

边界条件和风险

llm 必须至少支持顶层外推里访问 llm-2；源码没有对小层数做保护。
hybrid 分支依赖 pa/preff 的物理合理性；pa=0 或异常参考压力会破坏 sig_hybrid 分支判定。
z2sig.def 分支期望恰好读取 llm 行高度；文件过短会触发 Fortran 读错误。
esasig.def 优先级高于 z2sig.def；两者同时存在时不会读取 z2sig.def。

复现要点

准备 esasig.def 或 z2sig.def，并确认工作目录是 disvert_noterre 打开这些文件的目录。
先设置 pa/preff 和 hybrid，再调用 disvert_noterre。
用 ap/bp/aps/bps 构造压力层时，界面和层中公式分别为 ap+ps*bp 与 aps+ps*bps。
若要完全复现 dyn1d，继续执行 init_vertical_layers 把这些状态同步到物理侧 vertical_layers_mod。

待确认

当前 Windows 工作树把链接文件呈现为一行路径文本；在原始类 Unix 工作树中它可能是符号链接。
esasig.def 分支中的 s/dsig 权重当前是否被后续代码使用；本例程生成 s 后未见继续写出到模块状态。
dyn1d 运行目录中默认提供的是 z2sig.def 还是 esasig.def，需结合具体实验目录确认。