co2cloud_mod.F90

路径

LMDZ.MARS\libf\phymars\co2cloud_mod.F90

所属目录/模块

libf/phymars

文件定位

CO2 云形成的主调度模块。根据物理时间步驱动微时间子循环，每个子循环内依次调用 improvedco2clouds 微物理方案、重力沉降（newsedim），并做非负保护。包含 CLFvaryingCO2 次网格云覆盖分数计算（含饱和指数 SatIndex）、CO2 冰粒半径和密度更新、1 µm 光学厚度计算、CO2 饱和诊断输出。被 physiq 在每个物理时间步调用。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• SAVE 变量：imicroco2（微时间子循环步数）、sigma_iceco2、microtimestep、dev2、Qext1bins、Qextv1mic、radv、rb_cldco2、firstcall——全部 SAVE + THREADPRIVATE。
• firstcall：首次调用时读取 imicroco2，计算 CO2 冰粒半径网格和 1 µm 消光系数（从 optprop_co2ice_1mic.dat 读取）。
• 标准输出：firstcall 打印 bin-Qext 信息；NaN 检测时打印诊断并终止。
• 诊断输出（#ifndef MESOSCALE）：SatIndex、SatIndexmap、satuco2、precip_co2_ice_rate、co2ice_cond_rate、pdtcloudco2、riceco2、Tau3D1mic、tau1mic、co2cloudfrac。
• 数据文件依赖：optprop_co2ice_1mic.dat（CO2 冰 1 µm 光学属性，位于 datadir）。

核心逻辑

1. firstcall 初始化（行 300-406）：

   • 读取 imicroco2（默认 MESOSCALE=2，否则=30），计算 microtimestep = ptimestep/imicroco2。
   • 计算 CO2 冰粒半径网格 rad_cldco2（100 bins，对数等间距，rmin=1nm 到 rmax=5µm）和边界 rb_cldco2。
   • 从 optprop_co2ice_1mic.dat 读取 10000 点消光数据，插值到 rb_cldco2 网格得到 Qext1bins。
   • 将 rb_cldco2 转为对数空间存储。

2. 初始化倾向和默认值（行 410-428）：所有累积器清零；rhocloudco2 默认 rho_dust；riceco2 默认 0。

3. 计算层质量和厚度（行 432-441）：masse = (pplev(l)-pplev(l+1) + (bp(l)-bp(l+1))) / g（非 MESOSCALE 时含 bp 修正）。

4. CLFvaryingCO2 次网格云覆盖（行 447-565）：

   • 有效温度 zt = pt + pdt*ptimestep，有效 CO2 vap zq_co2vap = pq(co2) + pdq(co2)*ptimestep。
   • SatIndex（satindexco2=true 时，行 473-511）：在层 12-26（约 12-85 km）计算 Brunt-Väisälä 频率 N、风速 |u|、密度 ρ，按 Spiga et al. (2012) 公式 SatIndex = sqrt(Fo*λ_H*N/(ρ*|u|³)) 计算饱和指数；取层 12-26 最大值为 SatIndexmap。若 satindexco2=false，设 SatIndexmap=0.05。
   • tcond 计算（行 515）：调用 tcondco2(ngrid,nlay,pplay,zq_co2vap,tcond)。
   • cloud fraction（行 519-554）：仅 SatIndexmap ≤ 0.1 时计算。三段逻辑：若 tcond ≥ zt+zdelt（完全饱和），co2cloudfrac=1；若 tcond ≤ zt-zdelt（不饱和），co2cloudfrac=mincloud；否则按线性插值计算。mincloud=0.1 为下限。

5. 主微时间循环（行 582-951）：

   • 有效 tracer 调整（行 574-580）：CLFvaryingCO2 时 pqeff(ccnco2/co2_ice) /= co2cloudfrac。
   • 关键分支差异：improvedco2clouds 接收的是 pteff/pqeff——在 CLFvaryingCO2 路径中 pteff 被覆盖为 pt（行 570），但 pqeff 经过 cloudfrac 调整。这与 watercloud 中 improvedclouds 接收原始 pt/pq 不同。
   • 每个 microstep：
     • 阶梯输入累加 pdt/pdq 到 sum_subpdt/sum_subpdq（行 586-628），覆盖 co2/co2_ice/ccnco2/dust 及条件性 h2o 相关/meteor tracer。
     • 调用 improvedco2clouds（行 632-633）。
     • 非负保护（行 635-700）：co2_ice/co2 双向互反（阈值 1e-12）；ccnco2_number/mass 阈值 1/1e-20 与 dust 互反；meteor CCN 同理（meteo_flux）；ccnco2_h2o 相关三组 tracer 同理（co2useh2o）。
     • 累加 cloud tendency 到 sum_subpdq/sum_subpdt（行 704-746）。
   • 重力沉降（sedimentation=true 时，行 750-950）：
     • 计算沉降时温度和 tracer ztsed/zqsed，调用 updaterice_microco2 得到 riceco2 和 rhocloudco2t。
     • rsedcloudco2 = max(riceco2*(1+nuiceco2_sed)³, rdust)。
     • 对 co2_ice、ccnco2_mass/number、meteor、h2o 相关分别调用 newsedim 沉降。
     • 计算沉降倾向 subpdqsed 并累加回 sum_subpdq。

