动能、Bernoulli、位势涡度与位势高度算子

输入范围

dyn3d_common\bernoui.F
dyn3d_common\enercin.F90
dyn3d_common\tourpot.F90
dyn3d_common\geopot.F
dyn3dpar\bernoui_p.F
dyn3dpar\enercin_p.F
dyn3dpar\tourpot_p.F
dyn3dpar\geopot_p.F

Mars 运行参与度：条件经过。四个算子在 caldyn 动力倾向计算、COMMON leapfrog 耗散守恒修正和诊断路径中被活跃调用；MARS phymars/leapfrog_nogcm 仅调用 geopot。

页面定位

本页覆盖 Mars 动力核心中四个能量/涡度/位势短算子的签名、公式、Include 依赖、串并行对照和调用链路。这些算子共同支撑 caldyn 的动量倾向组装和 leapfrog 的能量守恒修正。算子本身的差分格式细节由 差分算子总览 覆盖；它们在耗守恒修正中的角色由 Mars 水平耗散链路 覆盖；caldyn 的完整调用顺序由 caldyn 串并行对照 覆盖。

算子概要

enercin — 动能

签名

SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin )
  REAL, INTENT(IN)  :: vcov   (ip1jm,   llm)
  REAL, INTENT(IN)  :: ucov   (ip1jmp1, llm)
  REAL, INTENT(IN)  :: vcont  (ip1jm,   llm)
  REAL, INTENT(IN)  :: ucont  (ip1jmp1, llm)
  REAL, INTENT(OUT) :: ecin   (ip1jmp1, llm)

公式

在 C-grid 上，标量点 P(i,j) 的动能由周围四个风场格点的 covariant × contravariant 乘积面积加权平均：

ecin(ij+1,l) = 0.5 * ( ucov(ij,l)  * ucont(ij,l)   * alpha3p4(ij+1)
                     + ucov(ij+1,l)* ucont(ij+1,l)  * alpha1p2(ij+1)
                     + vcov(ij-iim,l)*vcont(ij-iim,l)* alpha1p4(ij+1)
                     + vcov(ij+1,l)* vcont(ij+1,l)   * alpha2p3(ij+1) )

alpha1p2/alpha3p4/alpha1p4/alpha2p3 来自 comgeom.h，是由 inigeom/iniconst 初始化的面积权重对。ucov × ucont 给出 u² 的真实度量（协变 × 逆变 = 物理量的平方），v 同理。

极点处理：北极 ij=1..iip1 用 vcov*vcont*aire 面积加权求和除以 apoln；南极 ij=ip1jm+1..ip1jmp1 用对应南极面积除以 apols。周期性修正 ecin(ij,l) = ecin(ij+iim,l)。

并行变体

enercin_p.F 用 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 包裹外层 l 循环。子域范围 ijb=ij_begin, ije=ij_end+iip1，北极条件 pole_nord 时 ijb 上移 iip1（跳过北极点），南极条件 pole_sud 时 ije 下移 iip1（跳过南极点）。极点动能用 SSUM 替代串行 SUM 做行求和。

bernoui — Bernoulli 函数

签名

SUBROUTINE bernoui (ngrid, nlay, pphi, pecin, pbern)
  INTEGER nlay, ngrid
  REAL pphi(ngrid*nlay), pecin(ngrid*nlay), pbern(ngrid*nlay)

公式

Bernoulli 函数 = 位势 + 动能，再施加 filtreg 正则滤波：

pbern(ijl) = pphi(ijl) + pecin(ijl)
CALL filtreg( pbern, jjp1, llm, 2, 1, .true., 1 )

pphi 由 geopot 输出（位势高度），pecin 由 enercin 输出（动能）。滤波参数 (2,1,.true.,1) 表示 2 阶精度、1 次通过、带极点滤波。

Bernoulli 函数在 caldyn 中被 dudv2 用来计算动量倾向的梯度部分：dudv2(tsurpk, pkf, bern, du, dv)。

并行变体

bernoui_p.F 把加法循环改为 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 在 l 维度上，子域范围 ijb=ij_begin, ije=ij_end+iip1（pole_sud 时 ije=ij_end）。滤波用 filtreg_p(pbern, jjb, jje, jjp1, llm, 2, 1, .true., 1)，其中 jjb/jje 按 pole_sud 调整。

tourpot — 位势涡度

签名

SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot )
  REAL, INTENT(IN)  :: vcov    (ip1jm,   llm)
  REAL, INTENT(IN)  :: ucov    (ip1jmp1, llm)
  REAL, INTENT(IN)  :: massebxy(ip1jm,   llm)
  REAL, INTENT(OUT) :: vorpot  (ip1jm,   llm)

公式

