photochemistry_mod.F90

路径

LMDZ.MARS\libf\aeronomars\photochemistry_mod.F90

所属目录/模块

libf\aeronomars

文件定位

photochemistry_mod.F90 定义 photochemistry_mod，是 calchim_mod 逐列调用的主光化学积分器。它把 GCM tracer mole fraction zycol 映射到内部化学物种数组，计算在线或离线光解率、可选离子光电离率、气相/异相反应速率和反应索引矩阵，然后在 lswitch-1 以下逐层用自适应子时间步和 LAPACK dgesv 解隐式线性系统，最后把数密度结果写回 zycol。

该文件实现的是 ASIS 光化学求解框架的 Mars 反应网络，核心对象包括基础 C/O/H/N 光化学、可选氘化学、可选离子化学、冰/尘异相项以及 NO/O2 夜辉诊断。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• photochemistry 和 gcmtochim 各有 firstcall，均为 SAVE/THREADPRIVATE；首次调用后缓存内部物种索引并不再重复分配/检查。
• indice 会 allocate(indice_phot/indice_3/indice_4)，这些数组来自 types_asis，用于后续矩阵填充；页面未发现释放路径。
• indice 每次 firstcall 打印 nb_phot/nb_reaction_4/nb_reaction_3，并在计数不等于传入容量时中止。
• 在线离子路径写入并读取 param_v4_h::jion，依赖 jthermcalc_e107 和 phdisrate 的共享状态副作用。
• 未定义编译宏 LAPACK 时，主求解和 define_dt 的矩阵求解分支都会调用 abort_physic。
• reactionrates 内部硬编码 hetero_dust=.false.、hetero_ice=.true.，因此冰异相反应活跃，尘异相反应默认关闭。

核心逻辑

1. firstcall 设置内部物种索引：基础 17 个物种固定为 CO2、CO、O、O1D、O2、O3、H、H2、OH、HO2、H2O2、H2O、N、N2D、NO、NO2、N2；若 ionchem 为真追加 15 个离子/电子物种，若 deutchem 为真追加 6 个氘相关物种。随后调用 indice 建立反应化学计量表。
2. GCM tracer 到内部数组：gcmtochim 首次检查所有必需 tracer_mod::igcm_* 索引，之后把 zycol 的对应 tracer 拷贝到内部 rm，低于 1.e-30 的 mole fraction 清零，再乘 dens 得到数密度 c。
3. 光解和光电离率：jonline 时若 sza<=113 调 photolysis_online；夜间直接清零 v_phot。若 ionchem，再调 jthermcalc_e107 和逐层 phdisrate(ig,nlayer,2,sza,ilay)，把 jion 的 CO2/O2/O/NO/CO/N2/N/H 光电离通道放入 v_phot(:,nphot+1:nphot+18)。jonline=.false. 时改调 photolysis，并先把 tau 按 7 hPa 查找表口径重标定。
4. 保存诊断光解率：jo3=v_phot(:,5)，jh2o=v_phot(:,7)。
5. 反应速率表：reactionrates 先按 nphot+nphotion 初始化光解样反应数，再填 O、O(1D)、HOx、N/NOx、CO、氘、离子和异相反应。二体反应进 v_4，三体/同物种二次反应进 v_3，quench/异相项按光解样一阶反应进 v_phot。
6. 逐层隐式积分：对 ilev=1:lswitch-1，在 time<ptimestep 循环中先用 fill_matrix 形成矩阵、产生项和损失项，再由 define_dt 按误差准则修正子步，最后解 (I + J*dt)*cnew = c_old。低于 1.e-30*dens 的数密度置零。
7. 电荷中性修正：若 ionchem，每个子步后把电子数密度强制设为所有正离子数密度之和。
8. 夜辉和写回：积分结束后计算 em_no/em_o2，再调用 chimtogcm 将 1:lswitch-1 的内部数密度除以 dens 写回 zycol。

伪代码

