soil_settings.F

路径

LMDZ.MARS\libf\phymars\soil_settings.F

所属目录/模块

libf/phymars

文件定位

soil_settings.F 是土壤初始化设置例程。它从 startfi.nc 重启文件中读取土壤层深度几何、热惯量、土壤温度和地热通量等字段，并在层数不匹配时执行垂直插值。由 phyetat0_mod 在 GCM 初始化阶段调用一次。

定义的符号

依赖的模块

调用的关键例程

输入

输出

同时写入 comsoil_h 中的：layer/mlayer（层深度）、inertiedat/inertiesoil（热惯量）、volcapa（体积热容量）、flux_geo（地热通量）。

共享状态与副作用

• 写入 comsoil_h module 变量：layer/mlayer（层深度几何）、inertiedat/inertiesoil（热惯量）、volcapa（体积热容量）、flux_geo（地热通量）
• 写入 qsoil：土壤 tracer 数组（h2o_vap_soil, h2o_ice_soil, h2o_vap_ads）
• 日志输出：大量 write(*,*) 打印读取状态和字段范围

核心逻辑

1. 深度坐标（行 100–168）

1. 查询层数（行 106）：从 startfi.nc 读取 subsurface_layers 维度长度 dimlen。
2. 层数不匹配（行 108–123）：dimlen ≠ nsoil 时设 interpol=.true.，分配 oldmlayer。
3. 旧格式检测（行 128–153）：查询 inertiedat 维度数 ndims：
   • ndims=1：旧格式（热惯量仅地表），设 olddepthdef=.true.，用公式 oldmlayer(k) = sqrt(887.75/π) * (2^(k-0.5) - 1) 构建旧深度。
   • ndims≠1：新格式，读取 soildepth 坐标到 mlayer 或 oldmlayer。
4. 构建新深度（行 156–162）：interpol=.true. 时用公式 mlayer(k) = lay1*(1+k^2.9*(1-exp(-k/20)))，lay1=2e-4。
5. 构建层界面（行 165–168）：layer(k) = (mlayer(k)+mlayer(k-1))/2，最后一层外推。

2. 体积热容量（行 170–183）

• volcapa 由 tabfi.F 在读取 controle 表时设置。
• 若 volcapa ≤ 0，设为默认值 1e6 J/m³/K。

3. 热惯量（行 185–294）

3.1 当日气候热惯量 inertiedat（行 190–243）：

• ndims=1（旧格式）：读取地表值 surfinertia，复制到所有层。
• ndims≠1（新格式）：直接读取 3D 字段，或存入 oldinertiedat 待插值。

3.2 PEM 热惯量 inertiesoil（行 245–294）：

• 字段不存在时：从 inertiedat 复制。
• 字段存在时：直接读取，或存入 oldinertiesoil 待插值。

4. 土壤温度（行 298–349）

• 字段不存在时：从 tsurf 复制到所有层。
• 字段存在时：直接读取，或存入 oldtsoil 待插值。

5. 垂直插值（行 352–432）

• olddepthdef 路径（行 355–380）：旧格式深度用 inertiesoil/volcapa 缩放，插值 tsoil。
• interpol 路径（行 381–432）：新格式层数不匹配时，用 interp_line 对 inertiedat/inertiesoil/tsoil 分别插值到新网格。

6. 地热通量（行 436–445）

• 读取 flux_geo；不存在时设为 0。

7. 吸附 tracer（行 448–501）

• adsorption_soil=.true. 时读取 h2o_vap_soil、h2o_ice_soil、h2o_vap_ads；不存在时设为 0。

8. 报告（行 505–523）

• 打印 volcapa、inertiedat/inertiesoil/tsoil 的最小最大值。

伪代码

dimlen = inquire("subsurface_layers")
if dimlen ≠ nsoil: interpol = true

ndims = inquire_field_ndims("inertiedat")
if ndims == 1:   ! 旧格式
  olddepthdef = true
  oldmlayer = sqrt(887.75/π)*(2^(k-0.5)-1)
else:
  read "soildepth" → mlayer or oldmlayer

if interpol: mlayer = lay1*(1+k^2.9*(1-exp(-k/20)))
layer(k) = (mlayer(k)+mlayer(k-1))/2

if volcapa ≤ 0: volcapa = 1e6

read inertiedat → inertiedat or oldinertiedat
if inertiesoil exists: read → inertiesoil or oldinertiesoil
else: inertiesoil = inertiedat

if tsoil exists: read → tsoil or oldtsoil
else: tsoil = tsurf (all layers)

if olddepthdef or interpol:
  interp_line(oldgrid, oldval, mlayer, newval)  ! 对每个字段

read flux_geo (default 0)
if adsorption_soil: read h2o_vap_soil, h2o_ice_soil, h2o_vap_ads (default 0)

report min/max of volcapa, inertiedat, inertiesoil, tsoil

参与的主题流程

写法特点

• 固定格式 Fortran：使用 c 列注释、DO/ENDDO 大写。
• 旧格式兼容：ndims=1 检测旧格式（热惯量仅地表），用公式重建深度坐标。
• 垂直插值：层数不匹配时用 interp_line 线性插值，支持 olddepthdef（旧格式深度缩放）和 interpol（新格式直接插值）两种路径。
• 大量日志：每个字段读取都打印状态和范围。
• phyetat0_mod 唯一调用方：仅在 GCM 初始化时调用一次。

复现要点

• startfi.nc 必须包含 subsurface_layers 维度、soildepth 坐标、inertiedat 字段。
• inertiesoil 和 tsoil 可选；不存在时分别从 inertiedat 和 tsurf 初始化。
• flux_geo 可选；不存在时设为 0。
• volcapa 由 tabfi.F 读取 controle 表设置；若为 0 则用默认值 1e6。
• adsorption_soil 开关控制是否读取地下 h2o tracer。

待确认

• 行 138 oldmlayer(k) = sqrt(887.75/π) * (2^(k-0.5)-1) 中 887.75 的来源（待确认：推断为火星热扩散特征参数，但无文献引用）。
• 行 368 oldgrid(2:) = oldmlayer * (inertiesoil/volcapa) 中深度用 inertie/volcapa 缩放的物理含义（待确认：推断为扩散长度尺度 sqrt(κ*t) 的简化）。
• 日志中 phyetat0: 前缀（行 455 等）来自调用方名称，推断为复制粘贴遗留。

相关页面

• comsoil_h — 土壤共享参数文件，本文件写入层几何、热惯量等
• soil — 土壤热扩散求解器，使用本文件设置的 layer/mlayer/volcapa
• waterice_tifeedback_mod — 水冰热惯量反馈，修改 inertiesoil
• iniwritesoil — 土壤诊断 NetCDF 初始化
• soil_settings.md：本页面。
• water-cycle-config — adsorption_soil 开关说明
• water-cycle — 水循环主题页