位势涡度 = (滤波后的风场旋度 + 科氏参数) / Z 点空气质量：

vorpot = ( filtreg( d(vcov)/dx - d(ucov)/dy ) + fext ) / massebxy

三步计算：

1. 风场旋度：rot(ij,l) = vcov(ij+1,l) - vcov(ij,l) + ucov(ij+iip1,l) - ucov(ij,l)。这是 Z 点（角点）上的离散旋度。周期性修正 rot(iip1,j,l) = rot(1,j,l)。
2. 滤波：CALL filtreg(rot, jjm, llm, 2, 1, .FALSE., 1)。.FALSE. 表示不带极点滤波（旋度定义在 Z 点，不是 P 点）。
3. 除以质量：vorpot(ij,l) = (rot(ij,l) + fext(ij)) / massebxy(ij,l)。fext 来自 comgeom.h，是科氏参数加网格曲率项。massebxy 由 massbarxy 提供 Z 点面积加权空气质量。

vorpot 在 caldyn 中被 dudv1(vorpot, pbaru, pbarv, du, dv) 用来计算旋转诱导的动量倾向。

并行变体

tourpot_p.F 三段结构（旋度→滤波→除质量）分别用 OMP 和 filtreg_p：

• 旋度循环：c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 在 l 维度，子域 ijb=ij_begin-iip1（pole_nord 时 ijb=ij_begin），ije=ij_end（pole_sud 时 ije=ij_end-iip1-1）。
• 滤波：filtreg_p(rot, jjb, jje, jjm, llm, 2, 1, .FALSE., 1)，jjb=jj_begin-1（pole_nord 时 +1），jje=jj_end（pole_sud 时 -1）。
• 除质量：又一个 c$OMP DO 循环，边界同上。

geopot — 位势高度

签名

SUBROUTINE geopot ( ngrid, teta, pk, pks, phis, phi )
  INTEGER ngrid
  REAL teta(ngrid,llm), pks(ngrid), phis(ngrid), pk(ngrid,llm), phi(ngrid,llm)

公式

从地面向上逐层静力积分，计算每层中点的位势高度：

phi(ij,1) = phis(ij) + teta(ij,1) * (pks(ij) - pk(ij,1))     ! 近地面层
phi(ij,l) = phi(ij,l-1) + 0.5 * (teta(ij,l) + teta(ij,l-1))  ! 上层梯形积分
                      * (pk(ij,l-1) - pk(ij,l))

这里 teta 不是简单的位温 θ，而是经过 tpot2t 转换后的量 tsurpk = cpp * temp / pk，其中 cpp 和 pk 适配变比热 cp(T) 路径。注释原文："This computation (with teta = cp T / pk!) is identical to delta phi = R/RMD T/p delta p (r=R/RMD=cpp*kappa)"。

phis 是地表位势（地形），pks 是地表 Exner 函数，pk 是各层 Exner 函数。静力积分从地面（l=1）向上到模式顶（l=llm）。

并行变体

geopot_p.F 用 USE parallel_lmdz，子域范围 ijb=ij_begin, ije=ij_end+iip1（pole_sud 时 ije=ij_end），但没有 OMP 指令——地面层循环和逐层积分循环都是普通 DO。原因是垂直积分存在层间依赖 phi(ij,l) = phi(ij,l-1) + ...，无法并行化外层 l；而水平方向已在子域内顺序执行。

cp(T) 适配

活跃调用路径中 geopot 的第一个参数总是 tsurpk（经过 tpot2t 和 cpp*temp/pk 转换后的量），而非原始 teta。注释掉的旧代码 CALL geopot(..., teta, ...) 表明早期版本直接用位温。#ifdef 标记 ADAPTATION GCM POUR CP(T) 确认这是变比热改造。

caldyn 数据管线

caldyn（串行 L99-109 / 并行 L132-154）中这四个算子构成动量倾向组装的核心链：

tourpot(vcov, ucov, massebxy, vorpot)     → 位势涡度
  dudv1(vorpot, pbaru, pbarv, du, dv)     → 旋转诱导动量倾向
enercin(vcov, ucov, vcont, ucont, ecin)   → 动能
bernoui(ip1jmp1, llm, phi, ecin, bern)    → Bernoulli 函数
  dudv2(tsurpk, pkf, bern, du, dv)        → Bernoulli+压力动量倾向

前置依赖：tourpot 需要 massebxy（由 massbarxy 提供）；enercin 需要 ucont/vcont（由 covcont 提供）；bernoui 需要 phi（由 geopot 提供，已在 leapfrog 调用 caldyn 之前计算好）。

leapfrog 耗散守恒中的 enercin

leapfrog（串行 L699-728 / 并行 L1353-1401）的 apdiss 块中，enercin 被调用两次：

