photolysis_online.F

路径

LMDZ.MARS\libf\aeronomars\photolysis_online.F

所属目录/模块

libf\aeronomars

文件定位

photolysis_online.F 定义 photolysis_online_mod，是火星 GCM 在线光解率计算模块。它消费 photolysis_mod 在 init_photolysis 中预读的光谱网格、太阳通量、温度相关截面和量子产额，对每一列大气逐波长做单色辐射传输（两流近似 + 伪球面修正），累加得到 v_phot(nlayer, nb_phot_max) 的 15 个光解率。

本文件本身不读取任何外部文件，也不做化学矩阵积分；它只负责把 photolysis_mod 的共享谱表数据转化为逐层光解率，由 photochemistry_mod 在 jonline=.true. 时消费。

定义的符号

依赖的模块

调用的关键例程

输入

输出

光解反应编号

吸收气体编号

共享状态与副作用

• 本模块无自有 SAVE 或 THREADPRIVATE 变量；所有光谱数据通过 use photolysis_mod 访问。
• rtlink 内有 logical, save :: delta 并标记 !$OMP THREADPRIVATE(delta)，用于缓存 delta 缩放选择。
• 无文件 I/O、无 abort_physic、无 MPI 通信；所有输入数据已由 photolysis_mod::init_photolysis 预读并广播。
• v_phot 先初始化为零，波长循环累加后统一清除 <1e-30 的小值。

核心逻辑

1. 工作垂直网格：将 nlayer 层 GCM 网格复制到 nlev=nlayer+1 层界面网格（zpress/zalt/ztemp/zmmean），顶层设 zpress(nlev)=0，顶层高度和温度用外推。
2. 空气柱和 Rayleigh：setair 用静力方程 colinc = avo * 0.1 * dp / (mmean * g) 计算空气柱增量（g=3.72 m/s²），用 Ityaksov et al. 2008 CO₂ Rayleigh 截面 σ = 1.78e-26 * ν^(4+0.625) 计算 dtrl。
3. 温度相关截面：对 O₂（4 温度点线性）、CO₂（3 温度点对数）、O₃（2 温度点线性）、H₂O₂（JPL 多项式）和 NO₂（2 温度点+量子产额）调用各 set* 子例程，填充 dtgas 和 sj。
4. 温度无关截面：对 H₂O、HO₂、H₂、NO、N₂ 直接在循环中用 colinc * rm * xs 填充 dtgas 和 sj（含产额 yieldh2/yieldno/yieldn2）。
5. 氘化学：若 deutchem=.true.，HDO 用 0.5*xshdo 分两个通道（j_hdo_od 和 j_hdo_d）填充 sj。
6. 总气体光学厚度：dagas 对所有 10 种吸收气体的 dtgas 逐层逐波长求和。
7. 气溶胶和云：setaer 设置固定气溶胶廓线（scaleh=10km、gamma=0.03 波长依赖、omega=0.622、g=0.88）；setcld 设零云光学厚度。
8. 球面几何：sphers 按 Dahlback & Stamnes 1991 计算斜路径长度比 dsdh 和层间折射索引 nid，火星半径 3393 km。
9. 太阳通量缩放：fmars = f * (1/dist_sol)^2。
10. 波长循环：对 iw=1:nw-1，调用 rtlink 做单色辐射传输，得到直射/漫射上下行辐照度 edir/edn/eup 和作用量通量 fdir/fdn/fup；计算球面作用量通量 saflux = fmars * (fdir+fdn+fup)；逐反应逐层累加 v_phot += saflux * sj * (wu-wl)。
11. 小值清除：where (v_phot < 1e-30) v_phot = 0。

rtlink / ps2str 细节

• rtlink 做 delta 缩放（delta 逻辑变量 save 并 THREADPRIVATE）：当 g² > 0.2 时启用，缩放 tau/omega 后调用 ps2str。
• ps2str 实现 Toon et al. 1989 两流近似，选择半球平均法（hemispheric mean），用 tridiag 求解 2*nlev 维三对角系统。
• 伪球面修正：直射光束用 sphers 的 dsdh 斜路径衰减，而非平面平行近似。

伪代码

