topmons_mod.F90

路径

LMDZ.MARS\libf\phymars\topmons_mod.F90

所属目录/模块

libf/phymars

文件定位

源码依据（第 8-16 行注释及整体结构）：实现"山顶地形尘流"（Mountain Top Dust Flows）参数化方案，模拟火星上亚格点尺度山脉的坡风驱动尘埃抬升机制。

核心物理思想（第 211-214 行注释）：

• 亚格点尺度山脉在白天产生上坡气流（热力坡风）
• 山顶尖端的会聚气流将边界层（PBL）中的背景尘埃**裹挟（entrain）**至高空
• 被裹挟的尘埃作为 topdust 示踪物在山脉上方垂直输运
• topdust 在速度衰减或夜间时**脱轨（detrain）**回背景 dust

该模块与 rocketduststorm_mod 形成尘埃循环中两种集中垂直输运机制的互补：前者由辐射加热驱动（尘暴），后者由地形坡风驱动（山脉）。

定义的符号

依赖的模块

调用的关键例程

输入

topmons 主例程

topmons_setup

输出

topmons 主例程

topmons_setup

共享状态与副作用

• firstcall（SAVE, THREADPRIVATE，第 115 行）：声明但未在主逻辑中使用（推断为保留接口）。
• surfdat_h 全局修改：topmons_setup 写入 contains_mons、alpha_hmons、hsummit，供后续 topmons 和其他模块读取。
• 硬编码常量：
  • wmax = 10. m/s（第 164 行）：辐射驱动速度上限
  • secu = 3.（第 165 行）：子步安全系数
  • k0 = 0.25（第 166 行）：山顶会聚速度系数（对应最大约 3 m/s）
  • k1 = 5e-4（第 167 行）：动量衰减系数
  • k2 = 5e-3（第 168 行）：浮力加速度系数
  • mu0lim = 0.001（第 122 行）：白天判定阈值
• 内置火山列表（第 984-1001 行）：19 座火星火山的经纬度（Olympus Mons、Ascraeus Mons 等），ntop_max = 19。
• 1D 特殊值（第 1091 行）：hmax = 23162.1 m（推断为 Olympus Mons 高度），hsummit = 14000 m。
• 诊断输出：当前全部注释掉（第 647-695 行），包含 wup_top、wfin_top、entr 等 20+ 诊断量。

核心逻辑

按执行顺序（行号为源码行号）：

阶段 0：初始化（第 225-276 行）

1. 全部输出和中间变量清零。
2. 温度更新（第 265 行）：zt = pt + pdt * ptimestep（用上游倾向推进）。
3. 示踪物更新（第 271 行）：zq0 = pq + pdq * ptimestep（含火箭式尘暴后的背景尘埃）。

阶段 1：辐射加热（第 279-369 行）

1. 辐射传输（第 291 行）：调用 callradite 计算含 topdust 的加热率 zdtlw1/zdtsw1。
2. 对每列含山格点（mu0 > mu0lim 且 contains_mons(ig) 为真）：
   • 山脉温度廓线（第 310 行）：调用 t_topmons 计算山脉上方温度 t_top 和环境温度 t_env。
   • 加热率插值（第 322-346 行）：插值到界面，同 rocketduststorm。
   • 辐射垂直速度（第 356-368 行）： wrad = -deltahr / (g/cpp + max(zdtvert, -0.99*g/cpp)) 限幅 |wrad| ≤ wmax；向下分量强制置零（第 365-367 行：if wrad > 0 → wrad = 0），因为山脉坡风只产生向上输运。

阶段 2：垂直速度合成（第 372-474 行）

对每列含山格点：

1. 浮力驱动速度（第 382-429 行）：

   • 若 dt_top > 0（山脉上方比环境暖）：
     • 会聚速度（第 384 行）：w0 = -k0 * g * sqrt(dt_top / T_env(lsummit)) 负值表示向上。
     • 若 |w0| > |wrad|：从山顶层逐层向上推进，使用热羽流模型：

       wup²(l+1) = (1 - 2*k1*Δz) * wup²(l) + 2*k2*g*Δz * (newzt(l) - t_env(l)) / t_env(l)

       其中 k1 为动量衰减，k2 为浮力加速度。
     • 若 wup² > 0：wup = -sqrt(wup2)，wfin = wup。
     • 若 |wup| < |wrad| 或 wup² ≤ 0：退化到纯辐射方案 wfin = wrad。
   • 若 dt_top ≤ 0（山脉上方比环境冷）：退化到纯辐射方案。

