massflowrateco2.F90

路径

LMDZ.MARS\libf\phymars\massflowrateco2.F90

所属目录/模块

libf\phymars

文件定位

该文件定义 massflowrateco2_mod 模块，含 6 个 subroutine，构成一条计算链，目标是给出单个 CO2 冰粒的质量传输率 Ic（kg/s，>0 为增长/凝结，<0 为升华）。物理依据是 Listowski et al. (2014) 的扩散-热传导受限增长模型（eq.5/eq.6），2017 年由 C. Listowski 简化（用显式表面温度公式替代 Newton-Raphson 迭代）。

模块对外只暴露入口例程 massflowrateco2（唯一被外部 use 的符号，核验：improvedco2clouds_mod.F90:79）；其余 5 个例程是内部物性子程序，逐层计算混合气体热导率、扩散系数、N2 黏度等。

调用链（核验自源码 call 语句）：

massflowrateco2
└─ coefffunc                        (:61)
   ├─ KthMixNEW                     (:179)  CO2/N2 混合热导率
   │  ├─ KthCO2Scalab               (:348)  纯 CO2 热导率 (Scalabrin 2006)
   │  └─ KthN2LemJac                (:349)  纯 N2 热导率 (Lemmon & Jacobsen 2003)
   │     └─ viscoN2                 (:451)  N2 黏度
   └─ Diffcoeff                     (:181)  CO2/N2 扩散系数 (Reid 1987)

Ic 的唯一使用点是 improvedco2clouds（核验 :531）：把单粒速率乘以 CCNCO2 数浓度、时间步与 tauscaling 得到 dMice（每千克空气的 CO2 冰质量变化），再受“可用热能上限” facteurmax 与可用 CO2/CO2 冰量钳制，最后更新 co2_ice/co2 tracer 与潜热加热 subpdtcloudco2。

定义的符号

符号	类型	行号	作用
`massflowrateco2_mod`	module	1	容器模块
`massflowrateco2`	subroutine	7	入口：算单粒 CO2 冰质量传输率 `Ic`
`coefffunc`	subroutine	85	算 Listowski 模型系数 `kmix/Lsub/C0/C1/C2/Ak`
`Diffcoeff`	subroutine	215	CO2 在 N2 中的扩散系数（Reid et al. 1987）
`KthMixNEW`	subroutine	270	CO2/N2 混合气热导率（Mason & Saxena 1958，不用黏度）
`KthN2LemJac`	subroutine	365	纯 N2 热导率（Lemmon & Jacobsen 2003）
`viscoN2`	subroutine	481	纯 N2 黏度（Lemmon & Jacobsen 2003）
`KthCO2Scalab`	subroutine	541	纯 CO2 热导率（Scalabrin et al. 2006）

依赖的模块

use 模块	only 列表	用途（所在子程序）	待确认
`comcstfi_h`	`pi`	几何因子（`massflowrateco2`、`coefffunc`）	否
`tracer_mod`	`rho_ice_co2`	CO2 冰密度（kg·m⁻³），用于 Kelvin 因子 `Ak`（`coefffunc`）。注意：是常数版 `rho_ice_co2`，非温度依赖的 `density_co2_ice`，见复现风险	否
`microphys_h`	`kbz`	玻尔兹曼常数（1.381e-23）（`coefffunc`）	否
`microphys_h`	`mco2`	CO2 摩尔质量（44e-3 kg·mol⁻¹）（`coefffunc`、`Diffcoeff`、`KthMixNEW`）	否
`microphys_h`	`mn2`	N2 摩尔质量（28.01e-3 kg·mol⁻¹）（`Diffcoeff`、`KthMixNEW`、`KthN2LemJac`、`viscoN2`）	否
`microphys_h`	`nav`	阿伏伽德罗常数（6.023e23）（`coefffunc`）	否
`microphys_h`	`rgp`	理想气体常数（8.3143 J·mol⁻¹·K⁻¹）（`coefffunc`）	否
`microphys_h`	`sigco2`	CO2 冰/汽表面张力（0.08 J·m⁻²）（`coefffunc`）	否

各常数行号：microphys_h.F90:15(nav),17(rgp),19(kbz),25(mco2),27(mn2),69(sigco2)；rho_ice_co2 见 tracer_mod :27（real,save）。

调用的关键例程

被调用例程	所在模块/文件	调用位置	作用
`coefffunc`	本文件	`:61`	算所有中间系数
`KthMixNEW`	本文件	`:179`	混合气热导率 `kmix`
`Diffcoeff`	本文件	`:181`	扩散系数 `Dv`
`KthCO2Scalab`	本文件	`:348`	纯 CO2 热导率
`KthN2LemJac`	本文件	`:349`	纯 N2 热导率
`viscoN2`	本文件	`:451`	N2 黏度
`dlog/exp/sqrt/...`	Fortran 内建	多处	算术

