水循环主题

主题目标

本页串联 LMDZ.MARS 中水汽、云冰、地表霜、地下冰和诊断输出的主流程，目标是让读者能沿源码复现一个物理时间步中的水循环计算：

1. physiq 中何处调用 watercloud，调用前后哪些 tendency 被注入。
2. watercloud_mod 如何按 microphys 分派到 simpleclouds 或 improvedclouds。
3. watersat 与 tcondwater 如何提供水冰饱和混合比和凝结温度。
4. vdifc 如何处理水汽地表交换，并在 adsorption_soil 分支耦合 soilwater 的吸附水/地下冰储库。
5. waterice_tifeedback 如何把地表/孔隙水冰反馈到土壤热惯量，再影响地温。
6. 配置开关、tracer、XIOS 输出如何决定水循环是否可运行、可诊断。

源码依据：本页核验了 physiq_mod.F、watercloud_mod.F、watersat_mod.F90、tcondwater.F90、vdifc_mod.F、soilwater.F90、waterice_tifeedback_mod.F90、conf_phys.F、deftank/callphys.def.GCM6、traceur.def.watercycle.GCM5、deftank/xios/field_def_physics.xml、deftank/xios/file_def_physics.xml。更细的文件级证据见相关子页。

入口

参与文件

主调用链

physiq
  -> callradite(... rice, nuice, rdust ...)          ! 本步辐射先使用上一状态/上一步留下的半径
  -> IF (water):
       call watercloud(... pq,pdq,pdt, tau, tauscaling,
                       rdust, rice, nuice, rsedcloud, rhocloud, totcloudfrac)
         -> IF (microphys):
              improvedclouds(...)                   ! 成核、增长/升华、CCN/尘埃耦合
            ELSE:
              simpleclouds(...)                     ! 瞬时饱和调整
         -> watersat / tcondwater / updaterice_* / updaterdust
       inject zdtcloud/zdqcloud into pdt/pdq
  -> IF (sedimentation):
       callsedim(... rice, rsedcloud, rhocloud ...)  ! 水冰沉降
  -> diagnostics / XIOS write_output

vdifc
  -> IF (water .and. iq == h2o_vap):
       watersat(ztsrf, pplev(:,1)) -> qsat
       implicit lower boundary pulls near-surface vapor toward qsat
       optional soilwater(...) when adsorption_soil
       update pdqsdif(h2o_ice), dwatercap_dif, pdtsrf, flux_ssice

physiq
  -> IF (surfaceice_tifeedback .or. poreice_tifeedback):
       waterice_tifeedback(qsurf(h2o_ice), pore_icefraction) -> inertiesoil_tifeedback
       soil(..., inertiesoil_tifeedback, ...)

计算阶段

1. 配置与 tracer 准备

水循环必须从配置开始复现。conf_phys.F 中 water 的源码默认值为 .false.，样例 callphys.def.GCM6 显式设为 .true.；完整水云微物理通常还需要 microphys=.true.、temp_dependent_m=.true.、cloud_adapt_ts=.true.、scavenging=.true.。源码校验会中止若 microphys 没有同时满足 water 和 doubleq，也会中止若 scavenging 没有 microphys。

最小水循环 tracer 是 h2o_ice 和 h2o_vap。微物理路径还依赖 ccn_number、ccn_mass、dust_number、dust_mass；HDO 路径依赖 hdo_vap、hdo_ice。traceur.def.watercycle.GCM5 的样例包含 co2、两类尘埃、两类 CCN、h2o_ice、h2o_vap。

源码依据：conf_phys.F:715 读取 water；:852 读取 microphys；:858 读取 supersat；:879 读取 cloud_adapt_ts；:906-911 校验 microphys；:921-926 校验 scavenging；traceur.def.watercycle.GCM5:7-8 给出 h2o_ice/h2o_vap。

2. 辐射先使用当前 rice

physiq 在 watercloud 前调用 callradite。该辐射调用以 rdust、rice、nuice 为 INOUT，用于当前步辐射和 taucloudtes 等量。随后 watercloud 使用同一组半径数组并更新它们，更新后的 rice 再供本步 callsedim、诊断和下一步辐射使用。

复现风险：不能把 watercloud 提前到辐射前，否则 rice 的时间层会改变。physiq-water-cycle-callpoints.md 明确把时序写成“本步辐射用上一步留下的 rice，水云随后更新供后续沉降/诊断和下一步辐射”。

