drag_noro_mod.F90

路径

LMDZ.MARS\libf\phymars\drag_noro_mod.F90

所属目录 / 模块

libf\phymars

文件定位

drag_noro_mod.F90 定义 drag_noro_mod 模块和 drag_noro 子程序，是 calldrag_noro_mod.F90 与 orodrag_mod.F90::ORODRAG 之间的地形重力波包装层。它不读取配置，也不决定哪些全局格点启用方案；上游 calldrag_noro 已按 ndomainsz 切分子域，并用 zstd > oro_gwd_zstdthread 填好 kgwd/kdx/ktest。

本例程的职责是把当前子域的温度、风和压力垂直顺序翻转为 ORODRAG 需要的排列，按静力关系构造 zgeom，调用 ORODRAG 计算地形重力波 tendency，再把 tendency 翻回原始层序并乘以 ptimestep 形成每个物理步长的增量。文件名中的 noro 不应理解为非地形重力波；源码注释写明它是 sub-grid scale orographic scheme，下游调用的是 ORODRAG。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 本文件没有 save 变量、THREADPRIVATE 声明、文件 I/O、日志输出或配置读取。
• pgam 在接口上是 intent(inout)，本文件没有直接写它，但会把它传给 ORODRAG 及其下游链路。
• pulow/pvlow/pustr/pvstr/d_t/d_u/d_v 和局部 tendency 数组在进入主计算前清零；若 kgwd=0 或 ktest 全部关闭，下游返回值仍从零态开始。
• pustr/pvstr 在 ORODRAG 返回后由本文件逐层累加，使用翻转后的半层压力厚度 inv_pplev(i,k+1)-inv_pplev(i,k)。

核心逻辑

1. 对当前子域 i=1..ngrid 初始化低层风、stress、增量数组和 ORODRAG tendency 工作数组为 0。
2. 将温度、风、全层压力和半层压力按垂直方向翻转：
   • inv_pt(i,k)=t(i,nlayer-k+1)
   • inv_pu/inv_pv/inv_pplay 同理；
   • inv_pplev(i,k)=pplev(i,nlayer+1-k+1)，最后令 inv_pplev(i,nlayer+1)=pplev(i,1)。
3. 用 r*T*log(p ratio) 构造 zgeom。源码注释称其为 g*dz 和 geopotential height；由于公式未除以 g，这里更接近位势差量而不是几何高度米。
4. 调用 ORODRAG(ngrid,nlayer,kgwd,kgwdim,kdx,ktest,ptimestep,inv_pplev,inv_pplay,zgeom,inv_pt,inv_pu,inv_pv,pvar,psig,pgam,pthe,...)。
5. 将 pdudt/pdvdt/pdtdt 翻回原始层序，并乘以 ptimestep 得到 d_u/d_v/d_t；同时把风 tendency 乘以 g 和半层压力厚度累加到 pustr/pvstr。

伪代码

drag_noro(current_subdomain):
    zero pulow, pvlow, pustr, pvstr
    zero d_t, d_u, d_v and pdudt, pdvdt, pdtdt

    for k in 1..nlayer:
        source_k = nlayer - k + 1
        inv_pt(:,k) = t(:,source_k)
        inv_pu(:,k) = u(:,source_k)
        inv_pv(:,k) = v(:,source_k)
        inv_pplay(:,k) = pplay(:,source_k)
        inv_pplev(:,k) = pplev(:,nlayer + 2 - k)
    inv_pplev(:,nlayer+1) = pplev(:,1)

    zgeom(:,nlayer) = r * inv_pt(:,nlayer) * log(inv_pplev(:,nlayer+1) / inv_pplay(:,nlayer))
    for k from nlayer-1 down to 1:
        zgeom(:,k) = zgeom(:,k+1)
                   + r * average(inv_pt(:,k), inv_pt(:,k+1))
                     * log(inv_pplay(:,k+1) / inv_pplay(:,k))

    ORODRAG(... inverted fields, topography statistics, kgwd/kdx/ktest ...)

    for k in 1..nlayer:
        original_k = nlayer + 1 - k
        d_u(:,original_k) = ptimestep * pdudt(:,k)
        d_v(:,original_k) = ptimestep * pdvdt(:,k)
        d_t(:,original_k) = ptimestep * pdtdt(:,k)
        pustr(:) += g * pdudt(:,k) * (inv_pplev(:,k+1)-inv_pplev(:,k))
        pvstr(:) += g * pdvdt(:,k) * (inv_pplev(:,k+1)-inv_pplev(:,k))

参与的主题流程

写法特点

• .F90 自由格式模块，局部自动数组统一按 ndomainsz 和 nlayer 声明。
• 垂直层序翻转是硬性接口约定；ORODRAG 的输入数组、返回 tendency 和 stress 积分都使用翻转后的层序。
• 源码将 t 注释为 m/s 是注释错误；接口语义和下游公式显示它是全层温度，单位 K。
• ptimestep 在本文件中参与 tendency 到 increment 的转换；上游 calldrag_noro 随后再除回 tendency。

复现要点

• 复现 drag_noro 时必须保留垂直翻转和翻回映射；直接把原始层序送入 ORODRAG 会改变地形 GW 廓线计算。
• zgeom 的公式是 r*T*log(p ratio) 累加；不要在未检查 ORODRAG 期望单位前把它改写成几何高度米。
• d_t/d_u/d_v 是每步增量，不是 tendency；调用方 calldrag_noro 会除以 ptimestep。
• pvar/psig/pgam/pthe 必须由上游按子域起点预切片，并与 kdx/ktest 的局部编号对齐。
• ngrid 是当前子域有效点数，不一定等于 ndomainsz；最后一个子域尤其需要按 ngrid 限制循环。

待确认

• pthe 的角度单位需从地形预处理或 surfdat_h 填充路径继续确认；本文件只透传。
• zgeom 在注释中称为 geopotential height，但公式没有除以 g；本页按源码公式记录为位势相关量，未推断下游物理单位。
• pgam 被声明为 intent(inout)，但是否会被 ORODRAG/OROSETUP 修改需在 orodrag_mod.F90 和 orosetup.F90 专页中继续核验。

相关页面

• calldrag_noro_mod：上游调度层，读取 oro_gwd_zstdthread 并分子域调用本例程。
• phymars/index：phymars 目录总览。
• orodrag_mod：地形 GW 参数化核心，调用 OROSETUP/GWSTRESS/GWPROFIL 并返回风温 tendency。
• gwstress_mod：ORODRAG 内底层应力初始化步骤。
• gwprofil_mod：ORODRAG 内 stress 廓线重算步骤。
• sugwd：地形 GW 调参常量初始化。
• yoegwd_h：地形 GW 共享状态定义。