vdif_kc.F

路径

LMDZ.MARS\libf\phymars\vdif_kc.F

所属目录/模块

libf/phymars

文件定位

vdif_kc.F 是 vdifc_mod.F 在未启用 callyamada4 且未启用 callatke 时使用的垂直湍流扩散 K 系数方案。它根据层界面高度、层中心高度、风速、位温、地表拖曳系数和上一时步/输入的 q2，更新湍流速度尺度 q=sqrt(q2)，并输出动量扩散系数 km 与标量扩散系数 kn。在 vdifc_mod.F 中，km/kn 对应后续隐式扩散使用的 zkv/zkh。

该例程还会在首次调用时从 tracer_mod:noms 查找 "co2" tracer。如果存在 CO2 tracer，它用 CO2 质量混合比构造一个修正的稳定度变量 zhc，使 Brunt-Vaisala 型稳定度项包含平均分子量变化的影响。

定义的符号

依赖的模块

调用的关键例程

上游调用点

输入

输出

共享状态与副作用

• 读取共享 tracer 名称：noms(:) 来自 tracer_mod，只用于首次定位 "co2" tracer。
• 持久局部状态：ico2、A、B、firstcall 都是 SAVE + !$OMP THREADPRIVATE，每个 OpenMP 线程保留自己的首次调用状态。
• 原地更新：q2 是 INTENT(INOUT)，例程会裁剪下限并改写内部界面和顶界面值。
• 无文件 I/O/诊断输出：本文件不写 NetCDF、不读配置文件，也不直接修改其他模块变量。

核心逻辑

1. 首次 CO2 tracer 定位（行 236-255）：若 firstcall 为真，遍历 noms(1:nq) 查找名称等于 "co2" 的 tracer。找到后设置 ico2，并用 m_co2=44.01e-3 kg/mol 与 m_noco2=33.37e-3 kg/mol 预计算 A 和 B。之后 firstcall=.false.。
2. 构造稳定度用变量 zhc（行 270-279）：若 ico2 非 0，则 zhc = teta * (A*zq(:,:,ico2)+B)；否则 zhc=teta。后续 n2 使用 zhc 的层间差分而不是原始 teta。
3. 初始化 q2/q 与地表边界（行 285-297）：所有界面 q2 至少为 q2min=1e-3；最底界面 q2(:,1) 由 b1^(2/3) * cd * (u1^2+v1^2) 重新给定并同样裁剪下限。
4. 计算垂直几何倒数（行 313-327）：unsdz 是层厚倒数，unsdzdec 是相邻层中心或边界到层中心的距离倒数。源码注释特别提醒顶界面 zlev(:,nlev) 必须由调用方提供。
5. 计算剪切与稳定度（行 331-380）：底界面剪切由第一层风速除以近地距离得到；内部界面 n2 用 g * dz^-1 * 2*(zhc_l-zhc_{l-1})/(zhc_l+zhc_{l-1}) 计算，若为负则置 0。内部剪切 m2 来自相邻层风速差，顶界面复制倒数第二个界面。
6. Blackadar 混合长度与稳定度函数（行 387-449）：混合长度 long = kappa*z/(1+z/long0)，其中 kappa=0.4、long0=160 m。gn=-long^2*n2/q2 被限制到 [-10,0.0233]，再计算 sn 和 sm；若 gn 被截断或导数项符号不合理，snq2/smq2 被置 0。
7. 初始 kn/km（行 457-486）：底界面取 knmin=kmmin=1e-5；内部界面按 long*q*sn 与 long*q*sm 计算；顶界面复制倒数第二个界面。
8. 内部界面稳态耦合求解（行 497-605）：对 ilev=2:nlev-1，先把剪切生产、浮力项和耗散项拆成 q、q^3、q*m2、q^3*m2 贡献；然后求一个使局地 q2 方程与速度梯度方程同时达到稳态的 m2cstat/qcstat。ilev=2 使用地表拖曳 cd 的特殊边界项，其余界面使用上下相邻界面的 kmpre/mpre。
9. 最终裁剪与重算扩散系数（行 584-605）：若新 q2 小于 q2min，重置为下限；随后用新 q2 重新计算 gn/sn/sm 以及最终 kn/km。
10. 上下边界收尾（行 616-622）：底界面 kn/km 再次设为最小值；顶界面 q2/q/kn/km 复制倒数第二个界面。

伪代码