2. 最大速度层 lwmax（第 447-450 行）：lwmax = minloc(wup) — 最大向上速度所在层，作为 PBL 尘埃注入层。

3. 子步时间步（第 455-466 行）：同 rocketduststorm，dtmax = min(Δz / (secu * |wfin|))。

4. 空气质量通量（第 468-471 行）：w(l) = wfin(l) * pplev(l) / (r * ztlev(l)) * subtimestep。

阶段 3：PBL 裹挟 + 垂直输运（第 476-586 行）

1. 密度与空气质量（第 480-490 行）：rho = pplay / (r * pt)；masse = (pplev(l) - pplev(l+1)) / g。

2. PBL 总空气质量（第 494-500 行）：从地表到 lwmax-1 层的气柱质量。

3. 裹挟通量（第 516 行）：

   entr = alpha_hmons * wup(lwmax) * rhobarz(lwmax) / masse_pbl

   负值表示从 PBL 抽出空气质量。保护：|entr| ≤ 1/ptimestep（第 518-520 行）。

4. 每个子步（第 525-553 行）：

   • PBL 内混合比（第 527-534 行）：加权平均 qbar_pbl。
   • PBL 尘埃抽出（第 538-544 行）：指数衰减模型 dq = (1 - exp(entr*subtimestep)) * q * Δp/g，每层逐层扣除。
   • Van Leer 平流（第 549-552 行）：从 lwmax 到 nlayer，调用本地 van_leer；额外传入 masse_pbl、dqm_pbl、alpha_hmons、qbar_pbl 作为 PBL 注入边界条件。

5. 混合比反解（第 558-572 行）：同 rocketduststorm，用 alpha_hmons 做面积分数转换。

6. 倾向计算（第 576-583 行）：dqvl = (zq_after - zq0_before) / ptimestep。

阶段 4：脱轨（第 588-624 行）

1. 脱轨系数（第 599-601 行，白天且含山时）：

   coefdetrain = (rhobarz(l+1)*|wfin(l+1)| - rhobarz(l)*|wfin(l)|) / masse(l)

   即质量通量散度除以层质量：散度为正时（质量减少）产生脱轨。

2. 指数脱轨（第 603-612 行）： dqdet = -(1 - exp(coefdetrain * ptimestep)) * topdust / ptimestep 保护：不超过可用 topdust 量（第 608-612 行）。

3. 夜间/非含山格点（第 618-622 行）：完全脱轨，dqdet = -topdust / ptimestep。

阶段 5：最终倾向（第 626-641 行）

pdqtop(topdust) = dqdet_topdust + dqvl_topdust
pdqtop(dust) = -dqdet_topdust + dqvl_dust

t_topmons 子例程（第 705-785 行）

1. 山顶层定位（第 744-747 行）：hsummit ≥ zlay(lsummit) 的最高层。
2. 山脉上方温度（第 749-753 行）：山顶层设为地表温度 zt(1)，向上按绝热递减率递减。
3. 环境温度映射（第 756-774 行）：从 lmons（山坡高度对应层）开始，逐层向上映射环境温度。
4. t_top 下界保护（第 776-781 行）：t_top 不低于 t_env。
5. 温差输出（第 783 行）：dt_top = t_top(lsummit) - t_env(lsummit)。

topmons_setup 子例程（第 954-1100 行）

1. 19 座火山经纬度（第 984-1001 行）：硬编码列表，Olympus Mons 首位（lon=-134, lat=18.4）。
2. 含山格点检测（第 1020-1042 行）：逐格点检查 boundslon/lat 是否包络任一火山坐标。
3. alpha_hmons 计算（第 1062-1084 行）： alpha_hmons = 0.5 * (hmons - hmin) / (hmax - hmin)，范围 [0, 0.5]。 alpha_hmons = 0 时禁用（山脉已被格点分辨）。
4. hsummit 计算（第 1087 行）：hsummit = summit - phisfi/g。

van_leer 子例程（第 792-951 行）

基于 vlz_fi.F 的本地副本（第 790 行注释明确），关键差异：

• 从 lwmax 层开始（非从底层），第 823、828、838 行。
• wq(lwmax) = -dqm_pbl / alpha_hmons（第 848 行）：PBL 注入边界条件，从 lwmax 层注入尘埃通量。
• 扩展格式（多层穿越）被注释掉（第 867-893 行），仅保留常规格式。

伪代码

