hdo_surfex_mod.F

路径

LMDZ.MARS\libf\phymars\hdo_surfex_mod.F

所属目录 / 模块

libf\phymars

文件定位

hdo_surfex_mod.F 定义 hdo_surfex_mod 模块和 hdo_surfex 子程序，是水循环中 HDO（半重水）地表-大气通量计算的核心例程。它被 vdifc_mod.F 中的 vdifc 垂直湍流扩散例程调用，在 H2O 水汽地表交换段（行 1350）之后执行，基于已计算的 H2O 通量和地表冰储量计算 HDO 通量。

本例程的核心功能是：在升华/凝结过程中考虑同位素分馏效应，计算 HDO 的地表通量。它使用平衡分馏系数（Lamb 公式）和动力学分馏系数（Jouzel & Merlivat 1984），并处理永久冰盖（watercaptag）的特殊情况。

源码注释表明作者为 L. Rossi 和 M. Vals（2019）。

定义的符号

符号	类型	行号	作用
`hdo_surfex_mod`	module	1	封装 HDO 地表通量计算例程
`hdo_surfex`	subroutine	7	基于 H2O 通量计算 HDO 地表通量，含分馏效应

依赖的模块

use 模块	only 列表	用途	待确认
`tracer_mod`	`igcm_h2o_vap`, `igcm_h2o_ice`, `igcm_hdo_vap`, `igcm_hdo_ice`, `qparentmin`	tracer 索引和最小值阈值	-
`surfdat_h`	`watercaptag`	永久冰盖标记	-
`geometry_mod`	`longitude_deg`, `latitude_deg`	经纬度（本例程未使用）	-
`comcstfi_h`	`pi`	圆周率，用于扩散系数公式	-
`microphys_h`	`nav`, `kbz`, `mh2o`, `mco2`, `mhdo`, `molco2`, `molh2o`, `molhdo`	物理常数：阿伏伽德罗常数、玻尔兹曼常数、分子质量	-
`write_output_mod`	`write_output`	诊断输出（本例程已注释掉）	-
`callkeys_mod`	`hdofrac`	HDO 分馏开关	-

调用的关键例程

被调用例程	所在模块 / 文件	调用位置	作用
无	-	-	本例程只使用内在函数，不调用其他用户例程

输入

输入	来源	类型/维度	单位	含义
`ngrid`	`vdifc`	integer scalar	grid points	大气列数
`nlay`	`vdifc`	integer scalar	layers	大气层数（本例程只使用第 1 层）
`nq`	`vdifc`	integer scalar	-	tracer 总数
`ptimestep`	`vdifc`	real scalar	s	物理时间步长
`zt`	`vdifc`	real `(ngrid,nlay)`	K	温度（本例程只使用 `zt(:,1)`）
`pplay`	`vdifc`	real `(ngrid,nlay)`	Pa	层压力（本例程只使用 `pplay(:,1)`）
`zq`	`vdifc`	real `(ngrid,nlay,nq)`	kg/kg	tracer 混合比（本例程只使用 `zq(:,1,igcm_hdo_vap)`）
`pqsurf`	`vdifc`	real `(ngrid,nq)`	kg/m²	地表 tracer 储量
`old_h2o_vap`	`vdifc`	real `(ngrid)`	kg/kg	子时间步前的 H2O 水汽混合比
`qsat`	`vdifc`	real `(ngrid)`	kg/kg	饱和混合比
`pdqsdif`	`vdifc`	real `(ngrid,nq)`, intent(inout)	kg/kg/s	tracer 地表通量 tendency（输入 H2O 通量，输出 HDO 通量）
`dwatercap_dif`	`vdifc`	real `(ngrid)`	kg/kg/s	永久冰盖升华/凝结 tendency

输出

输出	去向	类型/维度	单位	含义
`hdoflux`	`vdifc`	real `(ngrid)`	kg/kg/s	HDO 水汽地表通量
`pdqsdif(igcm_hdo_ice)`	`vdifc`	real `(ngrid)`, intent(inout)	kg/kg/s	HDO 冰 tracer 地表通量 tendency

共享状态与副作用

本文件没有 save 变量、COMMON、THREADPRIVATE、文件 I/O 或配置读取。
副作用是原地改写 pdqsdif(:,igcm_hdo_ice)（输入时为 H2O 通量，输出时为 HDO 通量）。
诊断输出已注释掉（write_output 调用），不影响运行。
geometry_mod 中的 longitude_deg 和 latitude_deg 被 use 但未在例程中使用。

核心逻辑

初始化：
- 计算饱和比 satu = old_h2o_vap / qsat
- 初始化分馏系数 alpha = 1.0、alpha_c = 1.0
主循环（对每个格点 ig）：
- 计算总 H2O 通量：h2oflux = pdqsdif(igcm_h2o_ice) + dwatercap_dif
升华分支（h2oflux <= 0）：
- 若地表冰储量 pqsurf(igcm_h2o_ice) > qparentmin：
  - HDO 通量 = H2O 通量 × (HDO 冰储量 / H2O 冰储量)
- 否则：HDO 通量 = 0
- 约束：不能超过地表冰储量（max 操作）
- 若在永久冰盖上且升华超过地表冰储量：
  - 添加额外贡献：dwatercap_dif × 2 × 155.76e-6 × 5（D/H = 5 SMOW）
