vdifc_mod.F

路径

LMDZ.MARS\libf\phymars\vdifc_mod.F

所属目录 / 模块

libf/phymars
MODULE vdifc_mod

文件定位

vdifc_mod.F 是 phymars 物理包的垂直扩散主模块。主例程 vdifc 在 physiq_mod.F 中由 calldifv 分支调用，接收当前累计的动力、热力和 tracer 物理倾向，把它们先加到工作态，再用隐式垂直扩散求解动量、热量和 tracer 的新增倾向。

这个文件不是单纯的湍流 K 系数例程。它把地表层交换系数、PBL K 系数方案、地表能量平衡、CO2 凝结质量变化、尘埃地表注入、H2O 地表/地下交换、HDO 同位素通量、坡面投影和诊断输出放在同一条时间步流程里。因此它是边界层、地表、尘埃循环和水循环之间的强耦合入口。

定义的符号

SAVE / THREADPRIVATE 状态

本文件没有定义 common block。除上述 SAVE 状态外，主要共享状态来自 use 的模块变量，例如 ustar/tstar、zdqsdif_ssi_tot 和 tracer 索引。

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• firstcall 首次进入时打印 CO2 凝结常数，并扫描 noms(iq) 找 co2 tracer。后续调用复用 ico2/m_co2/m_noco2/A/B/acond/bcond/ccond。
• ustar、tstar 来自 turb_mod，在两次 vdif_cd 之后被重置并写入当前步近地层摩擦速度和温度尺度。
• pq2 是 intent(inout)，由 yamada4、call_atke 或 vdif_kc 更新，供后续时间步和其他湍流路径使用。
• qsoil 是 intent(inout)，仅在 adsorption_soil 路径由 soilwater 修改地下水汽/冰/吸附水状态。
• zdqsdif_ssi_tot 来自 paleoclimate_mod，每次 vdifc 开头清零，水汽路径累计地下冰与大气/霜的总通量，随后写诊断。
• write_output 写出 surf_h2o_lh、zdqsdif_ssi_frost_tot、zdqsdif_ssi_atm_tot、flux_ssice 和每个坡面的 flux_ssice_slopeNN。
• 本文件不读写 restart 文件；restart 相关状态由调用链中的物理状态初始化/写出模块处理。

初始化 / firstcall / restart / diagnostic 行为

1. firstcall 段在行 255-283 执行一次：设置 CO2 凝结公式常数，打印参数，扫描 tracer 名称并为 CO2 gas 建立平均分子量换算系数。
2. 每个时间步行 296-317 清零所有输出倾向和水/地下冰诊断缓存，避免沿用上一步工作数组。
3. 行 350-388 根据 callcond 和 ico2 计算 CO2 体积分数与凝结温度；未启用或无 CO2 tracer 时采用 vmr_co2=0.95。
4. 本文件没有直接 restart 读写，但会修改 pq2、qsoil 和 zdqsdif_ssi_tot 等可能由其他模块初始化或诊断的共享状态。
5. 诊断输出集中在行 1388-1407，只输出水潜热和地下冰交换相关字段；HDO flux 诊断在源码中保留为注释。

核心流程

1. 准备主坡面状态和隐式系数

行 285-288 从 major_slope(ig) 取主坡面的 ptsrf_tmp 和 pqsurf_tmp，供大气主列扩散、起尘和 K 系数路径使用。行 324-344 构造两类隐式系数：za=(pplev(l)-pplev(l+1))/g 表示层空气质量，zb0 是后续把 K 系数或地表导度转为隐式三对角系数的基量。地表层 zb0(:,1)=ptimestep*pplev(:,1)/(r*ptsrf_tmp)。

2. 把前序物理倾向加入工作态

行 395-405 计算

zu = pu + pdufi*dt
zv = pv + pdvfi*dt
zh = ph + pdhfi*dt
zq = pq + pdqfi*dt

后面求出的 pd*dif 都是相对于这个“已含前序物理过程”的工作态再增加的垂直扩散倾向。

3. 地表交换系数和 PBL K 系数

