LMDZ.COMMON Wiki

vlsplt_p 并行 halo 注册、边界带/内区计算和 band 依赖

输入范围

dyn3dpar\vlsplt_p.F: L1-201（vlsplt_p 主例程）
dyn3dpar\vlsplt_p.F: L204-588（vlx_p）
dyn3dpar\vlsplt_p.F: L591-1036（vly_p）
dyn3dpar\vlsplt_p.F: L1040-1260（vlz_p）
dyn3dpar\vlsplt_p.F: L1293-1335（minmaxq_p）
dyn3dpar\mod_hallo.F90（Register_Hallo/SendRequest/WaitRecvRequest/WaitSendRequest）
dyn3dpar\parallel_lmdz.F90（ij_begin/ij_end/pole_nord/pole_sud）
dyn3dpar\bands.F90（band 分布和 ij_begin/ij_end 来源）

Mars 运行参与度：条件经过。Mars 并行 advtrac_p 中 iadv=10 路径经过 vlsplt_p 或经 vlspltgen_p 间接调用 vlx_p/vly_p/vlz_p。

例程定位

vlsplt_p 是 vlsplt 的并行版本。它将本地 band 切分为"边界带"和"内区"，先算边界带、发起 halo 通信、再算内区，用计算-通信重叠掩盖 MPI 延迟。三个方向内核 vlx_p/vly_p/vlz_p 都接收 band 范围参数，并使用 OpenMP DO SCHEDULE(STATIC,OMP_CHUNK) 并行化。

vlsplt_p 主例程：halo 重叠调度

本地 band 范围

USE parallel_lmdz
ijb = ij_begin
ije = ij_end
if (pole_nord) ijb = ijb + iip1    ! 北极进程跳过极点行
if (pole_sud)  ije = ije - iip1    ! 南极进程跳过极点行

ij_begin/ij_end 来自 parallel_lmdz，由 bands.F90 的 SetDistrib 按负载均衡设定。极点进程缩进一行以避免在极点行做 E-W 输送。

mw 扩展维度

与串行 vlsplt 的 mw(ip1jmp1,llm+1) 不同，并行版的垂直通量数组为 mw(ip1jmp1,llm+1,nqtot)，为每个 tracer 保留独立的垂直通量。

五步 split 与 halo 重叠

并行版不是一次性算完整个方向再交换 halo，而是把每个 E-W 和垂直方向调用拆成"边界带 + 内区"两段，中间插入异步 halo 交换：

! 第一步 vlx_p：先算边界带
vlx_p(zq, zm, mu, ij_begin, ij_begin+2*iip1-1, iq)   ! 左边界 2 行
vlx_p(zq, zm, mu, ij_end-2*iip1+1, ij_end, iq)       ! 右边界 2 行

! 注册 halo 并异步发送
Register_Hallo(zq, ip1jmp1, llm, 2, 2, 2, 2, MyRequest1)
Register_Hallo(zm, ip1jmp1, llm, 1, 1, 1, 1, MyRequest1)
SendRequest(MyRequest1)

! 计算内区（与通信重叠）
vlx_p(zq, zm, mu, ij_begin+2*iip1, ij_end-2*iip1, iq)

! 等待 halo 完成
WaitRecvRequest(MyRequest1)

第二步 vly_p 不需要 halo 交换，因为 N-S 方向的数据依赖完全在 band 内部。

第三步 vlz_p 使用与第一步相同的"边界带 + halo + 内区"模式，使用 MyRequest2：

vlz_p(..., ij_begin, ij_begin+2*iip1-1, iq)
vlz_p(..., ij_end-2*iip1+1, ij_end, iq)
Register_Hallo(zq, ..., MyRequest2)
Register_Hallo(zm, ..., MyRequest2)
SendRequest(MyRequest2)
vlz_p(..., ij_begin+2*iip1, ij_end-2*iip1, iq)
WaitRecvRequest(MyRequest2)

第四步 vly_p 再次无 halo。

第五步（最终 vlx_p）直接用 ij_begin, ij_end 全范围调用，因为此时 halo 已到位：

vlx_p(zq, zm, mu, ij_begin, ij_end, iq)

最后回写 q 并调用 WaitSendRequest 确保发送完成。

halo 注册参数

Register_Hallo(field, dim1, dim2, halo_s, halo_n, halo_e, halo_w, request) 中的 halo 宽度参数：

