LMDZ.COMMON Wiki

vlsplt 串行 split 内核和几何/边界处理

输入范围

dyn3d\vlsplt.F: L1-139（vlsplt 主例程）
dyn3d\vlsplt.F: L140-519（vlx，E-W 输送）
dyn3d\vlsplt.F: L520-888（vly，N-S 输送）
dyn3d\vlsplt.F: L889-1060（vlz，垂直输送）
dyn3d\vlsplt.F: L1091-1135（minmaxq 诊断）

Mars 运行参与度：条件经过。vlsplt 是串行 advtrac 中 iadv=10（Van Leer）路径的核心输送例程。本页聚焦三个方向内核的数值格式、几何处理和边界条件。父子 tracer 比例输送和饱和路径详见 vlsplt-serial-family-tracers。

例程定位

vlsplt 的输送核心由三个 RECURSIVE 子例程 vlx、vly、vlz 组成，分别处理 E-W、N-S 和垂直方向。三者共享 Van Leer 斜率限制框架，但各自的几何处理和边界条件差异显著。主例程 vlsplt 通过 x-y-z-y-x 五步 Strang split 将它们组合成二阶时间精度的输送算子。

x-y-z-y-x Strang Split

通量缩放

zzpbar = 0.5 * pdt    ! 水平半时间步
zzw    = pdt           ! 垂直全时间步
mu = pbaru * zzpbar    ! E-W 缩放通量
mv = pbarv * zzpbar    ! N-S 缩放通量
mw = w * zzw           ! 垂直缩放通量
mw(:,llm+1) = 0.       ! 顶部无通量

水平通量使用半时间步 0.5*pdt，垂直通量使用全时间步 pdt。这是 Strang split 的标准做法：水平两步各半步，垂直一步全步。

Split 顺序

vlx(zq, pente_max, zm, mu, iq)   ! E-W 半步
vly(zq, pente_max, zm, mv, iq)   ! N-S 半步
vlz(zq, pente_max, zm, mw, iq)   ! 垂直全步
vly(zq, pente_max, zm, mv, iq)   ! N-S 半步
vlx(zq, pente_max, zm, mu, iq)   ! E-W 半步

这是对称的 Strang split，保证二阶时间精度。x-y-z-y-x 而非 x-y-z-x-y 的选择使得首尾都是 vlx，有利于 E-W 方向的周期边界一致性。

工作数组和回写

主例程将 q(:,:,iq) 和 masse 复制到工作数组 zq 和 zm，split 全部在工作数组上操作，完成后回写：

CALL SCOPY(ijp1llm, q(1,1,iq), 1, zq(1,1,iq), 1)
CALL SCOPY(ijp1llm, masse, 1, zm(1,1,iq), 1)
! ... split 操作 ...
q(:,:,iq) = zq(:,:,iq)
q(ij+iim,l,iq) = q(ij,l,iq)    ! 经向周期包裹

回写时做 E-W 周期包裹：iip1 位置的值从 1 位置复制。

vlx：E-W 输送

签名

RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
USE infotrac, ONLY: nqtot,nqfils,nqdesc,iqfils,qperemin,masseqmin

masse 维度为 (ip1jmp1,llm,nqtot)，每个 tracer 有独立的质量场。u_m 是缩放后的 E-W 质量通量。

Van Leer 斜率

pente_max > -1e-5（标准 Van Leer scheme I）

! 相邻差分
dxqu(ij) = q(ij+1,l,iq) - q(ij,l,iq)

! 最大允许斜率
dxqmax(ij,l) = pente_max * min(|dxqu(ij-1)|, |dxqu(ij)|)

! 中心斜率（相邻同号取平均，异号置零）
IF(dxqu(ij-1)*dxqu(ij) > 0) THEN
   dxq(ij,l) = dxqu(ij-1) + dxqu(ij)
ELSE
   dxq(ij,l) = 0.    ! 局部极值
ENDIF
dxq(ij,l) = 0.5 * dxq(ij,l)

! Van Leer 限制
dxq(ij,l) = sign(min(|dxq(ij,l)|, dxqmax(ij,l)), dxq(ij,l))

pente_max 通常取 2（标准 Van Leer 限制）或 0（退化为上游格式）。CRAY 平台使用 cvmgp 向量函数替代条件分支。

pente_max < -1e-5（乘积斜率）

zz(ij) = dxqu(ij-1) * dxqu(ij)
zz(ij) = zz(ij) + zz(ij)    ! 2*dxqu(ij-1)*dxqu(ij)
IF(zz(ij) > 0) THEN
   dxq(ij,l) = zz(ij) / (dxqu(ij-1) + dxqu(ij))   ! 调和平均形式
