LMDZ.COMMON Wiki

ppm3d 垂直输送、层序翻转和质量通量关系

输入范围

dyn3d_common\ppm3d.F: L420-930（主例程垂直段 + FZPPM）
dyn3d_common\ppm3d.F: L2058-2078（zflip）

Mars 运行参与度：条件经过。本页覆盖 ppm3d 内部垂直输送逻辑、FZPPM 子例程和 zflip 工具子例程。垂直输送是 ppm3d 三段 split（x-y-z）的最后一段。和 ppm3d-entry-contract 一样，当前 COMMON 并行构建中 PPM 路径不可达。

例程定位

ppm3d 的输送分三段：水平 E-W（xtp）、水平 N-S（ytp）、垂直（FZPPM）。垂直段在主例程 L676-756 中先计算垂直质量通量 W，再调用 FZPPM 执行 PPM 垂直输送。zflip 是一个独立的层序翻转工具子例程，ppm3d.F 内部不调用它，LMDZ 的 interpre/interpost 自行实现翻转逻辑。

垂直质量通量计算（主例程 L676-714）

FZPPM 被调用前，主例程先从水平质量散度反推垂直质量通量，过程分三步：

第一步：柱质量散度累积（L680-688）。 对每一列 (i,j) 从顶层开始累加水平质量散度 DPI：

CRY(i,j) = sum(DPI(i,j,k), k=1..NLAY)

CRY 在此处被复用为柱总质量散度（不是经向 Courant 数）。

第二步：地面气压更新（L696）。 利用静力假设，地面气压增量等于柱质量散度：

PS2(i,j) = PS1(i,j) + CRY(i,j)

第三步：垂直质量通量递推（L700-708）。 从质量守恒递推每层界面通量：

W(i,j,1)   = DPI(i,j,1) - DBK(1)*CRY(i,j)
W(i,j,k)   = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)   k=2..NLAY-1
W(i,j,NLAY)= 0

DBK(k) 是混合 sigma-P 坐标 B 系数层差（BP(k+1)-BP(k)），DBK(k)*CRY(i,j) 项反映 sigma 坐标下气压变化对各层质量的影响。地面通量为零，顶层通量由散度驱动。

递推完成后，更新 delp2（时间 n+1 的层厚度）为 DAP(k) + DBK(k)*PS2(i,j)。

FZPPM 签名

subroutine FZPPM(IMR, JNP, NLAY, j1, DQ, WZ, P, DC, DQDT,
     &           AR, AL, A6, flux, wk1, wk2, wz2, delp, KORD)
参数 方向 说明
IMR IN 经向格点数
JNP IN 纬向格点数（含极点）
NLAY IN 垂直层数
j1 IN 极冠参数（LMDZ 用 2）
DQ INOUT tracer 密度（= Q × delp），入口为水平输送后的值，出口为垂直输送后的值
WZ IN 垂直质量通量（即主例程计算的 W）
P IN tracer 混合比 Q
DC OUT PPM 单调斜率（复用为主例程传入的数组）
DQDT WORK 层间差值（复用为工作数组）
AR, AL WORK PPM 右/左界面值
A6 WORK PPM 二次项系数（= 6×(Q̄ - (AL+AR)/2)）
flux WORK 垂直界面通量
wk1, wk2 WORK 工作数组（分别暂存 P 和 delp）
wz2 WORK 暂存 WZ 的纬度切片
delp IN 时间 n 的层厚度
KORD IN 垂直方向阶数选项（≥3 时使用 PPM，<3 时退化为上游）

FZPPM 算法流程

FZPPM 按纬度切片循环处理（L808-928），每个纬度 j 独立执行以下步骤：

1. 单调斜率估计（L787-801）

对内部层 k=2..NLAY-1，用三层加权平均计算 PPM 单调斜率：

c0 = delp(k) / (delp(k-1) + delp(k) + delp(k+1))
c1 = (delp(k-1) + 0.5*delp(k)) / (delp(k+1) + delp(k))
c2 = (delp(k+1) + 0.5*delp(k)) / (delp(k-1) + delp(k))
tmp = c0 * (c1*DQDT(k) + c2*DQDT(k-1))
DC(k) = sign(min(|tmp|, Qmax, Qmin), tmp)

其中 Qmax = max(P(k-1),P(k),P(k+1)) - P(k)，Qmin = P(k) - min(P(k-1),P(k),P(k+1))。这是 Colella & Woodward 1984 的单调约束实现。

2. 界面值初始猜测（L825-878）

大气顶（k=1）：三格 PPM 抛物线拟合，用 DQDT 和 delp 构造二次函数系数 a,b，求顶界面值 AL(1)。如果 wk1(1)*AL(1) ≤ 0（符号改变），则 AL(1)=0，flux(1)=0。

地面（k=NLAY）：两格 PPM 零梯度条件。AR(NLAY) 和 AL(NLAY) 由 DQDT(NLAY-1) 和 delp 构造。如果 wk1(NLAY)*AR(NLAY) ≤ 0，则 AR(NLAY)=0。

