caldyn 动力 tendency 分段

输入范围

dyn3d\caldyn.F
dyn3dpar\caldyn_p.F

Mars 运行参与度：条件经过。串行动力主循环调用 caldyn，并行动力主循环调用 caldyn_p；Mars 是否走 COMMON 版本取决于构建入口和并行配置。#ifdef NODYN 编译选项会将全部 tendency 置零并跳过。

页面定位

本页深入 caldyn/caldyn_p 的内部调用管线，按功能段拆解每个算子的角色、数据依赖和输出写入点。caldyn 串并行对照 提供高层比较；本页补充每个段的算子链路、并行 SAVE 语义、convmas 拆分原因和 advect 命名差异。

签名

串行和并行签名完全一致：

SUBROUTINE caldyn[_p]
 $ (itau, ucov, vcov, teta, ps, masse, pk, pkf, tsurpk, phis,
 $  phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)

注意 phi 是 IN 参数——它在 leapfrog 调用 caldyn 之前由 geopot 预计算传入。

管线总览

caldyn 内部按 5 段顺序执行，每段依赖前段输出：

A. 风场转换和压力：   covcont + pression [+ psextbar]
         ↓
B. 质量初始化：       massdair → massbar → massbarxy → flumass
         ↓
C. 位温通量和 dp：    dteta1 → convmas → dp 计算 → vitvert
         ↓
D. 涡度和能量：       tourpot → dudv1 → enercin → bernoui → dudv2
         ↓
E. 垂直平流：         ang 构造 → advect/advect_new_p → 周期性修正

段 A：风场转换和压力

CALL covcont[_p]  (llm, ucov, vcov, ucont, vcont)
CALL pression[_p] (ip1jmp1, ap, bp, ps, p)

covcont 把协变风转为逆变风，供后续 flumass 使用。pression 从 ap/bp/ps 构造层界压力 p。

串行额外调用 psextbar(ps, psexbarxy) 计算地表压力 XY 面积加权平均。并行版本此调用被注释掉（L84 cym），因为 psexbarxy 在并行 caldyn_p 中未被后续使用。

Include 依赖：dimensions.h, paramet.h, comgeom.h（airesurg 等在后续段使用）。USE 依赖：comvert_mod（ap, bp）。

详细签名和公式参见 压力与质量算子 family 和 风场转换与质量通量算子组。

段 B：质量初始化

CALL massdair[_p]  (p, masse)
CALL massbar[_p]   (masse, massebx, masseby)
CALL massbarxy[_p] (masse, massebxy)
CALL flumass[_p]   (massebx, masseby, vcont, ucont, pbaru, pbarv)

四步从层界压力 p 构造出每层空气质量 masse、U/V/Z 点质量平均 massebx/masseby/massebxy，再与逆变风相乘得到质量通量 pbaru/pbarv。

massebxy 在 Z 点（角点）上计算，被段 D 的 tourpot 使用。pbaru/pbarv 同时被段 C 和段 D 使用。

段 C：位温通量和 dp

CALL dteta1[_p]  (teta, pbaru, pbarv, dteta)
CALL convmas[_p/1_p+2_p] (pbaru, pbarv, convm)

dteta1 计算位温水平通量散度：先构造 hbxu = pbaru * 0.5*(teta(ij) + teta(ij+1)) 和 hbyv = pbarv * 0.5*(teta(ij) + teta(ij+iip1))，再调用 convflu 求散度，最后 filtreg 滤波。输出 dteta 是位温的水平汇聚 tendency。

convmas 计算质量水平通量散度 convm。并行版本拆为两步：convmas1_p 计算水平散度部分，$OMP BARRIER 同步后 convmas2_p 执行垂直积分部分。拆分原因是两步之间存在数据依赖——垂直积分需要水平散度的完整结果。

dp 直接计算（串行 L94-96，并行 L120-128）：

DO ij = ijb, ije
   dp(ij) = convm(ij, 1) / airesurg(ij)
ENDDO

地表压力 tendency 等于第一层质量散度除以地表面积。并行版本在 OMP MASTER 区内执行，用 ij_begin:ij_end 子域范围。

垂直速度：

CALL vitvert[_p] (convm, w)

从质量散度垂直积分得到 sigma 坐标垂直速度 w，被段 E 的 advect 使用。

段 D：涡度和能量