ELSE
   dxq(ij,l) = 0.
ENDIF

这是 harmonic-mean 风格的斜率，不使用 pente_max 做额外限制。

周期包裹

E-W 方向的周期边界通过 iip1 间距复制实现：

DO ij = iip1+iip1, ip1jm, iip1
   dxqu(ij) = dxqu(ij-iim)
   dxq(ij-iim,l) = dxq(ij,l)
   u_mq(ij,l) = u_mq(ij-iim,l)
ENDDO

计算域从 iip2 到 ip1jm，iip1 步长的端点通过 iim 间距复制保持周期一致。

通量计算

IF (u_m(ij,l) > 0) THEN
   zdum(ij,l) = 1. - u_m(ij,l) / masse(ij,l,iq)
   u_mq(ij,l) = u_m(ij,l) * (q(ij,l,iq) + 0.5*zdum(ij,l)*dxq(ij,l))
ELSE
   zdum(ij,l) = 1. + u_m(ij,l) / masse(ij+1,l,iq)
   u_mq(ij,l) = u_m(ij,l) * (q(ij+1,l,iq) - 0.5*zdum(ij,l)*dxq(ij+1,l))
ENDIF

zdum 表示上游格点未被通量覆盖的质量比例。通量 = 质量通量 × (上游混合比 + 半斜率修正 × 未覆盖比例)。这是标准的 Van Leer 上游加权格式。

CFL > 1 Semi-Lagrangian 累积

当 zdum < 0 时，表示通量超过了格点质量（CFL > 1），此时标记 iadvplus(ij,l) = 1 并置通量为零，等待特殊处理。

Semi-Lagrangian 累积段逐格点追踪通量跨越的多个格点：

! 正向通量（u_m > 0）
do while (zu_m > masse(ijq,l,iq))
   u_mq(ij,l) = u_mq(ij,l) + q(ijq,l,iq) * masse(ijq,l,iq)
   zu_m = zu_m - masse(ijq,l,iq)
   i = mod(i-2+iim, iim) + 1        ! E-W 周期索引
   ijq = (j-1)*iip1 + i
enddo
! 最后部分格点
u_mq(ij,l) = u_mq(ij,l) + zu_m *
     (q(ijq,l,iq) + 0.5*(1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))

反向通量（u_m < 0）结构对称。索引 i = mod(i-2+iim, iim) + 1 实现 E-W 周期环绕。这段代码不向量化（注释：cette partie est mal vectorisee），只在 n0 > 0 时执行。

通量散度更新

new_m = max(masse(ij,l,iq) + u_m(ij-1,l) - u_m(ij,l), masseqmin)
q(ij,l,iq) = (q(ij,l,iq)*masse(ij,l,iq) + u_mq(ij-1,l) - u_m(ij,l,iq)) / new_m
masse(ij,l,iq) = new_m

守恒形式：新混合比 = (旧 tracer 质量 + 净通量) / 新空气质量。masseqmin 防止除零。

vly：N-S 输送

签名

RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
USE infotrac, ONLY: nqtot,nqfils,nqdesc,iqfils,qperemin,masseqmin
USE comconst_mod, ONLY: pi
include "comgeom.h"    ! aire, rlonu, rlonv, apoln, apols 等

vly 比 vlx 复杂得多，因为需要处理极点几何。

首次初始化（SAVE 语义）

vly 在首次调用时计算并缓存以下几何量：

IF(first) THEN
   do i = 2, iip1
      coslon(i) = cos(rlonv(i))
      sinlon(i) = sin(rlonv(i))
      coslondlon(i) = coslon(i) * (rlonu(i)-rlonu(i-1)) / pi
      sinlondlon(i) = sinlon(i) * (rlonu(i)-rlonu(i-1)) / pi
   ENDDO
   ! i=1 从 iip1 复制（周期包裹）
   airej2  = SSUM(iim, aire(iip2), 1)           ! j=2 纬圈面积和
   airejjm = SSUM(iim, aire(ip1jm-iim), 1)      ! j=jjm-1 纬圈面积和
ENDIF

sinlon/coslon 和带 dlon 权重的版本用于极点斜率滤波。airej2/airejjm 是极点相邻纬圈的参考面积，用于面积加权平均。

极点平均混合比

DO i = 1, iim
   airescb(i) = aire(i+iip1) * q(i+iip1,l,iq)           ! j=2 格点
   airesch(i) = aire(i+ip1jm-iip1) * q(i+ip1jm-iip1,l,iq) ! j=jjm-1 格点
ENDDO
qpns = SSUM(iim, airescb, 1) / airej2    ! 北极平均（面积加权）
qpsn = SSUM(iim, airesch, 1) / airejjm   ! 南极平均（面积加权）