photolysis_online(nlayer, ..., v_phot):
  nlev = nlayer + 1
  build interface-level grid zpress/zalt/ztemp/zmmean (top: p=0, extrapolate z/T)

  call setair → colinc, dtrl (Rayleigh, Ityaksov 2008 CO2)

  dtgas = 0
  call seto2  → dtgas(:,:  ,a_o2),  sj(:,:,:  )  [T-dep, 4-point linear, Schumann-Runge]
  call setco2 → dtgas(:,:  ,a_co2), sj(:,:,:  )  [T-dep, 3-point log, O(1D)<167nm]
  call seto3  → dtgas(:,:  ,a_o3),  sj(:,:,:  )  [T-dep, 2-point linear, JPL yield]
  call seth2o2→ dtgas(:,:  ,a_h2o2),sj(:,:,:  )  [T-dep, JPL polynomial]
  call setno2 → dtgas(:,:  ,a_no2), sj(:,:,:  )  [T-dep, 2-point + yield]

  for ilay, iw:  # T-independent
    dtgas(ilay,iw,a_h2o) = colinc * rm(i_h2o) * xsh2o
    dtgas(ilay,iw,a_ho2) = colinc * rm(i_ho2) * xsho2
    dtgas(ilay,iw,a_h2)  = colinc * rm(i_h2)  * xsh2  * yieldh2
    dtgas(ilay,iw,a_no)  = colinc * rm(i_no)  * xsno  * yieldno
    dtgas(ilay,iw,a_n2)  = colinc * rm(i_n2)  * xsn2  * yieldn2
    sj(ilay,iw,j_*) = corresponding xs * yield

  if deutchem: sj(:,:,j_hdo_od) = sj(:,:,j_hdo_d) = 0.5 * xshdo

  dagas = sum(dtgas over all absorbers)

  call setaer → dtaer, omaer, gaer (fixed aerosol profile)
  call setcld → dtcld=0, omcld, gcld (no clouds)
  call sphers → dsdh, nid (spherical geometry, R=3393km)

  fmars = f * (1/dist_sol)^2
  v_phot = 0

  for iw = 1..nw-1:
    call rtlink → edir, edn, eup, fdir, fdn, fup (two-stream, pseudo-spherical)
    saflux = fmars * (fdir + fdn + fup)
    for i = 1..nphot:
      for ilay = 1..nlayer:
        v_phot(ilay,i) += saflux(ilay) * sj(ilay,iw,i) * (wu(iw) - wl(iw))

  where (v_phot < 1e-30) v_phot = 0

参与的主题流程

写法特点

• 固定格式 Fortran（.F），保留大量 TUV 风格注释块和 ASCII 装饰线。
• 全部内部子例程通过 contains 封装在 photolysis_online 内部，共享主例程的局部变量（隐式宿主关联），不通过参数传递 nw/wl/wc 等——而是通过 use photolysis_mod 模块变量访问。
• rtlink 中 delta 用 save + THREADPRIVATE 缓存，首次调用计算后保持不变。
• 火星重力 g=3.72 m/s² 和 Avogadro 常数 avo=6.022e23 在 setair 中硬编码为 parameter。
• 火星半径 3393 km 在 sphers 中硬编码。
• setcld 当前设 dtcld=0，但已预留 omcld=0.99 和 gcld=0.85 参数，为未来云光学模块留出接口。
• ps2str 实现 4 种两流方法（hemispheric mean、Eddington、quadrature、P-I），但硬编码选择 1=hemispheric mean。
• seto2 中 Schumann-Runge 波段（180-200 nm）低分辨率修正因子 factor=0.8（mopt==2）。

复现要点

• 调用前必须已完成 photolysis_mod::init_photolysis，否则所有 xs*/yield*/f/albedo 数组未初始化。
• nphot 由 photolysis_mod 模块变量提供；基础值 13，deutchem 开启时 calchim_mod 会设为 15，nb_phot_max 需 ≥ nphot。
• rm 混合比的物种索引必须与 photochemistry_mod 中传入的 i_co2/i_o2/i_o3 等一致。
• sza 在线路径的日夜阈值为 113°（由 photochemistry_mod 控制），不同于离线 photolysis 的 95° 阈值。
• 输出 v_phot 的小值清除阈值为 1e-30；小于此值的光解率被置零。
• dist_sol 缩放以 1 AU 为基准（(1/dist_sol)^2），与离线路径的 (1.52/dist_sol)^2 基准不同。

待确认

• rtlink 的 delta 保存变量在多次调用间保持首次计算值，是否在不同 SZA/气溶胶条件下需要重算待确认。
• setcld 预留的 omcld=0.99 和 gcld=0.85 参数来源文献待确认。
• ps2str 硬编码选择 hemispheric mean 方法，是否曾使用或计划使用其他方法待确认。
• 离线 photolysis 的 sza 阈值为 95° 而在线路径为 113°（由调用方控制），差异是否为有意设计待确认。

相关页面

• photolysis_mod：上游共享谱表模块；提供全部光谱输入数据。
• photochemistry_mod：下游使用点，jonline=.true. 时调用本例程获取 v_phot。
• photolysis：离线光解路径；与本文件在线路径互斥。
• calchim_mod：高层化学调度，控制在线/离线光解选择。
• chemistrydata：离线光解查找表数据模块。
• aeronomars/index：aeronomars 目录总览。