improvedclouds_mod.F90

路径

LMDZ.MARS\libf\phymars\improvedclouds_mod.F90

所属目录/模块

libf/phymars

文件定位

完整水冰云微物理核心模块。当 microphys=.true. 时由 watercloud_mod 调用，取代简化方案 simpleclouds_mod。实现：（1）尘埃 CCN 上的水冰异质成核（连续率→离散概率）；（2）冰晶质量增长/升华（阻力增长方案 + 隐式格式）；（3）冰完全升华时的尘核释放；（4）可选 HDO 同位素分馏；（5）自适应子时间步（adapt_imicro，本文件内）。CO2 云版 improvedco2clouds_mod 即基于本文件改编。

定义的符号

依赖的模块

调用的关键例程

输入

输出

注：输出为微物理后的绝对值（已加上之前倾向并经子步积分），由 watercloud_mod 转回倾向供后续沉降/输送使用。

共享状态与副作用

• SAVE + THREADPRIVATE：firstcall（初始化标志）、rb_cld（尺度 bin 对数边界）、sigma_ice（冰/CCN 分布方差）、test_flag（调试开关）。
• 写 module 变量：vo1（水分子体积，行 182）、rad_cld（半径网格，行 148-152）——这两个来自 microphys_h，在 firstcall 内填充。rad_cld 在 firstcall 内被本例程写入，属跨例程共享状态。
• firstcall（行 130-192）：构造尺度网格、写 vo1、sigma_ice，并 write(*,*) 打印 bin 信息。
• 诊断输出：write_output 写 zpotcond、count_micro；test_flag=.true. 时额外打印成核格点占比。
• 无文件读取（与 CO2 版的 Meteo_flux_Plane.dat 不同——本文件无流星体通量源）。

核心逻辑

1. firstcall 初始化（行 130-192）：

   • 体积比 vrat_cld = exp(log(rmax/rmin)/(nbin_cld-1)*3)，构造对数等比半径网格 rad_cld 与边界 rb_cld（行 143-161）。
   • rb_cld(:) = log(rb_cld(:))——存对数值供后续 derf 直接用（行 174）。
   • vo1 = mh2o/rho_ice（水分子体积），sigma_ice = sqrt(log(1+nuice_sed))。

2. 初始化倾向合并（行 197-229）：zq = (上游 tracer) + pdq*ptimestep（dust 仅 scavenging 时加；ccn、h2o 总是加；hdo 仅 hdo 时加）；zt = pt + pdt*ptimestep；subpdtcloud=0；微小值钳到 1e-30。

3. 算饱和（行 234-238）：dev2 = 1/(sqrt(2)*sigma_ice)；watersat 全网格算 zqsat；zpotcond = h2o_vap - zqsat（可凝结量，供自适应子步）。

4. 主循环（行 245-450，对每个 l、ig）：

   • 自适应子步（行 248）：cloud_adapt_ts 时 adapt_imicro 按 zpotcond 定 zimicro。
   • 子步 while 循环（行 252-448，直到 ending_ts）：
     • 单点重算 zqsat、水汽分压 ph2o = h2o_vap*(mmean/18)*pplay、饱和比 satu、microtimestep = ptimestep/zimicro（行 253-258）。
     • 保存 zq0；若 spenttime+microtimestep >= ptimestep 则截断为剩余时间并置末步（行 264-268）。
     • 成核（satu >= 1，行 273-309）：updaterccn 算 rdust；尘埃 log-normal 经 derf 展开为 binned n_aer/m_aer；nuclea 算 rate；dN/dM = Σ aer*(exp(-rate*dt)-1)（负值，即被激活转为冰核的量）；从 dust 减、向 ccn 加。
     • 冰增长（ccn_number*tauscaling >= 1，行 317-414）：updaterice_micro 算 rice；平衡饱和 seq = exp(2*sig_h2o*mh2o/(rho_ice*rgp*T*rice))（开尔文效应）；growthrate 算阻力 res；隐式格式 dMice = (h2o_vap - seq*zqsat)/(res*zqsat/(cste*No*rice)+1)，cste=4*pi*rho_ice*dt；钳到 [-h2o_ice, h2o_vap]；更新冰/汽；潜热 lw = (2834.3 - 0.28*(T-To) - 0.004*(T-To)²)*1e3，subpdtcloud = dMice*lw/cpp/dt（限幅 ±5/dt）。
     • HDO 分馏（hdo，行 359-393）：凝结时按 alpha_c 分馏系数算 dMice_hdo，升华时按冰相比例；钳幅后更新 hdo 汽/冰。
     • 尘核释放（行 400-414）：h2o_ice <= 1e-28 时汽吸收残冰、dust 吸收 ccn、ccn 归零（hdo 同理）。
     • 温度更新（行 425-428）：.not.activice 时 subpdtcloud=0；zt += subpdtcloud*microtimestep。
     • 自适应反馈（行 433-444）：累计 dMicetot，估 zdq 速率，重算 zimicro。
     • 累加 spenttime、count_micro。