极点混合比是相邻纬圈的面积加权平均，用于计算极点处的斜率。

N-S 斜率

! v 点差分
dyqv(ij) = q(ij,l,iq) - q(ij+iip1,l,iq)

! 标量点斜率（相邻 v 点平均）
dyq(ij,l) = 0.5 * (dyqv(ij-iip1) + dyqv(ij))
dyqmax(ij) = pente_max * min(|dyqv(ij-iip1)|, |dyqv(ij)|)

! 极点斜率（与极点平均的差）
dyq(ij,l) = qpns - q(ij+iip1,l,iq)              ! 北极
dyq(ip1jm+ij,l) = q(ip1jm+ij-iip1,l,iq) - qpsn  ! 南极

极点斜率 sin/cos 滤波

极点处 dyq 需要滤波以保证物理一致性（极点没有确定的经度方向）：

! 投影到 sin/cos 分量
dyn1 = SUM(sinlondlon(ij) * dyq(ij,l))       ! 北极 sin 分量
dyn2 = SUM(coslondlon(ij) * dyq(ij,l))       ! 北极 cos 分量
dys1 = SUM(sinlondlon(ij) * dyq(ip1jm+ij,l)) ! 南极 sin 分量
dys2 = SUM(coslondlon(ij) * dyq(ip1jm+ij,l)) ! 南极 cos 分量

! 重构
dyq(ij,l) = dyn1*sinlon(ij) + dyn2*coslon(ij)              ! 北极
dyq(ip1jm+ij,l) = dys1*sinlon(ij) + dys2*coslon(ij)       ! 南极

这相当于将极点斜率投影到最低两个 Fourier 模态（sin 和 cos），滤掉高阶分量。权重 dlon/pi 确保积分归一化。

GOTO 8888 和极点斜率归零

goto 8888
! ... 被跳过的限制器代码（fn/fs 比例限制）...
8888 continue
DO ij = 1, iip1
   dyq(ij,l) = 0.
   dyq(ip1jm+ij,l) = 0.
ENDDO

GOTO 8888 无条件跳过了 fn/fs 比例限制器段（L662-L676），直接到 continue。之后立即将所有极点 dyq 置零。这意味着 LMDZ 实际运行中极点处不使用斜率，极点的 Van Leer 修正被完全禁用。

被跳过的代码（L663-L676）实现了 fn = min(pente_max*|dyqv|/|dyq|, 1.) 比例限制，曾经用于约束极点斜率幅度，但当前已被弃用。注释块（L683-L747）中还保留了多次历史测试的残迹。

非极点斜率限制

DO ij = iip2, ip1jm
   IF(dyqv(ij)*dyqv(ij-iip1) > 0.) THEN
      dyq(ij,l) = sign(min(|dyq(ij,l)|, dyqmax(ij)), dyq(ij,l))
   ELSE
      dyq(ij,l) = 0.    ! 极值点
   ENDIF
ENDDO

与 vlx 相同的 Van Leer 限制逻辑，但作用在 N-S 方向。

通量计算

IF(masse_adv_v(ij,l) > 0) THEN
   qbyv(ij,l) = q(ij+iip1,l,iq) + dyq(ij+iip1,l) *
                 0.5*(1. - masse_adv_v(ij,l)/masse(ij+iip1,l,iq))
ELSE
   qbyv(ij,l) = q(ij,l,iq) - dyq(ij,l) *
                 0.5*(1. + masse_adv_v(ij,l)/masse(ij,l,iq))
ENDIF
qbyv(ij,l) = masse_adv_v(ij,l) * qbyv(ij,l)

与 vlx 相同的上游加权格式，但索引沿 j 方向（步长 iip1）。

通量散度更新

newmasse = masse(ij,l,iq) + masse_adv_v(ij,l) - masse_adv_v(ij-iip1,l)
q(ij,l,iq) = (q(ij,l,iq)*masse(ij,l,iq) + qbyv(ij,l) - qbyv(ij-iip1,l)) / newmasse
masse(ij,l,iq) = newmasse

N-S 方向的通量散度。注意 masse_adv_v 的索引间距是 iip1（一行经度格点数），而 vlx 中 u_m 的间距是 1。

极点 SSUM 累积和质量分配

! 北极
convpn  = SSUM(iim, qbyv(1,l), 1)           ! 净 tracer 通量
convmpn = SSUM(iim, masse_adv_v(1,l), 1)    ! 净质量通量
massepn = SSUM(iim, masse(1,l,iq), 1)       ! 极点空气质量
qpn = SUM(masse(ij,l,iq)*q(ij,l,iq))        ! 极点 tracer 质量
qpn = (qpn + convpn) / (massepn + convmpn)  ! 新混合比
DO ij = 1, iip1
   q(ij,l,iq) = qpn                          ! 均匀赋值
