LMDZ.MARS Wiki

vdifc 水汽地表交换段

本页只追踪 vdifc（垂直湍流扩散例程）中 h2o_vap 与地表/地下冰交换的代码段：自适应子时间步、隐式三对角求解的地表边界条件、watersat 饱和约束、soilwater 吸附耦合、地下冰（SSI）通量，以及输出 pdqsdif(igcm_h2o_ice)/dwatercap_dif。不覆盖 vdifc 的动量/位温扩散、CO2 质量变化方案、尘埃抬升或 HDO 段。

所属文件

LMDZ.MARS\libf\phymars\vdifc_mod.F

源码 1530 行。本页核心段 h2o_vap 处理在行 959–1315；其前置的「非水汽示踪物」段（含 h2o_ice 地表通量注入）在行 900–957。

所属模块

vdifc_mod（SUBROUTINE vdifc / vdifc_mod 内）

相关 use 与声明

行 对象 用途
19–22 igcm_h2o_vap, igcm_h2o_ice, igcm_co2 tracer 索引
23–24 watercaptag, frost_albedo_threshold, dryness, old_wsublimation_scheme 永久冰盖标记、霜反照率阈值、地表干度、旧升华方案开关
26 watersat 计算饱和混合比 qsat(T,p)
35 To（来自 microphys_h） 三相点温度（潜热公式、融化封顶）
36–37 coef_ssdif, h2oice_depth, lag_layer, zdqsdif_ssi_tot（paleoclimate_mod） 地下冰扩散系数、冰埋深、滞后层开关、SSI 总通量诊断
38 adsorption_soil（comsoil_h） 是否启用完整吸附/解吸地下水模型
44–47 water, latentheat_surfwater（callkeys_mod） 水循环主开关、地表水潜热反馈开关

段内关键局部变量（声明见行 145–251）：

变量 含义
za(ngrid,nlay) ρ·dz = (pplev(l)-pplev(l+1))/g，层质量（行 326）
zb0(ngrid,nlay) 扩散系数基量；zb0(:,1)=ptimestep·pplev(:,1)/(r·ptsrf_tmp)（行 343），地表项
zb,zc,zd,z1 隐式三对角求解的工作数组（向下扫描系数）
zq_tmp_vap(ngrid,nlay,nq) 子时间步内的水汽工作副本
zq1temp(ig) 求解出的第 1 层（近地层）水汽混合比
qsat(ig) 地表温度下的饱和混合比
zdqsdif_surf / zdqsdif_tot 子时间步地表水汽通量（湍流 / 含地下贡献的总和）
zqsurf(ig) 该坡度的地表水冰储量（去坡度投影后，kg/m²）
nsubtimestep(ig) 自适应子步数（make_tsub 给出）
dtmax=0.5 子步温度判据（K，行 188）
Tice, qsat_ssi, resist 地下冰温度、其饱和混合比、通量衰减阻力系数

控制条件与总体结构

整段在示踪物循环 do iq=1,nq 内，由 if ((water).and.(iq.eq.igcm_h2o_vap))（行 965）守卫，到行 1313 END IF 结束。内层是坡度循环 DO islope=1,nslope（行 967）与列循环 DO ig=1,ngrid（行 984），再套自适应子时间步循环 DO tsub=1,nsubtimestep(ig)（行 990）。

do iq: if water .and. iq==h2o_vap:
  DO islope:
    zqsurf = pqsurf(h2o_ice,islope)/cos(slope)      ! 去坡度投影
    make_tsub -> nsubtimestep                        ! 自适应子步数
    saved_h2o_vap = zq(:,1,h2o_vap)
    DO ig:
      subtimestep = ptimestep/nsubtimestep(ig)
      zq_tmp_vap = zq
      DO tsub=1,nsubtimestep(ig):
        构造 zb（地表项用 zcdv 或 zcdh，乘 dryness）
        向下扫描 zc/zd（隐式三对角）
        watersat -> qsat(ztsrf)
        zq1temp = zc(1) + zd(1)*qsat            ! 地表边界 = qsat
        zdqsdif_surf = rho·dryness·zcd·(zq1temp-qsat)
        zdqsdif_tot = zdqsdif_surf
        [无霜判据] exchange = ...
        if adsorption_soil: call soilwater(...) 改写 zq1temp/zdqsdif_tot
        elif 无霜可升华超量: 钳制 + 可选地下冰<->大气通量
        [有霜+地下冰]: frost<->ssi 通量
        zq_tmp_vap(1,h2o_vap)=zq1temp
        if latentheat_surfwater: zdtsrf = zdqsdif_tot·lh/pcapcal
        向上回代 zq_tmp_vap(2..nlay)
        zqsurf += zdqsdif_tot·subtimestep；钳到 >=0（非 watercap）
        累计 SSI/潜热诊断；冰存在时温度封顶 To
      ENDDO tsub
      pdtsrf = (ztsrf-ptsrf)/ptimestep
      pdqsdif(h2o_ice,islope) = (zqsurf - pqsurf/cos)/ptimestep
      if watercaptag .and. 超量升华: 走 dwatercap_dif，pdqsdif 钳到 -qsurf/dt
      zq_slope_vap(...,islope)=zq_tmp_vap
    ENDDO ig
  ENDDO islope
  zq(h2o_vap) = Σ_islope zq_slope_vap·subslope_dist      ! 坡度平均回大气
  pdqsdif(h2o_ice,:) *= cos(slope)；dwatercap_dif *= cos(slope)  ! 回投影