CALL tourpot[_p]  (vcov, ucov, massebxy, vorpot)
CALL dudv1[_p]    (vorpot, pbaru, pbarv, du, dv)
CALL enercin[_p]  (vcov, ucov, vcont, ucont, ecin)
CALL bernoui[_p]  (ip1jmp1, llm, phi, ecin, bern)
CALL dudv2[_p]    (tsurpk, pkf, bern, du, dv)

这是动力 tendency 的核心组装段。数据管线：

1. tourpot 在 Z 点计算位势涡度 vorpot = (filtreg(curl) + fext) / massebxy。
2. dudv1 用位涡和质量通量计算旋转对 du/dv 的贡献（Coriolis-like 项）。
3. enercin 在 P 点计算 C-grid 动能 ecin，使用 covariant × contravariant 面积加权。
4. bernoui 计算 Bernoulli 函数 bern = phi + ecin，再 filtreg 滤波。
5. dudv2 把 Bernoulli 和压力梯度项加到 du/dv 上。

du/dv 在段 D 结束时已包含旋转、Bernoulli 和压力梯度三种贡献。段 E 再叠加垂直平流。

详细签名和公式参见 动能、Bernoulli、位势涡度与位势高度算子。

段 E：垂直平流

DO l = 1, llm
   DO ij = ijb, ije
      ang(ij, l) = ucov(ij, l) + constang(ij)
   ENDDO
ENDDO
CALL advect[_new_p] (ang, vcov, teta, w, massebx, masseby, du, dv, dteta)

ang 是绝对角速度（ucov + constang），constang 来自 comgeom.h 的行星自转常量。

串行调用 advect（171 行），并行调用 advect_new_p（290 行），旧并行版 advect_p（217 行）不再被调用。advect_new_p 显式处理更多并行边界条件和 halo 交换。

advect 把 ang/vcov/teta 和 w/massebx/masseby 转交给 advx/advy/advz 底层内核，对 du/dv/dteta 叠加垂直平流贡献。

并行 ang 构造的循环范围：ijb = ij_begin - iip1，ije = ij_end + iip1，向两侧各扩展一列 halo 宽度。极点条件：pole_nord 时 ijb = ij_begin（不向北扩展），pole_sud 时 ije = ij_end（不向南扩展）。使用 OMP DO SCHEDULE(STATIC, OMP_CHUNK)。

详细对照参见 advect 文件组对照 和 advxyz 串行内核。

后处理：周期性修正和守恒诊断

周期性修正（串行 L124-133，并行 L186-202）：

在每层的 ij = 1, ip1jm, iip1（即每行第一个格点）处，检查 dv(ij,l) 和 dv(ij+iim,l) 是否相等，若不等则强制 dv(ij+iim,l) = dv(ij,l)。这是因为 V 点网格在东西边界处应满足周期性，但浮点精度差异可能导致不一致。

并行版本范围：ijb = ij_begin，ije = ij_end，pole_sud 时 ije = ij_end - iip1（南极行不需修正，因为不存在对应的 wrap-around 行）。使用 OMP DO SCHEDULE(STATIC, OMP_CHUNK) + END DO NOWAIT。

守恒诊断（conser 触发）：

IF (conser) THEN
   CALL sortvarc (itau, ucov, tsurpk, ps, masse, pk, phis, vorpot, phi, bern, dp, time, vcov)
ENDIF

sortvarc 输出守恒诊断变量，用于调试和质量监控。并行注释提到需要 collective communication（也在 advect 中出现）。

调用链上下文

caldyn/caldyn_p 在 leapfrog 主循环中的位置：

leapfrog/leapfrog_p 主循环 (label 2):
  pression + exner_hyb/exner_milieu     ← 压力和 Exner 刷新
  tpot2t + geopot                        ← 温度→位势
  [NODYN: 全部 tendency 置零]
  ── caldyn/caldyn_p ──                  ← 本页
  [generic + ok_guide: du nudging]       ← 通用行星 guide 修正
  caladvtrac                             ← tracer advection
  addfi                                  ← 物理 tendency 加和
  integrd                                ← 状态积分写回

phi 在 caldyn 之前由 geopot 计算并作为 IN 参数传入。caldyn 输出的 tendency 随后被 integrd 使用，更新 prognostic state ucov/vcov/teta/ps。

串并行差异

SAVE 语义

并行版本给 vcont, ucont, ang, p, massebx, masseby, psexbarxy, vorpot, ecin, bern, massebxy, convm 等局部数组加了 SAVE 属性。这是因为 OpenMP 并行区中线程的栈空间有限，大数组需要放到静态存储区。串行版本无此需求。

