orosetup.F90

快速理解

它做什么： Lott/Miller 地形重力波拖曳链的预处理层，被 ORODRAG 调用。

基本过程： 把翻转后的压力/温度/风/位势/地形统计量转换成后续 GWSTRESS/GWPROFIL 所需工作量。

关键结果： 地形高度层索引、低层风、Brunt-Väisälä 频率、Richardson 数、投影风速和阻塞层参数。

路径

LMDZ.MARS\libf\phymars\orosetup.F90

所属目录 / 模块

libf\phymars

文件定位

orosetup.F90 定义 orosetup_mod 和 OROSETUP，是 Lott/Miller 地形重力波拖曳链中的预处理层。它不直接输出风温 tendency，而是在 orodrag_mod.F90 的 ORODRAG 内被调用，负责把已经由 drag_noro_mod 翻转到地形拖曳层序的压力、温度、风场、位势高度和次网格地形统计量转换成后续 GWSTRESS/GWPROFIL 与阻塞层 wake drag 所需的工作量。

核心输出包括地形高度层索引 IKNU/IKNU2/IKENVH/IKCRITH/ICRIT、低层平均风 PULOW/PVLOW、半层密度 ZRHO、Brunt-Vaisala 频率平方 BV、Richardson 数 PRI、投影风速 ZVPH、地形方向角 ZPSI、阻塞层 leakiness ZZDEP 以及各向异性组合量 ZD1/ZD2/ZDMOD。下游 ORODRAG 随后用这些量初始化应力、重算应力廓线，并把应力差分或阻塞层拖曳换算成 pdudt/pdvdt/pdtdt。

定义的符号

依赖的模块

调用的关键例程

OROSETUP 在 IKNUL=IKNUB 的保护分支中有一次 WRITE(*,*) 诊断输出，提示低层平均区间退化并用相邻层的 BV/ZRHO 兜底。

输入

输出

共享状态与副作用

• 本文件没有 SAVE 变量、THREADPRIVATE、文件 I/O 或配置读取。
• pgam 是 intent(inout)；行 144 将其改为 max(pgam,gtsec)，因此调用方传入的局部地形各向异性数组会被本例程下限保护。
• ZTAU 是输出工作数组，但本例程只在全层循环中初始化 ZTAU(:,1:nlayer)=0.0；底部半层 ZTAU(:,nlayer+1) 由后续 GWSTRESS 初始化。
• 当 IKNUL=IKNUB 时，行 362 会向标准输出写诊断信息，并把 BV(:,nlayer+1)、ZRHO(:,nlayer+1) 设为相邻全层值。这是唯一运行时输出副作用。
• 只有 ktest(JL)==1 的分支会计算密度、稳定度、投影风、Richardson 数和阻塞层积分；但部分初始化循环仍覆盖 kidia:kfdia 全部局地点。

核心逻辑

1. 初始化局部范围和常数：kidia=1、kfdia=ngrid，ILEVH=nlayer/3，并构造 1/r、g**2/cpp 和 1.5*pi。
2. 初始化地形高度索引。IKNU/IKNU2/IKNUb/IKNUl 先设为 nlayer，pgam 用 gtsec 下限保护，然后从近地层向上搜索 zgeom/g 何时越过 4*pvar、3*pvar、2*pvar 和 1*pvar。
3. 用 nktopg 修正 IKNU/IKNUb，并在 IKNUb 到达顶部限制或 IKNUl <= IKNUb 时调整 IKNUl，避免 1-2 倍 pvar 平均区间退化。
4. 初始化半层边界和输出数组：BV/PRI/ZPSI/ZVPH/PULOW/PVLOW/IKCRITH/IKENVH/ICRIT 等被设置为表面或顶部默认值。
5. 从近地层向上计算全层之间的半层密度 ZRHO 和 BV。密度使用 rho=p/(r*T) 的两层平均形式，BV 使用源码中的压力差表达式并用 gssec 下限保护。
6. 在 IKNUb..IKNUl 之间做压力加权平均，得到低层风 PULOW/PVLOW，并用 gvsec 给低层风速范数 ZNORM 下限。
7. 根据低层风向和地形主轴角得到底部半层 ZPSI(:,nlayer+1)，再由 pgam 计算 B=1-0.18*pgam-0.04*pgam^2、C=0.48*pgam+0.3*pgam^2，得到 ZD1/ZD2/ZDMOD。
8. 对每层把实际风投影到低层应力方向，形成 ZVPF；随后把 ZVPF 插值到半层得到 ZVPH。若半层投影风速低于 gvsec，则强制为 gvsec 并记录 ICRIT=JK。
9. 在 1-2 倍 pvar 区间累加底部半层 BV(:,nlayer+1) 和 ZRHO(:,nlayer+1)，再按压力厚度平均；退化区间走 WRITE 诊断和相邻层兜底。
10. 用 BV、ZRHO、沿应力方向风切变和压力差计算 PRI，并用 grcrit 下限保护。
11. 从 IKNU2 以下积分方程 9 的 ZNU，当它跨过 gfrcrit 且 IKENVH 仍为表面默认值时，记录阻塞层顶部。
12. 在阻塞层以上继续积分 ZNUP，当它跨过 1.5 且 IKCRITH 仍为默认值时，记录重力波破碎动力混合高度；最后 IKCRITH=min(IKCRITH,IKNU)。
13. 对 JK >= IKENVH 的阻塞层重新计算层内 ZPSI，并用 zgeom(IKENVH)、当前层 zgeom 和 pvar*g 计算 ZZDEP。