3. watercloud 分派云方案

physiq_mod.F:1831 在 IF (water) 内调用 watercloud。watercloud_mod.F:315-349 的关键分支是：

• microphys=.true.：调用 improvedclouds。传入整物理步长 ptimestep，自适应子步在 improvedclouds 内部处理。该路径使用原始网格平均 pt/pq，不使用 CLFvarying 改写后的 pteff/pqeff。
• microphys=.false.：在外层 imicro 子步循环里调用 simpleclouds。该路径可接收 CLFvarying 后的云内有效 pteff/pqeff，每个子步把已有 pdt/pdq 累加成阶梯输入。

watercloud 返回 pdqcloud、pdtcloud、rice、nuice、rsedcloud、rhocloud、totcloudfrac。随后它统一执行非负保护、watersat 过饱和限制（microphys 且 supersat=false）、updaterice_micro/updaterice_typ 半径更新、rsedcloud=max(rice*(1+nuice_sed)^3, rdust)，并在 CLFvarying 时把云内 tendency 乘回 cloudfrac。

源码依据：watercloud_mod.F:316 调用 improvedclouds；:349 调用 simpleclouds；:564 调用 watersat 做 supersat=false 检查；:618 调用 updaterice_micro；:638 调用 updaterice_typ；:665 计算 rsedcloud；:677 进入 CLFvarying 缩放。

4. 简化云与完整微物理的差别

simpleclouds 是瞬时饱和调整：先用 pq + pdq*dt 和 pt + pdt*dt 得到局地状态，调用 watersat，若水汽超过 zqsat 就把超量全部转成水冰；若不饱和，就让现有水冰升华但不超过冰量。潜热用温度二次式：

lw = (2834.3 - 0.28*(T-To) - 0.004*(T-To)^2) * 1000

improvedclouds 是完整微物理：先计算 zqsat 和可凝结量，再在每个格点/层内做子步循环。若 satu >= 1，用 updaterccn、尘埃 log-normal bin、nuclea 把可用尘埃激活为 CCN；若 CCN 数量足够，用 updaterice_micro、sig_h2o、growthrate 和隐式 dMice 公式处理冰增长/升华；冰耗尽时释放 CCN 回尘埃。cloud_adapt_ts 为真时，adapt_imicro 会按可凝结量和经验幂律更新子步数。

待确认：improvedclouds 中 seq 的单位匹配、nuice_sed 同时用于成核分布方差和沉降的物理含义仍以子页的待确认项为准；这里仍按待确认项处理。

5. 饱和与凝结温度公式

watersat 给定温度 T 和压力 P，使用硬编码分子量比：

epsi = 18 / 44
psat = 100 * 10^(2.07023 - 0.00320991*T - 2484.896/T + 3.56654*log10(T))
if psat > P:
  qsat = 1
else:
  qsat = epsi*psat/(P - (1-epsi)*psat)
  qsat = max(qsat, 1e-30)

tcondwater 给定压力和总水混合比 q，先算：

x = log10(max(q, 1e-16) * p / (100 * 0.41))
tcond = 252.846556
      + 23.9229653*x
      + 2.21333897*x^2
      + 0.179978042*x^3
      + 0.0100068018*x^4
      + 0.000255145044*x^5

复现风险：epsi=18/44 和 0.41 都把干大气近似为 CO2；若要复现非纯 CO2 平均分子量情形，不能只改外部 mmean，因为这两个助手没有读取 mmean。tcondwater 的拟合适用范围也未在源码注释中明确。

6. physiq 注入 tendency

watercloud 返回后，physiq 按固定顺序注入：

1. activice=true 时，pdt += zdtcloud。若 activice=false，水云仍改水汽/水冰 tracer，但不把潜热反馈进温度。
2. 无条件注入 h2o_vap 和 h2o_ice 的 zdqcloud。
3. hdo=true 时注入 hdo_vap 和 hdo_ice。
4. microphys=true 时注入 ccn_mass 和 ccn_number，并对成对质量/数浓度做负值保护。
5. scavenging=true 时注入 dust_mass 和 dust_number，并做同样的成对负值保护。