subroutine topmons(..., pq, pdqtop, wfin, ...):

  #--- 阶段 0: 初始化 ---
  zt = pt + pdt * ptimestep
  zq0 = pq + pdq * ptimestep

  #--- 阶段 1: 辐射加热 ---
  call callradite(..., zq, ..., zdtlw1, zdtsw1, ...)
  for ig where daytime and contains_mons(ig):
    call t_topmons(...) → t_top, dt_top, t_env, lsummit, lmons
    插值加热率到界面
    wrad = -deltahr / (g/cpp + max(zdtvert, -0.99*g/cpp))
    wrad = clamp(wrad, -wmax, wmax)
    if wrad > 0: wrad = 0   # 仅保留向上分量

  #--- 阶段 2: 垂直速度合成 ---
  for ig where daytime and contains_mons(ig):
    if dt_top > 0:
      w0 = -k0 * g * sqrt(dt_top / T_env(lsummit))   # 会聚速度
      if |w0| > |wrad|:
        从山顶向上逐层: wup²(l+1) = (1-2k1Δz)*wup²(l) + 2k2*g*Δz*(newzt-T_env)/T_env
        wup = -sqrt(wup2); wfin = wup
      else:
        wfin = wrad   # 退化为辐射方案
    else:
      wfin = wrad
    lwmax = minloc(wup)   # 最大速度层 = PBL 注入层
    dtmax = min(Δz / (secu*|wfin|)); nsubtimestep = ceil(ptimestep/dtmax)
    w(l) = wfin(l) * pplev(l) / (r*ztlev(l)) * subtimestep

  #--- 阶段 3: PBL 裹挟 + Van Leer 输运 ---
  entr = alpha_hmons * wup(lwmax) * rhobarz(lwmax) / masse_pbl
  for tsub=1..nsubtimestep:
    qbar_pbl = Σ(dust * Δp/g) / masse_pbl    # PBL 平均浓度
    for l=1..lwmax-1:                         # PBL 内逐层扣除
      dqm = (1 - exp(entr*subtimestep)) * dust(l) * Δp/g
      dust(l) = dust(l) - dqm
    call van_leer(topdust_mass, lwmax, ..., dqm_pbl, qbar_pbl)  # lwmax→top
    # van_leer 内: wq(lwmax) = -dqm_pbl / alpha_hmons
  # 反解混合比
  topdust = alpha_hmons * (mr_topdust - mr_dust)
  dqvl = (zq_after - zq0_before) / ptimestep

  #--- 阶段 4: 脱轨 ---
  if daytime and contains_mons and topdust > 0.01*dust:
    coefdetrain = (rho*|w|(l+1) - rho*|w|(l)) / masse(l)  # 质量通量散度
    dqdet = -(1 - exp(coefdetrain*ptimestep)) * topdust / ptimestep
  else:
    dqdet = -topdust / ptimestep   # 完全脱轨

  #--- 阶段 5: 最终倾向 ---
  pdqtop(topdust) = dqdet + dqvl_topdust
  pdqtop(dust) = -dqdet + dqvl_dust

subroutine t_topmons(nlayer, summit, hsummit, hmons, zt, zlay, t_top, dt_top, t_env, lsummit, lmons):
  lsummit = 最高 zlay < hsummit 的层
  t_top(lsummit) = zt(1)   # 山顶温度 = 地表温度
  向上: t_top(l+1) = t_top(l) - (zlay(l+1)-zlay(l)) * g/cpp   # 绝热递减
  lmons = 最高 zlay < hmons 的层
  t_env(lsummit) = zt(lmons)   # 环境温度映射
  dt_top = t_top(lsummit) - t_env(lsummit)

subroutine topmons_setup(ngrid):
  for each of 19 predefined volcanoes:
    检查格点边界是否包络 → contains_mons(ig) = .true.
  alpha_hmons = 0.5 * (hmons - hmin) / (hmax - hmin)
  hsummit = summit - phisfi/g

参与的主题流程

写法特点

• 19 座火山硬编码（第 984-1001 行）：包含 Olympus Mons、Ascraeus Mons、Arsia Mons 等，经纬度来自 Mars Global Surveyor 数据。添加新火山需修改源码。
• 坡风物理模型（第 384、406-409 行）：
  • w0 = -k0 * g * sqrt(dt_top / T_env) — 基于温差的会聚速度
  • wup² = (1-2k1Δz)*wup² + 2k2*g*Δz*(T_new-T_env)/T_env — 热羽流方程，含动量衰减和浮力加速度
  • 参数 k0=0.25, k1=5e-4, k2=5e-3（第 166-168 行注释）：调试为最大约 3 m/s。
