Exner 函数与 cp(T) 变比热路径

输入范围

dyn3d_common\exner_hyb_m.F90
dyn3d_common\exner_milieu_m.F90
dyn3d_common\cpdet_mod.F90
dyn3dpar\exner_hyb_p_m.F90
dyn3dpar\exner_milieu_p_m.F90

Mars 运行参与度：条件经过。Exner 函数在每个时间步由 COMMON leapfrog 和 MARS phymars/leapfrog_nogcm 调用；cpdet_mod 的 tpot2t/t2tpot 在动力和物理耦合中被广泛调用；cpofT 路径由 run.def 配置决定。

页面定位

本页覆盖 Mars 热力学转换的两个核心模块族：Exner 函数计算（exner_hyb vs exner_milieu）和变比热路径（cpdet_mod 中的 cpdet、tpot2t/t2tpot）。这些模块是 leapfrog 时间步中压力→Exner→温度→位势链路的中间环节。pression 的压力计算由 压力与质量算子 family 覆盖；Exner 的下游使用（geopot）由 动能、Bernoulli、位势涡度与位势高度算子 覆盖。

模块文件一览

Exner 函数选择：pressure_exner

leapfrog.F L264-268 按 pressure_exner 逻辑量选择路径：

if (pressure_exner) then
  CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
else
  CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
endif

pressure_exner 来自 comvert_mod，由 iniconst 初始化为 disvert_type == 1（标准 hybrid 坐标 → exner_hyb），可通过 getin('pressure_exner', ...) 覆盖。Mars 复杂垂直坐标通常使用 disvert_type != 1，走 exner_milieu 路径。

两个变体的签名完全一致：

SUBROUTINE exner_hyb(ngrid, ps, p, pks, pk, pkf)
SUBROUTINE exner_milieu(ngrid, ps, p, pks, pk, pkf)
  INTEGER ngrid
  REAL ps(ngrid), p(ngrid,llmp1), pks(ngrid), pk(ngrid,llm)
  REAL, optional :: pkf(ngrid,llm)

输入是地表压力 ps 和层界压力 p（由 pression 计算），输出是地表 Exner 函数 pks、各层 Exner 函数 pk 和可选的滤波版本 pkf。

exner_hyb — 递归 alpha/beta 法

公式

标准 Exner 函数定义：pk = Cp * (p/preff)^kappa，其中 preff 来自 comvert_mod（参考压力），kappa 和 cpp 来自 comconst_mod。

exner_hyb 使用 Fr. Hourdin 的递归方法保证能量守恒（总能量/内能/势能比例关系）：

1. Shallow Water 特例（llm==1）：pks = (cpp/preff)*ps，pk = 0.5*pks。首调用时检查 kappa==1 且 cpp==r。
2. 地表 Exner：pks(ij) = cpp * (ps(ij)/preff)^kappa。
3. alpha/beta 系数：从顶层（l=llm）向下递归计算。alpha(ij,llm)=0，beta(ij,llm)=1/(1+2*kappa)。逐层向下：dellta = p(l)*(1+2k) + p(l+1)*(beta(l+1)-(1+2k))，alpha(l) = -p(l+1)/dellta * alpha(l+1)，beta(l) = p(l)/dellta。
4. 底层 pk：pk(ij,1) = (p(1)*pks - 0.5*alpha(2)*p(2)) / (p(1)*(1+kappa) + 0.5*(beta(2)-(1+2k))*p(2))。
5. 向上递推：pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)，l=2..llm。
6. 可选滤波：若 present(pkf)，复制 pk 到 pkf 再 filtreg(pkf, jmp1, llm, 2, 1, .TRUE., 1)。

并行变体

exner_hyb_p_m.F90 使用 USE parallel_lmdz，子域 ijb=ij_begin, ije=ij_end。每个水平循环包裹 !$OMP DO SCHEDULE(STATIC)，垂直递归的层间依赖通过 !$OMP BARRIER 同步。firstcall 声明为 !$OMP THREADPRIVATE。滤波用 filtreg_p(pkf, jjb, jje, jmp1, llm, 2, 1, .TRUE., 1)。