第一次 vdif_cd 在行 416-419 用原始风 pu/pv 计算各坡面 zcdv_true/zcdh_true。随后行 423-456 根据 callrichsl 是否包含 gustiness 修正，得到实际地表动量/热量导度 zcdv/zcdh，并写 ustar/tstar。

行 461-479 三选一计算垂直 K 系数：

• callyamada4：调用 yamada4，输出 zkv/zkh/zkq。
• callatke：调用 call_atke，输出 zkv/zkh，并把第 nlay+1 层界面复制为第 nlay 值。
• 否则：调用 vdif_kc，按传统 PBL K 系数方案更新 pq2/zkv/zkh。

行 480-488 是 1D doubleq 保护：ngrid==1 时内部层 zkh/zkv 下限设为 80。

4. 动量 U/V 隐式扩散

行 503-529 解 U，行 544-567 解 V。两者都使用同一 Thomas 算法样式：先从顶层向下递推 zc/zd，再从地表向上回代。地表边界使用 zb(:,1)=zcdv_tmp*zb0(:,1)，内部层使用 zb(:,2:nlay)=zkv(:,2:nlay)*zb0(:,2:nlay)。求解后 zu/zv 已包含摩擦后的风廓线。

行 579-619 随后用摩擦后的 zu/zv 第二次调用 vdif_cd。源码注释说明这套风用于起尘和升华，因此尘埃和水汽通量使用的是已经受地表摩擦影响的近地风背景。

5. 热力变量、地表能量平衡和 CO2 凝结

行 637-777 用 zkh/zcdh 对热力变量做隐式扩散，并把大气第一层和地表能量平衡联立。地表温度解包括：

• 地表热容量 pcapcal*ptsrf；
• 湍流热通量项 cpp*zb(1)*zc(1)；
• 线性化 Planck 辐射项 zdplanck*ptsrf，除非 tke_heat_flux 非零；
• 土壤/辐射传入的 pfluxsrf*ptimestep。

如果 callcond 且非 turb_resolved，行 666-681 会判断是否有层低于 CO2 凝结位温，必要时把 llnt(ig) 设为 300。行 698-771 在收敛循环中更新 zhs/zt，并在 co2clouds=.false. 时根据 zt<ztcond 计算 pdtc/zcondicea/dmice。若 co2clouds=.true.，源码当前把 zcondicea/dmice/pdtc 置零。

多坡面情况下，行 779-801 对非主坡面只重算地表温度倾向 pdtsrf(ig,islope)，大气热力廓线仍来自主坡面隐式解。

6. 尘埃和普通 tracer

行 808-889 先计算近地密度 rho=zb0(:,1)/ptimestep，再按 lifting/doubleq/submicron/dustinjection/rdstorm 选择地表尘埃通量：

• lifting=.false. 时 pdqsdif_tmp=0。
• doubleq 或 submicron 的旧路径直接用 -alpha_lift(...) 或 -dustliftday*alpha_lift(...) 生成 mass/number/submicron/stormdust 通量。
• dustinjection=1 还要求 local_time 在 ti_injection_sol 与 tf_injection_sol 之间；CO2 地表冰 pqsurf_tmp(:,igcm_co2)>=1 时不抬升。
• 默认路径调用 dust_windstress_lift，由风应力阈值方案计算 pdqsdif_tmp。

行 900-957 对“非水汽”tracer 做隐式垂直混合。普通 tracer 在地表层加入 (-pdqsdif_tmp)*ptimestep，并把 pdqsdif_tmp*cos(slope) 写入每个坡面的 pdqsdif。h2o_ice 和 hdo_ice 在这里走无地表注入的特例；真正的 H2O 地表冰增减由后面的 h2o_vap 专段覆盖。

7. H2O 地表、地下和坡面交换

行 959-1315 是 H2O 水汽专段，已有更细的 section 页：vdifc-water-surface-exchange。主文件页保留复现主线：

