LMDZ.COMMON Wiki

advn 通量更新：advnx/advny/advnz 的三方向质量守恒更新

输入范围

dyn3d_common\advn.F: L470-733（advnx）
dyn3d_common\advn.F: L735-863（advny）
dyn3d_common\advn.F: L865-977（advnz）

Mars 运行参与度：条件经过。本页聚焦 advn 的三个通量更新子例程的 sigma 参数、mode 分支、极区处理和守恒更新机制。这些子例程接收重构子例程（advnqx/advnqy/advnqz）的界面值，计算 tracer 通量并更新质量场和混合比场。调度顺序和调用关系详见 advn-dispatcher-order。

例程定位

advnx/advny/advnz 共享相同的核心逻辑：基于界面重构值计算 sigma 参数，根据上游方向选择格点参数，使用分段线性或抛物线公式计算 tracer 通量，最后通过通量散度同时更新质量和 tracer 混合比。关键差异在于 advnx 有 mode=1/2 分支和 Semi-Lagrangian 极区处理，advny 有极区面积加权分配，advnz 有无通量边界。

advnx：E-W 通量更新

SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
参数 类型 方向 说明
q real(ip1jmp1,llm) INOUT tracer 混合比（被更新）
qg real(ip1jmp1,llm) IN 左界面重构值（来自 advnqx）
qd real(ip1jmp1,llm) IN 右界面重构值（来自 advnqx）
masse real(ip1jmp1,llm) INOUT 格点空气质量（被更新）
u_m real(ip1jmp1,llm) IN E-W 质量通量（已缩放 0.5*pdt）
mode integer IN 1=分段线性，2=抛物线修正

Sigma 参数

zdq = qd(ij,l) - qg(ij,l)
if (abs(zdq) > prec) then
   zsigd(ij,l) = (q(ij,l) - qg(ij,l)) / zdq   ! 格点值在界面间的相对位置
   zsigg(ij,l) = 1 - zsigd(ij,l)
else
   zsigd(ij,l) = 0.5
   zsigg(ij,l) = 0.5
   qd(ij,l) = q(ij,l)
   qg(ij,l) = q(ij,l)
endif

prec 精度阈值：CRAY 用 1e-24，其他平台用 1e-15。低于此阈值时退化为 1 阶（界面值等于格点值）。

上游方向选择

根据风向选择上游格点的 sigma 和界面值：

风向 上游格点 zsigp zsigm zqp zqm zq
u_m >= 0 ij zsigd(ij) zsigg(ij) qd(ij) qg(ij) q(ij)
u_m < 0 ij+1 zsigg(ij+1) zsigd(ij+1) qg(ij+1) qd(ij+1) q(ij+1)

zsigp 是从上游方向到格点值的相对距离，zsigm 是从格点值到下游方向的相对距离。zqp 是上游界面值，zqm 是下游界面值。

CFL 参数

zu = abs(u_m(ij,l))
ladvplus(ij,l) = zu > zm        ! CFL > 1 标记
zsig = zu / zm                   ! CFL 数
if (zsig == 0) zsigp = 0.1       ! 避免除零

mode=1 分段线性通量

if (zsig <= zsigp) then
   u_mq = u_m * zqp                              ! 未穿越上游界面
else
   u_mq = sign(zm,u_m) * (zsigp*zqp + (zsig-zsigp)*zqm)  ! 穿越上游界面
endif

mode=2 抛物线通量

if (zsig <= zsigp) then
   u_mq = u_m * (zqp - 0.5*zsig/zsigp*(zqp - zq))         ! 未穿越上游界面
else
   zz = 0.5*(zsig - zsigp)/zsigm
   u_mq = sign(zm,u_m) * (0.5*(zq+zqp)*zsigp               ! 穿越上游界面
      + (zsig-zsigp)*(zq + zz*(zqm-zq)))
endif

Semi-Lagrangian 极区处理（CFL>1）

当 ladvplus(ij,l)=.true.（即 |u_m| > masse）时，常规通量公式可能不守恒。advnx 对这些格点使用 Semi-Lagrangian 累加：