字段 halo 宽度 含义
zq 2, 2, 2, 2 tracer 场，四周各 2 行 halo
zm 1, 1, 1, 1 质量场，四周各 1 行 halo

zq 需要 2 行 halo 是因为 Van Leer 斜率计算需要左右各两个格点的信息；zm 只需 1 行因为质量更新只涉及相邻格点。

为什么边界带宽是 2*iip1

E-W 方向每行有 iip1 个格点（含周期端点），2 行 halo 对应 2*iip1 个连续元素。先算边界 2 行意味着计算覆盖 [ij_begin, ij_begin+2*iip1-1] 和 [ij_end-2*iip1+1, ij_end]，这些格点的 Van Leer 模板需要的邻居数据在本 rank 内已可用。通信完成后，内区 [ij_begin+2*iip1, ij_end-2*iip1] 可以使用从邻居收到的 halo 数据。

vlx_p：E-W 方向内核

签名差异

RECURSIVE SUBROUTINE vlx_p(q, pente_max, masse, u_m, ijb_x, ije_x, iq)
USE Parallel_lmdz

与串行 vlx 相比，增加 ijb_x/ije_x 参数限定计算范围，使主例程可以分段调用。

极点行排除

ijb = ijb_x
ije = ije_x
if (pole_nord .and. ijb==1) ijb = ijb + iip1
if (pole_sud .and. ije==ip1jmp1) ije = ije - iip1

北极进程（拥有全局第 1 行）和南极进程（拥有全局最后一行）跳过极点行，因为极点没有确定的经度方向。

OMP 并行化

所有主要循环使用 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 和 c$OMP END DO NOWAIT。OMP_CHUNK 来自 parallel_lmdz，控制静态调度的块大小。NOWAIT 允许循环间重叠。

数值计算

与串行 vlx 完全相同的 Van Leer 斜率限制、上游通量和 Semi-Lagrangian CFL > 1 累积逻辑。周期包裹使用相同的 mod(i-2+iim,iim)+1 公式。

父子 tracer

父子 tracer 的 masseq/Ratio 计算和递归调用与串行版相同，但所有循环限定在 [ijb, ije] 范围内，递归调用传递原始的 ijb_x/ije_x。

vly_p：N-S 方向内核

签名

RECURSIVE SUBROUTINE vly_p(q, pente_max, masse, masse_adv_v, iq)
USE parallel_lmdz
USE comconst_mod, ONLY: pi

vly_p 不接收 ijb_x/ije_x 参数，而是直接使用 parallel_lmdz 的 ij_begin/ij_end 和 pole_nord/pole_sud。

OMP THREADPRIVATE

vly_p 有大量 SAVE 变量需要 THREADPRIVATE：

SAVE temps0,...,temps5, first, testcpu
c$OMP THREADPRIVATE(temps0,...,temps5, first, testcpu)
SAVE sinlon, coslon, sinlondlon, coslondlon
c$OMP THREADPRIVATE(sinlon, coslon, sinlondlon, coslondlon)
SAVE airej2, airejjm
c$OMP THREADPRIVATE(airej2, airejjm)

这些几何量和计时变量每个线程需要独立副本。串行版 vly 只需 SAVE first, testcpu。

扩展 band 范围

vly_p 在计算斜率时需要扩展 band 范围以获取邻居数据：

! v 点斜率：扩展 2*iip1
ijb = ij_begin - 2*iip1
ije = ij_end + iip1

! 标量点斜率：扩展 iip1
ijb = ij_begin - iip1
ije = ij_end + iip1

这些扩展区域的数据来自之前 vlx_p 阶段的 halo 交换。

极点条件处理

与串行 vly 无条件处理两个极点不同，并行版按进程位置条件处理：

if (pole_nord) then
   ! 仅北极进程计算北极 qpns 和 dyq
endif
if (pole_sud) then
   ! 仅南极进程计算南极 qpsn 和 dyq
endif

sin/cos 投影和归零逻辑与串行版相同，但只在拥有对应极点的进程上执行。注释 c ym tout cela ne sert pas a grand chose 表示极点 dyq 最终仍被归零，投影计算是残留代码。

通量和散度更新范围

