nucleaco2.F90

路径

LMDZ.MARS\libf\phymars\nucleaco2.F90

所属目录/模块

libf\phymars

文件定位

该文件定义 nucleaco2_mod 模块，含一个 subroutine nucleaco2 和一个 module 级辅助 function fshapeco2。职责是：在给定 CO2 分压、温度、过饱和度和一组 CCN（凝结核）尺度分布的前提下，按经典异质成核理论计算 CO2 冰在固体基底上的异质成核率 nucrate(nbinco2_cld)（每个 CCN 尺度 bin 一个速率）。

文件头注释说明：算法来自 Pruppacher & Klett (1978) 水冰在固体基底成核的公式，由 Keese (JGR, 1989) 细化；水冰版本由 F. Montmessin 编写、J.-B. Madeleine 移植到 LMD/GCM（2011-10）、A. Spiga 优化（2012-02）；CO2 版本由 C. Listowski 与 J. Audouard 在水冰版本基础上改写（2016-2017）。

注意（核验自调用方）：本例程对“成核基底种类”不可知，由调用方通过形参 teta（接触参数 m=cosθ）和 n_ccn（核分布）区分。improvedco2clouds_mod.F90 用同一例程对三类基底各调用一次：尘埃 CCN、流星尘 CCN（均传 mtetaco2=0.95），以及 H2O 冰 CCN（传水的接触参数 mteta）。见“调用的关键例程 / 输入”。

定义的符号

依赖的模块

说明：fshapeco2 函数本身不 use 任何模块，仅做纯算术。

调用的关键例程

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

调用方（核验）：

输入

复现风险：n_ccn 的物理单位取决于调用方在各 bin 中放入的量（核验 improvedco2clouds_mod.F90 中由对数正态分布积分得到的数浓度 n_aer），本例程只把它作为线性权重乘进 nucrate，不做归一化。复现需保证 n_ccn 的定义与调用方一致。

输出

共享状态与副作用

• 无 common block，不写文件，无诊断打印。
• 读取 microphys_h 的 SAVE/THREADPRIVATE 变量 rad_cldco2（只读）及若干 PARAMETER 常数。
• 本例程内部无 SAVE 变量、无 firstcall；纯函数式（给定输入即给定输出，除依赖 rad_cldco2 当前值外）。
• fshapeco2 为纯函数。

核心逻辑

1. mtetalocal = dble(teta)：接触参数转双精度。
2. 计算环境量：
   • nco2 = pco2/(kbz*temp)：CO2 分子数密度（个/m³，理想气体）。
   • rstar = 2*sigco2*vo2co2/(kbz*temp*log(sat))：临界胚半径（Kelvin/开尔文成核半径，m）。log(sat) 要求 sat>1。
   • gstar = 4π*rstar³/(3*vo2co2)：临界胚所含分子数。
3. fshapeco2simple = (2+m)(1-m)²/4：平表面（CCN 半径 ≫ 临界半径）极限下的接触因子。
4. 遍历 i = 1..nbinco2_cld：
   • 若 n_ccn(i) < 1e-10：该 bin 无核，nucrate(i)=0，跳过。
   • 否则：
     • 选接触因子 zefshapeco2：若 rad_cldco2(i) > 3000*rstar（核远大于临界胚，近似平面）用 fshapeco2simple，否则调 fshapeco2(m, rad_cldco2(i)/rstar) 算曲率修正。
     • fistar = (4/3)π*sigco2*rstar²*zefshapeco2：形成临界胚的活化能（J）。
     • deltaf = (2*desorpco2 - surfdifco2 - fistar)/(kbz*temp)，并钳制到 [-100, 100]。
     • 若 deltaf == -100（被钳到下限，指数极小）：nucrate(i)=0。
     • 否则按经典成核率公式： nucrate(i) = sqrt(fistar/(3π*kbz*temp*gstar²)) * kbz*temp*rstar²*4π * (nco2*rad_cldco2(i))² / (zefshapeco2*nusco2*m0co2) * exp(deltaf)。

fshapeco2(cost, rap)（cost=m, rap=r_CCN/r*）：

• yeah = sqrt(1 - 2*cost*rap + rap²)（即 φ）。
• 三项相加再乘 1/2，得到经典异质成核几何因子 f(m,x)（Pruppacher & Klett 形式）。

伪代码