凝结分支（h2oflux > 0）：
- 若 hdofrac 为真（启用分馏）：
  - 计算 H2O 扩散系数 Dv（基于分子运动论）
  - 计算 HDO 扩散系数 Dv_hdo
  - 计算平衡分馏系数 alpha = exp(13525/T² - 0.0559)（Lamb 公式）
  - 计算动力学分馏系数 alpha_c（Jouzel & Merlivat 1984）
- 否则：alpha_c = 1.0
- 若 old_h2o_vap > qparentmin：
  - HDO 通量 = alpha_c × H2O 通量 × (HDO 水汽 / H2O 水汽)
- 否则：HDO 通量 = 0
- 若 hdofrac 为真：约束 HDO 通量不超过可用 HDO 水汽（min 操作）
输出：设置 hdoflux = pdqsdif(igcm_hdo_ice)

伪代码

hdo_surfex(ngrid, nlay, nq, ptimestep, zt, pplay, zq, pqsurf,
           old_h2o_vap, qsat, pdqsdif, dwatercap_dif, hdoflux):

    satu = old_h2o_vap / qsat  ! 饱和比
    alpha = 1.0                 ! 平衡分馏系数
    alpha_c = 1.0               ! 动力学分馏系数

    for ig in 1..ngrid:
        h2oflux = pdqsdif(igcm_h2o_ice) + dwatercap_dif

        if h2oflux <= 0:  ! 升华
            if pqsurf(igcm_h2o_ice) > qparentmin:
                pdqsdif(igcm_hdo_ice) = pdqsdif(igcm_h2o_ice) * (pqsurf(igcm_hdo_ice) / pqsurf(igcm_h2o_ice))
            else:
                pdqsdif(igcm_hdo_ice) = 0

            pdqsdif(igcm_hdo_ice) = max(pdqsdif(igcm_hdo_ice), -pqsurf(igcm_hdo_ice) / ptimestep)
            hdoflux = pdqsdif(igcm_hdo_ice)

            if watercaptag and (-h2oflux * ptimestep) > pqsurf(igcm_h2o_ice):
                ! 额外贡献：永久冰盖 D/H = 5 SMOW
                hdoflux += dwatercap_dif * 2 * 155.76e-6 * 5

        else:  ! 凝结
            if hdofrac:  ! 启用分馏
                Dv = f(zt, pplay, mh2o, molco2, molh2o)  ! H2O 扩散系数
                Dv_hdo = f(zt, pplay, mhdo, molco2, molhdo)  ! HDO 扩散系数
                alpha = exp(13525 / T² - 0.0559)  ! Lamb 公式
                alpha_c = (alpha * satu) / (alpha * (Dv/Dv_hdo) * (satu - 1) + 1)
            else:
                alpha_c = 1.0

            if old_h2o_vap > qparentmin:
                pdqsdif(igcm_hdo_ice) = alpha_c * pdqsdif(igcm_h2o_ice) * (zq(igcm_hdo_vap) / old_h2o_vap)
            else:
                pdqsdif(igcm_hdo_ice) = 0

            if hdofrac:
                pdqsdif(igcm_hdo_ice) = min(pdqsdif(igcm_hdo_ice), zq(igcm_hdo_vap) / ptimestep)

            hdoflux = pdqsdif(igcm_hdo_ice)

参与的主题流程

主题	参与方式
水循环 / HDO 同位素	`vdifc` 中 H2O 水汽地表交换段之后调用本例程，计算 HDO 的地表通量，用于追踪水循环中的同位素分馏
垂直湍流扩散	`vdifc_mod.F` 在处理 `h2o_vap` tracer 时调用本例程，作为水汽地表边界条件的一部分
永久冰盖	当 `watercaptag(ig)` 为真且升华超过地表冰储量时，从永久冰盖补充 HDO（D/H = 5 SMOW）

写法特点

固定格式 Fortran（.F 扩展名），使用传统 Fortran 77 风格的注释和续行符 &。
只使用大气第 1 层（近地层）的温度、压力和 HDO 水汽混合比。
pdqsdif 是 intent(inout)，输入时包含 H2O 通量，输出时被覆盖为 HDO 通量。
geometry_mod 中的经纬度被 use 但未在例程中使用，可能是历史遗留或调试用途。
诊断输出（write_output 调用）已被注释掉。
扩散系数公式基于分子运动论，考虑了 CO2-H2O 和 CO2-HDO 的二元扩散。
分馏系数公式使用 Lamb 参数化（替代了注释掉的 Merlivat 1984 公式）。
永久冰盖的 D/H 比值硬编码为 5 SMOW（2 × 155.76e-6 × 5）。

复现要点

必须确保 tracer_mod 中的 tracer 索引已正确初始化，特别是 igcm_hdo_vap 和 igcm_hdo_ice。
qsat 必须在调用前由 watersat 或等效例程计算。
pdqsdif(:,igcm_h2o_ice) 必须在调用前包含 H2O 冰通量信息。
qparentmin 是 tracer 最小值阈值，用于避免除零。
若 hdofrac 为假，所有分馏效应被禁用，HDO 通量直接正比于 H2O 通量和 HDO/H2O 比值。
永久冰盖逻辑依赖 watercaptag 标记和 dwatercap_dif 输入。

待确认

geometry_mod 中的经纬度被 use 但未使用，需确认是否为历史遗留或计划用于诊断输出。
注释掉的 write_output 调用是否在调试版本中启用。
永久冰盖 D/H = 5 SMOW 的假设是否适用于所有冰盖区域。
Lamb 分馏系数公式（exp(13525/T² - 0.0559)）的适用温度范围。
扩散系数公式中分子质量的组合方式（(molco2+molh2o)²）是否为标准二元扩散公式。