1. 每个坡面先把 pqsurf(h2o_ice) 和 watercap 除以 cos(slope)，转成真实坡面面积上的储量。
2. make_tsub 根据 pdtsrf、地表冰量、watercaptag 和 dtmax=0.5 K 选择子时间步数。
3. 每个子步用 zcdv 或 zcdh 作为地表导度，并乘 dryness，对水汽做隐式扩散，地表边界拉向 watersat(ztsrf,pplev(:,1)) 的 qsat。
4. adsorption_soil 时调用 soilwater，把地下水汽、吸附水、地下冰与大气边界条件耦合；否则在 lag_layer 和 h2oice_depth>0 时走地下冰-大气或地下冰-霜通量。
5. latentheat_surfwater 时用源码中的水冰升华潜热公式修正 pdtsrf；有冰时温度上升被封顶到 To。
6. 子步结束后写 pdqsdif(h2o_ice,islope)、dwatercap_dif、zdqsdif_ssi_atm_tot/frost_tot/ssi_tot，再把坡面水汽廓线按 subslope_dist 平均回大气主列。

这个段落是地表/水循环核心耦合点：它同时改变地表水冰、永久水冰盖、地下水相关 tracer、近地水汽、地表温度倾向和地下冰诊断通量。

8. HDO 地表交换

行 1317-1380 在 hdo 且 iq==igcm_hdo_vap 时执行。它先对 HDO vapor 建立与普通 tracer 类似的垂直混合系数，再逐坡面调用 hdo_surfex。调用输入包括 H2O 专段保存的 saved_h2o_vap、更新后地表温度的 qsat_tmp、H2O/HDO 地表通量数组以及 dwatercap_dif。hdo_surfex 返回坡面 hdoflux，随后按 subslope_dist 平均成 hdoflux_meshavg 并作为地表通量边界条件更新 HDO vapor 廓线。

9. 诊断和最终倾向

行 1388-1407 写水潜热和地下冰交换诊断。行 1424-1439 把工作态差分成输出倾向：

pdudif = (zu - (pu+pdufi*dt))/dt
pdvdif = (zv - (pv+pdvfi*dt))/dt
pdhdif = (zhs - (ph+pdhfi*dt + latcond*dmice/cpp/ppopsk))/dt
pdqdif = (zq - (pq+pdqfi*dt))/dt

physiq_mod.F 的调用点随后把这些倾向累加到全局物理倾向，并把 dwatercap_dif 加到 dwatercap。

伪代码

vdifc:
  if firstcall:
    initialize CO2 condensation constants
    find co2 tracer index in noms(:)

  select major-slope surface state
  zero all output tendencies and water/SSI diagnostics
  build za and zb0 implicit-diffusion coefficients
  compute CO2 vmr and condensation temperature if callcond

  zu,zv,zh,zq = state + incoming physical tendencies * dt

  call vdif_cd with original pu,pv
  derive zcdv/zcdh, ustar/tstar on major slope
  choose PBL K path:
    if callyamada4: call yamada4
    else if callatke: call call_atke
    else: call vdif_kc
  if doubleq and 1D: enforce zkv/zkh >= 80 inside the column

  solve implicit U diffusion using zkv and surface zcdv
  solve implicit V diffusion using zkv and surface zcdv

  call vdif_cd again with friction-modified zu,zv
  recompute conductances for lifting and sublimation

  solve implicit thermal diffusion plus surface energy balance
  if callcond: iterate CO2 condensation mass correction unless co2clouds path zeros it
  solve non-major slope surface temperature tendencies

  if lifting:
    compute dust/stormdust/submicron surface flux by configured branch
  else:
    dust flux = 0

  for each tracer:
    if normal tracer branch:
      solve implicit tracer diffusion
      apply dust-like surface source except h2o_ice/hdo_ice special cases
      write pdqsdif across slopes

    if water and tracer is h2o_vap:
      for each slope and grid point:
        choose H2O subtimestep count with make_tsub
        for each subtimestep:
          solve vapor diffusion with qsat surface boundary
          optionally call soilwater
          optionally compute lag-layer SSI atmosphere/frost fluxes
          apply latent-heat surface temperature feedback
          update vapor profile, surface ice and SSI diagnostics
        integrate pdtsrf, pdqsdif(h2o_ice), dwatercap_dif
      average slope vapor profiles back to grid-cell atmosphere

    if hdo and tracer is hdo_vap:
      call hdo_surfex per slope
      solve HDO vapor diffusion with mesh-averaged HDO surface flux

  write water/SSI diagnostics
  convert final zu,zv,zhs,zq work states to pdudif,pdvdif,pdhdif,pdqdif

