orodrag_mod.F90

路径

LMDZ.MARS\libf\phymars\orodrag_mod.F90

所属目录 / 模块

libf\phymars

文件定位

orodrag_mod.F90 定义 orodrag_mod 和核心例程 ORODRAG，是 Lott/Miller 地形重力波与次网格地形阻塞拖曳链的核心计算层。上游 calldrag_noro_mod 负责读取开关、按 zstd 筛选格点并分子域调用；drag_noro_mod 负责翻转垂直层序和构造 zgeom；本文件则调用 OROSETUP、GWSTRESS、GWPROFIL 生成半层重力波应力 ZTAU，再把应力垂直差分和阻塞层拖曳转换为 pdudt/pdvdt/pdtdt。

文件名中 orodrag 对应 orographic gravity wave drag。它和 nonoro_gwd_ran_mod / nonoro_gwd_mix_mod 的非地形 GW 路径不同；physiq_mod.F 中由 calllott 开关进入 calldrag_noro 链后才会走到本例程。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 本文件没有本地 SAVE 变量、THREADPRIVATE、文件 I/O、配置读取或诊断输出。
• pgam 是 intent(inout)；OROSETUP 中会执行 pgam=max(pgam,gtsec)，因此调用链会修改传入的地形各向异性局部数组。
• ZTAU、BV、ZRHO、ZVPH、PRI 等为本例程局部工作数组，由 OROSETUP/GWSTRESS/GWPROFIL 依次填充或改写。
• 只在 JL=1..kgwd 循环中写 pdudt/pdvdt/pdtdt；调用方 drag_noro 在调用前已把这些输出工作数组清零，因此未被 kgwd/kdx 覆盖的格点保持零 tendency。

核心逻辑

1. 调用 OROSETUP。该步骤按翻转后的压力、风、温度、zgeom 和地形统计量计算：
   • 低层风 PULOW/PVLOW；
   • 阻塞层顶 IKENVH、波破碎混合高度 IKCRITH、临界层 ICRIT；
   • 半层密度 ZRHO、Brunt-Vaisala 频率 BV、低层风速投影 ZVPH；
   • 地形方向组合量 ZD1/ZD2/ZDMOD、压力厚度 ZZDEP 和角度 ZPSI。
2. 调用 GWSTRESS。该步骤在底层半层初始化 ZTAU(:,nlayer+1)。
3. 调用 GWPROFIL。该步骤沿半层重算 ZTAU 垂直廓线。
4. 将 ZDUDT/ZDVDT/ZDTDT 初始化为 0。
5. 对每个层 JK=1..nlayer 和每个有效 GW 点 JL=1..kgwd：
   • 用 JI=kdx(JL) 找到当前子域局部格点；
   • 计算半层压力厚度 ZDELP=pplev(JI,JK+1)-pplev(JI,JK)；
   • 由 stress 差分计算公共因子 ZTEMP=-g*(ZTAU(JI,JK+1)-ZTAU(JI,JK))/(ZVPH(JI,nlayer+1)*ZDELP)；
   • 用 PULOW/PVLOW 和 ZD1/ZD2/ZDMOD 把公共因子投影为 ZDUDT/ZDVDT。
6. 如果 JK >= IKENVH(JI)，说明当前翻转层序下位于阻塞层内，则改用阻塞层 wake drag 公式：
   • 计算 zb=1-0.18*pgam-0.04*pgam^2 与 zc=0.48*pgam+0.3*pgam^2；
   • 计算 zconb=2*ptimestep*gkwake*psig/(4*pvar)；
   • 由 ZPSI 得到 zzd1、障碍物长宽比 ratio，再得 zbet=max(0,2-1/ratio)*zconb*ZZDEP*zzd1*zabsv；
   • 令 ZDUDT=-pu/ptimestep*zbet/(1+zbet)，ZDVDT=-pv/ptimestep*zbet/(1+zbet)。
7. 写出 pdudt/pdvdt。
8. 用 pu/pv + ptimestep*tendency 得到拖曳后的风 ZUST/ZVST，计算动能差 ZDIS，再以 ZDEDT=ZDIS/ptimestep、ZDTDT=ZDEDT/cpp 写出 pdtdt。

伪代码