6. 最终倾向计算（行 953-1026）：

   • pcondicea = sum_subpdq(co2_ice) / imicroco2。
   • pdqs_sedco2/pdqs_sedccn = sum_subpdqs / imicroco2。
   • pdtcloudco2 = sum_subpdt/imicroco2 - pdt。
   • pdqcloudco2(tracer) = sum_subpdq(tracer)/imicroco2 - pdq(tracer)。
   • 关键：pdqcloudco2(co2) = 0（行 982），注释说"this tendency is computed in co2condens"——CO2 气相倾向由 co2condens 处理。

7. 更新粒子半径和光学厚度（行 1028-1106）：

   • 用最终状态调用 updaterice_microco2 更新 riceco2 和 rhocloudco2。
   • 计算 1 µm 消光 Qext1bins2（基于 log-normal 分布和 Qext1bins 插值）。
   • co2useh2o 时调用 updaterice_micro 更新水冰 rice/rhocloud。
   • 调用 updaterdust 更新 rdust。

8. CO2 升华第一层修正（行 1111-1120）：若 pdpsrf*ptimestep > 0.9*(pplev(1)-pplev(2))，用上层 riceco2 覆盖；极端情况（>0.9*(pplev(1)-pplev(3))）用第 3 层覆盖第 2 层。

9. CO2 饱和诊断（行 1122-1135）：调用 co2sat 计算饱和蒸汽压 zqsatco2；计算饱和比 satuco2 = q(co2) * (mmean/mco2*1e3) * pplay / zqsatco2。

10. CLFvaryingCO2 均值化（行 1139-1178）：所有倾向乘以 co2cloudfrac，将云内值映射回网格平均。

11. 1 µm 柱光学厚度（行 1182-1187）：tau1mic = sum(Qext1bins2, dim=layer)。

12. 诊断输出（行 1188-1209，#ifndef MESOSCALE）。

伪代码

co2cloud(ngrid, nlay, ptimestep, ...):

  if firstcall:
    imicroco2 = getin("imicroco2")  // default: 30 (global), 2 (MESOSCALE)
    microtimestep = ptimestep / imicroco2
    构建 rad_cldco2 半径网格 (100 bins, 1nm~5µm)
    读取 optprop_co2ice_1mic.dat → Qext1bins (1µm 消光)
    rb_cldco2 转对数空间

  初始化: sum_subpdq=0, sum_subpdt=0, rhocloudco2=rho_dust, riceco2=0
  计算 masse, epaisseur

  // --- CLFvaryingCO2 次网格方案 ---
  if CLFvaryingCO2:
    zt = pt + pdt*ptimestep
    zq_co2vap = pq(co2) + pdq(co2)*ptimestep
    if satindexco2:
      SatIndex = sqrt(Fo*λ_H*N/(ρ*|u|³))  // 层 12-26
      SatIndexmap = maxval(SatIndex, 12:26)
    else:
      SatIndexmap = 0.05
    tcond = tcondco2(pplay, zq_co2vap)
    // cloudfrac ∈ [mincloud, 1] 由 tcond 与 zt±zdelt 的关系决定
    // 仅 SatIndexmap ≤ 0.1 时计算

  // --- 主微时间循环 ---
  pqeff = pq
  pteff = pt  // 注意：覆盖了 CLFvaryingCO2 中计算的 pteff
  if CLFvaryingCO2:
    pqeff(ccnco2/co2_ice) /= co2cloudfrac

  for microstep = 1..imicroco2:
    累加 pdt, pdq 到 sum_subpdt, sum_subpdq (阶梯输入)
    subpdqcloudco2, subpdtcloudco2 = improvedco2clouds(pteff, pqeff, ...)
    非负保护: co2↔co2_ice, ccnco2↔dust, ccnco2_h2o, meteor
    累加 cloud tendency 到 sum
    if sedimentation:
      计算 ztsed, zqsed, riceco2, rhocloudco2t
      rsedcloudco2 = max(riceco2*(1+nuiceco2_sed)³, rdust)
      对 co2_ice, ccnco2, meteor, h2o 相关分别调用 newsedim
      累加沉降倾向回 sum

  // --- 最终倾向 ---
  pcondicea = sum_subpdq(co2_ice) / imicroco2
  pdqs_sedco2 = sum_subpdqs_sedco2 / imicroco2
  pdtcloudco2 = sum_subpdt/imicroco2 - pdt
  pdqcloudco2(tracer) = sum_subpdq(tracer)/imicroco2 - pdq(tracer)
  pdqcloudco2(co2) = 0  // 由 co2condens 处理

  // --- 更新粒子半径和光学厚度 ---
  updaterice_microco2 → riceco2, rhocloudco2
  Qext1bins2 = 基于 log-normal 和 Qext1bins 的 1µm 消光
  if co2useh2o: updaterice_micro → rice, rhocloud
  updaterdust → rdust
  CO2 升华第一层修正

  // --- CO2 饱和诊断 ---
  zqsatco2 = co2sat(pteff + (pdt+pdtcloudco2)*ptimestep)
  satuco2 = q(co2) * (mmean/mco2*1e3) * pplay / zqsatco2

  // --- CLFvaryingCO2 均值化 ---
  if CLFvaryingCO2:
    所有倾向 *= co2cloudfrac

  // --- 1µm 柱光学厚度 ---
  tau1mic = sum(Qext1bins2, dim=layer)

