LMDZ.COMMON Wiki

ppm3d limiter、正值性检查和极区/边界滤波

输入范围

dyn3d_common\ppm3d.F: L1534-1611（lmtppm）
dyn3d_common\ppm3d.F: L1686-1791（qckxyz）
dyn3d_common\ppm3d.F: L1895-1979（filns）
dyn3d_common\ppm3d.F: L1981-2056（filew）
dyn3d_common\ppm3d.F: L1793-1893（filcr）
dyn3d_common\ppm3d.F: L716-726（主例程调用点）

Mars 运行参与度：条件经过。本页覆盖 ppm3d 内部的单调性约束、正值性检查和负值填充滤波。和 ppm3d-entry-contract 一样，当前 COMMON 并行构建中 PPM 路径不可达。

例程定位

ppm3d 使用五种约束机制保证输送结果的物理合理性：

lmtppm：PPM 单调性约束

subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)
参数 说明
DC 半斜率（0.5 × mismatch），由 xmist/ymist 计算
A6 抛物线曲率 3*(2*P - AL - AR)
AR 右界面值（INOUT，被修改）
AL 左界面值（INOUT，被修改）
P 格点平均值
IM 向量长度
LMT 约束选项（= _ORD - 3）

三档约束

LMT=0（Full monotonicity，_ORD=3）：

完全禁止新极值。对每个格点：

  1. 若 DC=0（局部平坦），直接展平：AR=AL=P, A6=0。
  2. 计算 da1=AR-AL, da2=da1², A6DA=A6*da1。
  3. 若 A6DA < -da2（左端有极值风险）：A6=3*(AL-P), AR=AL-A6。
  4. 若 A6DA > da2（右端有极值风险）：A6=3*(AR-P), AL=AR-A6。

LMT=1（Semi-monotonic，_ORD=4）：

允许过冲但禁止新局部极值。跳过条件：|AR-AL| >= -A6（抛物线已单调）。否则：

  1. 若 P 是局部极小（P<AR 且 P<AL）：展平。
  2. 若 AR>AL：用 A6=3*(AL-P) 修正右端，AR=AL-A6。
  3. 否则：用 A6=3*(AR-P) 修正左端，AL=AR-A6。

LMT=2（Positive-definite，_ORD=5）：

仅防止负值。在 LMT=1 的基础上增加检查：先计算抛物线最小值 fmin = P + 0.25*(AR-AL)²/A6 + A6/12。若 fmin >= 0，跳过（不会出负值）。只有 fmin < 0 时才执行与 LMT=1 相同的修正。

调用点

lmtppm 由 fxppm（L1091）、fyppm（L1336）和 FZPPM（L892-899）调用，每次传入 LMT = _ORD - 3。LMT>2（_ORD>=6）时不调用 lmtppm，即无约束 PPM。

qckxyz：正值性检查和垂直填充

subroutine qckxyz(Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,cross,IC,NSTEP)
参数 说明
Q tracer 密度 DQ（INOUT，被修复）
qtmp 工作数组
cosp/acosp 纬向 cos / 1/cos
cross 是否调用 filcr（LMDZ 传入 .false.）
IC tracer 编号（用于诊断输出）
NSTEP 子步编号（用于诊断输出）

处理流程

qckxyz 对每个垂直层执行三级滤波加垂直填充：

  1. filns（N-S 滤波）：从南北邻居借质量填充负值。返回 ipy：0=无负值，1=有修复。若 ipy=0 跳过该层。
  2. filew（E-W 滤波）：从东西邻居借质量填充负值。返回 ipx。若 ipx=0 跳过该层。
  3. filcr（对角滤波）：从对角邻居借质量填充负值。仅在 cross=.true. 时调用（顶层和内部层）。返回 icr。
  4. 垂直填充：若仍有负值，从上下层借质量。

分层处理

顶层（L=1）：负值直接加到 L=2，自身归零。

Q(i,j1,2) += Q(i,j1,1)    ! 负值下移
Q(i,j1,1)  = 0.

内部层（L=2..NLAY-1）：先从上层借，剩余负值传给下层。

qup = Q(L-1), qly = -Q(L)
dup = min(qly, qup)        ! 最多借完上层全部
Q(L-1) = qup - dup         ! 上层被借走
Q(L)   = dup - qly         ! 本层修复（可能仍有剩余负值）
Q(L+1) += Q(L)             ! 剩余负值传给下层
Q(L)    = 0.

底层（L=NLAY）：从 NLAY-1 借，不够部分记录为"地面质量源"。

qup = Q(NLAYM1), qly = -Q(NLAY)
dup = min(qly, qup)
Q(NLAYM1) = qup - dup
sum += qly - dup           ! 累加地面质量源
Q(NLAY) = 0.