核心计算步骤

1. 去坡度投影与自适应子步（行 970–979）

zqsurf = pqsurf(igcm_h2o_ice,islope)/cos(pi·def_slope_mean(islope)/180)，把坡度投影后的地表冰储量还原为单位真实面积上的量。make_tsub（定义在行 1445）按地表温度变化率估升华，子步数 ntsub = ceiling(|ptimestep·dtsurf/dtmax|)；只有在永久冰盖 watercaptag 或霜厚 ≥frost_albedo_threshold 时才细分，否则 ntsub=1。目的是把地表水冰升华潜热释放分摊到足够小的步长，控制每步温度变化 ≤ dtmax=0.5 K。

2. 隐式三对角向下扫描（行 992–1016）

地表交换系数 zb(:,1)：old_wsublimation_scheme 用动量整体系数 zcdv，否则用热整体系数 zcdh（行 994–1000），再乘地表干度 dryness（行 1001）。各层 zb=zkh·zb0。从顶层向下递推消元系数 zc/zd（标准 Thomas 算法前扫），到第 1 层得 zc(1)。

3. 饱和地表边界条件（行 1018–1028）

call watersat(1,ztsrf,pplev(1),qsat) 求地表温度下饱和混合比。地表边界取 zq1temp = zc(1) + zd(1)·qsat，即近地层水汽被拉向 qsat（地表是饱和源/汇）。湍流地表通量

zdqsdif_surf = rho · dryness · zcd · (zq1temp - qsat)

正值表示水汽向上（升华），负值表示向地表沉积（凝霜）。zdqsdif_tot 初始化为 zdqsdif_surf。

4. 地下/吸附耦合（行 1035–1204，三选一分支）

5. 子步水/温预算（行 1208–1248）

zq_tmp_vap(1,h2o_vap)=zq1temp 后，若 latentheat_surfwater：潜热 lh=(2834.3-0.28(T-To)-0.004(T-To)²)·1e3（升华焓，J/kg），地表温度倾向 zdtsrf = zdqsdif_tot·lh/pcapcal。向上回代得各层水汽。子步末更新 ztsrf += (pdtsrf+zdtsrf)·subtimestep、zqsurf += zdqsdif_tot·subtimestep，并把 zqsurf 钳到 ≥0（非 watercap）。冰仍存在时把温升封顶到 To（防止有冰却升温越过三相点，行 1243–1248）。

6. 子步积分回主步（行 1256–1311）

pdtsrf = (ztsrf-ptsrf)/ptimestep；地表水冰倾向

pdqsdif(ig,igcm_h2o_ice,islope) = (zqsurf - pqsurf(h2o_ice)/cos) / ptimestep

watercaptag 且升华超过现有霜时，多出的部分计入 dwatercap_dif（永久冰盖储库消耗，<0），并把 pdqsdif(h2o_ice) 钳到只消耗现有霜。坡度循环结束后，zq(h2o_vap) 用 subslope_dist 做坡度加权平均回到网格大气值（行 1290–1299），pdqsdif(h2o_ice) 和 dwatercap_dif 再乘 cos(slope) 回到坡度投影量（行 1301–1311）。

与「非水汽示踪物」段的关系（行 900–957）

h2o_vap 之前，所有其他示踪物（含 h2o_ice）在行 900–957 处理。其中地表系数 zb(:,1)=0（行 906，无湍流地表通量），但 h2o_ice/hdo_ice 走特例：地表边界只做层间混合（行 924–932），普通示踪物则把抬升通量 -pdqsdif_tmp(iq) 作为地表源注入（行 937–939）。pdqsdif(iq,islope)=pdqsdif_tmp(iq)·cos(slope)（行 951）。注意：h2o_ice 的地表倾向真正的物理值是在后面的 h2o_vap 段算出并写入 pdqsdif(igcm_h2o_ice,...) 的——这里的循环条件 (.not.water).or.(.not.iq==h2o_vap).or.(.not.iq==hdo_vap)（行 901–902）逻辑上恒为真（见待确认），但 h2o_vap 自身在 else 体外、由行 965 的专段处理。

边界条件与保护逻辑

输入数据来源

输出与副作用

输出 含义
pdqsdif(ig,igcm_h2o_ice,islope) 地表水冰升华/沉积倾向（kg/m²/s，坡度投影）
dwatercap_dif(ig,islope) 永久冰盖储库消耗倾向（<0，坡度投影）
pdtsrf(ig,islope) 地表温度倾向（含水冰潜热反馈）
zq(ig,:,igcm_h2o_vap) 近地及各层水汽混合比（坡度平均后写回，供后续 pdqdif 在行 1436 算大气倾向）
zdqsdif_ssi_tot(:, :) 地下冰总交换通量诊断（write_output 'flux_ssice'，行 1400）
surf_h2o_lh 地表水潜热通量诊断

zq 更新后在行 1436 经 pdqdif=(zq-(pq+pdqfi·dt))/dt 转成大气水汽扩散倾向输出。

复现注意事项

待确认