参与的主题流程

写法特点

• 固定格式 Fortran，主例程很长，靠注释编号组织流程；复现时更适合按执行段拆分而不是按变量声明阅读。
• 隐式扩散反复使用同一组 za/zb/zc/zd 工作数组；每个物理量在进入求解前都会重设地表或内部层的 zb。
• vdif_cd 调用两次：第一次服务于 K 系数和初始地表层诊断，第二次用摩擦后风服务于起尘和升华。
• 水汽段使用自适应子时间步，并在坡面真实面积和水平投影面积之间来回转换；这对质量收支很关键。
• firstcall、SAVE、THREADPRIVATE 只用于 CO2 凝结常数和 tracer 索引缓存，不用于水/尘埃状态累积。
• pdqsdif 的符号约定需按源码使用：尘埃注入路径使用负的 alpha_lift 或 dustliftday*alpha_lift 写入，H2O ice 倾向由 zqsurf 的差分得到。

复现要点

• 先把前序物理倾向加到 zu/zv/zh/zq，最后再相对同一工作基准差分输出；不要直接相对原始 pu/pv/ph/pq 计算扩散倾向。
• major_slope 决定大气主列的热力和 K 系数上下文；水冰和地表温度仍逐坡面处理。
• 动量和热量扩散使用不同地表导度：动量用 zcdv，热量/水汽新方案多处用 zcdh。
• callyamada4、callatke、fallback vdif_kc 是互斥路径；doubleq && ngrid==1 的 80 下限在三种路径之后统一施加。
• H2O 水汽地表交换必须保留 make_tsub 子步、dryness 缩放、old_wsublimation_scheme 的 zcdv/zcdh 选择、watercaptag 的 watercap 透支逻辑和坡面 cos 投影。
• adsorption_soil 与非吸附 lag-layer 地下冰路径是不同实现；不要同时套用两套地下通量公式。
• hdo_surfex 必须在 H2O 专段之后执行，因为它依赖 saved_h2o_vap、更新后的 qsat_tmp、H2O 通量和 dwatercap_dif。
• 诊断 zdqsdif_ssi_tot 每步开头清零，只有当前步水汽段累计值；它不是 restart 累积量。

待确认

• 行 900-902 的普通 tracer 守卫写成 (.not.water).or.(.not.iq.eq.igcm_h2o_vap).or.(.not.iq.eq.igcm_hdo_vap)。由于一个 iq 不可能同时等于 H2O vapor 和 HDO vapor，这个条件按字面几乎恒真。已有 vdifc-water-surface-exchange 标为复现风险；本页按源码行为记录，是否应为 .and. 需对照上游版本确认。
• co2clouds=.true. 分支行 763-768 将 zcondicea/dmice/pdtc 置零，注释中保留了 pcondicea_co2microp 的使用痕迹；该路径的物理意图需结合 CO2 云主流程确认。
• hfmax 作为 vdifc 入参声明但本文件未见实际使用；是否为历史接口保留待确认。
• compute_Tice 中 indexice 是 integer，但行 1503 赋值 0.；Fortran 会转换为 0，是否为笔误不影响当前源码字面复现。

相关页面

• physiq_mod - physiq 中 calldifv 调用 vdifc 并累加输出倾向。
• vdif_cd_mod - 地表层 Cd/Ch 交换系数。
• vdif_kc - fallback K 系数方案。
• yamada4 - Mellor-Yamada K 系数路径。
• lmdz_call_atke - ATKE K 系数路径。
• dust_windstress_lift - 风应力起尘默认路径。
• soilwater - 地下水汽/吸附水/地下冰求解器。
• hdo_surfex_mod - HDO 地表-大气通量。
• paleoclimate_mod - lag-layer 地下冰共享状态。
• vdifc-water-surface-exchange - H2O 地表交换段细页。
• dust-cycle, water-cycle, radiation