vdif_cd_mod.F90

路径

LMDZ.MARS\libf\phymars\vdif_cd_mod.F90

所属目录/模块

libf/phymars

文件定位

vdif_cd_mod.F90 是 LMDZ.MARS 的地表层交换系数计算模块。它提供唯一子程序 vdif_cd，根据 Monin-Obukhov 相似理论计算地表动量拖曳系数 Cd 和热量拖曳系数 Ch。callrichsl 控制是否进入 Richardson/Colaïtis 地表层路径：.false. 时直接使用简单中性对数公式，.true. 时迭代热粗糙度长度 z0t（Colaïtis et al. 2013）；在迭代路径内部，callatke 再选择 ATKE 或 D.E. England/Dyer 稳定性函数。被 vdifc_mod.F 调用于地表通量计算。

定义的符号

模块与子程序

关键局部变量

模块级常量

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 模块常量 karman, nu：SAVE + !$OMP THREADPRIVATE，仅初始化一次。
• 诊断输出：当 write_outputs=.true. 时，通过 write_output 写 7 个 XIOS 字段（rib_dry_vdif_cd, rib_vdif_cd, fm_vdif_cd, fh_vdif_cd, z0t_vdif_cd, z0_vdif_cd, Reynolds_vdif_cd），均取 iflat（平坦地形）切片。
• 无文件 I/O，无全局状态修改。

核心逻辑

1. 初始化（行 144-167）：设定迭代参数（itemax=10, tol_iter=0.01），Dyer 参数化常量（bm=bh=16, betam=betah=5, alphah=1），计算 lambda 和 ric_colaitis。
2. 水汽浮力修正（行 171-186）：若 include_waterbuoyancy，计算虚拟温度 tsurf_v（霜面用 qsat，非霜面用实际 pqvap）和 temp_v；否则直接用物理温度。
3. 简单路径（callrichsl=.false.，行 188-199）：经典对数风廓线公式 Cd = (κ/ln(1+z₁/z₀))²，Cd=Ch。
4. Colaïtis 迭代路径（callrichsl=.true.，行 200-279）：
   • 计算中性拖曳系数 cdn, chn（行 207-213）。
   • 迭代循环（最多 10 次）： a. 计算 zu2（近地表风速 + 对流浮力修正，行 215）。 b. 计算体 Richardson 数 rib（England et al. 1995 公式，行 218-219）。 c. 计算稳定性函数 fm, fh：
     • callatke=.true.：ATKE 方案，稳定用 cn*(1-Ri/ric)，不稳定用 atan 反正切（行 228-243）。
     • callatke=.false.：D.E. England/Dyer 路径；0<Ri<ric_colaitis 时稳定用 (ric-Ri)²/ric²，Ri>=ric_colaitis 时源码设 fm=fh=1，不稳定时用 sqrt(1-λ*bm*Ri) 等 Dyer 形式（行 245-263）。 d. 重算 Reynolds 数和 z0t（Zilitinkevich 型公式，行 266-267），检查收敛。
   • 最终系数：pcdv = cdn * fm, pcdh = chn * fh（行 273-274）。
5. 诊断输出（行 281-297）：写 rib_dry, rib, fm, fh, z0t, z0, Reynolds 到 XIOS。

伪代码

vdif_cd:
  初始化 Dyer 常量, 迭代参数

  if include_waterbuoyancy:
    计算虚拟温度 tsurf_v（霜面用 qsat, 非霜面用 pqvap）
    计算 temp_v（第一层虚拟位温）

  if .not. callrichsl:
    ! 简单对数风廓线
    Cd = Ch = (κ / ln(1 + z1/z0))²
  else:
    ! Colaïtis et al. 2013 迭代
    for each (ig, islope):
      cdn = (κ / ln(z1/z0))²
      迭代（最多 10 次）:
        chn = cdn * ln(z1/z0) / ln(z1/z0t)
        zu2 = u² + v² + (ln(1 + 0.7*wstar + 2.3*wstar²))²
        rib = f(Ri, 温度差, zu2)
        if callatke:
          ATKE 稳定性函数（sm, prandtl → fm, fh）
        else:
          Dyer/England 稳定性函数（0<Ri<Ric 用二次式；Ri>=Ric 时 fm=fh=1；Ri<=0 用不稳定 Dyer 形式）
        Re = κ * sqrt(fm) * sqrt(zu2) * z0 / (ln(z1/z0) * ν)
        z0t_new = z0 * exp(-κ*7.3*Re^0.25*Pr^0.5 + 5*κ)
        检查收敛
      pcdv = cdn * fm
      pcdh = chn * fh

  if write_outputs:
    写 7 个 XIOS 诊断字段（取 iflat 切片）

参与的主题流程

写法特点

• 两级路径开关：callrichsl 首先控制是否进入 Colaïtis 迭代；只有进入迭代后，callatke 才控制使用 ATKE 还是 D.E. England/Dyer 稳定性函数。
• 迭代热粗糙度长度：z0t 通过 Zilitinkevich 型公式迭代至收敛（容差 0.01*z0），最多 10 步。
• 水汽浮力修正：include_waterbuoyancy 启用时，虚拟温度考虑水汽分子量差异和地表霜。
• 模块常量 SAVE/THREADPRIVATE：karman, nu 仅初始化一次，OpenMP 安全。
• 诊断输出取 iflat 切片：所有 write_output 调用只输出平坦地形的诊断。

复现要点

• callrichsl=.true. 时才进入迭代路径；callrichsl=.false. 时直接用简单对数公式。
• callatke 仅在 callrichsl=.true. 时生效，控制稳定性函数类型。
• callatke=.false. 且 Ri>=ric_colaitis 时，源码设 fm=fh=1；复现时不能把稳定二次式外推到该区间。
• include_waterbuoyancy 影响虚拟温度计算，进而影响 Richardson 数。
• 迭代收敛容差为 tol_iter * pz0(ig)，即相对粗糙度的 1%。

待确认

• karman 和 nu 的 DATA 初始化 + SAVE + THREADPRIVATE 组合在 OpenMP 并行时是否会有初始化竞争（首次进入时 DATA 赋值）。
• zu2(ig) 在行 215-216 中只用了 ig 索引但 pu/pv 是二维数组——这里隐式取 pu(ig,1) 和 pv(ig,1)，即第一层风速，是否为设计意图。
• ite 声明为 REAL（行 121）但用作迭代计数器（ite = ite+1），应为 INTEGER（不影响功能但不规范）。

相关页面

• turb_mod — 提供 turb_resolved 与 wstar 等湍流/PBL 共享状态语义。
• pbl_parameters_mod — 边界层参数诊断模块
• comslope_mod — 坡面几何共享参数
• vdifc-water-surface-exchange — vdifc_mod.F 中地表交换段，记录 vdif_cd 调用方和 Cd/Ch 后续用于通量计算的上下文
• water-cycle — 水循环主题页