ORODRAG(inputs from drag_noro):
    OROSETUP(...) -> layer indices, low-level wind, BV, density, ZPSI, ZZDEP, ZD1/ZD2/ZDMOD
    GWSTRESS(...) -> initialize bottom half-level ZTAU
    GWPROFIL(...) -> recompute vertical ZTAU profile

    zero local tendency work arrays

    for each model layer JK:
        for each active GW point JI = kdx(JL):
            dp = pplev(JI,JK+1) - pplev(JI,JK)
            stress_factor = -g * (ZTAU(JI,JK+1)-ZTAU(JI,JK))
                            / (ZVPH(JI,nlayer+1) * dp)
            du = projection of stress_factor using PULOW, PVLOW, ZD1, ZD2, ZDMOD
            dv = projection of stress_factor using PULOW, PVLOW, ZD1, ZD2, ZDMOD

            if JK is inside blocked layer:
                compute B, C, obstacle aspect ratio and zbet
                du = -pu / ptimestep * zbet/(1+zbet)
                dv = -pv / ptimestep * zbet/(1+zbet)

            pdudt(JI,JK) = du
            pdvdt(JI,JK) = dv

            wind_after = wind_before + ptimestep * tendency
            kinetic_energy_loss = 0.5 * (|wind_before|^2 - |wind_after|^2)
            pdtdt(JI,JK) = kinetic_energy_loss / ptimestep / cpp

参与的主题流程

写法特点

• .F90 自由格式模块，但保留 ECMWF/IFS 历史注释和大量被注释掉的旧变量。
• kgwd/kdx 控制有效点循环；本例程不遍历所有 ngrid 点计算 tendency。
• ILEVP1=nlayer+1 只用于读取 ZVPH(:,nlayer+1)。
• ZVIDIS 被初始化和累加，但源码注释标明它是无效/未使用变量。
• 温度 tendency 来自拖曳前后水平风动能差，并除以 cpp；不是由热力方程单独求解。

复现要点

• 必须保留 OROSETUP -> GWSTRESS -> GWPROFIL -> tendency 的调用顺序；ZTAU 是顺序改写的 inout 工作数组。
• 上游 drag_noro 已翻转垂直层序。若直接用原始 physiq 层序调用本例程，会改变阻塞层和 stress 差分含义。
• pgam 会在 OROSETUP 内被下限保护；如果复用同一局部数组，需要考虑该 intent(inout) 副作用。
• 阻塞层分支会覆盖由 stress 差分得到的 ZDUDT/ZDVDT，而不是叠加。
• pvar 进入 zconb 分母，ZVPH(:,nlayer+1) 和 ZDELP 也在 stress tendency 分母中；源码没有额外零保护，依赖上游筛选和参数范围。
• pdtdt 对应动能损失换热，可能随阻塞层分支或 stress 分支的风 tendency 改变；复现时不要只复制动量 tendency。

待确认

• zgeom 在上游 drag_noro 页面中已标注为位势相关量；本文件按 OROSETUP 源码使用 zgeom/g 判定高度，未重新解释其单位。
• if(jk.ge.ikenvh) 的注释写“level lower than blocked layer”，该判断依赖翻转层序；需要在 orosetup.F90 文件页中继续核对所有 IK* 索引约定。
• ZTEMP 分母没有对 ZVPH(JI,nlayer+1) 或 ZDELP 做保护；本页只记录源码行为，不推断实际运行中是否可能为零。
• ZVIDIS 当前未被输出或读取，可能是历史诊断残留。

相关页面

• calldrag_noro_mod：physiq 调用层，读取开关和 zstd 阈值，分子域调用 drag_noro。
• drag_noro_mod：上游包装层，翻转层序、构造 zgeom 并把本例程 tendency 转为每步增量。
• gwstress_mod：本例程内部调用的底层应力初始化步骤。
• gwprofil_mod：本例程内部调用的 stress profile 重算步骤。
• phymars 目录：本模块所属目录。
• orosetup：本例程直接调用的预处理层，提供 IK* 层索引、BV/ZRHO/ZVPH/ZPSI/ZD* 和阻塞层 ZZDEP。
• sugwd：记录 gkwake/gkdrag/grcrit 等地形 GW 调参常量初始化来源。
• yoegwd_h：记录这些共享变量的定义位置和 THREADPRIVATE 属性。