convmas 拆分原因

convmas 并行拆为 convmas1_p + convmas2_p 是因为内部存在两步数据依赖：第一步计算水平质量通量散度（纯水平操作，可并行），第二步做垂直积分（需要第一步的完整结果）。中间的 BARRIER 保证所有线程完成第一步后再进入第二步。

数据流图

                    ucov, vcov, teta, ps, pk, pkf, tsurpk, phi
                              │
段 A:  covcont(ucov,vcov) → ucont, vcont
       pression(ap,bp,ps) → p
       [psextbar(ps) → psexbarxy]   (仅串行)
                              │
段 B:  massdair(p) → masse
       massbar(masse) → massebx, masseby
       massbarxy(masse) → massebxy
       flumass(massebx,masseby,vcont,ucont) → pbaru, pbarv
                              │
段 C:  dteta1(teta,pbaru,pbarv) → dteta     (convflu + filtreg)
       convmas(pbaru,pbarv) → convm
       dp(ij) = convm(ij,1) / airesurg(ij)
       vitvert(convm) → w
                              │
段 D:  tourpot(vcov,ucov,massebxy) → vorpot
       dudv1(vorpot,pbaru,pbarv) → du, dv   (旋转贡献)
       enercin(vcov,ucov,vcont,ucont) → ecin
       bernoui(phi,ecin) → bern             (phi+ecin, filtreg)
       dudv2(tsurpk,pkf,bern) → du, dv      (叠加 Bernoulli+压力)
                              │
段 E:  ang = ucov + constang
       advect[_new_p](ang,vcov,teta,w,massebx,masseby) → du,dv,dteta
                              │
后处理: dv 周期性修正
        [conser: sortvarc 诊断输出]
                              │
                    du, dv, dteta, dp, w, pbaru, pbarv, masse

复现要点

• caldyn 不是时间积分器，只计算 tendency。实际状态更新由 integrd 完成。
• phi 必须是 caldyn 调用前由 geopot 刷新过的值——过时的 phi 会导致 Bernoulli 项错误。
• 段 A 的 covcont 输出 ucont/vcont 同时被段 B 的 flumass 和段 D 的 enercin 使用，不能省略。
• 段 C 的 dp 计算在并行版本中由 OMP MASTER 独占执行，其余线程在后续 BARRIER 等待。
• 并行版本 convmas 拆分不可合并：两步之间的 BARRIER 是正确性保证。
• advect_new_p 替代了旧 advect_p；如果复现代码仍使用 advect_p，结果可能不同。
• 周期性修正是浮点精度 workaround，非物理修正。复现时如用双精度可能不需要。
• NODYN 编译选项下 caldyn 不被调用，全部 tendency 在 leapfrog 中被置零。

待确认

• psexbarxy 在并行版本中被注释掉后，是否有其他地方（如 leapfrog_p 初始化段）等价计算。
• dteta1_p 的滤波参数 (jjp1, llm, 2, 2, .true., 1) 与串行的 (jjp1, llm, 2, 2, .true., 1) 是否完全一致——需在 dteta1_p.F 中确认。
• advect_new_p 和 advect_p 的具体算法差异需在 advect 页面深入对比。
• conser 逻辑量在 leapfrog 中的设置条件需由 leapfrog 重点页 确认。
• Write_Field_p 的 #ifdef DEBUG_IO 默认是否编译开启。

相关页面

• caldyn 串并行对照 — 高层对照，不展开段内细节。
• 压力与质量算子 family — 段 A/B 算子详细签名。
• 风场转换与质量通量算子组 — 段 A/B 风场和通量算子。
• 动能、Bernoulli、位势涡度与位势高度算子 — 段 D 算子详细签名。
• Exner 函数与 cp(T) 变比热路径 — pk/pkf/tsurpk 的预计算来源。
• 水平微分/耗散 stencil family 深度页 — 段 D 算子的底层差分实现。
• advect 文件组对照 — 段 E advect/advect_new_p 包装层。
• advxyz 串行内核 — 段 E 底层方向内核。
• 差分算子总览 — 所有算子的分类索引。
• Mars 水平耗散链路 — caldyn 之后的耗散段。
• leapfrog 动力与输送 — caldyn 在时间推进中的调度位置。
• bilan_dyn — 使用 caldyn 输出的诊断例程。