5. 诊断输出（行 452-465）。

伪代码

improvedclouds(ngrid, nlay, ptimestep, ..., zt, zq):

  if firstcall:
    构造对数等比半径网格 rad_cld, rb_cld（存 log 值）
    vo1 = mh2o/rho_ice
    sigma_ice = sqrt(log(1+nuice_sed))

  zq = 上游tracer + pdq*ptimestep   (dust 仅 scavenging; hdo 仅 hdo)
  zt = pt + pdt*ptimestep
  钳 zq >= 1e-30

  dev2 = 1/(sqrt(2)*sigma_ice)
  zqsat = watersat(max(1,zt), pplay)
  zpotcond = h2o_vap - zqsat

  for l, ig:
    if cloud_adapt_ts: zimicro = adapt_imicro(ptimestep, zpotcond)
    spenttime = 0; ending_ts = false
    while not ending_ts:
      zqsat = watersat(zt, pplay)            // 单点
      ph2o = h2o_vap*(mmean/18)*pplay
      satu = h2o_vap/zqsat
      microtimestep = ptimestep/zimicro
      zq0 = zq
      if spenttime+microtimestep >= ptimestep:
        microtimestep = ptimestep-spenttime; ending_ts = true

      // --- 成核 (satu >= 1) ---
      rdust = updaterccn(dust_mass, dust_number)
      n_aer/m_aer = 尘埃 log-normal 经 derf 展开
      rate = nuclea(ph2o, T, satu, n_aer)       // 接触参数 mteta
      dN/dM = Σ aer*(exp(-rate*dt)-1)
      dust += dM/tau; ccn -= dM/tau (数同理)

      // --- 冰增长 (ccn_number*tau >= 1) ---
      rice = updaterice_micro(h2o_ice, ccn_mass, ccn_number)
      seq = exp(2*sig_h2o(T)*mh2o/(rho_ice*rgp*T*rice))
      res = growthrate(T, pplay, ph2o/satu, rice)
      cste = 4*pi*rho_ice*microtimestep
      dMice = (h2o_vap - seq*zqsat)/(res*zqsat/(cste*No*rice)+1)
      dMice = clip(dMice, [-h2o_ice, h2o_vap])
      h2o_ice += dMice; h2o_vap -= dMice
      lw = (2834.3 - 0.28*(T-To) - 0.004*(T-To)^2)*1e3
      subpdtcloud = clip(dMice*lw/cpp/dt, ±5/dt)

      if hdo: dMice_hdo 按分馏系数/冰相比例, 更新 hdo

      // --- 尘核释放 ---
      if h2o_ice <= 1e-28:
        h2o_vap += h2o_ice; h2o_ice = 0
        dust += ccn; ccn = 0  (hdo 同理)

      if not activice: subpdtcloud = 0
      zt += subpdtcloud*microtimestep
      钳 zq >= 1e-30
      if cloud_adapt_ts: zimicro = adapt_imicro(ptimestep, |dMicetot/spenttime|)
      spenttime += microtimestep

参与的主题流程

写法特点

