ppm3d 水平输送和经纬向通量计算

输入范围

dyn3d_common\ppm3d.F: L584-672（主例程水平段 + cross term 调用）
dyn3d_common\ppm3d.F: L932-1047（xtp）
dyn3d_common\ppm3d.F: L1049-1107（fxppm）
dyn3d_common\ppm3d.F: L1109-1123（xmist）
dyn3d_common\ppm3d.F: L1125-1189（ytp）
dyn3d_common\ppm3d.F: L1191-1271（ymist）
dyn3d_common\ppm3d.F: L1273-1350（fyppm）
dyn3d_common\ppm3d.F: L1352-1440（yadv）
dyn3d_common\ppm3d.F: L1442-1532（xadv）

Mars 运行参与度：条件经过。本页覆盖 ppm3d 内部水平输送的八个子例程和主例程调用模式。和 ppm3d-entry-contract 一样，当前 COMMON 并行构建中 PPM 路径不可达。

例程定位

ppm3d 水平输送采用 x-y directional splitting。每个垂直层 k 独立处理：先算交叉项（cross terms），再做 E-W 输送（xtp），最后做 N-S 输送（ytp）。E-W 和 N-S 各有一个 PPM 通量计算子例程（fxppm/fyppm）和一个斜率初始化子例程（xmist/ymist）。交叉项由 xadv/yadv 处理。

主例程调用模式

ppm3d 主例程在 L584-672 对每个 tracer IC 和每个垂直层 k 执行以下序列：

交叉项初始化（L586-614）：内联计算 E-W 和 N-S 方向的 1 阶上游交叉项存入 wk1(:,:,:,1) 和 wk1(:,:,:,2)。E-W 部分区分高 CFL 区（j<=JS 或 j>=JN，用 Semi-Lagrangian 周期包裹）和 Eulerian 区。
交叉项修正（L645-665，如果 cross=.true.）：
- xadv(wk1(:,:,2), UA, ..., DC2, iad)：对 N-S advected 场做 E-W 交叉修正
- yadv(wk1(:,:,1), VA, ..., PV, jad)：对 E-W advected 场做 N-S 交叉修正
- q(i,j,k,IC) += DC2(i,j) + PV(i,j)：叠加交叉修正到 tracer
E-W 输送（L667-668）：xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,wk1(:,:,2),CRX,fx1,xmass,IORD)
N-S 输送（L670-671）：ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,wk1(:,:,1),CRY,DC2,ymass,WK1(:,:,3),WK1(:,:,4),WK1(:,:,5),WK1(:,:,6),JORD)

iad 和 jad 根据 IORD/JORD 决定：>=2 用 2 阶，否则用 1 阶。

cross 在 L300 被 data cross /.true./ 初始化，LMDZ 调用时始终为 .true.。

xtp：E-W 输送

subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC,fx1,xmass,IORD)

参数	说明
`IMR`	经向格点数
`JNP`	纬向格点数（含极点）
`IML`	周期包裹半宽（= max integer CFL）
`j1,j2`	活动纬度范围
`JN,JS`	Eulerian 区边界（由 CFL>1 确定）
`PU`	E-W 界面密度（= 0.5*(delp2(i)+delp2(i-1))）
`DQ`	tracer 密度（INOUT，被更新）
`Q`	混合比（经交叉修正后的值）
`UC`	Courant 数（= CRX）
`fx1`	工作数组，最终存放质量加权通量
`xmass`	E-W 质量通量（= PU * CRX）
`IORD`	E-W 方案阶数

两条路径

xtp 根据纬度位置选择 Eulerian 或 Semi-Lagrangian（FFSL）路径：

Eulerian 路径（JS < j < JN，CFL<1 区域）：

IORD	公式	说明
`1` 或极区边界	`fx1(i) = qtmp(iu)`	1 阶上游，`iu = i - uc(i,j)`
`2` 或高纬度	`fx1(i) = qtmp(iu) + DC(iu)*(sign(1,uc)-uc)`	2 阶 van Leer，DC 由 `xmist` 计算
`>=3`	`call fxppm(...)`	PPM 抛物线重构

最终通量：fx1(i) = fx1(i) * xmass(i,j)。

Semi-Lagrangian 路径（j<=JS 或 j>=JN，CFL>1 区域）：

两条路径都先填充周期包裹的 ghost cells（qtmp(-IML:0) 和 qtmp(IMP:IMP+IML)），确保大 CFL 数时上游取值不越界。

IORD=1 使用分数部分 rut = uc - int(uc) 和整数偏移 ISAVE = i - int(uc) 做分段累加。IORD>=2 使用 van Leer 斜率加整数格累积。最终通量：fx1(i) = PU(i,j) * fx1(i)。

DQ 更新

fx1(IMP) = fx1(1)                     ! 经向周期闭合
DQ(i,j) += fx1(i) - fx1(i+1)          ! 通量散度

高纬度 van Leer 切换