参与的主题流程

写法特点

• 自由格式 Fortran 90（.F90 后缀）。
• SAVE + THREADPRIVATE：imicroco2、sigma_iceco2、microtimestep、dev2、Qext1bins、Qextv1mic、radv、rb_cldco2、firstcall。
• firstcall 模式：首次调用时读取参数、构建半径网格、读取光学属性文件、打印诊断。
• 阶梯输入法：与 watercloud 相同，每个 microstep 累加 pdt/pdq。
• 非负保护逻辑：co2↔︎co2_ice 双向互反（阈值 1e-12），ccnco2↔︎dust 互反（阈值 1/1e-20），meteor 和 h2o 相关 CCN 同理。
• pdqcloudco2(co2) = 0（行 982）：CO2 气相倾向故意设为零，由 co2condens 负责计算——这是与 watercloud 的关键差异。
• pteff 覆盖（行 570）：CLFvaryingCO2 路径中先计算 pteff（行 449-539），但行 570 将其覆盖为 pt。实际传给 improvedco2clouds 的是 pteff=pt + cloudfrac 调整后的 pqeff——与 watercloud 中 improvedclouds 接收原始 pt/pq 不同。
• SatIndex（Spiga et al. 2012）：仅在 satindexco2=true 时计算，用于调制云覆盖分数。
• 编译宏：MESOSCALE 影响 imicroco2 默认值、层质量计算（bp 修正）和诊断输出。
• 硬编码常数：mincloud=0.1、beta=0.85（不规则粒子形状修正）、threshold=1e-30、threshold_2=1e-13、rmin_cld=1e-9 m、rmax_cld=5e-6 m、Fo=7.5e-7 J/m³、lambdaH=150 km。
• 数据文件依赖：optprop_co2ice_1mic.dat（10000 点 CO2 冰 1 µm 消光系数）。

复现要点

• imicroco2 默认值依赖编译宏 MESOSCALE（2 vs 30），可通过 imicroco2 覆盖。
• pdqcloudco2(co2) 始终为零——CO2 气相倾向由 co2condens 在外部计算，调用者必须知道这一点。
• 非负保护的执行顺序（co2↔︎co2_ice → ccnco2↔︎dust → meteor → h2o 相关）影响最终倾向值。
• CO2 升华第一层修正阈值 0.9 与 watercloud 相同。
• optprop_co2ice_1mic.dat 必须存在于 datadir，否则 abort_physic 终止。
• 沉降半径 rsedcloudco2 = max(riceco2*(1+nuiceco2_sed)³, rdust)，使用 nuiceco2_sed（来自 tracer_mod）。
• SatIndex 仅在层 12-26（约 12-85 km）计算，反映 CO2 云主要形成在高层大气。

待确认

• zdelt（行 250 声明、行 524-536 使用）：Delta T for temperature distribution，声明但未赋值——可能是未初始化变量的 bug，或依赖编译器默认零初始化。复现时需确认其值。
• density_co2_ice（行 107 use）：在 co2cloud 本体中未被直接调用，是否为 improvedco2clouds 内部所需——Fortran use 不自动传递，需确认是否多余。
• nuice（输入参数）：声明为输入但代码中未使用，可能是接口遗留。
• r3n_q（从 tracer_mod use）：代码中未直接使用，可能为 improvedco2clouds 或其他模块所需。

复现风险

• zdelt 未初始化：CLFvaryingCO2 云覆盖分数计算依赖 zdelt，若其值不确定则整个次网格方案不可复现。
• pdqcloudco2(co2) = 0 是硬编码行为：调用者必须确保 co2condens 正确计算 CO2 气相倾向，否则 CO2 质量不守恒。
• optprop_co2ice_1mic.dat 缺失将导致运行终止。
• SatIndex 中 NN = sqrt(g/zt(iq,l) * ...) 使用 iq 而非 ig（行 479）——iq 是 tracer 循环变量，此处应为 ig（grid point），可能是 bug。

相关页面

• co2-saturation-helpers — CO2 饱和与凝结温度
• updaterad — 冰粒/尘埃半径更新
• newsedim_mod — 单示踪物重力沉降
• improvedco2clouds_mod- CO2 云微物理核心方案
• co2condens- CO2 地表/大气凝结（计算 pdqcloudco2(co2)）
• co2snow- CO2 雪/地表沉积
• CO2 循环主题页