rocketduststorm_mod.F90

路径

LMDZ.MARS\libf\phymars\rocketduststorm_mod.F90

所属目录/模块

libf/phymars

文件定位

源码依据（第 12-20 行注释及整体结构）：实现"火箭式尘暴"（Rocket Dust Storm, RDS）参数化方案，用于模拟火星尘暴季节中分离尘埃层（detached dust layers）的形成机制。

核心物理思想（第 188-194 行注释）：

• 风暴尘埃吸收太阳辐射产生加热
• 加热驱动垂直对流（上升气流）
• 上升气流通过绝热冷却平衡辐射加热
• stormdust 示踪物随垂直气流输运
• 在高空或弱风区，stormdust 逐渐"脱轨"（detrain）为背景 dust

该模块是尘埃循环主题页中处理尘暴期间垂直输运的专用机制，与背景尘埃的重力沉降（callsedim_mod.F）形成互补。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 模块级变量：dustliftday(ngrid)（SAVE, ALLOCATABLE, THREADPRIVATE），由 compute_dtau_mod.F90 填充，在本模块中仅分配/释放，不直接读写。
• 运行时参数：coeff_detrainment（来自 callkeys_mod），默认 0.02（conf_phys.F:382-383），可通过 callphys.def 覆盖。
• 硬编码常量：
  • coefmin = 0.025（第 151 行）：最小脱轨分数
  • wmin = 0.25 m/s（第 152 行）：触发完全脱轨的低速阈值
  • wmax = 10. m/s（第 153 行）：垂直速度上限
  • secu = 3.（第 160 行）：子步安全系数，限制尘埃每子步穿越层厚
• 诊断输出：write_output('rds_lapserate', ..., 'K/m', lapserate)（第 530-532 行）。
• 无文件 I/O。

核心逻辑

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

阶段 0：初始化与风暴检测（第 208-253 行）

1. 初始化（第 208-228 行）：将所有输出数组和中间变量清零，detrain(:,:) = 1.（默认完全脱轨）。
2. 提取示踪物（第 232-239 行）：从 pq 复制出 zq_dust_mass/number 和 zq_stormdust_mass/number。
3. 风暴检测（第 244-253 行）：对每列 ig，逐层检查 if (stormdust_mass > dust_mass * 1e-4) 若任一层满足，标记 storm(ig) = .true.，否则跳过该列。

阶段 1：辐射传输计算（第 259-267 行）

调用 callradite 计算含风暴尘埃的辐射加热率 zdtlw1（长波）和 zdtsw1（短波）。同时输出 dsodust、dsords、dsotop、tau_pref_scenario、tau_pref_gcm 等诊断量。

阶段 2：垂直速度计算（第 272-343 行）

对每列 ig，若 storm(ig) 为真：

1. 加热率插值到界面（第 292-318 行）：对风暴尘埃和背景尘埃分别计算层间界面的加热率 zdtlw1_lev(l+1)、zdtsw1_lev(l+1)（风暴）和 zdtlw_lev(l+1)、zdtsw_lev(l+1)（背景）。 使用线性插值：val_lev = (val_l * (zlay_l+1 - zlev_l+1) + val_l+1 * (zlev_l+1 - zlay_l)) / (zlay_l+1 - zlay_l)。

2. 环境温度递减率（第 322-326 行）： zdtvert(l+1) = (T_lev(l+1) - T_lev(l)) / (zlay_l+1 - zlay_l)

3. 额外加热 vs 绝热冷却平衡（第 332-340 行）：

   • 额外加热：deltahr(l) = (zdtlw1_lev + zdtsw1_lev) - (zdtlw_lev + zdtsw_lev)
   • 垂直速度：wrad(l) = -deltahr(l) / (g/cpp + max(zdtvert(l), -0.99*g/cpp)) 分母中 g/cpp 为干绝热递减率，max 防止环境递减率接近绝热时除零。
   • 限幅：wrad(l) = max(-wmax, min(wmax, wrad(l)))，即 |wrad| ≤ 10 m/s。

阶段 3：垂直输运（第 349-446 行）

对每列 ig，若 storm(ig) 为真：

1. 浓度混合比（第 352-359 行）：将风暴 fraction 内的示踪物浓度转换：

   mr_stormdust_mass(l) = dust_mass(l) + stormdust_mass(l) / totstormfract(ig)

   推断：totstormfract 为风暴占据的面积分数，此步将风暴区域内的浓度从平均浓度反推为风暴内部浓度。

2. 子步时间步（第 367-378 行）：

   • 计算最大允许时间步 dtmax：dtmax = min over layers of (Δzlev / (secu * |wrad|))，限制每子步穿越 ≤ 1 层。
   • 子步数：nsubtimestep = int(ptimestep / dtmax)
   • 子步长：subtimestep = ptimestep / nsubtimestep

