gwstress_mod.F90

路径

LMDZ.MARS\libf\phymars\gwstress_mod.F90

所属目录 / 模块

libf\phymars

文件定位

gwstress_mod.F90 定义 gwstress_mod 模块和 GWSTRESS 子程序，是 orodrag_mod.F90::ORODRAG 中地形重力波参数化链的低层应力初始化步骤。本例程在 OROSETUP 计算出 Brunt-Väisälä 频率、低层风速和地形统计量之后、GWPROFIL 重算应力廓线之前被调用。

本例程的职责是按 Lott & Miller (1997) 方程 (17) 计算地形重力波在大气底层（半层 nlayer+1）的初始应力 ZTAU(:,nlayer+1)。它只设置底层一个半层的应力值，不计算垂直廓线，也不生成风或温度 tendency。源码注释表明它源自 ECMWF IFS 地形重力波拖曳方案，由 F. Lott 于 1993 年移植到 IFS，2022 年由 J. Liu 重写。

定义的符号

符号	类型	行号	作用
`gwstress_mod`	module	1	封装地形重力波低层应力计算例程
`GWSTRESS`	subroutine	7	按方程 (17) 计算底层半层 `ZTAU(:,nlayer+1)` 的地形重力波应力

依赖的模块

use 模块	only 列表	用途	待确认
`dimradmars_mod`	`ndomainsz`	固定局部数组第一维；接口参数 `ktest`、`ICRIT`、`IKNU`、`ZVPH` 等均按 `ndomainsz` 声明	-
`yoegwd_h`	`gkdrag`, `gtsec`, `gvcrit`	`gkdrag` 进入方程 (17) 的应力公式；`gtsec` 是小应力截断阈值（当前未生效）；`gvcrit` 是低层风速截断阈值（当前未生效）	-

调用的关键例程

被调用例程	所在模块 / 文件	调用位置	作用
无	-	-	本例程只使用内在函数 `AMAX1`、`sqrt` 和数组循环，不调用其他用户例程

输入

输入	来源	类型/维度	单位	含义
`ngrid`	`ORODRAG`	integer scalar	grid points	大气列数
`nlayer`	`ORODRAG`	integer scalar	layers	大气全层数；半层数组使用 `nlayer+1`
`ktest`	`calldrag_noro` 筛选结果	integer `(ndomainsz)`	-	GW 点开关图；`ktest(jl)==1` 时该点参与计算
`ZRHO`	`OROSETUP`	real `(ndomainsz,nlayer+1)`	kg/m³	半层密度；只读取 `ZRHO(:,nlayer+1)`
`BV`	`OROSETUP`	real `(ndomainsz,nlayer+1)`	s⁻¹	Brunt-Väisälä 频率；只读取 `BV(:,nlayer+1)`
`pvar`	`surfdat_h:zstd` 经上游传入	real `(ndomainsz)`	m	次网格地形标准差
`psig`	`surfdat_h:zsig` 经上游传入	real `(ndomainsz)`	-	次网格坡度统计量
`zgeom`	`drag_noro` 构造的翻转位势	real `(ndomainsz,nlayer)`	-	全层位势高度；本例程读取但实际未参与当前可执行路径
`zdmod`	`OROSETUP`	real `(ndomainsz)`	-	地形各向异性组合幅度 √(τ₁²+τ₂²)；参见方程 (17)/(18)
`ICRIT`	`OROSETUP`	integer `(ndomainsz)`	level index	临界层位置；接口保留，当前可执行路径中未使用
`IKNU`	`OROSETUP`	integer `(ndomainsz)`	level index	`4*pvar` 对应层；接口保留，当前可执行路径中未使用
`ZVPH`	`OROSETUP`	real `(ndomainsz,nlayer+1)`	m/s	半层低层风速 U_H；只读取 `ZVPH(:,nlayer+1)`

输出

输出	去向	类型/维度	单位	含义
`ZTAU`	返回 `ORODRAG`，随后传给 `GWPROFIL`	real `(ndomainsz,nlayer+1)`, intent(inout)	stress	地形重力波半层应力；本例程只设置 `ZTAU(:,nlayer+1)` 即底层半层的初始应力值

共享状态与副作用

本文件没有 save 变量、COMMON、THREADPRIVATE、文件 I/O、配置读取或诊断输出。
唯一副作用是原地改写 ZTAU(:,nlayer+1)。ZTAU 是 intent(inout)，调用前由 OROSETUP 初始化过；本例程覆盖底层半层的值。
lo 逻辑变量在当前可执行路径中被赋值但从未读取；对应的 IF(LO) ZTAU=0.0 已被注释掉。
ICRIT、IKNU、zgeom 作为接口参数存在但不参与当前可执行计算路径。

核心逻辑

