growthrate.F

路径

LMDZ.MARS\libf\phymars\growthrate.F

所属目录 / 模块

libf\phymars

文件定位

growthrate.F 定义 growthrate_mod，提供水冰晶增长/升华阻力的局部解析辅助例程。它不直接推进水冰质量，而是根据温度、层中压强、饱和水汽分压和当前冰晶半径，计算水汽扩散系数 Dv 以及总阻力 res = Rk + Rd*(1+lambda*knudsen)。

当前源码搜索中，唯一外部调用点在 improvedclouds_mod.F90 的水冰增长段：improvedclouds 先用 sig_h2o 算 Kelvin 平衡饱和修正 seq，再调用 growthrate 得到 res，最后把 res 放入隐式 dMice 公式，决定本微物理子步内水汽和水冰之间交换多少质量。因此本文件是完整水云微物理中“单颗冰晶增长阻力/扩散系数”的计算部件，而不是云方案主调度器。

定义的符号

依赖的模块

调用的关键例程

上游调用点

输入

输出

共享状态与副作用

growthrate_mod 不定义保存型 module 变量，不读写文件，不打印日志，也不读取运行时配置。它只读取三个共享常数来源：tracer_mod::rho_ice、comcstfi_h::pi、microphys_h 的分子物理常数。副作用限于写入实参 res 和 Dv。

需要注意的是，rho_ice 是 tracer_mod 中的 save 变量，并声明为 OpenMP THREADPRIVATE；本文件不初始化它。复现 growthrate 时必须保证上游 tracer 初始化已经给 rho_ice 赋值，否则阻力量级会错误。

核心逻辑

1. 由温度计算 CO2 热导率 k = (0.17913*temp - 13.9789) * 4.184e-4。
2. 由温度和参考温度 To 计算水冰升华潜热 Lv = (2834.3 - 0.28*(temp-To) - 0.004*(temp-To)^2) * 1.e3。
3. 构造平均自由程常数 afactor = 0.707*rgp/(4*pi*molco2^2*nav)。
4. 用 H2O 热速度、H2O/CO2 碰撞截面和压强计算水汽扩散系数 Dv。
5. 用 temp/pmid*afactor/rcrystal 得到 Knudsen 数 knudsen。
6. 用 (1.333+0.71/knudsen)/(1+1/knudsen) 得到非连续修正系数 lambda。
7. 计算潜热释放热阻 Rk = Lv^2*rho_ice*mh2o/(k*rgp*temp^2)。
8. 计算水汽扩散阻力 Rd = rgp*temp*rho_ice/(Dv*psat*mh2o)。
9. 输出 res = Rk + Rd*(1+lambda*knudsen)；源码注释明确没有直接修正 Dv，而是把修正放入阻力项。

伪代码

growthrate(temp, pmid, psat, rcrystal):
  k  = thermal_conductivity_CO2(temp)
  Lv = latent_heat_H2O(temp)

  afactor = 0.707 * R / (4*pi*molco2^2*Avogadro)

  Dv = one_third * thermal_speed_H2O(temp)
       * k_B * temp
       / (pi * pmid * (molco2+molh2o)^2 * sqrt(1+mh2o/mco2))

  knudsen = temp / pmid * afactor / rcrystal
  lambda  = (1.333 + 0.71/knudsen) / (1 + 1/knudsen)

  Rk = Lv^2 * rho_ice * mh2o / (k * R * temp^2)
  Rd = R * temp * rho_ice / (Dv * psat * mh2o)

  res = Rk + Rd * (1 + lambda * knudsen)
  return res, Dv

in improvedclouds:
  seq = Kelvin_equilibrium_saturation(...)
  growthrate(T, pplay, ph2o/satu, rice, -> res, -> Dv)
  dMice = (h2o_vap - seq*zqsat) / (res*zqsat/(4*pi*rho_ice*dt*No*rice) + 1)

参与的主题流程

写法特点

• 文件扩展名为 .F，但内部使用 module/subroutine 结构；源码采用固定格式续行 &。
• Kelvin 平衡饱和项 seq 不在本文件计算，只在注释中保留公式；当前调用方 improvedclouds 已在调用前计算 seq。
• 源码注释保留了旧方案的显式增长率公式 growth = ... 和半径更新公式，但当前实际运行路径只输出阻力 res。
• Dv 的非连续修正行被注释掉；实际修正体现在 res = Rk + Rd*(1+lambda*knudsen)。

复现要点

• pmid 必须按压力传入。虽然参数注释写了 (K)，Dv 和 knudsen 公式都把它作为 Pa 量级的压强。
• psat 在当前主调用中来自 ph2o/satu，不是 Kelvin 修正后的 seq*zqsat；Kelvin 修正由 improvedclouds 单独处理。
• 不要把 Dv 再除以 (1+lambda*knudsen)；源码明确把这一行注释掉，并把修正放在 Rd 的阻力乘子里。
• rcrystal 不在本文件内部做下限保护；调用方在计算 seq 时对 rice 用了 max(rice,1.e-7)，但传给 growthrate 的仍是原始 rice(ig,l)。
• rho_ice 由 tracer_mod 运行时状态提供，不是本文件常数；复现实验需先确认 tracer 物性初始化路径。

待确认

• res 的严格物理量纲没有在源码注释中给出；本页只记录它在 improvedclouds 隐式分母中的实际用法。
• rcrystal 若极小或为零时本文件没有保护，是否完全依赖 updaterice_micro 的输出保证正半径，需要在 updaterad.F90 专页或后续源码核验中确认。

相关页面

• improvedclouds_mod
• watercloud_mod
• water-cycle
• water-cycle-config
• tracer_mod
• comcstfi_h
• microphys_h：提供本例程读取的分子常数、参考温度和气体常数。