本文件不调用其他源文件中的例程。

输入

入口 massflowrateco2(P,T,Sat,Radius,Matm,Ic)：

输入	来源	类型	单位	含义
`P`	调用方	real, intent(in)	Pa	层压力（`improvedco2clouds` 传 `pplay(ig,l)`）
`T`	调用方	real, intent(in)	K	温度（传 `zt(ig,l)`）
`Sat`	调用方	real(kind=8), intent(in)	—	CO2 过饱和比 S（传 `satu`）
`Radius`	调用方	double, intent(in)	m	CO2 冰粒半径（传 `riceco2(ig,l)`）
`Matm`	调用方	real, intent(in)	g·mol⁻¹	大气平均摩尔质量（传 `mmean(ig,l)`），注意单位是 g·mol⁻¹（见 `coefffunc` 内多处 ×1e-3 换算）

输出

输出	去向	类型	单位	含义
`Ic`	调用方	double, intent(out)	kg·s⁻¹	单个 CO2 冰粒的质量传输率，>0 增长、<0 升华

调用方对 Ic 的处理（核验 improvedco2clouds_mod.F90:534-553）：先判 isnan(Ic) 或 Ic==0 直接置 0；否则 dMice = Nccnco2 * Ic * microtimestep * tauscaling，再用 facteurmax（热能上限）和可用 CO2/CO2 冰量双向钳制，更新 tracer 和潜热加热率。

共享状态与副作用

无 common block，不写文件，无诊断打印。
全部依赖为只读 PARAMETER 常数 + 只读 rho_ice_co2（module SAVE 变量）。
各子程序无 SAVE 变量、无 firstcall，纯函数式（给定输入即给定输出）。
KthN2LemJac/KthCO2Scalab/viscoN2 把大量经验拟合系数硬编码为局部变量/PARAMETER。

核心逻辑

massflowrateco2（入口，:61-74）：

调 coefffunc 得 kmix（混合热导率）、Lsub（升华潜热）、C0/C1/C2、Ak（Kelvin 因子）。
显式表面温度：Tsurf = (1/C1)*ln(Sat/Ak) + T（Listowski 2017 简化，省去 Newton-Raphson，源码注释称误差 <0.6%）。
cond = 4π*Radius*kmix（热传导导率几何项）。
Ic = cond*(Tsurf - T)/Lsub：表面与环境温差驱动的质量流。

coefffunc（:152-204）：

平表面饱和压 psat = 1.382e12 * exp(-3182.48/T) (Pa)；pco2 = psat*S。
升华潜热 4 阶多项式 Lsub = l0 + l1*T + l2*T² + l3*T³ + l4*T⁴（Azreg-Aïnou 形式，J/kg）。
大气密度 rhoatm = P*Matm/(rgp*T)（先 g·m⁻³，再 ×1e-3 转 kg·m⁻³）。
调 KthMixNEW 得 kmix，调 Diffcoeff 得 Dv（cm²/s→×1e-4 转 m²/s）。
Fuchs-Sutugin（FS）非连续修正：对 Dv 和 kmix 各做一次。分别算 CO2 热速度 vthco2、空气热速度 vthatm，由此得 Knudsen 数与 lambda（Dahneke/Monchick & Mason 形式），把 Dv 和 kmix 除以 (1+lambda*Kn)。
装配系数：Ak = exp(2*sigco2*mco2/(rgp*rho_ice_co2*T*rc))（Kelvin 曲率因子）；C0/C1/C2 按 Listowski 2014 eq.6 给出。入口只用到 C1 和 Ak（C0/C2 计算但未被入口使用，见待确认）。

Diffcoeff（:249-257）：Reid et al. 1987 经验式，Diff = 0.00143*T^1.75/(Pbar*sqrt(Mab)*(dva^(1/3)+dvb^(1/3))²)，扩散体积 CO2=26.9、N2=18.5，输出 cm²/s。

KthMixNEW（:320-353）：Mason & Saxena (1958)/Wassiljeva (1904) 混合规则，不用黏度。算 CO2/N2 摩尔分数 x1/x2，按临界参数算 Gamma、平移热导率 lambda_trans、交叉系数 A12/A21，调 KthCO2Scalab、KthN2LemJac 得纯组分热导率，按 Wassiljeva 加权求和，末尾 ×1e-3（mW→W）。

KthN2LemJac（:406-469）：Lemmon & Jacobsen 2003 纯 N2 热导率，先调 viscoN2 得黏度，再用稀薄项 k1 + 残余项 k2（多项式 × delta^d × exp(-gamma*delta^l)）。