设定水平循环范围 kidia=1、kfdia=ngrid。
对每个水平格点 jl，若 ktest(jl)==1（地形 GW 活跃点）：
- 将阻塞层高度 ZBLOCK 设为 0.0（IKENVH 相关分支已注释掉）。
- 有效山高 ZEFF = max(0, 2*pvar - ZBLOCK)。
- 底层应力公式（方程 17）：
```
ZTAU(jl,nlayer+1) = ZRHO * GKDRAG * psig * ZEFF² / 4 / pvar
                    * ZVPH * zdmod * sqrt(BV)
```
- 评估截断条件（当前被注释掉，不生效）：若应力太小 (<GTSEC)、或临界层太低 (ICRIT >= IKNU)、或低层风速太小 (<GVCRIT)，则置零。
若 ktest(jl) != 1，直接把 ZTAU(jl,nlayer+1) 设为 0.0。

伪代码

GWSTRESS(ngrid, nlayer, ktest, ZRHO, BV, pvar, psig, zgeom, zdmod,
         ICRIT, IKNU, ZVPH, ZTAU):
    kidia = 1
    kfdia = ngrid

    for jl in kidia..kfdia:
        if ktest(jl) == 1:
            ZBLOCK = 0.0       ! IKENVH branch commented out
            ZVAR = pvar(jl)
            ZEFF = max(0, 2*ZVAR - ZBLOCK)

            ! Equation (17): low-level gravity wave stress
            ZTAU(jl,nlayer+1) = ZRHO(jl,nlayer+1) * GKDRAG * psig(jl)
                                * ZEFF**2 / 4 / ZVAR
                                * ZVPH(jl,nlayer+1) * zdmod(jl)
                                * sqrt(BV(jl,nlayer+1))

            ! Truncation check (currently commented out, not active):
            ! LO = (ZTAU < GTSEC) or (ICRIT >= IKNU) or (ZVPH < GVCRIT)
            ! if LO: ZTAU(jl,nlayer+1) = 0.0
        else:
            ZTAU(jl,nlayer+1) = 0.0

参与的主题流程

主题	参与方式
地形重力波 / 次网格地形拖曳	`calldrag_noro` 筛选地形格点，`drag_noro` 调用 `ORODRAG`，`ORODRAG` 先调用 `OROSETUP`，再调用本例程初始化底层应力，最后调用 `GWPROFIL` 重算廓线
主物理时间步 `physiq`	`physiq_mod.F` 在 `calllott` 为真时进入 `calldrag_noro` 链；本例程是该链内部的底层应力初始化环节
垂直坐标 / 层序转换	调用前 `drag_noro` 已翻转垂直层序；因此本文中的 `nlayer+1` 对应翻转后的最低半层（即物理大气的底层）

写法特点

.F90 自由格式模块，但保留 ECMWF/IFS 历史标签式注释（100 CONTINUE、300 CONTINUE）和未使用的 CONTINUE 标签。
接口包含大量注释掉的历史参数（IKCRIT、ISECT、IKHLIM、IKCRITH、IKENVH、PVAR1、pgam、zd1、zd2、znu、ZTFR），表明接口经历过精简。
ICRIT、IKNU、ZVPH 在接口注释中被标注为 "actually not used"，当前可执行路径确实未读取它们；zgeom 也未参与计算。
lo 逻辑变量和对应的 IF(LO) ZTAU=0.0 已被注释掉，但变量声明保留。这意味着当前版本不做底层应力截断。
AMAX1 是旧式 Fortran 内在函数，等价于现代的 MAX；保留它可能是为了与 IFS 源码风格一致。
方程 (17) 中 ZEFF² / 4 / pvar 的写法等价于 (2*pvar)² / 4 / pvar = pvar（当 ZBLOCK=0 时），但源码保留了通用形式以支持 ZBLOCK 非零的情况。

复现要点

复现时必须先让 OROSETUP 或等效步骤初始化 ZTAU 和 BV/ZVPH/ZRHO，否则本例程会读取未定义输入。
ZBLOCK 当前硬编码为 0.0；如果恢复 IKENVH 分支，需要确保 zgeom 和 RG（重力加速度）已正确传入。
ktest 必须由上游按子域筛选填入；对 ktest(jl)!=1 的格点，本例程直接把底层 ZTAU 设为 0.0，会覆盖 OROSETUP 的初始化值。
当前版本不做底层应力截断（IF(LO) 已注释）；如果 GWPROFIL 依赖截断行为，复现时需要确认是否需要恢复。
zdmod 来自 OROSETUP 中 zd1/zd2 的合成，代表地形各向异性组合幅度；它在方程 (17) 中是线性乘子。

待确认

注释掉的 IF(LO) ZTAU=0.0 截断逻辑是否在旧 IFS 版本中启用过、是否计划恢复，需结合 sugwd 和 IFS 历史版本确认。
IKENVH 阻塞层分支注释中引用了 RG（重力加速度），但当前 use 列表中没有引入；如果恢复该分支，需要额外 use comcstfi_h, only: g 或类似依赖。
ICRIT >= IKNU 截断条件的物理含义（临界层低于 4*pvar 层时截断）需结合 OROSETUP 的层索引定义继续核验。
AMAX1(0., 2*ZVAR-ZBLOCK) 中第一个参数 0. 的类型是默认 real；在 ZEFF 赋值后未做负值保护以外的其他检查。