LMDZ.MARS Wiki

soilwater.F90

路径

LMDZ.MARS\libf\phymars\soilwater.F90

所属目录/模块

libf/phymars

文件定位

地下水/吸附水/地下冰储库的一维垂直求解器。计算 regolith（风化层）内从地表到约 18 m 深处（layer(nsoil)）的水汽、吸附水、地下冰三相剖面，并返回地下与大气之间的水汽通量 zdqsdifrego，作为 vdifc 计算大气水汽垂直廓线的下边界条件。同时计算地表冰与地下的交换。

只在 callphys.def 中 adsorption_soil=.true. 且 water=.true. 时由 vdifc 调用。需要在 traceur.def 定义三个地下 tracer：h2o_vap_soil（地下水汽）、h2o_ice_soil（地下冰）、h2o_ads_vap（吸附水）。

本文件是单个自由格式的 subroutine（无 module 封装），作者 Pierre-Yves Meslin，Lucas Lange 清理。

定义的符号

符号 类型 行号 作用
soilwater subroutine 1 地下水汽/吸附水/冰的吸附-扩散-相变一维求解器

依赖的模块

use 模块 only 列表 用途 待确认
comsoil_h igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg, ads_const_D, ads_massive_ice 地下 tracer 索引、层深、吸附方案开关、孔隙率
comcstfi_h pi, r 圆周率、干空气气体常数
tracer_mod igcm_h2o_vap 大气水汽 tracer 索引
surfdat_h watercaptag 永久极冠标记（极冠点跳过 regolith 计算）
geometry_mod cell_area, latitude_deg 网格面积（全球质量统计）、纬度 待确认：latitude_deg 在主体未见使用
write_output_mod write_output 诊断输出（#ifndef MESOSCALE）

调用的关键例程

被调用例程 所在模块/文件 调用位置 作用
write_output write_output_mod 行 1218-1226 输出地下通量、冰饱和度、水/冰柱质量等诊断字段

本例程不调用其他物理子程序，是一个自包含的数值求解器。

输入

输入 来源 类型/维度 单位 含义
ngrid vdifc INT - 格点数（调用时为 1，逐点调用）
nlayer vdifc INT - 大气层数
nq vdifc INT - 大气 tracer 数
nsoil vdifc INT - 地下层数
nqsoil vdifc INT - 地下 tracer 数（=3）
ptsrf vdifc (ztsrf) REAL (ngrid) K 地表温度
ptsoil vdifc (ptsoil(ig,:,islope)) REAL (ngrid,nsoil) K 地下各层温度
ptimestep vdifc (subtimestep) REAL s 子时间步长
exchange vdifc LOGICAL (ngrid) - 当前步是否与大气交换水汽
qsat_surf vdifc (qsat) REAL (ngrid) kg/kg 地表温度下的饱和混合比
pq vdifc (zq_tmp_vap) REAL (ngrid,nlayer,nq) kg/kg 大气 tracer 混合比
pa,pb,pc,pd vdifc (za/zb/zc/zd) REAL (ngrid,nlayer) - vdifc 隐式扩散三对角系数
pdqsdifpot vdifc (zdqsdif_surf) REAL (ngrid) kg/m²/s 不含地下交换的潜在地表通量
pplev vdifc REAL (ngrid,nlayer+1) Pa 大气压力层界
rhoatmo vdifc (rho) REAL (ngrid) kg/m³ 第一层大气密度（当前未实际使用）
writeoutput vdifc LOGICAL - 是否为最后子步（触发输出）
pqsurf vdifc (zqsurf) REAL (ngrid) kg/m² 地表水冰量（intent(in)）

输出

输出 去向 类型/维度 单位 含义
pqsoil vdifc (qsoil) REAL (ngrid,nsoil,nqsoil) 见下 地下三相 tracer（inout）：水汽（kg/m³孔隙气）、冰、吸附水（kg/m³regolith）
zdqsdifrego vdifc (zdqsdif_regolith) REAL (ngrid) kg/m²/s 地下→大气水汽通量（正=向外）
zq1temp2 vdifc (zq1temp_regolith) REAL (ngrid) kg/kg 交换后地表上方临时水汽混合比
saturation_water_ice vdifc (pore_icefraction) REAL*8 (ngrid,nsoil) - 孔隙水冰饱和度（孔隙被冰填充的体积分数）

共享状态与副作用