• PBL 指数裹挟（第 540-544 行）：dq = (1 - exp(entr*subtimestep)) * q * Δp/g，比线性近似更保守（防止 PBL 内浓度变负）。
• qbar_number_pbl 赋值错误（第 534 行）：qbar_number_pbl(ig) = masse_pbl_dust_mass(ig) / masse_pbl(ig) — 此处使用了 masse_pbl_dust_mass 而非 masse_pbl_dust_number，推断为代码 bug，导致 number 分量的 PBL 平均混合比被错误地设为 mass 的值。
• Van Leer 中 wq(lwmax) 边界条件（第 848 行）：wq(lwmax) = -dqm_pbl / alpha_hmons，从 PBL 注入的尘埃通量被注入到 lwmax 层；除以 alpha_hmons 是因为被输运变量已转换为山脉内部浓度。
• 辐射驱动速度向下分量被截断（第 365-367 行）：wrad > 0 → wrad = 0，使辐射加热仅产生向上速度，这与 rocketduststorm_mod 不同（后者保留双向分量）。
• 夜间完全脱轨（第 618-622 行）：太阳落山后所有 topdust 立即脱轨回背景尘埃，不保留残余。
• van_leer 注释掉的扩展格式（第 867-893 行）：多层穿越格式被完全注释，仅保留常规格式，因为子步时间步已保证单步不穿越多层。

复现要点

• 19 座火山列表（第 984-1001 行）：必须保持一致。任何遗漏或坐标偏差都会导致 contains_mons 判定不同，进而使 topdust 方案在错误格点触发或不触发。
• alpha_hmons 公式（第 1065 行）：0.5 * (hmons - hmin) / (hmax - hmin)，因子 0.5 是关键——范围被限制在 [0, 0.5]。重新实现时不可改为 1.0。
• k0, k1, k2 参数（第 166-168 行）：分别控制山顶会聚速度、动量衰减和浮力加速度。注释标注了不同速度上限的对应值（max 10 m/s vs max 3 m/s），当前使用 max 3 m/s 的配置。
• wq(lwmax) = -dqm_pbl / alpha_hmons（第 848 行）：PBL 注入边界条件中的 alpha_hmons 除法是必须的，否则输运变量量纲不匹配。
• PBL 指数裹挟（第 540 行）：1 - exp(entr*subtimestep) 而非线性 entr*subtimestep，确保不产生负浓度。
• 夜间完全脱轨（第 618-622 行）：日落后 topdust 在下一个时间步全部转回 dust，不能保留。
• dt_top > 0 才产生坡风（第 382 行）：温度差为负时（山脉上方比环境冷），无浮力驱动。
• 复现风险：第 534 行 qbar_number_pbl(ig) = masse_pbl_dust_mass(ig)/masse_pbl(ig) 疑似 bug——应为 masse_pbl_dust_number。重新实现时需决定是按源码复制（保留 bug）还是修正。
• 复现风险：van_leer 的扩展格式（多层穿越）被注释掉。若子步数不足（例如 nsubtimestep 计算有误），单步穿越多层会因缺少扩展格式而产生不可预测结果。
• 复现风险：topmons_setup 在介观尺度模式（MESOSCALE 宏）下 hmin=hmax=0，导致全部 contains_mons = .false.，即 topflows 方案被完全禁用。重新实现时需注意此编译选项。

待确认

• 待确认：第 534 行 qbar_number_pbl(ig) = masse_pbl_dust_mass(ig)/masse_pbl(ig) 是否为 bug（应为 masse_pbl_dust_number）。此 bug 影响 number 分量的 PBL 平均浓度的计算。
• 待确认：firstcall（第 115 行）的原始用途——在当前主逻辑中未被使用，推断为残留代码或保留接口。
• 待确认：1D 模式下 hmax = 23162.1 m 的来源——推断为 Olympus Mons 海拔，但未找到文档确认。
• 待确认：k0=0.25, k1=5e-4, k2=5e-3 的物理推导来源——注释仅提供"max 3 m/s"的目标速度信息，未引用原始文献。
• 待确认：t_topmons 中 t_env 的映射逻辑（第 756-774 行）——环境温度取自更远层的温度再映射回 lsummit，该映射方法的物理依据和文献引用不明。

相关页面

• rocketduststorm_mod.md：辐射驱动尘暴方案，与 topmons 互补——前者不需要山脉，后者需要亚格点地形。
• vlz_fi.md：Van Leer 垂直输运格式的本源，本模块的 van_leer 从此复制但增加了 PBL 注入边界条件。
• callsedim_mod.md：重力沉降调度器，与 topmons 形成"沉降 vs 集中抬升"的互补。
• newsedim_mod.md：单示踪物沉降核心。
• updaterad.md：尘埃/水冰有效半径反演模块。
• dust-cycle.md：尘埃循环主题页。