正风向（u_m > 0）：

ijq = ij
do while (zu_m > masse(ijq,l))
   u_mq(ij,l) = u_mq(ij,l) + q(ijq,l)*masse(ijq,l)  ! 累加完全穿越格点
   zu_m = zu_m - masse(ijq,l)
   i = mod(i-2+iim,iim) + 1                           ! 经向周期包裹
   ijq = (j-1)*iip1 + i
enddo
! 最后部分穿越格点使用 mode=2 插值
zsig = zu_m / masse(ijq,l)
if (zsig <= zsigd(ijq,l)) then
   u_mq = u_mq + zu_m*(qd(ijq) - 0.5*zsig/zsigd*(qd(ijq)-q(ijq)))
else
   zz = 0.5*(zsig-zsigd)/zsigg
   u_mq = u_mq + masse(ijq)*(0.5*(q(ijq)+qd(ijq))*zsigd
      + (zsig-zsigd)*(q(ijq)+zz*(qg(ijq)-q(ijq))))
endif

负风向（u_m < 0）：类似逻辑，但向 ij+1 方向累加，使用 zsigg 和 qg。

进入条件：n0 > 1（至少有 2 个 CFL>1 格点）才执行 Semi-Lagrangian 累加。若 prt_level > 9，打印穿越格点总数。

排除极点行：mod(ij,iip1) /= 0 条件确保极点行不进入 Semi-Lagrangian 处理。

经向周期包裹

通量计算后，对每行最后一个格点复制首格点的通量和标记：

do l=1,llm
   do ij = iip1+iip1, ip1jm, iip1
      u_mq(ij,l) = u_mq(ij-iim,l)
      ladvplus(ij,l) = ladvplus(ij-iim,l)
   enddo
enddo

通量散度更新

do l=1,llm
   do ij = iip2+1, ip1jm
      new_m = masse(ij,l) + u_m(ij-1,l) - u_m(ij,l)
      q(ij,l) = (q(ij,l)*masse(ij,l) + u_mq(ij-1,l) - u_mq(ij,l)) / new_m
      masse(ij,l) = new_m
   enddo
   do ij = iip1+iip1, ip1jm, iip1
      q(ij-iim,l) = q(ij,l)
      masse(ij-iim,l) = masse(ij,l)
   enddo
enddo

质量和 tracer 同时被更新，保证下一步 split 使用一致的场。

advny：N-S 通量更新

SUBROUTINE advny(q,qs,qn,masse,v_m)
参数 类型 方向 说明
q real(ip1jmp1,llm) INOUT tracer 混合比（被更新）
qs real(ip1jmp1,llm) IN 南界面重构值（来自 advnqy）
qn real(ip1jmp1,llm) IN 北界面重构值（来自 advnqy）
masse real(ip1jmp1,llm) INOUT 格点空气质量（被更新）
v_m real(ip1jm,llm) IN N-S 质量通量（已缩放 0.5*pdt）

Sigma 参数

zdq = qn(ij,l) - qs(ij,l)
if (abs(zdq) > prec) then
   zsign(ij) = (q(ij,l) - qs(ij,l)) / zdq
   zsigs(ij) = 1 - zsign(ij)
else
   zsign(ij) = 0.5
   zsigs(ij) = 0.5
endif

prec 使用 1e-15（与 advnx 相同）。zsign/zsigs 是 1D 数组（按 ij 索引），因为每层独立处理。

上游方向选择

风向 上游格点 zsigp zsigm zqp zqm zq
v_m >= 0 ij+iip1 zsign(ij+iip1) zsigs(ij+iip1) qn(ij+iip1) qs(ij+iip1) q(ij+iip1)
v_m < 0 ij zsigs(ij) zsign(ij) qs(ij) qn(ij) q(ij)

通量公式（仅抛物线）

advny 没有 mode 分支，始终使用抛物线公式：