photochemistry(...):
  if firstcall:
    assign fixed species indices
    append ion species if ionchem
    append deuterium species if deutchem
    call indice(...) to allocate/fill reaction stoichiometry tables
    firstcall = false

  call gcmtochim(...)  # zycol -> rm -> c

  if jonline:
    if sza <= 113:
      call photolysis_online(..., v_phot)
      if ionchem:
        call jthermcalc_e107(...)
        for ilay = 1..lswitch-1: call phdisrate(..., chemthermod=2, ilay)
        copy jion channels into v_phot after nphot
    else:
      v_phot = 0
  else:
    tau = tau * 7 / press(1)
    call photolysis(..., v_phot)

  jo3 = v_phot(:,5)
  jh2o = v_phot(:,7)
  hetero_dust = false
  hetero_ice = true
  call reactionrates(..., v_phot, v_3, v_4, ind_norec, ind_orec)

  for ilev = 1..lswitch-1:
    time = 0
    dt_guess = ptimestep
    while time < ptimestep:
      call fill_matrix(ilev, mat1, prod, loss, c, ...)
      call define_dt(..., dt_corrected, ...)
      cap dt_corrected at remaining physical timestep
      mat = identity + mat1*dt_corrected
      solve mat*cnew = c(ilev,:) with dgesv
      zero very small cnew
      c(ilev,:) = cnew
      if ionchem: electron = sum(positive ions)
      time += dt_corrected
      dt_guess = dt_corrected

  em_no = O*N*NO_recombination_rate
  em_o2 = 0.75*O*O*CO2*O_recombination_rate
  call chimtogcm(...)  # c -> zycol

参与的主题流程

写法特点

• 主例程包含 6 个内部子例程，indice_* 类型来自 types_asis，反应化学计量和速率计算分离。
• 编译必须有 LAPACK 宏和 dgesv 链接，否则运行时会在求解处中止。
• jionos 在主例程内硬编码为 .true.，没有外部配置输入。
• hetero_ice 硬编码 .true.，hetero_dust 硬编码 .false.。
• phychemrat=1，名义化学步长等于物理时间步；实际子步由 define_dt 自适应切分。
• chimtogcm 只写回 1:lswitch-1，因此 .not.unichim 时 lswitch:nlayer 留给热层化学路径处理。

复现要点

• nesp 必须与 ionchem/deutchem 对应的内部物种数一致：基础 17，离子化学加 15，氘化学加 6。
• nb_phot_max/nb_reaction_3_max/nb_reaction_4_max 必须与 indice 实际填充计数一致，否则 indice 会中止。
• 在线光解夜间阈值是 sza>113，而 phdisrate 自身还有 zenit>140 夜间判断；本文件进入在线夜间分支时已经把全部 v_phot 清零。
• 离线路径会原地修改局部参数 tau=tau*7/press(1)，这是查找表 7 hPa 口径换算。
• 若 ionchem 为真，param_read_e107/fill_data_thermos 等上游热层参数初始化必须已使 jthermcalc_e107 和 phdisrate 可用。

待确认

• reactionrates 中 i027 和 i056 速率分别被连续写入两次（1894-1898、2148-2152），未见 i028 对应赋值；是否为编号遗漏或有意复用待确认。
• phdisrate 在本文件中固定 chemthermod=2，因此不进入 chemthermod==3 的扩展通道；离子光化学是否需要这些通道需与物理设计确认。
• indice 分配的 types_asis 共享索引数组未在本文件释放；是否由进程生命周期管理待确认。
• comcstfi_h 在 reactionrates 中全模块 use，需进一步逐项反应确认具体使用的常数。

相关页面

• calchim_mod：上游调度方，准备列状态并调用 photochemistry。
• chemistrydata：离线光解查找表读取方。
• photolysis：离线光解查找表插值器，jonline=false 时生成基础 v_phot。
• photolysis_online.F：在线光解率计算。
• photolysis_mod：在线光解初始化和共享光解参数；本文件读取其 nphot 来布局光解槽位。
• types_asis：ASIS 反应索引类型和 indice_phot/indice_3/indice_4 共享数组定义。
• paramfoto_compact：提供 phdisrate，本文件离子路径读取其 jion 输出。
• param_read_e107：提供 E107 光化学表和分支比。
• iono_h：电子温度来自 calchim_mod 对 temp_elect 的调用，本文件参数 temp_elect(nlayer) 是数组而非函数。