pbl_parameters_mod.F90

路径

LMDZ.MARS\libf\phymars\pbl_parameters_mod.F90

所属目录/模块

libf/phymars

文件定位

pbl_parameters_mod.F90 是行星边界层（PBL）诊断参数计算模块。它从温度、风场和热羽流输出中计算摩擦速度 ustar、摩擦温度 tstar 和 Monin-Obukhov 长度，并基于 Monin-Obukhov 理论在地表层内插值位温和风速。这些是纯诊断量，不影响模型的物理倾向或时间推进。

参考文献：Colaïtis et al. (2013), J. Geophys. Res. Planets, 118, 1468–1487, doi:10.1002/jgre.20104。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• karman（=0.41）和 nu（=0.001）：DATA 初始化 + SAVE + !$OMP THREADPRIVATE（行 111–115），每个 OpenMP 线程独立副本。
• 诊断输出：非 MESOSCALE 模式下，通过 write_output 输出 10 个诊断字段（tke_pbl、rib_pbl、cdn_pbl、fm_pbl、uv、uvtrue、chn_pbl、fh_pbl、B_pbl、zot_pbl、zz1），以及 physiq 中追加的 T_out_*、u_out_*、u_star、teta_star、vvv、vhf。
• 不影响物理倾向：本模块是纯诊断，不修改任何物理状态变量。

核心逻辑

Part I：Richardson/Reynolds/热粗糙度/稳定度系数（行 162–290）

1. 初始化（行 170–184）：清零 ustar、tstar、reynolds、rib、pcdv、pcdh；pz0tcomp 初始为 0.1*pz0；迭代参数 itemax=10、tol_iter=0.01。
2. 稳定度函数参数（行 188–194）：Dyer 参数集（bm=bh=16、alphah=1、betam=betah=5）；ric_colaitis = betah/(betam²)。
3. ATKE vs Colaïtis 切换（行 196–200）：callatke=.true. 时用 ATKE 的 ric，否则用 ric_colaitis。
4. 水浮力修正（行 204–217）：include_waterbuoyancy=.true. 时计算虚拟温度 tsurf_v 和 temp_v，考虑水汽浮力效应；地表有霜时用 watersat 计算饱和混合比。
5. 迭代求 z0t（行 219–290）：对每个网格点：
   • 计算中性拖曳系数 cdn = (karman/ln(z1/z0))²
   • 近地表风速 zu2 = u²+v² + (ln(1+0.7w*+2.3w*²))²（含浮力修正）
   • 迭代循环（最多 10 次，收敛阈值 0.01*z0）：
     • 计算中性热拖曳系数 chn
     • 计算 Bulk Richardson 数 rib（England et al. 1995 公式）
     • 稳定度函数双路径：
       • ATKE 路径（callatke=.true.）：稳定用 sm = max(smmin, cn*(1-rib/ric))，不稳定用 sm = 2/π*(cinf-cn)*atan(-rib/ri0)+cn；Prandtl 数从 Venayagamoorthy & Stretch (2010)
       • Colaïtis 路径：稳定用 fm = ((ric-rib)/ric)²（rib<ric 时），不稳定用 Dyer 公式 fm = sqrt(1-λbm*rib)
     • Reynolds 数 reynolds = karman*sqrt(fm)*sqrt(zu2)*z0/(ln(z1/z0)*ν)
     • 热粗糙度 z0t = z0*exp(-karman*7.3*Re^0.25*Pr^0.5+5*karman)
   • 拖曳系数 pcdv = cdn*fm、pcdh = chn*fh

Part II：ustar/tstar/内插值（行 294–338）

1. 摩擦量（行 299–308）：rib ≥ ric 时 ustar=tstar=0（稳定无湍流）；否则 ustar = sqrt(cdv)*sqrt(zu2)、tstar = -cdh*(pts-ph[:,1])/sqrt(cdv)
2. MO 内插（行 310–326）：
   • zout < z0t：u=0、theta=pts
   • z0t ≤ zout < z0：u=0
   • zout ≥ z0：u = ustar*ln(zout/z0)/(karman*sqrt(fm))、theta = pts + tstar*sqrt(fm)*ln(zout/z0t)/(karman*fh)
3. 无热羽流 + 对流调整特例（行 333–336）：calltherm=.false. && calladj=.true. 时，位温直接取第一层值。
4. 静力修正（行 337）：T_out = theta * exp((zout/z1)*ln(p1/ps))^rcp

Part III：垂直速度方差/热通量廓线（行 341–378）

仅 calltherm=.true. 时执行，按 Spiga et al. (2010) QJRMS：

• x = zout/zmax
• 无量纲热通量 dvhf：x≤0.3 用对数公式，0.3<x≤1 用线性公式，x>1 为 0
• 无量纲速度方差 dvvv：x≤1 用 2.05*x^(2/3)*(1-0.64x)²，x>1 为 0
• vhf = dvhf*hfmax、vvv = dvvv*wstar²

Part IV：诊断输出（行 381–400）