subroutine nucleaco2(pco2, temp, sat, n_ccn, nucrate, vo2co2, teta):
    m   = dble(teta)
    nco2  = pco2 / (kbz*temp)
    rstar = 2*sigco2*vo2co2 / (kbz*temp*log(sat))
    gstar = 4*pi*rstar^3 / (3*vo2co2)
    fsimple = (2+m)*(1-m)^2 / 4

    for i in 1..nbinco2_cld:
        if n_ccn(i) < 1e-10:
            nucrate(i) = 0
        else:
            if rad_cldco2(i) > 3000*rstar:
                fshape = fsimple                      # 平面近似
            else:
                fshape = fshapeco2(m, rad_cldco2(i)/rstar)

            fistar = (4/3)*pi*sigco2*rstar^2*fshape
            deltaf = (2*desorpco2 - surfdifco2 - fistar)/(kbz*temp)
            deltaf = clamp(deltaf, -100, 100)

            if deltaf == -100:
                nucrate(i) = 0
            else:
                nucrate(i) = sqrt(fistar/(3*pi*kbz*temp*gstar^2))
                           * kbz*temp*rstar^2*4*pi
                           * (nco2*rad_cldco2(i))^2
                           / (fshape*nusco2*m0co2)
                           * exp(deltaf)

function fshapeco2(cost, rap):           # f(m,x) 几何因子
    phi  = sqrt(1 - 2*cost*rap + rap^2)
    A    = 1 + ((1-cost*rap)/phi)^3
    y    = (rap-cost)/phi
    B    = rap^3*(2 - 3*y + y^3)
    C    = 3*cost*rap^2*(y-1)
    return 0.5*(A + B + C)

参与的主题流程

写法特点

• 自由格式 Fortran（.F90），续行用行尾 &。
• 主要算术用 DOUBLE PRECISION；temp、teta 为 REAL 输入，内部用 dble() 提升。
• 硬编码数值：1e-10（空核阈值）、3000.*rstar（平面近似阈值）、deltaf 钳制区间 [-100,100]、fshapeco2simple 的 /4 与 fshapeco2 末尾 0.5。
• 物理常数（sigco2、desorpco2、nusco2、surfdifco2 等）全部来自 microphys_h 的 PARAMETER，本文件不重复定义。
• deltaf 钳到 -100 即视作成核率为 0，是防 exp 下溢的数值保护；钳到 +100 防上溢但仍按公式算。
• 接触参数 teta 由调用方按基底种类传入（CO2 基底 0.95 vs 水冰基底 mteta），本例程不区分。

复现要点

• 临界半径公式要求 sat>1，否则 log(sat)<=0 使 rstar 非物理；调用方需保证仅在过饱和列调用（核验：improvedco2clouds 在 satu 相关分支内调用）。
• 三类基底（尘埃/流星尘/水冰）共用本例程，差异只在传入的 n_ccn 分布和 teta；复现时勿把 CO2 接触参数误用于水冰基底。
• vo2co2 在调用方用随温度变的 CO2 冰密度算出（m0co2/rho_ice_co2T），不是常数。
• rad_cldco2 是 THREADPRIVATE 的 bin 半径数组，由 CO2 云方案初始化（见 co2cloud_mod / microphys_h），本例程只读取。
• 成核率单位与“每核每秒”一致；调用方再乘 microtimestep 经泊松概率转成本步成核数。
• 注意 deltaf == -100. 的相等判断是浮点比较，但由于上一行用 min(max(...),100.d0) 钳制，下限恰好取到字面值 -100.，比较成立；复现需保持同样钳制+比较写法。

待确认

• n_ccn 各 bin 的精确单位（数浓度 vs 每 bin 绝对数）：本文件只作线性权重；需结合 improvedco2clouds_mod 的 n_aer 构造确认。
• 2*desorpco2 - surfdifco2 这一组合（Keese 1989 形式）中各项符号约定的物理出处，注释只给了数值来源（Bachnar 2016）。

相关页面

• improvedco2clouds_mod：CO2 云微物理核心，本例程唯一调用方。
• co2cloud_mod：CO2 云总调度。
• co2-saturation-helpers：CO2 饱和/凝结温度，提供过饱和度上下文。
• density_co2_ice：CO2 冰密度，vo2co2 计算的输入。
• CO2 循环父级主题页。
• microphys_h：定义本例程全部物理常数与 rad_cldco2。