1. 耗散前基准：covcont → enercin(vcov,ucov,vcont,ucont,ecin0)，记录耗散前动能。
2. 耗散后守恒修正：ucov += dudis; vcov += dvdis 后重新 covcont → enercin → ecin，计算动能差 ecin0 - ecin 作为热力修正注入温度场。

详见 Mars 水平耗散链路。

leapfrog 中 geopot 的多点调用

geopot 在 COMMON leapfrog.F 中被调用 6 次（串行），在 leapfrog_p.F 中被调用 4 次（并行），在 MARS phymars/leapfrog_nogcm.F 中被调用 5 次。每次都在 tpot2t → tsurpk = cpp*temp/pk 之后调用，确保用 cp(T) 适配后的温度。

主要调用位置：

• 初始时间步（L426）：caldyn 之前的首次位势计算。
• Matsuno 步（L545）：Matsuno 前向步的位势刷新。
• leapfrog 步（L838, L872）：leapfrog 时间步的位势刷新。
• forward 步（L1006, L1031）：forward 模式的位势刷新。

MARS phymars/leapfrog_nogcm.F 中保留全部 geopot 调用（5 次），但不调用 enercin/bernoui/tourpot——nogcm 路径无动力倾向计算。

诊断调用方

enercin 在诊断路径中有额外调用：

• bilan_dyn.F L422 / bilan_dyn_p.F L362：能量收支诊断。
• diagedyn.F L163：动力诊断。

tourpot/bernoui/geopot 在诊断路径中无额外调用（仅 geopot 在 iniacademic.F90 L249 有学术测试初始化调用）。

串并行差异汇总

输入输出数据流

leapfrog:
  tpot2t(teta) → temp; tsurpk = cpp*temp/pk
  geopot(ip1jmp1, tsurpk, pk, pks, phis) → phi       ─┐
                                                       │
caldyn:                                                │
  massbarxy(masse) → massebxy ─┐                       │
  covcont(ucov,vcov) → ucont,vcont ─┐                  │
  tourpot(vcov, ucov, massebxy) → vorpot              │
    dudv1(vorpot, pbaru, pbarv) → du,dv               │
  enercin(vcov, ucov, vcont, ucont) → ecin            │
  bernoui(ip1jmp1, llm, phi←──────┘, ecin) → bern     │
    dudv2(tsurpk, pkf, bern) → du,dv (叠加)           │
                                                       │
leapfrog dissip block:                                 │
  covcont(ucov,vcov) → ucont,vcont                     │
  enercin(vcov,ucov,vcont,ucont) → ecin0  (耗散前)    │
  dissip(vcov,ucov,teta,p) → dudis,dvdis,dtetadis     │
  ucov += dudis; vcov += dvdis                         │
  covcont → ucont,vcont (更新后)                       │
  enercin → ecin (耗散后)                              │
  temp += (ecin0 - ecin) / cpdet(temp)                 │

复现要点

• enercin 需要 covariant 和 contravariant 两种风场；如果只有 covariant，需先通过 covcont 转换。
• bernoui 的 filtreg 调用改变了输出 pbern 的值，不能跳过——未滤波的 phi+ecin 不满足 dudv2 的输入假设。
• tourpot 在 Z 点（角点）上计算，因此滤波参数 .FALSE. 表示不带极点处理；而 bernoui 在 P 点上用 .true.。
• geopot 的 teta 参数必须是 tsurpk = cpp*temp/pk 而不是原始位温，否则会丢失 cp(T) 变比热修正。
• 并行 tourpot_p 的旋度循环需要 ijb=ij_begin-iip1 扩展，因为旋度 stencil 需要左侧一列 halo 数据。

待确认

• filtreg 的参数 (2,1,.true.,1) 和 (2,1,.FALSE.,1) 的精确语义需在 filtreg 页面中补齐（当前只有参数位置推断）。
• SSUM 外部函数在 enercin_p.F 中声明为 EXTERNAL，具体实现在哪个文件中需追踪（可能在 parallel_lmdz 或独立 helper 中）。
• geopot_p 没有 OMP 指令是否影响并行性能——其子域内水平循环可能已被外层 OMP 区域覆盖，需在 leapfrog_p 上下文中确认。

相关页面

• caldyn 串并行对照
• caldyn 动力 tendency 分段
• Mars 水平耗散链路
• leapfrog 耗散、sponge 和守恒修正
• 差分算子总览
• 水平微分/耗散 stencil family 深度页
• 压力与质量算子 family
• 风场转换与质量通量算子组
• 质量、压力、位势和能量例程总览
• Exner 函数与 cp(T) 变比热路径
• leapfrog 动力与输送
• bilan_dyn