3. 空气质量通量（第 380-383 行）：

   w(l) = wrad(l) * pplev(l) / (r * ztlev(l)) * subtimestep

   即静力学下 w = ρ * v * Δt = (P / (rT)) * v * Δt。

4. Van Leer 平流（第 396-404 行）：

   • 对 mr_stormdust_mass 和 mr_stormdust_number 分别调用本地 van_leer 子例程。
   • 循环 nsubtimestep 次。

5. 混合比反解（第 411-431 行）：

   • 常规情形（mr_stormdust ≥ mr_dust，第 416-420 行）：

     stormdust = totstormfract * (mr_stormdust - mr_dust)
     dust = mr_dust

   • 异常情形（mr_stormdust < mr_dust，第 422-427 行）：

     dust = (1 - totstormfract) * mr_dust + totstormfract * mr_stormdust
     stormdust = 0

     推断：输运后风暴浓度低于背景，说明风暴"耗尽"，将风暴质量按比例合并回背景。

6. 倾向计算（第 435-446 行）：

   dqvl_stormdust = (zq_stormdust - pq_stormdust) / ptimestep
   dqvl_dust = (zq_dust - pq_dust) / ptimestep

阶段 4：脱轨（第 449-482 行）

1. 脱轨系数（第 453-471 行）：

   • 若 |wrad| < wmin（0.25 m/s）或 dust_mass > 10000 * stormdust_mass：detrain = 1.（完全脱轨）
   • 若 wmin ≤ |wrad| ≤ wmax：二次多项式

     detrain = coeff_detrainment * ((1-coefmin)/(wmin-wmax)² * (|wrad|-wmax)² + coefmin)

   • 若 |wrad| > wmax：detrain = coefmin（最小脱轨，保持高速上升）

2. 脱轨倾向（第 475-482 行）：

   dqdet_stormdust = -detrain * stormdust_mass / ptimestep

阶段 5：最终倾向（第 484-498 行）

对每列每层：

pdqrds(stormdust) = dqdet_stormdust + dqvl_stormdust
pdqrds(dust) = -dqdet_stormdust + dqvl_dust

即 stormdust 减少（脱轨 + 输运），dust 增加（脱轨贡献，输运也可能改变背景）。

伪代码

subroutine rocketduststorm(..., pq, pdqrds, wrad, ...):

  #--- 阶段 0: 初始化与风暴检测 ---
  storm(:) = false
  for ig=1..ngrid:
    for l=1..nlayer:
      if stormdust_mass(ig,l) > dust_mass(ig,l) * 1e-4:
        storm(ig) = true; break

  #--- 阶段 1: 辐射传输 ---
  call callradite(..., pq, ..., zdtlw1, zdtsw1, dsodust, dsords, dsotop, ...)

  #--- 阶段 2: 垂直速度 ---
  for ig where storm(ig):
    for l=1..nlayer:
      插值加热率到界面：zdtlw1_lev, zdtsw1_lev（风暴）；zdtlw_lev, zdtsw_lev（背景）
      deltahr(l) = (zdtlw1_lev + zdtsw1_lev) - (zdtlw_lev + zdtsw_lev)
      wrad(l) = -deltahr(l) / (g/cpp + max(zdtvert(l), -0.99*g/cpp))
      wrad(l) = clamp(wrad, -wmax, wmax)

  #--- 阶段 3: 垂直输运 ---
  for ig where storm(ig):
    mr_stormdust = dust + stormdust / totstormfract    # 风暴内部浓度
    dtmax = min over layers of (Δz / (secu * |wrad|))  # 限制穿越层数
    nsubtimestep = int(ptimestep / dtmax)
    subtimestep = ptimestep / nsubtimestep
    for tsub=1..nsubtimestep:
      w(l) = wrad(l) * pplev(l) / (r * ztlev(l)) * subtimestep  # 空气质量通量
      call van_leer(stormdust_mass, w, wq)
      call van_leer(stormdust_number, w, wq)
    # 反解混合比
    if mr_stormdust ≥ mr_dust:
      stormdust = totstormfract * (mr_stormdust - mr_dust)
    else:
      dust = (1-totstormfract)*mr_dust + totstormfract*mr_stormdust; stormdust=0
    dqvl = (zq_after - pq_before) / ptimestep

  #--- 阶段 4: 脱轨 ---
  for ig, l:
    if |wrad| < wmin or dust > 10000*stormdust:
      detrain = 1
    elif |wrad| ∈ [wmin, wmax]:
      detrain = coeff_detrainment * ((|wrad|-wmax)² * (1-coefmin)/(wmin-wmax)² + coefmin)
    else:
      detrain = coefmin
    dqdet_stormdust = -detrain * stormdust / ptimestep

  #--- 阶段 5: 最终倾向 ---
  pdqrds(stormdust) = dqdet_stormdust + dqvl_stormdust
  pdqrds(dust) = -dqdet_stormdust + dqvl_dust