• 自由格式 Fortran 90（.F90），但用 REAL*8/DOUBLE PRECISION 混合精度——成核相关量（ph2o、satu、n_aer、m_aer、Mo、No 等）用双精度，避免 derf 与 exp(-rate*dt) 的精度损失。
• derf 误差函数：声明为 REAL*8 :: derf（行 70），用于把 log-normal 分布积分到各 bin。
• SAVE + THREADPRIVATE：firstcall、rb_cld、sigma_ice、test_flag——OpenMP 下每线程独立。
• 概率成核：exp(-rate*microtimestep)-1 把连续成核率转为该子步内被激活的份额（负号表示从可用气溶胶中扣除）。
• 隐式增长格式：dMice 解析解（行 335）避免显式格式的稳定性约束。
• 开尔文效应：seq 用 sig_h2o(T) 表面张力修正小冰粒的平衡饱和。
• 潜热限幅：subpdtcloud*microtimestep 限在 ±5 K（行 351-356），稳定数值。
• 自适应子步：adapt_imicro 用幂律 zimicro = ceiling(coef*min(max(alpha*|potcond|^beta,5),7000))，coef = ptimestep/defstep，defstep = 88775*5/960 ≈ 462 s（iphysiq=5 的 7.5 min）。
• 硬编码常数：bin 半径范围 rmin_cld=0.1e-6/rmax_cld=10e-6/rbmin_cld=1e-10/rbmax_cld=1e-2；潜热系数 2834.3/0.28/0.004；adapt_imicro 的 alpha=1.88e5/beta=0.458（高自转倾角值，代码内另有注释掉的保守值/当代火星值）。

复现要点

• 入口 ptimestep 是整物理步；子时间步在本例程内部经 adapt_imicro（或固定 imicro）细分——这与简化方案、以及 watercloud_mod ELSE 分支（外部子步循环）的关键区别。唯一调用点 watercloud_mod.F:316，仅 microphys=.true. 时进入。
• tauscaling 把相对尘埃量转绝对量（成核用绝对量），增量再除回 tauscaling 转相对量写回 zq——成核/释放两侧都要除。
• 平衡饱和 seq 用 max(rice,1e-7) 下限保护，防止极小冰粒导致 seq 非物理放大。
• 隐式 dMice 分母 res*zqsat/(cste*No*rice)+1：No = ccn_number*tauscaling + 1e-30 防零除。
• HDO 分馏系数 alpha = exp(16288/T² - 9.34e-2)（行 367）；alpha_c 为含动力学修正的"真实"分馏系数（行 369），仅 hdofrac=.true. 时计算，否则 alpha_c=1。
• 潜热 lw 为水冰升华潜热的温度依赖式（行 347），与 CO2 版 Azreg-Aïnou 多项式不同。

待确认

• microphys_h::mco2、molco2 在 use 列表中（行 19），仅 HDO 分馏的 Dv_hdo 计算（行 365）使用 molco2、mco2——非 HDO 运行时为未使用依赖。
• seq 公式中 rgp（J/mol/K）与 mh2o（kg/molecule? 还是 kg/mol）的单位匹配需核对 microphys_h 定义——推断：mh2o 此处应为摩尔质量量纲才能与 rgp 配平，但行 182 vo1=mh2o/rho_ice 又把 mh2o 当单分子质量用。待确认：mh2o 在两处是否同一量纲，或 seq 公式隐含单位换算。
• nuice_sed 同时用于 sigma_ice（成核分布方差）与沉降——本文件只读不写。

复现风险

• adapt_imicro 的幂律系数 alpha/beta 当前为"高自转倾角"标定值（行 490-491），当代火星模拟应核对是否需改用注释中的当代值（行 492-493），否则子步数偏差影响微物理积分精度。
• 概率成核 exp(-rate*dt) 在 microtimestep 很大时（cloud_adapt_ts=.false. 且 imicro 小）可能高估单步激活量。
• 潜热限幅 ±5 K/步是数值稳定措施，强凝结事件下会人为削弱温度反馈。
• derf 为编译器内置/外部双精度误差函数，不同平台实现差异可能影响成核分布尾部。

相关页面

• watercloud_mod — 水云主调度 - 唯一调用方（microphys=.true. 分支）
• simpleclouds_mod — 简化水冰云方案 - microphys=.false. 时的替代方案
• water-saturation-helpers — 水饱和与凝结温度 - watersat
• improvedco2clouds_mod — CO2 云微物理 - 基于本文件改编的 CO2 版
• updaterad — 冰粒/尘埃/CCN 半径更新 - updaterice_micro/updaterccn 半径反演
• nuclea.md：水冰异质成核率。
• growthrate - 冰晶增长阻力与水汽扩散系数
• sig_h2o.md：水冰表面张力。
• 水循环主题页