exner_milieu — 简化中点法（Mars 专用）

公式

源码 WARNING 注释明确：这是 Mars 复杂垂直坐标的专用版本，不满足总能量/内能/势能的比例关系（F. Forget, 2001）。

1. Shallow Water 特例（llm==1）：与 exner_hyb 相同。
2. 地表 Exner：与 exner_hyb 相同，pks(ij) = cpp * (ps(ij)/preff)^kappa。
3. 各层 pk（l=1..llm-1）：直接中点公式 pk(ij,l) = cpp * (2*preff)^(-kappa) * (p(l) + p(l+1))^kappa。
4. 顶层外推（l=llm）：对数等距外推 pk(ij,llm) = pk(ij,llm-1)^2 / pk(ij,llm-2)。
5. 可选滤波：同 exner_hyb。

并行变体

exner_milieu_p_m.F90 结构与 exner_hyb_p 类似，使用 USE parallel_lmdz，子域 + OMP STATIC 循环 + !$OMP THREADPRIVATE(firstcall)。由于 exner_milieu 没有层间递归（每层独立计算），并行化更简单。

cpdet_mod — 变比热模块

ini_cpdet

初始化 nu_venus 和 t0_venus（在 comconst_mod 中声明）：

if (cpofT) then
    nu_venus = 0.35
    t0_venus = 460.
else
    nu_venus = 0.
    t0_venus = 0.
endif

由 gcm.F90 L192 / gcm.F L220 / nogcm.F90 L184 在初始化阶段调用。cpofT 是 control_mod 中的逻辑量，由 conf_gcm 从 run.def 读取（默认 .False.）。

cpdet(t) — 变比热函数

FUNCTION cpdet(t)
  if (cpofT) then
      cpdet = cpp * (t / t0_venus)**nu_venus
  else
      cpdet = cpp
  endif

当 cpofT=.true. 时，比热随温度变化：cp(T) = cpp * (T/460)^0.35。这是 Venus 厚大气的适配，Mars 通常 cpofT=.false.（常比热）。

在 leapfrog.F L722 / leapfrog_p.F L1412 中被内联调用，用于 conservative dissipation 修正：dtec = (ecin0-ecin)/cpdet(temp)。在 diagedyn.F L229 中用于诊断能量计算。

tpot2t — 位温→温度

SUBROUTINE tpot2t(npoints, yteta, yt, ypk)
  if (cpofT) then
      yt = yteta**nu_venus + nu_venus * t0_venus**nu_venus * log(ypk/cpp)
      yt = yt**(1./nu_venus)
  else
      yt = yteta * ypk / cpp
  endif

常比热路径：T = theta * pk / cpp（标准 Exner 反变换）。 变比热路径：非线性反演 T = (theta^nu + nu * t0^nu * ln(pk/cpp))^(1/nu)。

t2tpot — 温度→位温

SUBROUTINE t2tpot(npoints, yt, yteta, ypk)
  if (cpofT) then
      yteta = yt**nu_venus - nu_venus * t0_venus**nu_venus * log(ypk/cpp)
      yteta = yteta**(1./nu_venus)
  else
      yteta = yt * cpp / ypk
  endif

常比热路径：theta = T * cpp / pk。变比热路径是非线性正变换。

代码中包含一段 "ATMOSPHERE PROFONDE" 深大气修正（ratio_mod = mmm0/mmm(p) 分子量修正），但被 if (1 .EQ. 0) 永久禁用，当前不影响任何活跃路径。

并行变体

cpdet_mod.F90 在 #ifdef CPP_PARA 条件下包含四个并行子例程：

_glo_p 变体使用 #include "dimensions.h" 和 #include "paramet.h" 获取网格维度，直接在完整动力网格上操作，比多次调用 _p 切片更高效。