通量计算使用 ij_begin-iip1 到 ij_end（考虑 v 点偏移），散度更新使用 ij_begin 到 ij_end，均受极点条件调整。极点 SSUM 累积也条件化：

if (pole_nord) then
   convpn = SSUM(iim, qbyv(1,l), 1)
   ! ... 北极质量更新 ...
endif
if (pole_sud) then
   convps = -SSUM(iim, qbyv(ip1jm-iim,l), 1)
   ! ... 南极质量更新 ...
endif

vlz_p：垂直方向内核

签名

RECURSIVE SUBROUTINE vlz_p(q, pente_max, masse, w, ijb_x, ije_x, iq)
USE Parallel_lmdz

与 vlx_p 一样接收 ijb_x/ije_x 参数。

动态分配的工作数组

REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: wq
REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: Ratio
LOGICAL, SAVE :: first = .TRUE.
!$OMP THREADPRIVATE(first)

IF (first) THEN
   first = .FALSE.
   !$OMP MASTER
   ALLOCATE(wq(ip1jmp1, llm+1, nqtot))
   ALLOCATE(Ratio(ip1jmp1, llm, nqtot))
   !$OMP END MASTER
   !$OMP BARRIER
END IF

与串行版的栈上固定数组不同，并行版使用 ALLOCATABLE, SAVE 数组，在首次调用时由 OMP MASTER 线程分配。OMP BARRIER 确保所有线程等待分配完成。注释解释这样做是因为这些变量需要被所有线程可见。

其他工作数组

REAL, SAVE :: dzq(ip1jmp1,llm), dzqw(ip1jmp1,llm), adzqw(ip1jmp1,llm)

这些数组是共享的 SAVE 工作区，不是隐式 THREADPRIVATE。并行正确性依赖每次 vlz_p 调用只在当前 [ijb, ije] band 上读写，并由 OMP DO 把同一数组的不同 ij/l 片段分给线程；dzq 的垂直边界由 OMP MASTER 写入后用 BARRIER 同步，内部层由前面的 OMP 循环覆盖后再读。

OMP MASTER + BARRIER 模式

垂直边界条件由主线程设置：

!$OMP MASTER
DO ij=ijb,ije
   dzq(ij,1) = 0.
   dzq(ij,llm) = 0.
ENDDO
!$OMP END MASTER
!$OMP BARRIER

零通量边界 wq(ij,llm+1,iq)=0 和 wq(ij,1,iq)=0 也用相同模式。

父子 tracer 中的 BARRIER

! 父子 Ratio/masseq 计算
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
DO l=1,llm
   DO ij=ijb,ije
      masseq(ij,l,iq2) = max(masse(ij,l,iq)*q(ij,l,iq), masseqmin)
      ...
      w(ij,l,iq2) = wq(ij,l,iq)   ! 注意：写入 w 而非 wq
   enddo
enddo
!$OMP END DO NOWAIT
...
!$OMP BARRIER

这里有一个值得注意的修正：w(ij,l,iq2) = wq(ij,l,iq)，将父 tracer 的垂直通量 wq 赋给子 tracer 的 w 输入。注释 correction bug le 15mai2015 标记了这个修正。递归调用 vlz_p 使用 w（而非 wq）作为通量输入，因此需要将父代的 wq 转写到子代的 w。

递归调用后的通量散度更新前有额外 !$OMP BARRIER，确保所有线程的 wq 已就绪。

与串行 vlsplt 的关键差异

特征 串行 vlsplt 并行 vlsplt_p
计算范围 全域 iip2 到 ip1jm band 范围 ij_begin 到 ij_end
halo 交换 无 vlx_p/vlz_p 前后 Register_Hallo + SendRequest
通信重叠 无 边界带先算 → halo 发送 → 内区计算 → 等待
极点处理 无条件 pole_nord/pole_sud 条件分支
OMP 无 DO SCHEDULE(STATIC,OMP_CHUNK)
THREADPRIVATE 无 vly_p 几何量和计时变量、vlz_p 的 first/testcpu/temps*
工作数组分配 栈上固定数组 vlz_p 的 wq/Ratio 动态分配
E-W split 调用 单次全域 三次分段（边界+内区）
极点 dyq GOTO 8888 跳过 注释掉 GOTO + 条件处理
垂直通量维度 mw(llm+1) mw(llm+1,nqtot) 每 tracer

复现要点

相关页面

待确认