核心逻辑

  1. 分配 + 首次初始化（行 242-364）：分配 SAVE 数组；逐点逐层设置 midlayer_dz/interlayer_dz、porosity_ice_free（ads_massive_ice 时冰超 tol_massiveice=50 kg/m³ 的格点设孔隙率 0.999999）、tortuosity=1.5、rho_soil=1300、meshsize=5e-6、D0=porosity/tortuosity；从 pqsoil 读入三相初值；算 saturation_water_ice 和 porosity。

  2. massive ice 孔隙率更新（行 366-377）：每步根据当前冰量重置 porosity_ice_free。

  3. 扩散系数（行 382-458，跳过 watercaptag 极冠点）：

    • saturation_water_ice = ice/(rho_H2O_ice*porosity_ice_free)，封顶 0.999。
    • 热速度 vth、H2O-CO2 碰撞积分 omega（Mellon & Jakosky 1993）、分子扩散 Dm、Knudsen 扩散 Dk。
    • 中层 D_mid = D0*(1-Sw)²/(1/Dm+1/Dk)（Meslin 2010 饱和度依赖）；ads_const_D 时用 default_diffcoeff=4e-4。
    • 层间 D_inter 由中层值插值；闭合标志 close_* 时层间饱和度设 0.999。

  4. 吸附/解吸系数（行 463-549）：按 choice_ads（1=固定 D0，2=2020 修正 DeltaQ，0=无吸附）算平衡吸附系数 k_ads_eq、吸附/解吸时间常数 Ka/Kd，及单层饱和后的第二段 Ka2/Kd2（双线性吸附等温线）。算 C/E/F/E2/F2 系数；保存上一步 znsoilprev/adswprev/iceprev；算总量 ztot1。

  5. 主求解循环（行 554-1098，逐点，跳过极冠）：

    • 算地表层 delta0/alpha0/beta0（耦合 vdifc 系数和地表冰 pqsurf）。
    • 升华外循环（while sublimation_flag，行 604）：
      • 单层饱和内循环（while recompute_all_cells_ads_flag，行 618）：自上而下算三对角系数 zdelta/zalpha/zbeta（有冰层固定为饱和 nsatsoil），底层显式解 znsoil(nsoil)，再自下而上回代解各层 znsoil 和 adswater_temp；检查每层是否越过单层饱和 adswater_sat_mono，越界则切换 A/B 系数段并标记整列重算。
      • 算层间通量 flux、dztot1；做凝结/升华：有冰层用 ice=ztot1+dztot1*dt-porosity*nsatsoil，冰耗尽则清零并重启升华循环；过饱和无冰层则置冰标志重启。
    • 地表通量（行 839-852）：exchange 时 zdqsdifrego 由 vdifc 系数耦合算出；否则按地表冰用 D_mid 梯度算。
    • 特例（行 860-1085）：无交换且地表冰不足以供应整个时间步时，改用通量边界条件，把剩余地表冰全部升华转入地下（重跑一遍带 znsoilprev2 的求解）。
    • 更新 ztot1/h2otot。

  6. choke 闭合检测（行 1100-1138）：冰饱和度超 choke_fraction=0.8 时，按层间通量方向设 close_top/close_bottom（阻止进一步扩散）；冰量回落则重开。

  7. 质量统计 + 写回（行 1140-1186）：算每点水柱/冰柱质量；把 znsoil/ice/adswater 写回 pqsoil；极冠点 saturation_water_ice=-1；exch 数值化。

  8. 全球总量 + 输出（行 1188-1228）：全球水/冰质量；n+1；重排 flux 后写诊断。

伪代码

soilwater(... pqsoil[inout], zdqsdifrego[out], zq1temp2[out], saturation_water_ice[out]):

  if .not.allocated: allocate SAVE arrays
  if firstcall_soil:
    设层间距、孔隙率、密度、meshsize、D0
    adswater_sat_mono = 2.6998e-7 * S * rho_soil
    从 pqsoil 读 znsoil/ice/adswater 初值

  if ads_massive_ice: 按冰量更新 porosity_ice_free

  for ig (跳过 watercaptag 极冠):
    for ik: Sw = ice/(rho_ice*porosity_free); porosity = porosity_free*(1-Sw)
    算 vth, omega, Dm, Dk, D_mid, D_inter (饱和度依赖扩散)
    按 choice_ads 算 Ka/Kd/Ka2/Kd2, k_ads_eq, C/E/F/E2/F2
    保存 znsoilprev/adswprev/iceprev; ztot1 = porosity*znsoil + ice

    算地表层 delta0/alpha0/beta0 (耦合 vdifc 的 pa..pd 和 pqsurf)

    while sublimation_flag:
      while recompute_all_cells_ads_flag:   // 单层饱和迭代
        自上而下: zdelta/zalpha/zbeta (有冰层固定 nsatsoil)
        底层显式 znsoil(nsoil)
        自下而上回代: znsoil(ik), adswater_temp(ik)
        if adswater_temp 越过/回落单层饱和: 切 A/B 段, 标记整列重算
      算 flux, dztot1
      凝结/升华: 有冰 ice = ztot1+dztot1*dt - porosity*nsatsoil
                 冰耗尽 -> 清零 + 重启; 过饱和无冰 -> 置冰标志 + 重启

    if exchange: zq1temp2 = beta0 + alpha0*znsoil(1)/rho
                 zdqsdifrego = porosity_free*pb(1)/dt*(zq1temp2 - znsoil(1)/rho)
    else:        zdqsdifrego = D_mid(1)/midlayer_dz(0)*(znsoil(1) - qsat_surf*rhosurf)

    if 无交换 且 地表冰不足供应整步:   // 特例
      zdqsdifrego = -(pqsurf + pdqsdifpot*dt)/dt
      重跑求解 (用 znsoilprev2), 把剩余地表冰转入地下

    更新 ztot1, h2otot

  choke 检测: Sw > 0.8 -> close_top/close_bottom 按通量方向
  写回 pqsoil = (znsoil, ice, adswater)
  算质量, 写诊断输出

参与的主题流程

主题 参与方式
水循环 (water-cycle) 地下储库：regolith 吸附水/地下冰与大气水汽交换，提供 vdifc 下边界通量
水循环主题页 主题入口

写法特点

复现要点

待确认

复现风险

相关页面