lwflux.F

快速理解

它做什么： 汇总长波通量、冷却率和净辐射收支。被 lwmain 调用（lwmain_mod.F:160）。

基本过程： 合并 xi 与带外散射通量 lwdiff → 构造 ksidb → 调 lwdiff 算散射通量 → 汇入 netrad。

关键结果： 每层冷却率 coolrate、地表下行通量 fluxground、大气顶出射通量 fluxtop，供 lwi 半隐式修正。

路径

LMDZ.MARS\libf\phymars\lwflux.F

所属目录/模块

libf\phymars

文件定位

lwflux.F 定义 lwflux_mod 模块，提供长波 LTE 辐射主流程中的通量、冷却率和净辐射收支汇总例程 lwflux。它位于 lwmain_mod.F 的长波调用链后段：lwmain 先调用 lwb 生成谱带 Planck 函数和半层差分，再准备气体/气溶胶总光学性质，随后在第 160 行调用 lwflux。

本例程把 CO2 15 微米带内的净交换率 xi 与带外散射通量 lwdiff 合并。源码先按层-层、地面-层、层-太空、地面-太空组合构造 ksidb，把谱带 1 和 2 的 CO2 贡献合并到索引 3；再调用 lwdiff 计算 iir=3:nir 的扩散散射通量；最后把两部分汇入 netrad，输出每层冷却率 coolrate、地表下行长波通量 fluxground 和大气顶出射长波通量 fluxtop。

定义的符号

符号	类型	行号	作用
`lwflux_mod`	module	1	封装长波 LTE 通量、冷却率和辐射收支汇总例程。
`lwflux`	subroutine	7	合并 CO2 15 微米带内净交换率和带外扩散通量，输出 `coolrate/fluxground/fluxtop/netrad`。

依赖的模块

use 模块	only 列表	用途	待确认
`dimradmars_mod`	`ndlo2`, `nir`, `ndlon`, `nuco2`, `nflev`	提供辐射列数、红外谱带数、CO2 谱带数和垂直层维度；用于输入输出和局部工作数组声明。	否
`yomlw_h`	`nlaylte`, `xi`, `xi_ground`, `gcp`	提供 LTE 长波层数、净交换率表、地面邻层修正系数和 `g/cp` 换算因子。	否
`lwdiff_mod`	`lwdiff`	计算 CO2 15 微米带外 `iir=3:nir` 的红外散射通量贡献。	否

调用的关键例程

被调用例程	所在模块/文件	调用位置	作用
`lwdiff`	`lwdiff_mod` / `lwdiff.F`	`lwflux.F` 第 201 行	用 `bsurf/btop/dbsublay` 和总光学性质求带外散射通量，写入 `fluxdiff`。

输入

输入	来源	类型/维度	单位	含义
`ig0`	`lwmain`	`integer`	grid index offset	当前批量列在全局物理网格中的起始偏移；访问 `xi(ig0+jl,...)` 和 `xi_ground(ig0+jl,...)`。
`kdlon`	`lwmain`	`integer`	column	当前向量化处理的水平列数，循环范围为 `1:kdlon`。
`kflev`	`lwmain`	`integer`	layer	调用方传入的垂直层数；形参数组按它声明，核心循环用 `nlaylte`。
`dp`	`lwmain`	`real(ndlo2,kflev)`	Pa	层压强厚度；用于 `coolrate = gcp * netrad / dp`。
`bsurf`	`lwb` -> `lwmain`	`real(ndlo2,nir)`	Planck 函数值	地表谱带 Planck 函数；用于地面交换、地面到太空交换和 `lwdiff`。
`btop`	`lwb` -> `lwmain`	`real(ndlo2,nir)`	Planck 函数值	大气顶谱带 Planck 函数；传给 `lwdiff` 构造带外扩散顶边界。
`blev`	`lwb` -> `lwmain`	`real(ndlo2,nir,kflev+1)`	Planck 函数值	层界谱带 Planck 函数；本例程只在地面邻层修正中读 `blev(:,ja,1)`。
`blay`	`lwb` -> `lwmain`	`real(ndlo2,nir,kflev)`	Planck 函数值	层中心谱带 Planck 函数；用于所有净交换率乘积。
`dbsublay`	`lwb` -> `lwmain`	`real(ndlo2,nir,2*kflev)`	Planck 差分量	半层 Planck 差分；原样传给 `lwdiff`。
`tlay`, `tlev`, `dt0`	`lwmain`	`real(...)`	K	注释说明用于 G2D 输出；当前 `lwflux.F` 主体未读取这些参数。
`emis`	`lwmain`	`real(ndlo2)`	-	地表发射率；原样传给 `lwdiff`。
`tautotal`, `omegtotal`, `gtotal`	`lwmain`	`real(ndlo2,kflev,nir)`	optical depth / -	每层每红外谱带的总光学厚度、单次散射反照率和不对称因子；原样传给 `lwdiff`。

输出

输出	去向	类型/维度	单位	含义
`coolrate`	`lwmain` -> 下游加热率流程	`real(ndlo2,kflev)`	K/s	每层长波冷却率，按 `gcp * netrad / dp` 从净辐射收支换算。
`netrad`	`lwmain` -> `lwi`	`real(ndlo2,kflev)`	W/m2	每层净辐射收支，包含 CO2 带内净交换和 `lwdiff` 带外扩散通量散度。
`fluxground`	`lwmain`	`real(ndlo2)`	W/m2	到达地面的下行长波通量；由层到地面交换项减去 `fluxdiff(:,2,1)` 得到。
`fluxtop`	`lwmain`	`real(ndlo2)`	W/m2	大气顶出射长波通量；由到太空交换项加上 `fluxdiff(:,1,nlaylte+1)` 得到。