zsig = abs(v_m(ij,l)) / zm
if (zsig == 0) zsigp = 0.1
if (zsig <= zsigp) then
   v_mq = v_m * (zqp - 0.5*zsig/zsigp*(zqp - zq))
else
   zz = 0.5*(zsig - zsigp)/zsigm
   v_mq = sign(zm,v_m) * (0.5*(zq+zqp)*zsigp
      + (zsig-zsigp)*(zq + zz*(zqm-zq)))
endif

通量散度更新

do ij = iip2, ip1jm
   new_m = masse(ij,l) + v_m(ij,l) - v_m(ij-iip1,l)
   q(ij,l) = (q(ij,l)*masse(ij,l) + v_mq(ij,l) - v_mq(ij-iip1,l)) / new_m
   masse(ij,l) = new_m
enddo

极区面积加权分配

advny 在更新非极区格点后，对极点做特殊处理：

北极（j=1 行）：

convpn  = SSUM(iim, v_mq(1,l), 1)      ! 进入北极的总 tracer 通量
convmpn = SSUM(iim, v_m(1,l), 1)        ! 进入北极的总质量通量
massen  = SSUM(iim, masse(1,l), 1)       ! 北极当前总质量
new_m   = massen + convmpn
q(1,l)  = (q(1,l)*massen + convpn) / new_m
do ij = 1, iip1
   q(ij,l) = q(1,l)                      ! 极点各经度格 tracer 均匀
   masse(ij,l) = new_m * aire(ij)/apoln  ! 按面积比例分配质量
enddo

南极（j=jjp1 行）：

convps  = -SSUM(iim, v_mq(ip1jm-iim,l), 1)
convmps = -SSUM(iim, v_m(ip1jm-iim,l), 1)
masses  = SSUM(iim, masse(ip1jm+1,l), 1)
new_m   = masses + convmps
q(ip1jm+1,l) = (q(ip1jm+1,l)*masses + convps) / new_m
do ij = ip1jm+1, ip1jmp1
   q(ij,l) = q(ip1jm+1,l)
   masse(ij,l) = new_m * aire(ij)/apols
enddo

aire(ij)/apoln 和 aire(ij)/apols 是格点面积与极冠总面积之比，保证极点各经度格的质量按面积比例分配。SSUM 是对 iim 个经度格求和。

advnz：垂直通量更新

SUBROUTINE advnz(q,qh,qb,masse,w_m)
参数 类型 方向 说明
q real(ip1jmp1,llm) INOUT tracer 混合比（被更新）
qh real(ip1jmp1,llm) IN 上界面重构值（来自 advnqz）
qb real(ip1jmp1,llm) IN 下界面重构值（来自 advnqz）
masse real(ip1jmp1,llm) INOUT 格点空气质量（被更新）
w_m real(ip1jmp1,llm+1) IN 垂直质量通量（已缩放 pdt）

Sigma 参数

zdq = qb(ij,l) - qh(ij,l)
if (abs(zdq) > prec) then
   zsigb(ij,l) = (q(ij,l) - qh(ij,l)) / zdq
   zsigh(ij,l) = 1 - zsigb(ij,l)
   zsigb(ij,l) = min(max(zsigb(ij,l), 0.), 1.)  ! 额外的 [0,1] clip
else
   zsigb(ij,l) = 0.5
   zsigh(ij,l) = 0.5
endif

prec 使用 1e-13（比 advnx/advny 的 1e-15 宽松）。zsigb 额外做 [0,1] clip，防止数值误差导致 sigma 越界。

上游方向选择

风向 上游格点 zsigp zsigm zqp zqm zq
w_m >= 0 l（当前层） zsigb(ij,l) zsigh(ij,l) qb(ij,l) qh(ij,l) q(ij,l)
w_m < 0 l-1（下层） zsigh(ij,l-1) zsigb(ij,l-1) qh(ij,l-1) qb(ij,l-1) q(ij,l-1)

通量公式（仅抛物线）

advnz 没有 mode 分支，始终使用抛物线公式：