xtp 在纬度 j <= j1vl 或 j >= j2vl 时强制使用 van Leer（即使 IORD>=3），其中 jvan = max(1, JNP/18)，j1vl = j1+jvan，j2vl = j2-jvan。这是为了避免极区附近 PPM 抛物线重构不稳定。

fxppm：E-W PPM 通量

subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD)

PPM 抛物线重构和界面通量计算，分四步：

左界面值：AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1)-DC(i))/3
右界面值：AR(i) = AL(i+1)，周期包裹 AR(IMR) = AL(1)
曲率：A6(i) = 3*(2*p(i) - AL(i) - AR(i))
单调约束：如果 LMT = IORD-3 <= 2，调用 lmtppm 限制抛物线

界面通量（对每个格点界面 i）：

UT(i) > 0:  flux(i) = AR(i-1) + 0.5*UT*(AL(i-1)-AR(i-1) + A6(i-1)*(1-2/3*UT))
UT(i) <= 0: flux(i) = AL(i)   - 0.5*UT*(AR(i)-AL(i)     + A6(i)*(1+2/3*UT))

这是标准 PPM 积分公式：对抛物线在 [0, UT] 或 [UT, 0] 区间积分得到通过界面的质量。

周期包裹：AL(0)=AL(IMR), AR(0)=AR(IMR), A6(0)=A6(IMR)。

xmist：E-W 斜率初始化

subroutine xmist(IMR,IML,P,DC)

4 阶中心差分加单调约束：

tmp  = (8*(p(i+1)-p(i-1)) + p(i-2)-p(i+2)) / 24
Pmax = max(p(i-1), p(i), p(i+1)) - p(i)
Pmin = p(i) - min(p(i-1), p(i), p(i+1))
DC(i)= sign(min(|tmp|, Pmax, Pmin), tmp)

DC 是半斜率（0.5 × mismatch），保证重构值不超出邻域极值。

ytp：N-S 输送

subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL,JORD)

参数	说明
`acosp`	纬向格面积倒数（1/cos(lat) 归一化）
`RCAP`	极冠面积倒数
`VC`	N-S Courant 数（= CRY）
`ymass`	N-S 质量通量
`fx`	工作数组，存放 N-S 通量
`A6,AR,AL`	PPM 工作数组（传入 fyppm）
`JORD`	N-S 方案阶数

三条路径

JORD	公式	说明
`1`	`fx(i,j1) = p(i,JT)`	1 阶上游，`JT = j1 - VC(i,j1)`
`2`	`fx = p(i,JT) + (sign(1,VC)-VC)*DC2(i,JT)`	2 阶 van Leer，DC2 由 `ymist` 计算
`<=0` 或 `>=3`	`call fyppm(...)`	PPM 抛物线重构

最终通量：fx(i,j1) = fx(i,j1) * ymass(i,j1)。

DQ 更新

DQ(i,j) += (fx(i,j) - fx(i,j+1)) * acosp(j)    ! 面积加权通量散度

极区处理

ytp 在极区做特殊质量守恒处理：

sum1 = Σ fx(i,j1)      ! 南极边界总通量
sum2 = Σ fx(i,J2+1)    ! 北极边界总通量
DQ(i,1)   = DQ(1,1)   - sum1 * RCAP    ! 南极：均匀分配
DQ(i,JNP) = DQ(1,JNP) + sum2 * RCAP    ! 北极：均匀分配

如果 j1 ≠ 2（扩大极冠），还会把极冠第二行（j=2 和 j=JMR）也设为极区平均值。

fyppm：N-S PPM 通量

subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)

N-S 方向的 PPM 重构，结构与 fxppm 类似但有极区特殊处理：

左界面：AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1)-DC(i,2))/3
右界面：AR(i,1) = AL(i,2)
极区交叉（L1317-1323）：跨极点半宽 IMH = IMR/2 的交叉赋值
- AL(i,1) = AL(i+IMH,2)，AL(i+IMH,1) = AL(i,2)
- AR(i,JNP) = AR(i+IMH,JMR)，AR(i+IMH,JNP) = AR(i,JMR)
LMDZ.3.3 补丁（L1328-1329）：AR(IMR,1)=AL(1,1)，AR(IMR,JNP)=AL(1,JNP)，修复经向周期格点的极区界面一致性。
曲率：A6(i,j11) = 3*(2*p(i,j11) - AL(i,j11) - AR(i,j11))
单调约束：LMT = JORD-3，若 <=2 则调用 lmtppm
通量计算：同 fxppm 公式，使用 VC 符号选择上游方向

ymist：N-S 斜率初始化

subroutine ymist(IMR,JNP,j1,P,DC,ID)

ytp 以 ID=4 调用 ymist，执行 4 阶斜率计算：

j=2（南边界）：tmp = (8*(p(i,3)-p(i,1)) + p(i+IMH,2)-p(i,4))/24，使用跨经度半宽 IMH 的邻居。前 IMH 个格点用 +IMH 偏移，后 IMH 个用 -IMH 偏移。
j=JMR（北边界）：类似，使用 p(i,JNP) 和跨经度邻居。
内部（j=3..JMR-1）：标准 4 阶 (8*(p(i,j+1)-p(i,j-1)) + p(i,j-2)-p(i,j+2))/24。