非 MESOSCALE 模式下输出 10+ 个诊断字段。

伪代码

初始化 ustar=tstar=reynolds=rib=pcdv=pcdh=0
z0tcomp = 0.1*z0
设置 Dyer 参数 (bm=bh=16, alphah=1, betam=betah=5)
ric_4interp = callatke ? ric_ATKE : ric_colaitis

if include_waterbuoyancy:
  计算虚拟温度 tsurf_v, temp_v（含水汽浮力）

for each grid point:
  cdn = (karman/ln(z1/z0))²
  zu2 = u²+v² + (ln(1+0.7w*+2.3w*²))²   ! 浮力修正风速
  while (z0t not converged, max 10 iter):
    chn = cdn * ln(z1/z0) / ln(z1/z0t)
    rib = (g/tsurf_v) * sqrt(z1*z0) * (ln(z1/z0))² / ln(z1/z0t) * (Tv-Tsurf)/zu2
    if callatke:
      稳定: sm=max(smmin, cn*(1-rib/ric)), Prandtl=...
      不稳定: sm=2/π*(cinf-cn)*atan(-rib/ri0)+cn, Prandtl=...
    else:
      稳定: fm=((ric-rib)/ric)², fh=same
      不稳定: fm=sqrt(1-λbm*rib), fh=...
    Re = karman*sqrt(fm)*sqrt(zu2)*z0/(ln(z1/z0)*ν)
    z0t = z0*exp(-karman*7.3*Re^0.25*Pr^0.5+5*karman)
  cdv = cdn*fm; cdh = chn*fh

for each z_out:
  ustar = sqrt(cdv)*sqrt(zu2)
  tstar = -cdh*(Ts-T1)/sqrt(cdv)
  if zout < z0t: u=0, theta=Ts
  elif zout < z0: u=0
  else: u=ustar*ln(zout/z0)/(k*sqrt(fm)), theta=Ts+tstar*sqrt(fm)*ln(zout/z0t)/(k*fh)
  T_out = theta * hydrostatic_correction

if calltherm:
  x = zout/zmax
  dvhf = piecewise(x)    ! Spiga et al. 2010
  dvvv = piecewise(x)
  vhf = dvhf*hfmax; vvv = dvvv*wstar²

参与的主题流程

写法特点

• 固定格式风格：虽为 .F90，但使用大写关键字、DATA 语句、列对齐注释等旧式风格。
• !$OMP THREADPRIVATE：karman 和 nu 为线程私有（行 113）。
• #ifndef MESOSCALE：诊断输出仅在非 MESOSCALE 模式下编译（行 388–400）。
• 双稳定度函数路径：callatke 开关切换 ATKE 方案与 Colaïtis (2013) 方案。
• 迭代求 z0t：热粗糙度长度通过 Reynolds 数和 Prandtl 数迭代求解（最多 10 次）。
• 硬编码插值高度：z_out = [3., 2., 1., 0.5, 0.1] m（行 2661 in physiq_mod.F）。
• 行 335 疑点：u_out(:,n)=(sqrt(cdn(:))*sqrt(zu2(ig))/karman)*log(zout/pz0(:)) 中 zu2(ig) 应为 zu2(:)（循环变量 ig 在外层循环已结束，此处用标量索引可能为 bug）。

复现要点

• 本模块是纯诊断，不影响物理倾向或模型状态。
• callatke 开关决定稳定度函数的计算路径：ATKE（湍流动能方案）vs Colaïtis (2013)。
• include_waterbuoyancy 控制是否考虑水汽浮力（虚拟温度修正）。
• calltherm 控制是否计算热羽流廓线（Part III）。
• 迭代求 z0t 时，收敛阈值为 0.01*z0，最多 10 次迭代。

待确认

• 行 335 zu2(ig) 中 ig 来自外层 DO ig=1,ngrid 循环，循环结束后 ig 的值取决于编译器（推断：Fortran 标准规定循环变量在循环结束后值为 ngrid+1，此处 zu2(ig) 越界或取最后一个元素，可能是 bug；待确认）。
• 行 226 IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.) 被注释掉（待确认：是否因 turb_resolved 不再适用或效果不佳）。
• 行 113 karman 和 nu 的 THREADPRIVATE 与 SAVE 共存（待确认：SAVE 在 OpenMP 并行时是否多余）。

相关页面

• vdifc-water-surface-exchange — vdifc 中水汽地表交换段，也涉及 PBL 参数
• water-saturation-helpers — watersat 模块，本文件调用以计算地表饱和混合比
• water-cycle — 水循环主题页
• soilwater — 土壤水求解器，与 PBL 地表交换相关
• vdif_cd_mod — 垂直扩散交换系数；callatke 与 callrichsl 决定地表层交换系数路径
• vdif_kc — 垂直扩散 K 系数方案；在 vdifc fallback 路径中更新 q2/km/kn
• turb_mod — 湍流/PBL 共享状态模块，提供 turb_resolved、q2 和 wstar 等数组语义