waterice_tifeedback_mod.F90

路径

LMDZ.MARS\libf\phymars\waterice_tifeedback_mod.F90

所属目录/模块

libf\phymars

文件定位

该文件定义 waterice_tifeedback_mod 模块，含唯一 subroutine waterice_tifeedback。职责是：根据地表水冰覆盖量或土壤孔隙冰填充率，重算一套土壤热惯量廓线 newtherm_i(ngrid,nsoil,nslope)，供 soil 例程在求解地温时替代默认热惯量 inertiedat。这是水冰对地表/地下热状态的“热惯量反馈”。

源码头注释说明两种机制（互斥，当前不能同时开启）：

• surfaceice_tifeedback：地表存在足够厚的块状水冰时，构造三层热惯量（冰层 + 过渡层 + 下伏风化层）。物理依据 Madeleine et al. (2009, Icarus)。
• poreice_tifeedback：孔隙填充冰，按 Siegler et al. (2012, JGR) 公式 sqrt(I_reg² + porosity*pore_filling*I_pureice²) 抬高热惯量。

历史：JBM 为纯冰编写（2008-2012），LL 于 2024 移入模块/改 .F90 并加入孔隙冰路径。

定义的符号

依赖的模块

调用的关键例程

调用方（核验）：

调用实参（核验 physiq_mod.F:737-739）：icecover = qsurf(:,igcm_h2o_ice,:)、poreice = pore_icefraction，输出 newtherm_i = inertiesoil_tifeedback，随后该热惯量传入 soil(...)（:741）。

输入

注意：ngrid/nsoil/nslope 及 icecover/poreice/newtherm_i 在声明区未写 intent（见“写法特点”），但按用途 icecover、poreice 为输入、newtherm_i 为输出。

输出

共享状态与副作用

• 无 common block，不写文件，无诊断打印。
• 所有 module 依赖均为只读（rho_ice、layer、inertiedat、porosity_reg、watercaptag、inert_h2o_ice、两个开关）。
• 无 SAVE 变量、无 firstcall；纯计算，给定输入即给定输出。
• 局部变量 inertie_purewaterice = 2100 在声明时初始化（纯水冰热惯量，仅孔隙冰路径用）。注意此值（2100）与地表冰 inert_h2o_ice（默认 2400）不同，前者硬编码、后者可配置。

核心逻辑

1. 先把 newtherm_i 全置 0。
2. 若 surfaceice_tifeedback 为真（地表块状冰路径），对每个坡面、每个格点：
   • 冰厚 icedepth = icecover/rho_ice（m）。
   • 若 icedepth >= layer(nsoil) 或 watercaptag(ig) 为真（冰极厚或永久冰盖）：整列 newtherm_i = inert_h2o_ice。
   • 否则若 icedepth < layer(1)（冰极薄，忽略）：整列 newtherm_i = inertiedat。
   • 否则（冰介于其间）：
     • 找冰/风化层边界 iref：第一个满足 layer(ik) <= icedepth < layer(ik+1) 的 ik+1。
     • 边界以上各层 ik=1..iref-1 设为 inert_h2o_ice（纯冰）。
     • 边界层 iref 用热传导串联公式取过渡值：sqrt( (layer(iref)-layer(iref-1)) / ( (icedepth-layer(iref-1))/inert_h2o_ice² + (layer(iref)-icedepth)/inertiedat(ig,ik)² ) )。
     • 边界以下各层 ik=iref+1..nsoil 保留 inertiedat。
3. 否则若 poreice_tifeedback 为真（孔隙冰路径），对每个坡面整层向量化： newtherm_i(:,:,islope) = sqrt(inertiedat(:,:)² + porosity_reg*poreice(:,:,islope)*inertie_purewaterice²)。
4. 两个开关都为假时，newtherm_i 保持全 0（调用方此时不会用它，见复现风险）。

伪代码