leapfrog 中的调用链路

每个 leapfrog 时间步中，Exner 和 cp(T) 转换的典型调用顺序：

pression(ip1jmp1, ap, bp, ps, p)           → p 层界压力

if pressure_exner:
  exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)  → pks, pk, pkf
else:
  exner_milieu(ip1jmp1, ps, p, pks, pk, pkf) → pks, pk, pkf

tpot2t(ijp1llm, teta, temp, pk)            → temp 温度
tsurpk = cpp * temp / pk                    → 位势计算用的 teta 替代量
geopot(ip1jmp1, tsurpk, pk, pks, phis, phi) → phi 位势高度

caldyn(ucov,vcov,teta,ps,...,pk,pkf,tsurpk,phis,phi,...)

在耗散守恒修正段中：

tpot2t(ijp1llm, teta, temp, pk)            → temp
covcont → enercin → ecin0                   → 耗散前动能
dissip → ucov+=dudis; vcov+=dvdis
covcont → enercin → ecin                    → 耗散后动能
dtec = (ecin0 - ecin) / cpdet(temp)         → 动能差/变比热
temp += dtec                                 → 温度修正
t2tpot(ijp1llm, temp, ztetaec, pk)          → 位温修正

调用方总览

Exner 函数

cpdet_mod 使用方

串并行差异汇总

输入输出数据流

conf_gcm:
  run.def → cpofT (logical)
  iniconst: disvert_type → pressure_exner (logical)

gcm initialization:
  ini_cpdet → nu_venus, t0_venus (comconst_mod)

leapfrog each timestep:
  pression(ap, bp, ps) → p
  if pressure_exner:
    exner_hyb(ps, p) → pks, pk, pkf
  else:
    exner_milieu(ps, p) → pks, pk, pkf
  tpot2t(teta, pk) → temp
  tsurpk = cpp * temp / pk
  geopot(tsurpk, pk, pks, phis) → phi

dissipative correction (if apdiss):
  tpot2t(teta, pk) → temp
  cpdet(temp) → local cp value
  t2tpot(temp, pk) → ztetaec

复现要点

• pressure_exner 的选择由 disvert_type 决定：Mars 通常 disvert_type != 1，走 exner_milieu，但该变体不保证能量守恒——这意味着 Mars 运行的 energy budget 可能有系统性偏差。
• cpofT 默认为 .False.，Mars 通常使用常比热 cpp。Venus 厚大气设置 cpofT=.True.。
• tpot2t 和 t2tpot 不是简单互逆——当 cpofT=.true. 时存在数值精度差异。
• exner_hyb 和 exner_milieu 的 pkf 是 optional 参数；guide_mod 调用时不传 pkf（5 参数形式），leapfrog 调用时传 pkf（6 参数形式）。
• 并行 _glo_p 变体比多次 _p 切片调用更高效，leapfrog_p 统一使用 _glo_p，calfis_p 使用 _p 切片。

待确认

• Mars 当前 run.def/conf_gcm 中 disvert_type 的实际值和 pressure_exner 的覆盖值需追踪确认。
• exner_milieu 的能量不守恒对 Mars long-run energy budget 的影响需在 bilan_dyn 页面中量化。
• cpdet_mod 中深大气修正（ratio_mod）的 Venus 条件和未来激活路径需在 Venus 专题页中覆盖。
• exner_hyb 的 alpha/beta 递归方法的具体数学推导（Fr. Hourdin 的 note）未包含在源码中。

相关页面

• 压力与质量算子 family
• 动能、Bernoulli、位势涡度与位势高度算子
• 质量、压力、位势和能量例程总览
• caldyn 串并行对照
• Mars 水平耗散链路
• leapfrog 耗散、sponge 和守恒修正
• 垂直离散构造
• conf_gcm
• COMMON 状态/控制模块约定
• guide 配置与初始化约定