LMDZ.MARS Wiki

gwprofil_mod.F90

路径

LMDZ.MARS\libf\phymars\gwprofil_mod.F90

所属目录 / 模块

libf\phymars

文件定位

gwprofil_mod.F90 定义 gwprofil_mod 模块和 GWPROFIL 子程序，是 orodrag_mod.F90::ORODRAG 中地形重力波参数化链的应力廓线重算步骤。上游 GWSTRESS 已给出低层或初始的 ZTAU，本例程沿垂直半层从高层向低层扫描，按 IKNU2、临界层 ICRIT、波 Richardson 数和 GRCRIT 破碎阈值改写 ZTAU(:,jk)，最后在阻塞层顶 IKENVH 与动力混合高度 IKCRITH 之间做压力线性插值。

它不计算风、温度 tendency，也不读取配置；下游 ORODRAG 用本例程返回的 ZTAU 垂直差分来计算 pdudt/pdvdt/pdtdt。源文件注释把本例程称为 stress profile for gravity waves，并说明其源自 Lott/IFS 地形重力波拖曳方案。

定义的符号

符号 类型 行号 作用
gwprofil_mod module 1 封装地形重力波 stress profile 计算例程。
GWPROFIL subroutine 7 根据阻塞层、临界层、稳定度和波破碎判据重算 ZTAU 半层应力廓线。

依赖的模块

use 模块 only 列表 用途 待确认
dimradmars_mod ndomainsz 固定局部子域数组第一维；kdx 把 kgwd 个有效地形 GW 点映射回 1..ndomainsz。 -
yoegwd_h gkdrag, grcrit, gssec, gtsec gkdrag 进入 ZNORM stress 归一化；grcrit 是波 Richardson 数破碎阈值；gssec/gtsec 是除零和小 stress 保护阈值。 -

调用的关键例程

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

输入

输入 来源 类型/维度 单位 含义
ngrid ORODRAG integer scalar grid points 大气列数；本例程接口保留但源码内未使用。
nlayer ORODRAG integer scalar layers 大气全层数；半层数组使用 nlayer+1。
kgwd calldrag_noro -> drag_noro -> ORODRAG integer scalar points 当前子域中进入地形 GW 计算的格点数。
kdx calldrag_noro 筛选结果 integer (ndomainsz) - 第 ji 个 GW 点在局部子域中的格点编号 jl。
ktest calldrag_noro 筛选结果 integer (ndomainsz) - GW 点开关图；接口保留但源码内未读取。
IKCRITH OROSETUP integer (ndomainsz) level index 重力波破碎的动力混合高度。
ICRIT OROSETUP integer (ndomainsz) level index 地形重力波遇到临界层的位置。
IKENVH OROSETUP integer (ndomainsz) level index 阻塞层顶。
IKNU OROSETUP integer (ndomainsz) level index 4*pvar 对应层；本例程只用来缓存 ZTAU(:,IKNU+1)。
IKNU2 OROSETUP integer (ndomainsz) level index 3*pvar 对应层；控制从高层常值段转入递推段。
pplev drag_noro 传入的翻转半层压力 real (ndomainsz,nlayer+1) Pa 半层压力；末段用于在 IKCRITH 到 IKENVH 之间压力线性插值 ZTAU。
ZRHO OROSETUP real (ndomainsz,nlayer+1) kg/m3，按上游语义 半层密度，用于 ZNORM。
BV OROSETUP real (ndomainsz,nlayer+1) s-1，按上游语义 Brunt-Vaisala 频率，用于 ZNORM 和波位移。
ZVPH OROSETUP real (ndomainsz,nlayer+1) m/s，按上游语义 半层低层风速或相位速度尺度，用于 ZALFA。
PRI OROSETUP real (ndomainsz,nlayer+1) - 半层平均流 Richardson 数。
zdmod OROSETUP real (ndomainsz) - 地形各向异性组合幅度；参与 zoro。
psig surfdat_h:zsig 经上游传入 real (ndomainsz) - 次网格坡度。
pvar surfdat_h:zstd 经上游传入 real (ndomainsz) m 次网格地形标准差；参与 zoro 并用 max(pvar,1.0) 保护。

输出

输出 去向 类型/维度 单位 含义
ZTAU 返回 ORODRAG real (ndomainsz,nlayer+1), intent(inout) stress，按上游 GW 链语义 地形重力波半层应力廓线；输入来自 GWSTRESS，输出被 ORODRAG 用于风和温度 tendency。