极冠斜率（j1=2）：在南极和北极用跨经度差分 tmp = 0.25*(p(i,2)-p(i+IMH,2))，后半经度取反 DC(i+IMH) = -DC(i)。

如果 j1 ≠ 2（扩大极冠），极点斜率直接设为零。

所有斜率都做单调约束 DC = sign(min(|tmp|, Pmax, Pmin), tmp)。

xadv：E-W 交叉项

subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)

计算 E-W 方向的 advective cross term。输入 p 是已经做过 N-S 平流的场（wk1(:,:,2)），输出 adx 是 E-W 交叉修正量。

只在高 CFL 区域（j<=JS 或 j>=JN）做 Semi-Lagrangian 交叉修正。Eulerian 区域（JS+1 到 JN-1）做简单 upwind：

IAD=2: adx(i,j) = qtmp(ip)-p(i,j) + ru*(a1*ru+b1)    ! 2阶
IAD=1: adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1))       ! 1阶

极点处交叉项设为零：adx(i,1) = 0, adx(i,JNP) = 0。如果 j1 ≠ 2，极冠第二行也设为零。

yadv：N-S 交叉项

subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)

计算 N-S 方向的 advective cross term。输入 p 是已经做过 E-W 平流的场（wk1(:,:,1)），输出 ady 是 N-S 交叉修正量。

yadv 准备跨极区 ghost cells（wk(i,-1:0) 和 wk(i,JNP+1:JNP+2)），使用半宽 IMH 跨经度取值以匹配极点物理位置。

IAD=2: ady(i,j) = wk(i,jp) + rv*(a1*rv+b1) - wk(i,j)   ! 2阶
IAD=1: ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))         ! 1阶

极点处做纬向平均：sum1 = Σ ady(i,1) / IMR，然后所有经度格点赋相同值，保证极区质量守恒。

通量与质量守恒

水平输送的质量守恒由以下机制保证：

E-W 周期闭合：fx1(IMP) = fx1(1)，通量散度 DQ += fx1(i)-fx1(i+1) 在求和时自动消去内部项。
N-S 极区分配：ytp 把所有经度格点的极区边界通量求和，乘以 RCAP（极冠面积倒数）后均匀分配给极点所有经度格。这保证进入/离开极区的总质量精确守恒。
面积加权：N-S 通量乘以 acosp(j)（纬向格面积倒数），补偿球面格点面积随纬度的变化。
Semi-Lagrangian 整数格累积：当 CFL>1 时，xtp 的 FFSL 路径显式累加穿越的整格质量（fx1(i) += Σ qtmp(ist)），保证大 Courant 数下的质量精确。

复现要点

xtp 的 Eulerian/FFSL 切换由 JS/JN 控制，这两个值在主例程 L470-488 通过扫描 CRX 找到首个 |CRX|>1 的纬度确定。复现时必须先正确计算 JS/JN。
fxppm 和 fyppm 的 LMT 参数每次调用重新计算（LMT = _ORD - 3），不做 SAVE。注释掉的代码（L1058-1079）显示早期版本按分辨率自动选择 LMT，现已废弃。
ymist 的 ID=4 路径使用跨经度半宽（IMH=IMR/2）邻居，这是因为极点附近经度方向物理距离缩短，4 阶差分需要跨越更大经度范围取样。
fyppm 的 LMDZ.3.3 补丁（L1328-1329）修复了 AR(IMR) 在极区的周期一致性。如果没有这个补丁，最后一个经度格点的 N-S PPM 界面值会和第一个格点不匹配。
cross terms（xadv/yadv）的结果叠加到 q 后，xtp/ytp 使用的是已修正的 wk1(:,:,1) 和 wk1(:,:,2)，不是原始混合比。这意味着 cross term 对水平输送有直接影响。
极点处 xadv 的交叉项被强制设为零（L1527-1530），这是因为极点无经向方向定义，E-W 交叉修正无物理意义。

待确认

cross 标志对 Mars 干大气 tracer 的实际影响。cross=.true. 启用交叉项修正，但交叉项的大小和符号依赖具体风场和 tracer 分布。
ymist 的 ID=2 路径（2 阶简单斜率）在当前 LMDZ 调用中不被使用（ytp 固定传 ID=4），可能在其他 PPM 调用场景中存在。

ppm3d 水平输送和经纬向通量计算

输入范围

例程定位

主例程调用模式

xtp：E-W 输送

两条路径

DQ 更新

高纬度 van Leer 切换

fxppm：E-W PPM 通量

xmist：E-W 斜率初始化

ytp：N-S 输送

三条路径

DQ 更新

极区处理

fyppm：N-S PPM 通量

ymist：N-S 斜率初始化

xadv：E-W 交叉项

yadv：N-S 交叉项

通量与质量守恒

复现要点

相关页面

待确认