伪代码

OROSETUP(pressure, wind, temperature, zgeom, subgrid-orography):
    set constants and vertical search limits

    for each local column:
        initialize IKNU*, IKENVH, IKCRITH, ICRIT and work arrays
        pgam = max(pgam, gtsec)

    for thresholds 4*pvar, 3*pvar, 2*pvar, 1*pvar:
        scan from surface upward
        compare zgeom/g against threshold
        record the layer where the logical test changes

    constrain terrain-height indices with nktopg

    for active columns and half-levels:
        compute density and BV from pressure, temperature and constants
        apply gssec lower bound to BV

    average pu/pv between 1*pvar and 2*pvar levels:
        PULOW/PVLOW = pressure-weighted low-level wind
        ZVPH(:,nlayer+1) = max(low-level wind speed, gvsec)

    compute terrain orientation terms:
        ZPSI = terrain-axis angle - low-level flow angle
        ZD1/ZD2/ZDMOD = anisotropic mountain coefficients from pgam and ZPSI

    project each layer wind into the low-level stress plane:
        ZVPF -> half-level ZVPH
        if ZVPH < gvsec:
            set ZVPH = gvsec
            record ICRIT

    average bottom half-level BV/ZRHO between 1*pvar and 2*pvar:
        if interval is degenerate:
            print diagnostic and copy adjacent layer values

    compute PRI with a grcrit lower bound

    integrate ZNU below 3*pvar:
        first crossing of gfrcrit defines IKENVH

    integrate ZNUP above blocked layer:
        first crossing of 1.5 defines IKCRITH
        IKCRITH = min(IKCRITH, IKNU)

    for blocked-layer levels:
        recompute ZPSI from local layer wind
        compute ZZDEP leakiness from envelope height and current height

参与的主题流程

写法特点

• .F90 自由格式模块，但保留 ECMWF/IFS 和 Lott/Miller 方案的历史注释。
• 本例程大量使用 nlayer 侧作为表面、从 nlayer 向小层号扫描的层序约定；这依赖上游 drag_noro 的翻转。
• BV 注释说是 Brunt-Vaisala frequency，但表达式和下限变量按 N^2 使用；页面按源码计算量记录为频率平方。
• kentp、ncount、zulow、zvlow 等局部变量只初始化或声明，未参与后续计算。
• 注释中保留多处开发者疑问，例如 ILEVH=nlayer/3 “maybe not enough for Mars”、BV 公式待替换、PRI 中 dp 是否应为 dp^2。

复现要点

• 必须先按 drag_noro 的层序翻转输入；直接用原始物理层序调用会改变所有 JK >= IKENVH、IKNU* 和 pressure-thickness 逻辑。
• pgam 会被原地改写为不小于 gtsec，复现实验若比较原始地形统计量，需要在调用前后区分。
• IKNUb..IKNUl 的压力平均区间决定 PULOW/PVLOW、底部 BV/ZRHO 和后续应力尺度；nktopg、pvar 和 zgeom/g 的一致性会直接影响拖曳强度。
• ZVPH 和风切变用 gvsec 下限保护，BV 用 gssec 下限保护，PRI 用 grcrit 下限保护；这些 yoegwd_h 参数必须与运行配置一致。
• 源码只对底部 BV/ZRHO 平均的 IKNUL=IKNUB 情况显式兜底；低层风平均和其他分母仍依赖前置索引修正与输入范围。
• ZZDEP 只在 JK >= IKENVH 的阻塞层写入，其他层保持 0；下游 wake drag 分支不能把 0 当作已计算 leakiness。

待确认

• 源码注释称 ILEVH=nlayer/3 “maybe not enough for Mars”，本页只记录字面行为，未判断该上限对不同垂直分辨率是否充分。
• BV 公式旁注释说 “Use N^2=g/T[1/(cpp*T)+dT/dz] to replace in the future”；是否应替换属于方案演进问题，非本页结论。
• PRI 公式旁注释说 “Here dp maybe dp^2 ? Need ask Francois lott later”；本页按源码记录，未确认理论公式。
• zgeom 在上游页面中按位势相关量记录；本页按源码 zgeom/g 和 pvar*g 用法说明，不重新定义其单位来源。

相关页面

• orodrag_mod：直接调用本例程，并消费 IK*、BV/ZRHO/ZVPH/ZPSI/ZZDEP/ZD* 输出。
• drag_noro_mod：上游包装层，翻转垂直层序并构造 zgeom。
• calldrag_noro_mod：读取地形 GW 开关并筛选进入拖曳计算的格点。
• gwstress_mod：使用本例程输出的低层风、稳定度、密度和地形方向量初始化底部应力。
• gwprofil_mod：使用本例程输出的临界层、阻塞层和稳定度信息重算垂直应力廓线。
• phymars 目录：本模块所属目录。
• sugwd：记录 gfrcrit/grcrit/gssec/gvsec/nktopg 等 GWD 调参常量初始化来源。
• yoegwd_h：记录这些共享变量的定义位置和 THREADPRIVATE 属性。