ENDDO

! 南极（符号相反）
convps  = -SSUM(iim, qbyv(ip1jm-iim,l), 1)
convmps = -SSUM(iim, masse_adv_v(ip1jm-iim,l), 1)
masseps = SSUM(iim, masse(ip1jm+1,l,iq), 1)
qps = (SUM(masse*q) + convps) / (masseps + convmps)
DO ij = ip1jm+1, ip1jmp1
   q(ij,l,iq) = qps
ENDDO

极点处理的关键特征：

vlz：垂直输送

签名

RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)

w 维度为 (ip1jmp1,llm+1)，包含 llm+1 个界面（顶部和底部）。

垂直斜率

! 界面差分（从 l=2 开始）
dzqw(ij,l) = q(ij,l-1,iq) - q(ij,l,iq)

! 标量点斜率
IF(dzqw(ij,l)*dzqw(ij,l+1) > 0.) THEN
   dzq(ij,l) = 0.5*(dzqw(ij,l) + dzqw(ij,l+1))
ELSE
   dzq(ij,l) = 0.
ENDIF
dzqmax = pente_max * min(|dzqw(ij,l)|, |dzqw(ij,l+1)|)
dzq(ij,l) = sign(min(|dzq(ij,l)|, dzqmax), dzq(ij,l))

与 vlx 相同的 Van Leer 限制框架，但作用在垂直方向。

顶底边界斜率归零

DO ij = 1, ip1jmp1
   dzq(ij,1)  = 0.    ! 顶层
   dzq(ij,llm) = 0.    ! 底层
ENDDO

顶层和底层的斜率被强制置零，防止边界处出现非物理的斜率外推。

通量计算

DO l = 1, llm-1
   IF(w(ij,l+1) > 0.) THEN
      sigw = w(ij,l+1) / masse(ij,l+1,iq)
      wq(ij,l+1) = w(ij,l+1) * (q(ij,l+1,iq) + 0.5*(1.-sigw)*dzq(ij,l+1))
   ELSE
      sigw = w(ij,l+1) / masse(ij,l,iq)
      wq(ij,l+1) = w(ij,l+1) * (q(ij,l,iq) - 0.5*(1.+sigw)*dzq(ij,l))
   ENDIF
ENDDO

sigw = w/masse 是 sigma 坐标的垂直 Courant 数。通量公式与 vlx/vly 结构相同。

无通量边界

DO ij = 1, ip1jmp1
   wq(ij,llm+1) = 0.    ! 底部（地面）
   wq(ij,1)     = 0.    ! 顶部（模式顶）
ENDDO

顶底两个界面都强制零通量。这与 vlsplt 主例程中 mw(:,llm+1) = 0. 一致，但 vlz 内部再次确保边界为零。

通量散度更新

newmasse = masse(ij,l,iq) + w(ij,l+1) - w(ij,l)
q(ij,l,iq) = (q(ij,l,iq)*masse(ij,l,iq) + wq(ij,l+1) - wq(ij,l)) / newmasse
masse(ij,l,iq) = newmasse

垂直方向不需要极点处理或周期包裹，结构最简单。

minmaxq 诊断

subroutine minmaxq(zq,qmin,qmax,comment)

受 #ifdef isminismax 编译保护。使用 ismin/ismax 外部函数找全局极值位置，如果超出 [qmin, qmax] 范围则打印位置和值。在 vlsplt 主例程中以注释形式存在（c call minmaxq(...)），开发调试时取消注释可追踪每步 split 后的极值变化。

默认范围：qmin = 0., qmax = 1.e33（在主例程 DATA 语句中定义）。

三方向对比

特征 vlx (E-W) vly (N-S) vlz (垂直)
行数 L140-519（380 行） L520-888（369 行） L889-1060（172 行）
斜率限制 Van Leer scheme I 或乘积 同左 同左
周期/边界 E-W 周期包裹（iim 间距） 极点特殊处理 顶底零通量
极点处理 无 SSUM 累积+面积加权分配 无
极点斜率 不适用 sin/cos 投影后归零（GOTO 8888） 顶底归零
CFL > 1 Semi-Lagrangian 累积（逐格点追踪） 无 无
通量索引间距 1（连续 i） iip1（跨一行） 1（连续 l）
通量散度 标准守恒 标准守恒 + 极点 SSUM 标准守恒
RECURSIVE 是（父子 tracer） 是（父子 tracer） 是（父子 tracer）
几何缓存 无 sinlon/coslon/airej2/airejjm（SAVE） 无

复现要点

相关页面

待确认