yamada4.F

路径

LMDZ.MARS\libf\phymars\yamada4.F

所属目录 / 模块

libf/phymars

文件定位

yamada4.F 定义 Mars 版本的 Mellor-Yamada 湍流闭合模块 yamada4_mod。主例程 yamada4 在 vdifc_mod 的 callyamada4 分支中调用，用风切变、位温梯度、层界面高度、地表拖曳系数和上一时步 q2 更新湍流动能尺度 q2、动量扩散系数 km、标量扩散系数 kn 和 TKE 扩散系数 kq。这些系数随后供垂直扩散求解使用。

源码注释说明该文件由 Earth 版本改成 Mars 版本，并按 iflag_pbl 选择 MY 2.0、MY 2.0 Fournier、MY 2.5、以及带 TKE 垂直扩散的 MY 2.5 分支。模块还包含 vdif_q2 和 vdif_q2e 两个 TKE 垂直扩散辅助例程；当前主例程只在 iflag_pbl==9 时调用隐式 vdif_q2，vdif_q2e 未在本文件内被调用。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• first、ipas、ric/rifc/b1/kap、firstcall、q2min/q2max/knmin/kmmin、ico2/A/B 均为 SAVE / THREADPRIVATE 状态。
• firstcall 只在首次调用扫描 tracer_mod:noms 查找 co2，并计算 CO2 分子量修正系数；但实际修正 teta=phc*(A*pq+B) 代码被注释，当前总是 teta=phc。
• l0 来自 turb_mod，本例程会按 sqrt(q2) 加权高度积分更新它。
• 当 iflag_pbl 不在 6..10 时，例程调用 abort_physic 中止。

核心逻辑

1. 首次调用扫描 tracer 名称，若存在 co2 则准备 A/B 分子量系数；当前实际稳定度变量直接取 teta=phc。
2. 用 ndt=ceiling(3840./(3699.*24./dt)) 计算 MY 2.5 分支的内部子步数。
3. 校验 iflag_pbl 在 6..10。
4. 计算层厚倒数、层中心间距、风切变平方 m2、位温梯度、Brunt-Vaisala 频率平方 n2、Richardson 数 ri、rif/alpha/sm 和 zz。
5. 初始或 iflag_pbl==6 时迭代 10 次估计混合长度 l、q2 和积分尺度 l0。
6. 每次调用根据当前 q2 更新 l0，再计算界面混合长度 l。
7. iflag_pbl==6：MY 2.0，直接令 q2=l**2*zz。
8. iflag_pbl==7：Fournier 分支，用前一步 km/mpre 估计准静态风切变并更新 q2。
9. iflag_pbl>=8：MY 2.5 分支，在 ndt 个内部子步中半隐式更新 q2；iflag_pbl==9 时额外调用 vdif_q2 做 TKE 垂直扩散。
10. 对非 iflag_pbl==9 分支，根据 km=l*sqrt(q2)*sm、kn=km*alpha、kq=l*sqrt(q2)*0.2 计算扩散系数并设置顶/底边界。
11. 最后对稳定边界层施加 Holtslag-Boville 风格的最小扩散约束，并据此回写 q2/km/kn/kq。

伪代码

yamada4(...):
  on first call:
    find CO2 tracer index and prepare molecular-mass coefficients
  teta = phc
  require 6 <= iflag_pbl <= 10
  compute vertical differences, shear m2, stability n2
  compute rif, alpha, sm, zz
  if first or iflag_pbl == 6:
    initialize l0, iterate l and q2 10 times
  update l0 from sqrt(q2)-weighted height
  compute l at interfaces
  if iflag_pbl == 6:
    q2 = l**2 * zz
  else if iflag_pbl == 7:
    update q2 from Fournier quasi-static branch
  else if iflag_pbl >= 8:
    for each internal substep:
      update q2 with MY 2.5 formula
      clamp q2 to q2min/q2max
      if iflag_pbl == 9:
        compute km/kn/kq and call vdif_q2()
  if iflag_pbl != 9:
    compute km/kn/kq from l, sqrt(q2), sm, alpha
  apply stable-layer minimum K constraint

参与的主题流程

写法特点

• 文件为固定格式 Fortran 风格，使用 statement functions frif/falpha/fsm/fl。
• fl 形式上有参数，但实际表达式直接使用外层 ig/k 和 zlev/q2/n2/l0，阅读时应按内联函数理解。
• iflag_pbl=10 通过合法性检查，并进入 iflag_pbl>=8 分支；但垂直扩散只在 iflag_pbl==9 执行，注释也只说明到 9。
• 诊断数组 rino/smyam/styam/... 位于 if(1.eq.0) 块内，当前不会执行。

复现要点

• iflag_pbl 必须在 6..10；不同值对应不同 q2 时间推进方式。
• zlev/zlay 必须严格单调并避免零层厚；源码直接用差值取倒数。
• q2 会被夹在 q2min=1e-10 与 q2max=1e2 之间，影响强不稳定或数值爆发场景。
• iflag_pbl==9 会对子步后的 q2 调用 vdif_q2；其他分支不会执行该隐式扩散。
• 稳定层最小扩散约束依赖 ustar 和 ngrid：1D 用 kminfact=0.3，3D 用 0.45。

待确认

• iflag_pbl=10 的预期物理含义需确认：它通过检查并走 >=8 的 q2 更新，但不触发 ==9 的 vdif_q2 垂直扩散。
• CO2 分子量修正代码被注释后，pq 和首次扫描出的 ico2/A/B 对当前计算没有实际影响；是否为有意停用需结合版本历史确认。
• knmin/kmmin 定义了 Mars 下限但当前未见用于最终夹值。

相关页面

• vdifc_mod
• vdif_kc
• turb_mod
• calltherm_interface_mod
• config-system