zsig = abs(w_m(ij,l)) / zm
if (zsig == 0) zsigp = 0.1
if (zsig <= zsigp) then
   w_mq = w_m * (zqp - 0.5*zsig/zsigp*(zqp - zq))
else
   zz = 0.5*(zsig - zsigp)/zsigm
   w_mq = sign(zm,w_m) * (0.5*(zq+zqp)*zsigp
      + (zsig-zsigp)*(zq + zz*(zqm-zq)))
endif

通量计算范围 l=2..llm，因为 l=1 和 l=llm+1 是边界。

无通量边界

do ij = 1, ip1jmp1
   w_mq(ij,llm+1) = 0.   ! 大气顶无通量
   w_mq(ij,1) = 0.        ! 地面无通量
enddo

通量散度更新

do l = 1, llm
   do ij = 1, ip1jmp1
      new_m = masse(ij,l) + w_m(ij,l+1) - w_m(ij,l)
      q(ij,l) = (q(ij,l)*masse(ij,l) + w_mq(ij,l+1) - w_mq(ij,l)) / new_m
      masse(ij,l) = new_m
   enddo
enddo

垂直方向没有周期包裹或极点处理，所有格点统一更新。

三方向通量更新对比

特征 advnx (E-W) advny (N-S) advnz (垂直)
mode 分支 有（1=线性/2=抛物线） 无（仅抛物线） 无（仅抛物线）
prec 阈值 1e-15 (CRAY 1e-24) 1e-15 (CRAY 1e-24) 1e-13 (CRAY 1e-24)
sigma clip 无 无 [0,1] min/max
Semi-Lagrangian 有（CFL>1 累加） 无 无
极区处理 极点行排除 SL 面积加权分配 无
周期包裹 经向 mod 循环 无 无
边界条件 无 极点 SSUM 顶底零通量
通量数组范围 ip1jmp1 ip1jm ip1jmp1×(llm+1)
更新范围 iip2+1..ip1jm iip2..ip1jm + 极点 1..ip1jmp1

守恒分析

advn 的通量更新同时修改 q 和 masse，这保证了质量和 tracer 的耦合守恒：

new_m = masse + flux_mass_in - flux_mass_out
q_new = (q*masse + flux_tracer_in - flux_tracer_out) / new_m

这意味着 q_new * new_m = q*masse + flux_tracer_in - flux_tracer_out，tracer 总量变化等于净通量。

守恒风险点：

  1. advnx Semi-Lagrangian 累加中，zu_m 逐格点减去 masse，最后的 zu_m 是剩余质量通量。若累加精度不足，可能导致 u_mq 不完全等于原始通量。
  2. advny 极区 SSUM 求和可能因浮点累加顺序导致微小不守恒。
  3. advnz 的 zsigb clip 到 [0,1] 会截断极端 sigma 值，可能破坏守恒。

复现要点

  1. advny 和 advnz 没有 mode 分支，始终使用抛物线公式。这意味着即使 mode=1（advn 主例程传入），N-S 和垂直方向仍使用抛物线通量。mode 参数只在 advnx 的 E-W 方向生效。
  2. advnz 的 prec=1e-13 比 advnx/advny 的 prec=1e-15 宽松两个量级。这意味着垂直方向更容易退化为 1 阶（当 |qb-qh| < 1e-13 时）。
  3. advnx 的 Semi-Lagrangian 累加只在 n0 > 1 时执行。n0=0 或 n0=1 时不执行，使用常规通量公式。这可能导致单格点 CFL>1 时的不守恒。
  4. zsigp=0.1 的 fallback（当 zsig==0 时）避免了除零，但引入了非物理的截断。这个 0.1 是经验值。
  5. advny 的 aire(ij)/apoln 面积加权分配要求 apoln = Σ aire(ij)（极冠总面积），否则极点质量总和不守恒。
  6. advnx 的 mod(ij,iip1) /= 0 条件确保极点行不进入 Semi-Lagrangian 处理，因为极点行的经向索引无法做有意义的 mod 操作。

相关页面

待确认