源码依据：physiq_mod.F:1837-1908。复现时不能打乱这个顺序，尤其是 CCN 和尘埃的成对保护，否则单精度舍入导致的负值修正会不同。

7. 与沉降、辐射和诊断的关系

水冰沉降不在 watercloud 内完成，而在 physiq 第 9b 段调用 callsedim。callsedim 使用 watercloud 更新出的 rice、rsedcloud、rhocloud，返回 zdqsed 和 zdqssed，再注入大气 pdq 和地表 dqsurf。因此顺序是“水云先凝结/升华并更新半径，沉降再用这些半径沉降”。

诊断通常在 zq = pq + pdq*ptimestep 后进行，使用包含水云、沉降等全部 tendency 的更新态。physiq 里可见 watersat 用于诊断过饱和度，write_output 输出 h2o_ice、h2o_vap、h2o_ice_s、rice、zdqcloud_ice、zdqcloud_vap 和二者之和。XIOS XML 还定义 mtot、icetot、reffice、Nccntot、Mccntot、flux_ssice、watercap、各坡面地表水冰等字段；是否进入文件取决于 file_def_physics.xml 是否引用。

8. vdifc 地表水汽交换

vdifc 的水汽段由 if ((water).and.(iq.eq.igcm_h2o_vap)) 控制。它对每个坡面把地表水冰 pqsurf(h2o_ice,islope) 除以 cos(slope) 转成真实坡面面积上的储量，然后用 make_tsub 在有霜/永久冰盖时按地表温度变化细分子步。

每个子步中，vdifc 构造隐式三对角扩散系数。地表边界调用 watersat(1,ztsrf,pplev(:,1),qsat)，把近地层水汽拉向地表温度下的饱和混合比：

zq1temp = zc(1) + zd(1)*qsat
zdqsdif_surf = rho * dryness * zcd * (zq1temp - qsat)

这里正值表示向上升华，负值表示向地表沉积。old_wsublimation_scheme 决定地表交换系数用动量导度 zcdv 还是热导度 zcdh；dryness 线性缩放地表水汽导度。若 latentheat_surfwater 为真，地表潜热通量会反馈到 pdtsrf。

输出侧，vdifc 更新 pdqsdif(h2o_ice)、dwatercap_dif、pdtsrf、zq(h2o_vap)，并对坡面量成对做 /cos 和 *cos 投影。watercaptag 永久冰盖上超出当前霜量的升华会走 dwatercap_dif。

9. soilwater 吸附和地下冰耦合

当 adsorption_soil=true，vdifc_mod.F:1052 逐点逐坡面调用 soilwater。soilwater 维护 regolith 中三类地下 tracer：h2o_vap_soil、h2o_ice_soil、h2o_ads_vap。

它在首次调用时分配 SAVE 数组并读取 pqsoil 初值；每步根据孔隙冰饱和度更新扩散系数，用 choice_ads 选择吸附/解吸方案，用嵌套迭代处理单层吸附饱和、地下凝结/升华和冰层关闭。返回给 vdifc 的关键量是 zdqsdifrego（地下到大气水汽通量）和 zq1temp2（耦合后的近地水汽边界值）。当孔隙冰饱和度超过 choke_fraction=0.8 时，soilwater 会按通量方向关闭相邻层间扩散。

待确认：choice_ads 在 comsoil_h 中声明为实数但代码用 ==0/1/2 比较；注释里“3 means no adsorption”与实际 choice_ads==0 不一致。latitude_deg、H2O/H2O_depth 等地图变量有遗留痕迹。复现地下水储库时应优先参考 soilwater.md 的待确认和风险列表。

10. 水冰热惯量反馈

waterice_tifeedback 不直接改 tracer，也不直接改地温；它只计算 newtherm_i(ngrid,nsoil,nslope)，随后由 soil 地温求解使用。physiq_mod.F:737 和 :2487 在 surfaceice_tifeedback.or.poreice_tifeedback 时调用该例程。

两条路径互斥：

• surfaceice_tifeedback=true：把地表水冰质量 qsurf(h2o_ice) 除以 rho_ice 得到厚度。冰厚超过土壤总深度或 watercaptag=true 时，整列热惯量设为 inert_h2o_ice；冰薄于第一层时保留 inertiedat；中间厚度时，上方设冰热惯量，边界层用串联导热公式过渡，下方保留风化层热惯量。
• poreice_tifeedback=true：用 sqrt(inertiedat^2 + porosity_reg*poreice*inertie_purewaterice^2) 抬高热惯量。