subroutine waterice_tifeedback(ngrid, nsoil, nslope, icecover, poreice, newtherm_i):
    newtherm_i = 0

    if surfaceice_tifeedback:
        for islope, ig:
            icedepth = icecover(ig,islope) / rho_ice
            if icedepth >= layer(nsoil) or watercaptag(ig):
                newtherm_i(ig,:,islope) = inert_h2o_ice           # 整列纯冰
            else if icedepth < layer(1):
                newtherm_i(ig,:,islope) = inertiedat(ig,:)        # 冰太薄，忽略
            else:
                iref = first ik+1 where layer(ik) <= icedepth < layer(ik+1)
                newtherm_i(ig, 1:iref-1, islope) = inert_h2o_ice  # 冰层
                newtherm_i(ig, iref, islope) =                    # 过渡层（串联）
                    sqrt( (layer(iref)-layer(iref-1)) /
                          ( (icedepth-layer(iref-1))/inert_h2o_ice^2
                          + (layer(iref)-icedepth)/inertiedat^2 ) )
                newtherm_i(ig, iref+1:nsoil, islope) = inertiedat # 下伏风化层

    else if poreice_tifeedback:
        for islope:
            newtherm_i(:,:,islope) =
                sqrt(inertiedat^2 + porosity_reg*poreice(:,:,islope)*inertie_purewaterice^2)

参与的主题流程

写法特点

• 自由格式 Fortran（.F90），续行行尾 &。
• 形参未声明 intent：ngrid/nsoil/nslope/icecover/poreice/newtherm_i 只用维度声明，无 intent(in/out)。复现时按用途判定方向。
• 局部 inertie_purewaterice=2100 与可配置的 inert_h2o_ice（默认 2400）是两个不同的热惯量值，分属孔隙冰路径与地表冰路径，勿混用。
• 过渡层公式是基于热传导的“串联电阻”式平均（按层内冰段/风化段的 1/I² 加权）。
• 两个反馈机制互斥（IF ... ELSE IF），源码注释明确“can not run with both options”。

复现要点

• 例程只算热惯量，不改任何 tracer/温度/质量；其输出 newtherm_i 经 soil 才影响地温。
• 触发条件：surfaceice_tifeedback 或 poreice_tifeedback 至少一个为真；二者互斥（同真时只走 surface 分支，因为 IF 在前）。调用方 physiq_mod.F:735 用 .or. 守卫，因此两个都假时不会调用本例程（newtherm_i 全 0 不会被使用）。
• surfaceice_tifeedback 还要求 water 与 callsoil 同时为真（核验 conf_phys.F:786-797 的中止检查）。
• 边界层过渡公式中 inertiedat(ig,ik) 的 ik 是紧邻其上的 DO ik=1,iref-1 循环（:115-117）正常结束后的残留值。按 Fortran 标准，正常结束的 DO 循环索引变量取 (iref-1)+1 = iref，因此该处实际等价于 inertiedat(ig,iref)，即边界层处的风化层热惯量。复现需照源码字面用 ik，结果与 iref 一致。
• icedepth 用 rho_ice（水冰密度）换算；layer 是土壤层深度数组。

待确认

• 过渡层公式 newtherm_i(ig,iref,islope)=sqrt(...inertiedat(ig,ik)...)（:119-120）中的 ik：此处 ik 的来源是紧邻其上的 DO ik=1,iref-1 循环（:115-117），不是边界查找循环。按 Fortran 标准，DO ik=1,iref-1 正常结束（未 EXIT）后循环变量取终止值 ik=iref。因此 inertiedat(ig,ik) 实际等于 inertiedat(ig,iref)，即边界层处的风化层热惯量，物理上合理。这是依赖"DO 循环结束后索引变量值"的隐式写法，可读性差但结果正确。推断：等价于显式写 inertiedat(ig,iref)；复现时可直接用 iref，与源码字面行为一致。
• 形参缺 intent 是否在某些编译器下有副作用；当前按惯例理解为 in/in/out。

相关页面

• soil.md：土壤热扩散求解，本例程输出的使用点（physiq_mod.F:741）。
• physiq 水循环调用点：本例程的调用上下文。
• water-saturation-helpers：水饱和基础公式。
• 水循环父级主题页。
• 待确认页：comsoil_h（定义 layer、inertiedat、porosity_reg）。
• surfdat_h：定义 inert_h2o_ice 和 watercaptag。