共享状态与副作用

核心逻辑

  1. 令 ilevh=nlayer/3，但该变量只出现在注释掉的打印片段中，不影响计算。
  2. 对每个 GW 点 jl=kdx(ji)，计算地形尺度：
    • zoro = psig * zdmod / 4 / max(pvar,1.0)
    • 缓存 ZTAU(jl,IKNU(jl)+1) 与顶半层 ZTAU(jl,nlayer+1) 到 Z_TAU。
  3. 从 jk=nlayer 向 2 递减扫描半层。
  4. 当 jk >= IKNU2(jl) 时，把 ZTAU(jl,jk) 设为顶半层缓存 Z_TAU(jl,nlayer+1)，形成上部常值 stress 段。
  5. 当 jk < IKNU2(jl) 时，计算局地归一化量与上一半层波位移：
    • ZNORM = gkdrag * ZRHO * sqrt(BV) * ZVPH * zoro
    • ZDZ2 = ZTAU(jl,jk+1) / max(ZNORM, gssec)
  6. 对 jk < IKNU2(jl) 的层，如果上一层 stress 小于 gtsec 或已到/低于临界层 jk <= ICRIT(jl)，直接把 ZTAU(jl,jk)=0。
  7. 否则用 PRI、BV、ZDZ2 和 ZVPH 计算波 Richardson 数 ZRIW。若 ZRIW < GRCRIT，按二次方程根得到新的 ZALPHA 和 ZDZ2N，用 ZNORM*ZDZ2N 限制破碎后的 stress；否则沿用 ZNORM*ZDZ2。
  8. 每层递推末尾执行 ZTAU(jl,jk)=min(ZTAU(jl,jk), ZTAU(jl,jk+1))，保证向下递推时 stress 不超过上一半层。
  9. 扫描完成后，缓存 IKENVH 与 IKCRITH 两处的 ZTAU。
  10. 对 jk=1..nlayer，若 IKCRITH(jl) < jk < IKENVH(jl)，按压力差在两个缓存 stress 之间线性插值，重新组织低层破碎时的 stress profile。

伪代码

GWPROFIL(..., ZTAU):
    for each active GW point jl = kdx(ji):
        zoro(jl) = psig(jl) * zdmod(jl) / 4 / max(pvar(jl), 1)
        cache top stress Z_TAU(jl,nlayer+1)

    for jk from nlayer down to 2:
        for each active point jl:
            if jk >= IKNU2(jl):
                ZTAU(jl,jk) = cached_top_stress

        for each active point jl:
            if jk < IKNU2(jl):
                ZNORM = gkdrag * ZRHO(jl,jk) * sqrt(BV(jl,jk)) * ZVPH(jl,jk) * zoro(jl)
                ZDZ2 = ZTAU(jl,jk+1) / max(ZNORM, gssec)

        for each active point jl:
            if jk < IKNU2(jl):
                if ZTAU(jl,jk+1) < gtsec or jk <= ICRIT(jl):
                    ZTAU(jl,jk) = 0
                else:
                    compute wave Richardson number ZRIW
                    if ZRIW < grcrit:
                        solve limited displacement ZDZ2N
                        ZTAU(jl,jk) = ZNORM * ZDZ2N
                    else:
                        ZTAU(jl,jk) = ZNORM * ZDZ2
                    ZTAU(jl,jk) = min(ZTAU(jl,jk), ZTAU(jl,jk+1))

    for each active point jl:
        cache ZTAU at IKENVH and IKCRITH

    for jk in 1..nlayer:
        for each active point jl:
            if IKCRITH(jl) < jk < IKENVH(jl):
                interpolate ZTAU linearly in pressure between IKENVH and IKCRITH

参与的主题流程

主题 参与方式
地形重力波 / 次网格地形拖曳 calldrag_noro 筛选地形格点，drag_noro 调用 ORODRAG，ORODRAG 调用本例程重算应力廓线，然后用 ZTAU 垂直差分生成风和温度 tendency。
主物理时间步 physiq physiq_mod.F 在 calllott 为真时进入 calldrag_noro 链；本例程是该链内部的 stress profile 环节。
垂直坐标 / 层序转换 调用前 drag_noro 已翻转垂直层序；因此本文中的 jk 与注释里的高低层含义要按 ORODRAG 接口层序理解。

写法特点

复现要点

待确认

相关页面