共享状态与副作用

lwflux.F 不定义保存型 module 变量，不读写文件，不读取运行时配置，也不打印诊断。它读取 yomlw_h 中的共享表 xi、xi_ground 和换算因子 gcp，并写入调用方传入的四个输出数组。

源码中有一个诊断分支 computeflux = .false.。当保持默认值时，flw_up/flw_dn/coefu/coefd 仅为局部未输出工作数组；如果源码维护者手动改成 .true.，该分支会计算各层宽带上下行通量，但当前没有对外输出或文件副作用。

核心逻辑

清零局部净交换率数组 ksidb。
对 i=1:nlaylte、j=i+1:nlaylte、ja=1:nuco2 构造层-层交换：xi(ig0+jl,ja,i,j) * (blay(j)-blay(i))，并写入反对称项 ksidb(j,i)=-ksidb(i,j)。
构造地面-层交换：xi(...,0,i) * (bsurf-blay(i))，再用 xi_ground 对地面和第一层之间的邻层贡献做修正。
构造层-太空交换和地面-太空交换，太空端 Planck 值按源码写法等效为零。
把 CO2 15 微米带索引 1 和 2 的交换率相加到 ksidb(:,3,:,:)。
先把 fluxdiff(:,1:2,nlaylte+1) 置零，再调用 lwdiff 计算带外散射通量。这个预清零是 lwdiff 回填顶层界面的接口前提。
对每层 i 汇总 ksidb(:,3,i,j) 得到带内 netrad，再叠加扩散通量散度：-fluxdiff(1,i+1)-fluxdiff(2,i+1)+fluxdiff(1,i)+fluxdiff(2,i)。
用 gcp * netrad / dp 得到冷却率。
对地面下行通量，累加所有层到地面的 xi * blay 贡献，再减去底界面向下扩散通量 fluxdiff(:,2,1)。
对顶层出射通量，累加地面和各层到太空的净交换贡献，再加上顶界面向上扩散通量 fluxdiff(:,1,nlaylte+1)。

伪代码

ksidb = 0

for each layer pair i < j and CO2 band ja:
  ksidb(i,j) = xi(i,j) * (blay(j) - blay(i))
  ksidb(j,i) = -ksidb(i,j)

for each layer i and CO2 band ja:
  ksidb(i,0) = xi(0,i) * (bsurf - blay(i))
  ksidb(0,i) = -ksidb(i,0)

adjust ground-first-layer exchange with xi_ground
build layer-space and ground-space exchange terms
ksidb(:,3,:,:) = ksidb(:,1,:,:) + ksidb(:,2,:,:)

zero fluxdiff at top interface
call lwdiff to fill diffusion flux interfaces

for each layer i:
  netrad(i) = sum_j ksidb(3,i,j)
  netrad(i) += diffusion_flux_divergence(i)
  coolrate(i) = gcp * netrad(i) / dp(i)

fluxground = layer_to_ground_exchange - bottom_downward_diffusion
fluxtop = exchange_to_space + top_upward_diffusion

参与的主题流程

主题	参与方式
辐射计算	长波 LTE 主流程的通量和冷却率汇总阶段；位于 `lwb/lwu/lwxd/lwxn/lwxb` 准备量之后，`lwi` 的温度导数/线性化步骤之前。
phymars 核心物理模块目录	属于 `libf\phymars` 的长波辐射辅助文件。

写法特点

固定格式 Fortran 文件，模块内只有一个子例程。
ksidb 使用下标 0:nflev+1 同时表示地面、物理层和太空边界；其中 0 是地面，nlaylte+1 是太空。
ksidb 对层-层、地面-层、层-太空、地面-太空交换都显式写反对称项，后续对 ksidb(:,3,i,j) 求和时只读取目标层对应的行。
ja=3 是源码用来保存 ja=1 与 ja=2 合计的工作槽；ksidb 因此声明为 nuco2+1。
computeflux 诊断分支在源码内硬编码为 .false.，不是运行时配置开关。

复现要点

ig0+jl 必须能落在 xi 和 xi_ground 的全局列索引范围内。
kdlon 不应超过 ndlo2/ndlon 静态数组容量。
kflev 的形参数组维度必须覆盖 nlaylte 层和 nlaylte+1 层界；源码循环上界使用 nlaylte。
dp(:,i) 不能为零或异常小，否则 coolrate = gcp * netrad / dp 会放大或产生除零风险。
复现风险：lwflux 只预清零 fluxdiff(:,:,nlaylte+1)，其余界面依赖 lwdiff 填写；若替换或跳过 lwdiff，必须明确初始化 fluxdiff。
复现风险：xi/xi_ground、bsurf/blay 的符号和单位约定直接决定 netrad/fluxground/fluxtop 的正负号；不要把 ksidb 的反对称项当作重复累加误删。

待确认

tlay、tlev、dt0 在当前源码主体中未读取，注释称仅用于 G2D 输出；是否仍需保留在接口中属于源码维护问题。
computeflux 分支计算的 flw_up/flw_dn 没有对外输出；是否曾用于历史诊断版本需要维护者确认。