内部（k=3..NLAY-1）：四阶插值，使用相邻层的单调斜率和层厚度加权：

AL(k) = wk1(k-1) + c1 + c2*(wk2(k)*(c1*(A1-A2) + A2*flux(k-1)) - wk2(k-1)*A1*flux(k))

然后 AR(k) = AL(k+1)（连续性），A6(k) = 3*(2*wk1(k) - AL(k) - AR(k))（二次项系数）。

3. 单调性约束施加（L892-899）

call lmtppm(flux_top,   A6_top,   AR_top,   AL_top,   wk1_top,   IMR, 0)
call lmtppm(flux_bottom,A6_bottom,AR_bottom,AL_bottom,wk1_bottom,IMR, 0)
if (LMT <= 2) call lmtppm(flux_int, A6_int, AR_int, AL_int, wk1_int, IMR*(NLAY-2), LMT)

顶/底边界始终施加单调约束（LMT=0）。内部层的约束取决于 KORD：KORD=3 时 LMT=0（全单调），KORD=4 时 LMT=1（半单调），KORD=5 时 LMT=2（正定）。

4. 垂直通量计算（L903-918）

对每个界面 k（实际索引 2 对应 k=1 和 k=2 之间），根据通量方向选择上/下界面值：

if wz2(i,k) > 0:     flux = AR(k) + 0.5*CM*(AL(k) - AR(k) + A6(k)*(1 - 2/3*CM))
else:                 flux = AL(k+1) + 0.5*CP*(AL(k+1) - AR(k+1) - A6(k+1)*(1 + 2/3*CP))

其中 CM = wz2/wk2(k) 和 CP = wz2/wk2(k+1) 是 Courant 数。最终通量 = wz2 × flux（界面值 × 质量通量）。

5. tracer 密度更新（L920-927）

DQ(i,j,1)   = DQ(i,j,1)   - flux(i,2)         ! 顶层：减去下界面流出
DQ(i,j,NLAY)= DQ(i,j,NLAY)+ flux(i,NLAY)       ! 底层：加上上界面流入
DQ(i,j,k)   = DQ(i,j,k)   + flux(i,k) - flux(i,k+1)   ! 内部：净通量

6. 混合比恢复（主例程 L731-738）

FZPPM 返回后，主例程用更新后的密度除以新层厚度恢复混合比：

Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k)

zflip 工具子例程

subroutine zflip(q, im, km, nc)

zflip 把三维数组 q(im, km, nc) 在垂直方向做翻转：q(i,k,IC) ← q(i, km+1-k, IC)。使用临时数组 qtmp(im, km) 暂存翻转后的值，然后写回。

ppm3d.F 内部不调用 zflip。LMDZ 的 interpre/interpost 自行实现垂直层序翻转（field_ppm(l) = field(llm-l+1)），因为 LMDZ 的翻转还涉及经向风取反和地面气压计算，不能简单调用 zflip。

zflip 是 PPM Transport Core 提供给外部使用者的通用工具。当外部代码的垂直层序与 PPM 内部约定（k=1 为大气顶）不一致时使用。

垂直输送在 split 顺序中的位置

ppm3d 主例程的输送顺序（L447-756）：

  1. 水平段（L447-674，每个垂直层 k 循环）：
    • 计算水平 Courant 数 CRX/CRY
    • 计算水平质量通量 xmass/ymass 和散度 DPI
    • 对每个 tracer IC：计算交叉项 → xtp（E-W）→ ytp（N-S）
  2. 垂直质量通量计算（L676-714）：从水平散度递推 W
  3. 垂直段（L716-748，每个 tracer IC 循环）：
    • FZPPM 垂直 PPM 输送
    • qckxyz 正值性检查（如果 fill=.true.）
    • 混合比恢复 Q = DQ / delp2

水平段按层循环（外层 k），垂直段按 tracer 循环（外层 IC），这是因为垂直输送需要所有层的完整柱信息。

复现要点

  1. FZPPM 的 DQ 入口值是水平输送后的 tracer 密度（= Q × delp1），不是原始混合比。复现结果必须从水平段输出开始。
  2. 垂直质量通量 W 由水平散度递推得出，不依赖外部输入。这意味着垂直和水平输送是耦合的，不能独立测试。
  3. KORD=max(3, KORD) 在主例程 L716 强制执行，即使用户传入 KORD<3，FZPPM 实际使用的值至少为 3。
  4. qckxyz 正值性检查在 FZPPM 之后执行，如果 fill=.true.（LMDZ 调用值为 .true.），会对所有负值做填充修正。
  5. LMDZ 不使用 zflip，自行在 interpre/interpost 中实现翻转。如果要用外部代码调用 ppm3d 并依赖 zflip，需注意 zflip 只做纯垂直翻转，不处理经向风取反和气压积分。

相关页面

待确认