vdif_kc:
  if firstcall:
    ico2 = index of "co2" in noms, or 0
    if ico2 exists:
      A = 1/m_co2 - 1/m_noco2
      B = 1/m_noco2
    firstcall = false

  if ico2 exists:
    zhc = teta * (A*zq(:,:,ico2) + B)
  else:
    zhc = teta

  q2 = max(q2, q2min)
  q = sqrt(q2)
  q2(:,1) = max(b1^(2/3) * cd * (u(:,1)^2 + v(:,1)^2), q2min)

  compute layer-thickness inverses unsdz and center-distance inverses unsdzdec
  compute shear m2/m and stable stratification n2
  if n2 < 0: n2 = 0

  for internal interfaces:
    long = Blackadar(kappa, long0, height_above_surface)
    gn = clamp(-long^2*n2/q2, gnmin, gnmax)
    sn, sm = stability_functions(gn)
    snq2, smq2 = Taylor q2 derivative terms, or 0 if clipped/sign mismatch
    kn = long*q*sn
    km = long*q*sm

  for internal interfaces:
    split q2 source terms into q, q^3, q*m2, q^3*m2
    solve stationary coupled q2/shear estimate:
      m2cstat = m2 - source_sum/(dt*2*km)
      mcstat = sqrt(m2cstat)
      if lowest internal interface:
        include surface-drag boundary term with cd
      else:
        use neighboring kmpre/mpre
      qcstat = (kmcstat / (sm/q2) / long)^(1/3)
    q2 = max(qcstat^2, q2min)
    recompute sn/sm and final kn/km

  bottom kn/km = minimum values
  top q2/q/kn/km = values from previous interface

参与的主题流程

写法特点

• 固定格式 Fortran：文件扩展名 .F，使用固定格式续行和数字标签 10001/10002。
• 只有一个外部模块依赖：除 tracer_mod:noms 外不依赖其他模块，核心算法完全在本例程内展开。
• SAVE + THREADPRIVATE 首次初始化：CO2 tracer 索引和系数每个 OpenMP 线程缓存一次。
• 稳定度只处理稳定层结：内部 n2<0 时直接置 0；源码注释说明对流调整应保证进入本例程的是稳定构型。
• 下限保护不等于上下限保护：q2min、knmin、kmmin 被实际使用；q2max 虽声明为参数，但源码没有对 q2 做上限裁剪。

复现要点

• vdif_kc 只在 vdifc 的 fallback 路径执行：callyamada4=.false. 且 callatke=.false.。若启用 Yamada4 或 ATKE，vdifc 会走其他例程。
• 传入的 zlev(:,nlay+1) 必须是真实顶界面；源码中原先补顶界面的代码已被注释掉。
• CO2 修正取决于 tracer_mod:noms 首次调用时是否存在 "co2"；如果没有，稳定度直接用 teta。
• n2 的负值会被钳到 0，因此复现不稳定层结时不能期待本例程显式产生对流增强。
• 底界面 q2 每次都会由 cd*(u1^2+v1^2) 重置，不只是沿用输入 q2。
• vdifc 在 vdif_kc 返回后，如果 doubleq 且 ngrid==1，还会把 zkh/zkv 的内部层最低值抬到 80；该保护不在本文件内。

待确认

• m2cstat 在行 539-541 直接进入 sqrt(m2cstat)，源码没有非负保护。复现风险：若源项组合使 m2cstat<0，可能产生 NaN；需开发者确认输入范围是否保证其非负。
• ilev=2 的地表边界分支在行 552 使用 u(igrid,3)/v(igrid,3)。复现风险：该写法假设 nlay>=3，常规 3D/1D 大气列应满足，但极简测试网格需要避免少于 3 层。
• q2max=1e2 在行 191 声明但未使用；不应把它当作实际上限。
• firstcall 缓存意味着若一次运行中动态改变 noms 或 tracer 顺序，本例程不会重新扫描 CO2 索引；当前流程通常在初始化后固定 tracer 表。

相关页面

• vdif_cd_mod — 地表层动量/热量交换系数，向 vdifc 提供 cd/zcdv_true_tmp 等输入背景。
• vdifc-water-surface-exchange — vdifc 中水汽地表交换段；本页覆盖的是同一 vdifc 流程更早的 K 系数计算。
• pbl_parameters_mod — PBL 诊断参数页，和本页共享边界层/稳定度背景，但该页不改变物理倾向。
• tracer_mod — noms 来源，决定 CO2 tracer 是否触发稳定度修正。
• phymars index — phymars 目录总览。