LMDZ 调用特殊性

主例程 L725-726 调用 qckxyz 时 cross 传 .false.：

if(fill) call qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2,
     &               cosp,acosp,.false.,IC,NSTEP)

这意味着顶层和内部层的 filcr 被跳过。但底层（L=NLAY）的 filcr 没有 if(cross) 保护，始终被调用。这可能是底层缺乏下层可借质量，需要更积极的填充。

诊断输出

filns：N-S 负值填充

subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)

从南北邻居借质量消除负值。面积加权使用 cosp（纬向 cos）和 acosp（1/cos）。

对内部格点 q(i,j) < 0：

dq = -q(i,j) * cosp(j)           ! 需要借的总量（面积加权）
dn = q(i,j+1) * cosp(j+1)        ! 北邻居可用量
d1 = min(dq, max(0, dn))          ! 从北借
q(i,j+1) = (dn-d1) * acosp(j+1)  ! 更新北邻居
ds = q(i,j-1) * cosp(j-1)        ! 南邻居可用量
d2 = min(dq-d1, max(0, ds))       ! 从南借
q(i,j-1) = (ds-d2) * acosp(j-1)  ! 更新南邻居
q(i,j) = (d2-dq) * acosp(j) + tiny  ! 修复本层（tiny 防止精确零）

边界处理：

返回 ipy=1 若有任何负值被发现。

filew：E-W 负值填充

subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)

从东西邻居借质量消除负值。先转置 q(i,j) → qtmp(j,i) 以改善向量化。

对内部经度 i=2..IMR-1，qtmp(j,i) < 0：

d1 = min(-qtmp(j,i), max(0, qtmp(j,i-1)))   ! 从西借
qtmp(j,i-1) -= d1
qtmp(j,i)   += d1
d2 = min(-qtmp(j,i), max(0, qtmp(j,i+1)))   ! 从东借
qtmp(j,i+1) -= d2
qtmp(j,i)   += d2 + tiny

经向周期包裹：i=1 从 i=IMR 借，i=IMR 从 i=1 借。

修复完成后转置回 q(i,j)。若无需修复（ipx=0），仍检查极点是否有负值并设置 ipx=1。

filcr：对角负值填充

subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)

从对角邻居借质量消除负值。做两遍扫描：

前向扫描（i=1..IMR-1）：对 q(i,j) < 0，从 NE（i+1,j+1）和 SE（i+1,j-1）借。

后向扫描（i=2..IMR）：对 q(i,j) < 0，从 NW（i-1,j+1）和 SW（i-1,j-1）借。

周期包裹：

借用逻辑与 filns 相同：面积加权 cosp/acosp，按可用量上限 max(0, dn) 借用，剩余量 dq 逐步减少。

最终检查 j1/j2 边界和极点是否有负值，设置 icr=1。

滤波顺序和调用关系

qckxyz 对每层执行三级滤波，形成递进式负值消除：

filns (N-S)     → 从最近邻居借，代价最小
  ↓ if ipy≠0
filew (E-W)     → 从经向邻居借，处理 filns 残留
  ↓ if ipx≠0
filcr (diagonal) → 从对角邻居借，最激进的水平填充
  ↓ if icr≠0
vertical fill    → 从上下层借，最后手段

每级只在上一级未能完全修复时才执行下一级。tiny = 1e-60 加到修复后的格点值上，防止精确零导致后续计算中出现除零。

复现要点

  1. lmtppm 的 LMT = _ORD - 3 每次调用重新计算（不做 SAVE），允许不同 tracer 使用不同 PPM 选项。注释掉的代码（L1058-1079）显示早期版本按分辨率自动选择 LMT。
  2. qckxyz 在 LMDZ 中被传入 cross=.false.，因此 filcr 只在底层被调用。若要在外部复现 qckxyz 行为，需要注意这个参数差异。
  3. filns 的极点处理使用 CAP1 = IMR*(1-cos((j1-1.5)*DP))/DP，这个值每次调用重新计算（注释掉的 SAVE 代码说明早期版本做过缓存）。
  4. filew 的转置（q(i,j) → qtmp(j,i)）是为向量化优化，不影响数值结果。若无向量化需求可以省略转置。
  5. 垂直填充的"地面质量源"是 tracer 不守恒的诊断信号。对 Mars 干大气 tracer，非零 sum 表明负值填充无法从上层完全补偿，可能需要检查上游风场或时间步长。
  6. tiny = 1e-60 的量级远小于任何物理 tracer 混合比，不影响数值精度但保证后续除法安全。

相关页面

待确认