参与的主题流程

写法特点

• 模块级变量 dustliftday 的跨模块共享：本模块分配但不填充，由 compute_dtau_mod.F90 计算，physiq_mod.F 消费。这是 Fortran 90 模块作为"共享状态容器"的典型模式。
• 本地 van_leer 与 vlz_fi.F 的复制（第 544 行注释明确）：本模块的 van_leer 是 vlz_fi 的简化版本（单列 vs 多列），但包含向上输运分支（第 601-641 行），而 vlz_fi.F 通过 goto 99 跳过了向上分支。推断原因：火箭式尘暴需要向上输运，而重力沉降只需向下。
• goto 控制流：van_leer 中保留了 goto 77（第 622、628 行）和 goto 88（第 662、667 行），与 vlz_fi.F 一致。
• 子步循环：为保持质量守恒，Van Leer 输运在多个子步内执行（第 396-404 行），子步数由最大垂直速度决定。
• totstormfract 的面积分数处理（第 355-358 行）：将格点平均浓度转换为风暴内部浓度，输运后再转换回格点平均。这是亚格点尺度参数化的标准做法。
• 脱轨系数的二次多项式（第 461-464 行）：在 wmin 和 wmax 之间平滑过渡，避免不连续。coeff_detrainment 为全局缩放因子。
• 注释掉的非负保护（第 500-525 行）：原计划防止负值，但被注释。推断：Van Leer 格式本身有非负保护（第 687-689 行），外层保护冗余。

复现要点

• 风暴检测阈值 1e-4（第 248 行）：stormdust_mass > dust_mass * 1e-4。这是硬编码的经验阈值，重新实现时必须保留，否则可能误触发或漏触发风暴。
• 垂直速度公式（第 335-336 行）：wrad = -deltahr / (g/cpp + max(zdtvert, -0.99*g/cpp))。分母中的 -0.99*g/cpp 是关键保护：若环境递减率接近绝热，分母趋零，max 防止发散。重新实现时必须保留此数值保护。
• totstormfract 的面积分数转换（第 355-358 行及第 416-427 行）：风暴内部浓度 = 格点平均 / totstormfract；输运后反解时，若 mr_stormdust < mr_dust，说明风暴"耗尽"，需特殊处理。重新实现时必须保留此逻辑分支。
• 子步时间步（第 367-378 行）：dtmax = min over layers of (Δz / (secu * |wrad|))，secu=3 为安全系数。重新实现时必须保留，否则单步穿越多层会导致 Van Leer 格式不稳定。
• 脱轨系数公式（第 461-464 行）：二次多项式 (1-coefmin)/(wmin-wmax)² * (|wrad|-wmax)² + coefmin。系数 coeff_detrainment（默认 0.02）为全局缩放，重新实现时需从 callphys.def 读取。
• van_leer 本地副本包含向上输运：与 vlz_fi.F 不同（后者 goto 99 跳过向上分支）。重新实现时若使用 vlz_fi 模块接口，需确认其支持向上输运（当前版本不支持，第 124 行 goto 99）。
• 复现风险：dustliftday 由 compute_dtau_mod 填充，本模块仅分配。若调用顺序错误（先 rocketduststorm 后 compute_dtau），dustliftday 未初始化，可能导致不确定行为。
• 复现风险：IRtoVIScoef 声明为 inout 但本例程未修改（第 86-87 行注释说明）。重新实现时可简化为 in，但需确认 callradite 接口兼容。

待确认

• 待确认：totstormfract 的计算逻辑（来源未在本文件中追踪，推断由上游 physiq_mod.F 或更高层传入）。
• 待确认：dustliftday 在 physiq_mod.F 中的具体使用方式（本文件仅分配，未追踪 compute_dtau_mod 的填充时序）。
• 待确认：van_leer 本地副本与 vlz_fi.F 的版本差异（第 544 行注释说"copied"，但向上分支是否完全一致未逐行比对）。
• 待确认：注释掉的非负保护（第 500-525 行）是否曾被启用，以及为何被注释（推断为冗余，但未找到历史记录）。

相关页面

• vlz_fi.md：本模块的 van_leer 从 vlz_fi.F 复制，但包含向上输运分支。
• callsedim_mod.md：背景尘埃的重力沉降，与本模块的火箭式尘暴垂直输运形成互补。
• newsedim_mod.md：单示踪物沉降核心，调用 vlz_fi（无向上输运）。
• topmons_mod.md：山顶地形尘流输运，同样有独立 Van Leer 实现。
• dust-cycle.md：尘埃循环主题页，本模块为其中的火箭式尘暴分支。