viscoN2（:522-528）：Lemmon & Jacobsen 2003 的 N2 黏度，碰撞积分 RGCS 用 ln(T*) 多项式拟合，输出 microPa·s。

KthCO2Scalab（:610-622）：Scalabrin et al. 2006 纯 CO2 热导率，对比态 Tr=T/Tc、rhor=rho/rhoc，多项式 k1 + 残余项 k2*exp(-5*rhor²)，× Lambdac，输出 mW·m⁻¹·K⁻¹。

伪代码

subroutine massflowrateco2(P, T, Sat, Radius, Matm, Ic):
    call coefffunc(P,T,Sat,Radius,Matm -> kmix,Lsub,C0,C1,C2,Ak)
    Tsurf = (1/C1)*ln(Sat/Ak) + T          # 显式表面温度 (Listowski 2017)
    cond  = 4*pi*Radius*kmix
    Ic    = cond*(Tsurf - T)/Lsub           # kg/s, >0 增长

subroutine coefffunc(P,T,S,rc,Matm -> kmix,Lsub,C0,C1,C2,Ak):
    psat = 1.382e12*exp(-3182.48/T);  pco2 = psat*S
    Lsub = poly4(T)                          # 升华潜热
    rhoatm = P*Matm/(rgp*T) * 1e-3           # kg/m3
    call KthMixNEW(-> kmix);  call Diffcoeff(-> Dv [m2/s])
    # Fuchs-Sutugin 修正 (对 Dv 与 kmix 各一次)
    Kn = ...;  lambda = (1.333+0.71/Kn)/(1+1/Kn)
    Dv   /= (1+lambda*Kn);  kmix /= (1+lambda*Kn)
    Ak = exp(2*sigco2*mco2/(rgp*rho_ice_co2*T*rc))
    C1 = Lsub*mco2/(rgp*T^2)                  # 入口实际使用 C1, Ak
    C0, C2 = ... (Listowski 2014 eq.6)

KthMixNEW: 调 KthCO2Scalab + KthN2LemJac，Wassiljeva 混合
KthN2LemJac: 调 viscoN2，Lemmon&Jacobsen 2003
KthCO2Scalab: Scalabrin 2006 对比态多项式
Diffcoeff: Reid 1987 经验式

参与的主题流程

主题	参与方式
CO2 循环 / CO2 云微物理	CO2 冰晶增长/升华的核心动力学。被 `improvedco2clouds` 在成核之后调用，决定每个 CCNCO2 粒子的质量变化速率，进而驱动 `co2_ice`↔︎`co2` 气相转移与潜热加热

写法特点

自由格式 Fortran（.F90），续行行尾 &，分隔注释块用大段 !===/!___。
混合精度：P、T、Matm 为 REAL，Sat、Radius、Ic 及全部中间量为 double/real(kind=8)，内部大量 dble() 提升。
Cpco2、Cpn2 用过时的 data 语句初始化（:147-148）。
大量硬编码经验拟合系数（Lemmon-Jacobsen N2、Scalabrin CO2、Reid 扩散体积等），来源在注释中标注了文献。
psat 系数 1.382e12、-3182.48 为平表面 CO2 饱和压拟合常数。
FS 修正中的 lambda 公式注释里作者自己存疑（“pas adaptee, Dahneke 1983? en fait si (Monchick&Black)”），保留原样。

复现要点

Ic 是单粒 kg/s，正负号即增长/升华；调用方再乘数浓度与时间步并钳制，复现勿混淆单粒速率与每千克空气速率。
Tsurf 用的是 Listowski 2017 显式简化式（非迭代），与 Listowski 2014 原式有 <0.6% 误差；复现时若用迭代解会有微小差异。
Kelvin 因子 Ak 用的是 tracer_mod 的常数 rho_ice_co2，不是温度依赖的 density_co2_ice。这是文件内的实际写法，复现需注意两者区别。
Matm 入口单位是 g·mol⁻¹（不是 kg·mol⁻¹）；coefffunc 内多处 ×1e-3 转 kg。
扩散系数与热导率都做了 FS 非连续修正，复现时不能省略，否则小粒子结果偏大。
各物性拟合式有各自适用温度/密度范围（火星条件），文件未做边界检查。

待确认

coefffunc 计算了 C0 和 C2 但入口 massflowrateco2 只用 C1 与 Ak。C0/C2 是否为历史遗留（旧 Newton-Raphson 路径）或供其他调用方使用——当前模块内无其他使用点。推断：2017 简化后遗留。
Ak 在 coefffunc 中既是 intent(out) 形参又在局部声明区出现于注释（:135），实际为输出形参；已确认通过参数表回传给入口。