复现风险：surfaceice_tifeedback 还依赖 water 与 callsoil；poreice_tifeedback 依赖 adsorption_soil。inert_h2o_ice 可由配置覆盖，样例 GCM6 用 800，而源码默认是 2400；孔隙冰路径的 inertie_purewaterice=2100 是例程内部硬编码值，不能混用。

状态变量

配置开关

数据依赖

与其他主题的耦合

输出结果

伪代码总览

for each physics timestep:
  call radiation with current rice/nuice/rdust

  if water:
    call watercloud:
      if CLFvarying: tcondwater(total water) -> cloudfrac
      if microphys:
        improvedclouds:
          watersat -> saturation
          if saturated: nuclea/updaterccn -> CCN activation
          if CCN present: growthrate/sig_h2o/updaterice_micro -> ice growth or sublimation
          if ice exhausted: release CCN back to dust
      else:
        for imicro substeps:
          simpleclouds:
            watersat -> qsat
            condense excess vapor or sublime existing ice
      enforce non-negative H2O/HDO/CCN/dust
      if not supersat: watersat -> remove excess vapor
      update rice, rsedcloud, rhocloud
      if CLFvarying: multiply cloud tendency by cloudfrac

    if activice: pdt += zdtcloud
    pdq(h2o_vap/h2o_ice) += zdqcloud
    if hdo: pdq(hdo_vap/hdo_ice) += zdqcloud
    if microphys: pdq(ccn_mass/number) += zdqcloud with pairwise guards
    if scavenging: pdq(dust_mass/number) += zdqcloud with pairwise guards

  if sedimentation:
    callsedim using rice/rsedcloud/rhocloud
    pdq += zdqsed
    dqsurf += zdqssed

  in vdifc for h2o_vap:
    for each slope:
      convert surface ice by /cos(slope)
      make_tsub if frost or watercap
      for each substep:
        watersat(surface T, surface pressure) -> qsat
        solve implicit vapor diffusion with lower boundary qsat
        if adsorption_soil: soilwater -> regolith flux and boundary vapor
        else if subsurface ice active: SSI flux branch
        update surface ice, watercap tendency, surface temperature
      project output back by *cos(slope)

  if surfaceice_tifeedback or poreice_tifeedback:
    waterice_tifeedback -> new soil thermal inertia
    soil uses new thermal inertia

  diagnostics:
    zq = pq + pdq*ptimestep
    write water vapor/ice, surface ice, radius, cloud tendencies, subsurface fluxes

源码疑点与待确认项

• vdifc_mod.F:901-902 的非水汽 tracer 段守卫写成 (.not. water).or.(.not. iq.eq.igcm_h2o_vap).or.(.not. iq.eq.igcm_hdo_vap)。按布尔逻辑，对任意 iq 至少一个“不等于”条件为真，整体近似恒真。推断它可能想排除 h2o_vap/hdo_vap，但源码实际未这样写；因为 h2o_vap 后续由专段覆盖，净效应仍需对照历史版本确认。复现时应按源码字面执行。
• watersat/tcondwater 硬编码 18/44 或 0.41，没有读取实时平均分子量；非纯 CO2 平均分子量影响未在这两个助手中体现。
• improvedclouds_mod.F90 中 Kelvin 项 seq 的单位量纲和 mh2o 用法需对照 microphys_h 定义继续核验；本页只按子页记录其源码行为。
• soilwater.F90 中 choice_ads 用浮点等值比较，并且注释与实际“无吸附”取值存在不一致；地下水复现应记录实际配置和编译器。
• soilwater 的非收敛路径会 stop，长积分极端条件下可能直接中止模型。
• waterice_tifeedback 两个反馈机制用 IF ... ELSE IF 互斥；若两个开关都为真，实际走 surfaceice_tifeedback，不应写成二者叠加。
• calchim 接收 zdqcloud/zdqscloud 的具体化学用途尚未核验，需化学模块页确认。

相关页面

• physiq 水循环调用点
• watercloud_mod
• simpleclouds_mod
• improvedclouds_mod
• water-saturation-helpers
• vdifc 水汽地表交换段
• soilwater
• waterice_tifeedback_mod
• 水循环配置